From 0dac5f84ef15259016f66a4d0f1e9e7ce14d3dfb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Dec 2020 12:00:47 +0100 Subject: [PATCH 01/21] dummy data layout --- src/constitutive.f90 | 7 +++++++ src/constitutive_thermal.f90 | 7 +++++++ src/homogenization.f90 | 5 +++++ src/homogenization_thermal.f90 | 37 ++++++++++++++++++++++++++++++++++ 4 files changed, 56 insertions(+) create mode 100644 src/homogenization_thermal.f90 diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 667a23127..36400fca7 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -202,11 +202,17 @@ module constitutive real(pReal) :: T end function constitutive_thermal_T + module subroutine constitutive_mech_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F integer, intent(in) :: co, ip, el end subroutine constitutive_mech_setF + module subroutine constitutive_thermal_setT(T,co,ip,el) + real(pReal), intent(in) :: T + integer, intent(in) :: co, ip, el + end subroutine constitutive_thermal_setT + ! == cleaned:end =================================================================================== module function crystallite_stress(dt,co,ip,el) result(converged_) @@ -414,6 +420,7 @@ module constitutive constitutive_restartRead, & integrateSourceState, & constitutive_mech_setF, & + constitutive_thermal_setT, & constitutive_mech_getP, & constitutive_mech_getLp, & constitutive_mech_getF, & diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 1a05b983f..01d517124 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -139,4 +139,11 @@ module function constitutive_thermal_T(co,ip,el) result(T) end function constitutive_thermal_T +! setter for homogenization +module subroutine constitutive_thermal_setT(T,co,ip,el) + real(pReal), intent(in) :: T + integer, intent(in) :: co, ip, el +end subroutine constitutive_thermal_setT + + end submodule constitutive_thermal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 686bb9885..df7369096 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -27,6 +27,8 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point + real(pReal), dimension(:), allocatable, public :: & + homogenization_T real(pReal), dimension(:,:,:), allocatable, public :: & homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment @@ -56,6 +58,9 @@ module homogenization num_homog !< pointer to mechanical homogenization numerics data end subroutine mech_init + module subroutine thermal_init + end subroutine thermal_init + module subroutine mech_partition(subF,ip,el) real(pReal), intent(in), dimension(3,3) :: & subF diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 new file mode 100644 index 000000000..59e7357b6 --- /dev/null +++ b/src/homogenization_thermal.f90 @@ -0,0 +1,37 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!-------------------------------------------------------------------------------------------------- +submodule(homogenization) homogenization_thermal + + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief Allocate variables and set parameters. +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_init() + + print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' + + allocate(homogenization_T(discretization_nIPs*discretization_Nelems), source=0.0_pReal) + +end subroutine thermal_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Partition T onto the individual constituents. +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_partition(T,ip,el) + + real(pReal), intent(in) :: T + integer, intent(in) :: & + ip, & !< integration point + el !< element number + + + call constitutive_thermal_setT(T,1,ip,el) + +end subroutine thermal_partition + + +end submodule homogenization_thermal From bc12ac44c3ccda8616ea255a2ab173b17f989cbc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Dec 2020 12:34:00 +0100 Subject: [PATCH 02/21] basic functionality for thermal homogenization --- src/constitutive.f90 | 27 ++++++++++++----------- src/constitutive_mech.f90 | 7 +++--- src/constitutive_thermal.f90 | 42 +++++++++++++++++++++++++++--------- 3 files changed, 51 insertions(+), 25 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 36400fca7..e68d90de6 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -112,13 +112,15 @@ module constitutive interface ! == cleaned:begin ================================================================================= - module subroutine mech_init + module subroutine mech_init(phases) + class(tNode), pointer :: phases end subroutine mech_init module subroutine damage_init end subroutine damage_init - module subroutine thermal_init + module subroutine thermal_init(phases) + class(tNode), pointer :: phases end subroutine thermal_init @@ -197,10 +199,10 @@ module constitutive real(pReal), dimension(3,3) :: P end function constitutive_mech_getP - module function constitutive_thermal_T(co,ip,el) result(T) - integer, intent(in) :: co, ip, el + module function thermal_T(ph,me) result(T) + integer, intent(in) :: ph,me real(pReal) :: T - end function constitutive_thermal_T + end function thermal_T module subroutine constitutive_mech_setF(F,co,ip,el) @@ -463,6 +465,8 @@ subroutine constitutive_init phases + print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) + debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList) debugConstitutive%basic = debug_constitutive%contains('basic') debugConstitutive%extensive = debug_constitutive%contains('extensive') @@ -471,15 +475,14 @@ subroutine constitutive_init debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1) debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1) -!-------------------------------------------------------------------------------------------------- -! initialize constitutive laws - print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) - call mech_init - call damage_init - call thermal_init - phases => config_material%get('phase') + + call mech_init(phases) + call damage_init + call thermal_init(phases) + + constitutive_source_maxSizeDotState = 0 PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index f5a5cd0a2..97bf7d853 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -319,7 +319,10 @@ contains !> @brief Initialize mechanical field related constitutive models !> @details Initialize elasticity, plasticity and stiffness degradation models. !-------------------------------------------------------------------------------------------------- -module subroutine mech_init +module subroutine mech_init(phases) + + class(tNode), pointer :: & + phases integer :: & el, & @@ -331,7 +334,6 @@ module subroutine mech_init Nconstituents class(tNode), pointer :: & num_crystallite, & - phases, & phase, & mech, & elastic, & @@ -341,7 +343,6 @@ module subroutine mech_init !------------------------------------------------------------------------------------------------- ! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? - phases => config_material%get('phase') allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) allocate(phase_elasticityInstance(phases%length), source = 0) allocate(phase_NstiffnessDegradations(phases%length),source=0) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 01d517124..bdcd0bc26 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -2,6 +2,12 @@ !> @brief internal microstructure state for all thermal sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_thermal + + type :: tDataContainer + real(pReal), dimension(:), allocatable :: T + end type tDataContainer + + type(tDataContainer), dimension(:), allocatable :: current interface @@ -49,8 +55,29 @@ contains !---------------------------------------------------------------------------------------------- !< @brief initializes thermal sources and kinematics mechanism !---------------------------------------------------------------------------------------------- -module subroutine thermal_init +module subroutine thermal_init(phases) + + class(tNode), pointer :: & + phases + + integer :: & + ph, & + Nconstituents + + print'(/,a)', ' <<<+- constitutive_mech init -+>>>' + + allocate(current(phases%length)) + + + do ph = 1, phases%length + + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + + allocate(current(ph)%T(Nconstituents)) + + enddo + ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID @@ -122,21 +149,16 @@ end subroutine constitutive_thermal_getRateAndItsTangents - ! getter for non-thermal (e.g. mech) -module function constitutive_thermal_T(co,ip,el) result(T) +module function thermal_T(ph,me) result(T) - integer, intent(in) :: co, ip, el + integer, intent(in) :: ph, me real(pReal) :: T - integer :: ho, tme - ho = material_homogenizationAt(el) - tme = material_homogenizationMemberAt(ip,el) + T = current(ph)%T(me) - T = temperature(ho)%p(tme) - -end function constitutive_thermal_T +end function thermal_T ! setter for homogenization From 8c6d759b557b8fed05b7b53c01ad7a28d349afba Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Dec 2020 12:45:08 +0100 Subject: [PATCH 03/21] consistent naming --- src/constitutive.f90 | 31 +++++++++++++++---------------- src/constitutive_mech.f90 | 26 +++++++++++++------------- src/constitutive_thermal.f90 | 25 +++++++++++-------------- 3 files changed, 39 insertions(+), 43 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e68d90de6..05653e3b4 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -174,25 +174,25 @@ module constitutive end subroutine mech_restartRead - module function constitutive_mech_getS(co,ip,el) result(S) - integer, intent(in) :: co, ip, el + module function mech_S(ph,me) result(S) + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: S - end function constitutive_mech_getS + end function mech_S - module function constitutive_mech_getLp(co,ip,el) result(Lp) - integer, intent(in) :: co, ip, el - real(pReal), dimension(3,3) :: Lp - end function constitutive_mech_getLp + module function mech_L_p(ph,me) result(L_p) + integer, intent(in) :: ph,me + real(pReal), dimension(3,3) :: L_p + end function mech_L_p module function constitutive_mech_getF(co,ip,el) result(F) integer, intent(in) :: co, ip, el real(pReal), dimension(3,3) :: F end function constitutive_mech_getF - module function constitutive_mech_getF_e(co,ip,el) result(F_e) - integer, intent(in) :: co, ip, el + module function mech_F_e(ph,me) result(F_e) + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: F_e - end function constitutive_mech_getF_e + end function mech_F_e module function constitutive_mech_getP(co,ip,el) result(P) integer, intent(in) :: co, ip, el @@ -421,12 +421,10 @@ module constitutive constitutive_restartWrite, & constitutive_restartRead, & integrateSourceState, & - constitutive_mech_setF, & constitutive_thermal_setT, & constitutive_mech_getP, & - constitutive_mech_getLp, & + constitutive_mech_setF, & constitutive_mech_getF, & - constitutive_mech_getS, & constitutive_initializeRestorationPoints, & constitutive_windForward, & PLASTICITY_UNDEFINED_ID, & @@ -667,7 +665,8 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) sourceType: select case (phase_source(so,ph)) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState(constitutive_mech_getS(co,ip,el), co, ip, el) ! correct stress? + call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& + co, ip, el) ! correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState(co, ip, el) @@ -1020,7 +1019,7 @@ subroutine crystallite_orientations(co,ip,el) call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& - constitutive_mech_getF_e(co,ip,el)))) + mech_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))))) if (plasticState(material_phaseAt(1,el))%nonlocal) & call plastic_nonlocal_updateCompatibility(crystallite_orientation, & @@ -1123,7 +1122,7 @@ function integrateSourceState(dt,co,ip,el) result(broken) enddo if(converged_) then - broken = constitutive_damage_deltaState(constitutive_mech_getF_e(co,ip,el),co,ip,el,ph,me) + broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) exit iteration endif diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 97bf7d853..f5be4863a 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1847,27 +1847,27 @@ end subroutine mech_restartRead ! getter for non-mech (e.g. thermal) -module function constitutive_mech_getS(co,ip,el) result(S) +module function mech_S(ph,me) result(S) - integer, intent(in) :: co, ip, el + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: S - S = constitutive_mech_S(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + S = constitutive_mech_S(ph)%data(1:3,1:3,me) -end function constitutive_mech_getS +end function mech_S ! getter for non-mech (e.g. thermal) -module function constitutive_mech_getLp(co,ip,el) result(Lp) +module function mech_L_p(ph,me) result(L_p) - integer, intent(in) :: co, ip, el - real(pReal), dimension(3,3) :: Lp + integer, intent(in) :: ph,me + real(pReal), dimension(3,3) :: L_p - Lp = constitutive_mech_Lp(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + L_p = constitutive_mech_Lp(ph)%data(1:3,1:3,me) -end function constitutive_mech_getLp +end function mech_L_p ! getter for non-mech (e.g. thermal) @@ -1883,15 +1883,15 @@ end function constitutive_mech_getF ! getter for non-mech (e.g. thermal) -module function constitutive_mech_getF_e(co,ip,el) result(F_e) +module function mech_F_e(ph,me) result(F_e) - integer, intent(in) :: co, ip, el + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: F_e - F_e = constitutive_mech_Fe(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + F_e = constitutive_mech_Fe(ph)%data(1:3,1:3,me) -end function constitutive_mech_getF_e +end function mech_F_e diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index bdcd0bc26..9e2807d42 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -109,32 +109,29 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, real(pReal) :: & my_Tdot, & my_dTdot_dT - real(pReal), dimension(3,3) :: Lp, S integer :: & - phase, & + ph, & homog, & instance, & - grain, & - source, & - constituent + me, & + so, & + co homog = material_homogenizationAt(el) instance = thermal_typeInstance(homog) - do grain = 1, homogenization_Nconstituents(homog) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) + do co = 1, homogenization_Nconstituents(homog) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + do so = 1, phase_Nsources(ph) + select case(phase_source(so,ph)) case (SOURCE_thermal_dissipation_ID) - Lp = constitutive_mech_getLp(grain,ip,el) - S = constitutive_mech_getS(grain,ip,el) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - S, Lp, phase) + mech_S(ph,me),mech_L_p(ph,me), ph) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - phase, constituent) + ph, me) case default my_Tdot = 0.0_pReal From f9f56a1755c6b64251d357c31047c98d34ce10cd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Dec 2020 13:57:37 +0100 Subject: [PATCH 04/21] documenting --- src/constitutive.f90 | 34 ---------------------------------- src/constitutive_mech.f90 | 24 +++++++++++++++++------- src/constitutive_thermal.f90 | 26 +++++++++++++++++--------- 3 files changed, 34 insertions(+), 50 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 05653e3b4..6348c18d6 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -352,23 +352,6 @@ module constitutive end subroutine source_damage_isoBrittle_deltaState - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - end subroutine constitutive_plastic_LpAndItsTangents - - module subroutine constitutive_plastic_dependentState(co,ip,el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -376,23 +359,6 @@ module constitutive el !< element end subroutine constitutive_plastic_dependentState - - - module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress tensor - real(pReal), intent(out), dimension(3,3,3,3) :: & - dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - end subroutine constitutive_hooke_SandItsTangents - end interface diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index f5be4863a..fedec379f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -530,7 +530,7 @@ end function plastic_active !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & +subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi, co, ip, el) integer, intent(in) :: & @@ -616,7 +616,7 @@ end subroutine constitutive_plastic_dependentState ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & +subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -1846,7 +1846,9 @@ module subroutine mech_restartRead(groupHandle,ph) end subroutine mech_restartRead -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) +!---------------------------------------------------------------------------------------------- module function mech_S(ph,me) result(S) integer, intent(in) :: ph,me @@ -1858,7 +1860,9 @@ module function mech_S(ph,me) result(S) end function mech_S -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get plastic velocity gradient (for use by non-mech physics) +!---------------------------------------------------------------------------------------------- module function mech_L_p(ph,me) result(L_p) integer, intent(in) :: ph,me @@ -1870,7 +1874,9 @@ module function mech_L_p(ph,me) result(L_p) end function mech_L_p -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get deformation gradient (for use by homogenization) +!---------------------------------------------------------------------------------------------- module function constitutive_mech_getF(co,ip,el) result(F) integer, intent(in) :: co, ip, el @@ -1882,7 +1888,9 @@ module function constitutive_mech_getF(co,ip,el) result(F) end function constitutive_mech_getF -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get elastic deformation gradient (for use by non-mech physics) +!---------------------------------------------------------------------------------------------- module function mech_F_e(ph,me) result(F_e) integer, intent(in) :: ph,me @@ -1895,7 +1903,9 @@ end function mech_F_e -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get second Piola-Kichhoff stress (for use by homogenization) +!---------------------------------------------------------------------------------------------- module function constitutive_mech_getP(co,ip,el) result(P) integer, intent(in) :: co, ip, el diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 9e2807d42..f2b61fb26 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -2,11 +2,11 @@ !> @brief internal microstructure state for all thermal sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_thermal - + type :: tDataContainer real(pReal), dimension(:), allocatable :: T end type tDataContainer - + type(tDataContainer), dimension(:), allocatable :: current interface @@ -56,10 +56,10 @@ contains !< @brief initializes thermal sources and kinematics mechanism !---------------------------------------------------------------------------------------------- module subroutine thermal_init(phases) - + class(tNode), pointer :: & phases - + integer :: & ph, & Nconstituents @@ -71,13 +71,13 @@ module subroutine thermal_init(phases) do ph = 1, phases%length - + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs allocate(current(ph)%T(Nconstituents)) enddo - + ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID @@ -145,8 +145,9 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, end subroutine constitutive_thermal_getRateAndItsTangents - -! getter for non-thermal (e.g. mech) +!---------------------------------------------------------------------------------------------- +!< @brief Get temperature (for use by non-thermal physics) +!---------------------------------------------------------------------------------------------- module function thermal_T(ph,me) result(T) integer, intent(in) :: ph, me @@ -158,10 +159,17 @@ module function thermal_T(ph,me) result(T) end function thermal_T -! setter for homogenization +!---------------------------------------------------------------------------------------------- +!< @brief Set temperature +!---------------------------------------------------------------------------------------------- module subroutine constitutive_thermal_setT(T,co,ip,el) + real(pReal), intent(in) :: T integer, intent(in) :: co, ip, el + + + current(material_phaseAt(co,el))%T(material_phaseMemberAt(co,ip,el)) = T + end subroutine constitutive_thermal_setT From a1facadf3fb9451dbc3402e2e8c18dd975019cab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Dec 2020 18:08:19 +0100 Subject: [PATCH 05/21] needed for MSC.Marc --- src/commercialFEM_fileList.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index d8ab6390d..371c85fd3 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -52,4 +52,5 @@ #include "homogenization_mech_none.f90" #include "homogenization_mech_isostrain.f90" #include "homogenization_mech_RGC.f90" +#include "homogenization_thermal.f90" #include "CPFEM.f90" From 92ec10b2518c488ec01e3246c373005c816ed74f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Dec 2020 07:46:26 +0100 Subject: [PATCH 06/21] consistent names --- src/homogenization_mech.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 783d08dd1..dd3e47c59 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -113,24 +113,24 @@ module subroutine mech_partition(subF,ip,el) el !< element number integer :: co - real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt(el))) :: Fs chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - F(1:3,1:3,1) = subF + Fs(1:3,1:3,1) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mech_isostrain_partitionDeformation(F,subF) + call mech_isostrain_partitionDeformation(Fs,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mech_RGC_partitionDeformation(F,subF,ip,el) + call mech_RGC_partitionDeformation(Fs,subF,ip,el) end select chosenHomogenization do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - call constitutive_mech_setF(F(1:3,1:3,co),co,ip,el) + call constitutive_mech_setF(Fs(1:3,1:3,co),co,ip,el) enddo From ebc4f671c88d3b3d2a12eb73464622d6ac4f63e7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Dec 2020 08:00:20 +0100 Subject: [PATCH 07/21] names follow structure --- src/commercialFEM_fileList.f90 | 4 ++-- ...l_dissipation.f90 => constitutive_thermal_dissipation.f90} | 4 ++-- ...externalheat.f90 => constitutive_thermal_externalheat.f90} | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) rename src/{source_thermal_dissipation.f90 => constitutive_thermal_dissipation.f90} (97%) rename src/{source_thermal_externalheat.f90 => constitutive_thermal_externalheat.f90} (98%) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 371c85fd3..a27abec79 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -33,8 +33,8 @@ #include "constitutive_plastic_disloTungsten.f90" #include "constitutive_plastic_nonlocal.f90" #include "constitutive_thermal.f90" -#include "source_thermal_dissipation.f90" -#include "source_thermal_externalheat.f90" +#include "constitutive_thermal_dissipation.f90" +#include "constitutive_thermal_externalheat.f90" #include "kinematics_thermal_expansion.f90" #include "constitutive_damage.f90" #include "source_damage_isoBrittle.f90" diff --git a/src/source_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 similarity index 97% rename from src/source_thermal_dissipation.f90 rename to src/constitutive_thermal_dissipation.f90 index f28567aa7..27653a9ef 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) source_thermal_dissipation +submodule(constitutive:constitutive_thermal) source_dissipation integer, dimension(:), allocatable :: & source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? @@ -96,4 +96,4 @@ module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT end subroutine source_thermal_dissipation_getRateAndItsTangent -end submodule source_thermal_dissipation +end submodule source_dissipation diff --git a/src/source_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 similarity index 98% rename from src/source_thermal_externalheat.f90 rename to src/constitutive_thermal_externalheat.f90 index 9ba4a051b..3ef96790e 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) source_thermal_externalheat +submodule(constitutive:constitutive_thermal) source_externalheat integer, dimension(:), allocatable :: & @@ -135,4 +135,4 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d end subroutine source_thermal_externalheat_getRateAndItsTangent -end submodule source_thermal_externalheat +end submodule source_externalheat From 228398e78709d9cbbbb9bf37af9a5aa40c20d3a5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Dec 2020 09:10:59 +0100 Subject: [PATCH 08/21] config follows structure --- src/lattice.f90 | 83 ++++++++++++++++++++++++++----------------------- 1 file changed, 44 insertions(+), 39 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 6af135e4e..c9b6d99ef 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -453,12 +453,13 @@ contains !-------------------------------------------------------------------------------------------------- subroutine lattice_init - integer :: Nphases, p,i + integer :: Nphases, ph,i class(tNode), pointer :: & phases, & phase, & mech, & - elasticity + elasticity, & + thermal print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) @@ -476,67 +477,71 @@ subroutine lattice_init lattice_mu, lattice_nu,& source=[(0.0_pReal,i=1,Nphases)]) - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length + phase => phases%get(ph) mech => phase%get('mechanics') elasticity => mech%get('elasticity') - lattice_C66(1,1,p) = elasticity%get_asFloat('C_11') - lattice_C66(1,2,p) = elasticity%get_asFloat('C_12') + lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11') + lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12') - lattice_C66(1,3,p) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal) - lattice_C66(2,2,p) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal) - lattice_C66(2,3,p) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal) - lattice_C66(3,3,p) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) - lattice_C66(4,4,p) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal) - lattice_C66(5,5,p) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal) - lattice_C66(6,6,p) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) + lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal) + lattice_C66(2,2,ph) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal) + lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal) + lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) + lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal) + lattice_C66(5,5,ph) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal) + lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) select case(phase%get_asString('lattice')) case('cF') - lattice_structure(p) = lattice_FCC_ID + lattice_structure(ph) = lattice_FCC_ID case('cI') - lattice_structure(p) = lattice_BCC_ID + lattice_structure(ph) = lattice_BCC_ID case('hP') - lattice_structure(p) = lattice_HEX_ID + lattice_structure(ph) = lattice_HEX_ID case('tI') - lattice_structure(p) = lattice_BCT_ID + lattice_structure(ph) = lattice_BCT_ID case('oP') - lattice_structure(p) = lattice_ORT_ID + lattice_structure(ph) = lattice_ORT_ID case('aP') - lattice_structure(p) = lattice_ISO_ID + lattice_structure(ph) = lattice_ISO_ID case default call IO_error(130,ext_msg='lattice_init: '//phase%get_asString('lattice')) end select - lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice')) + lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice')) - lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt') - lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt') + lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt') + lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt') - lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation + lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation do i = 1, 6 - if (abs(lattice_C66(i,i,p)) phase%get('thermal') + lattice_K(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) + lattice_K(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal) + lattice_K(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) + lattice_K(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,ph), & + phase%get_asString('lattice')) + lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal) + endif + + + lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) + lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) + lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) + lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), & phase%get_asString('lattice')) - lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal) - lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal) - - lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) - lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) - lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) - lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), & - phase%get_asString('lattice')) - - lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal) + lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal) ! SHOULD NOT BE PART OF LATTICE END call selfTest From a2d0a9e51152b302c50f670c4c77a39a30b81683 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Dec 2020 09:54:13 +0100 Subject: [PATCH 09/21] WIP: separating states --- src/constitutive.f90 | 166 ++++++++++++++++++++++++++++++----- src/constitutive_mech.f90 | 11 ++- src/constitutive_thermal.f90 | 64 ++++++++++++++ 3 files changed, 215 insertions(+), 26 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6348c18d6..a13904697 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -102,7 +102,7 @@ module constitutive type(tPlasticState), allocatable, dimension(:), public :: & plasticState type(tSourceState), allocatable, dimension(:), public :: & - sourceState + sourceState, thermalState integer, public, protected :: & @@ -139,21 +139,37 @@ module constitutive integer, intent(in) :: ph, me end subroutine mech_initializeRestorationPoints - module subroutine constitutive_mech_windForward(ph,me) + module subroutine thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me - end subroutine constitutive_mech_windForward + end subroutine thermal_initializeRestorationPoints + + + module subroutine mech_windForward(ph,me) + integer, intent(in) :: ph, me + end subroutine mech_windForward + + module subroutine thermal_windForward(ph,me) + integer, intent(in) :: ph, me + end subroutine thermal_windForward + + + module subroutine mech_forward() + end subroutine mech_forward + + module subroutine thermal_forward() + end subroutine thermal_forward - module subroutine constitutive_mech_forward - end subroutine constitutive_mech_forward module subroutine mech_restore(ip,el,includeL) - integer, intent(in) :: & - ip, & - el - logical, intent(in) :: & - includeL + integer, intent(in) :: ip, el + logical, intent(in) :: includeL end subroutine mech_restore + module subroutine thermal_restore(ip,el) + integer, intent(in) :: ip, el + end subroutine thermal_restore + + module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -776,6 +792,7 @@ subroutine constitutive_restore(ip,el,includeL) enddo call mech_restore(ip,el,includeL) + call thermal_restore(ip,el) end subroutine constitutive_restore @@ -784,12 +801,13 @@ end subroutine constitutive_restore !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -subroutine constitutive_forward +subroutine constitutive_forward() integer :: ph, so - call constitutive_mech_forward() + call mech_forward() + call thermal_forward() do ph = 1, size(sourceState) do so = 1,phase_Nsources(ph) @@ -802,7 +820,7 @@ end subroutine constitutive_forward !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine constitutive_results +subroutine constitutive_results() integer :: ph character(len=:), allocatable :: group @@ -826,7 +844,7 @@ end subroutine constitutive_results !-------------------------------------------------------------------------------------------------- !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- -subroutine crystallite_init +subroutine crystallite_init() integer :: & ph, & @@ -937,10 +955,12 @@ subroutine constitutive_initializeRestorationPoints(ip,el) me = material_phaseMemberAt(co,ip,el) call mech_initializeRestorationPoints(ph,me) + call thermal_initializeRestorationPoints(ph,me) - do so = 1, phase_Nsources(material_phaseAt(co,el)) + do so = 1, size(sourceState(ph)%p) sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me) enddo + enddo end subroutine constitutive_initializeRestorationPoints @@ -964,10 +984,13 @@ subroutine constitutive_windForward(ip,el) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - call constitutive_mech_windForward(ph,me) + call mech_windForward(ph,me) + call thermal_windForward(ph,me) + do so = 1, phase_Nsources(material_phaseAt(co,el)) sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state(:,me) enddo + enddo end subroutine constitutive_windForward @@ -1049,8 +1072,7 @@ function integrateSourceState(dt,co,ip,el) result(broken) me = material_phaseMemberAt(co,ip,el) converged_ = .true. - broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) if(broken) return do so = 1, phase_Nsources(ph) @@ -1067,8 +1089,7 @@ function integrateSourceState(dt,co,ip,el) result(broken) source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) enddo - broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) if(broken) exit iteration do so = 1, phase_Nsources(ph) @@ -1122,6 +1143,111 @@ function integrateSourceState(dt,co,ip,el) result(broken) end function integrateSourceState + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +function integrateThermalState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + so + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + broken, converged_ + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_thermal_collectDotState(ph,me) + if(broken) return + + do so = 1, phase_Nsources(ph) + size_so(so) = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) & + + thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do so = 1, phase_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me) + enddo + + broken = constitutive_thermal_collectDotState(ph,me) + broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) exit iteration + + do so = 1, phase_Nsources(ph) + zeta = damper(thermalState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) & + - thermalState(ph)%p(so)%subState0(1:size_so(so),me) & + - thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt + thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + thermalState(ph)%p(so)%state(1:size_so(so),me), & + thermalState(ph)%p(so)%atol(1:size_so(so))) + enddo + + if(converged_) then + broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) + exit iteration + endif + + enddo iteration + + broken = broken .or. .not. converged_ + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end function integrateThermalState + + !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index fedec379f..67b6f9fbe 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1485,7 +1485,7 @@ end subroutine mech_initializeRestorationPoints !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_mech_windForward(ph,me) +module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me @@ -1499,14 +1499,14 @@ module subroutine constitutive_mech_windForward(ph,me) plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me) -end subroutine constitutive_mech_windForward +end subroutine mech_windForward !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_mech_forward() +module subroutine mech_forward() integer :: ph @@ -1521,7 +1521,7 @@ module subroutine constitutive_mech_forward() plasticState(ph)%state0 = plasticState(ph)%state enddo -end subroutine constitutive_mech_forward +end subroutine mech_forward @@ -1678,8 +1678,7 @@ module subroutine mech_restore(ip,el,includeL) constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) - plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = & - plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el)) + plasticState(ph)%state(:,me) = plasticState(ph)%partitionedState0(:,me) enddo end subroutine mech_restore diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index f2b61fb26..f1675f0a1 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -145,6 +145,70 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, end subroutine constitutive_thermal_getRateAndItsTangents + +module subroutine thermal_initializeRestorationPoints(ph,me) + + integer, intent(in) :: ph, me + + integer :: so + + + do so = 1, size(sourceState(ph)%p) + thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) + enddo + +end subroutine thermal_initializeRestorationPoints + + + +module subroutine thermal_windForward(ph,me) + + integer, intent(in) :: ph, me + + integer :: so + + + do so = 1, size(sourceState(ph)%p) + thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) + enddo + +end subroutine thermal_windForward + + +module subroutine thermal_forward() + + integer :: ph, so + + + do ph = 1, size(thermalState) + do so = 1, size(sourceState(ph)%p) + thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state + enddo + enddo + +end subroutine thermal_forward + + +module subroutine thermal_restore(ip,el) + + integer, intent(in) :: ip, el + + integer :: co, ph, me, so + + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + do so = 1, size(sourceState(ph)%p) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) + enddo + + enddo + +end subroutine thermal_restore + + !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- From 27f4e4ce2abb980028f38509e92127ca32196593 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Jan 2021 22:15:18 +0100 Subject: [PATCH 10/21] separate state for thermal --- src/constitutive.f90 | 146 ++------------- src/constitutive_mech.f90 | 10 + src/constitutive_thermal.f90 | 212 ++++++++++++++++++++-- src/constitutive_thermal_dissipation.f90 | 2 +- src/constitutive_thermal_externalheat.f90 | 24 +-- src/homogenization.f90 | 12 +- src/homogenization_thermal.f90 | 4 +- src/thermal_conduction.f90 | 43 +++-- src/thermal_isothermal.f90 | 26 ++- 9 files changed, 299 insertions(+), 180 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index f9537a136..6a023022f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -90,6 +90,7 @@ module constitutive phase_kinematics !< active kinematic mechanisms of each phase integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) + thermal_Nsources, & phase_Nsources, & !< number of source mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase @@ -233,6 +234,16 @@ module constitutive ! == cleaned:end =================================================================================== + module function integrateThermalState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + end function + module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: co, ip, el @@ -665,31 +676,6 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) end function constitutive_damage_collectDotState -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_thermal_collectDotState(ph,me) result(broken) - - integer, intent(in) :: ph, me - logical :: broken - - integer :: i - - - broken = .false. - - SourceLoop: do i = 1, phase_Nsources(ph) - - if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) & - call source_thermal_externalheat_dotState(ph,me) - - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me))) - - enddo SourceLoop - -end function constitutive_thermal_collectDotState - - !-------------------------------------------------------------------------------------------------- !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model @@ -856,7 +842,7 @@ subroutine crystallite_init() cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements - + class(tNode), pointer :: & num_crystallite, & @@ -914,6 +900,9 @@ subroutine crystallite_init() do so = 1, phase_Nsources(ph) allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack enddo + do so = 1, thermal_Nsources(ph) + allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack + enddo enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -1144,111 +1133,6 @@ function integrateSourceState(dt,co,ip,el) result(broken) end function integrateSourceState - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -function integrateThermalState(dt,co,ip,el) result(broken) - - real(pReal), intent(in) :: dt - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - - integer :: & - NiterationState, & !< number of iterations in state loop - ph, & - me, & - so - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken, converged_ - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - converged_ = .true. - broken = constitutive_thermal_collectDotState(ph,me) - if(broken) return - - do so = 1, phase_Nsources(ph) - size_so(so) = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) & - + thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo - - iteration: do NiterationState = 1, num%nState - - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me) - enddo - - broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) exit iteration - - do so = 1, phase_Nsources(ph) - zeta = damper(thermalState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) & - - thermalState(ph)%p(so)%subState0(1:size_so(so),me) & - - thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt - thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - thermalState(ph)%p(so)%state(1:size_so(so),me), & - thermalState(ph)%p(so)%atol(1:size_so(so))) - enddo - - if(converged_) then - broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) - exit iteration - endif - - enddo iteration - - broken = broken .or. .not. converged_ - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end function integrateThermalState - - !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index fccccf00c..9a065a829 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1588,6 +1588,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me) enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) + enddo subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) @@ -1616,6 +1619,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me) + enddo endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1632,6 +1638,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) + enddo todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif @@ -1645,6 +1654,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif enddo cutbackLooping diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index f1675f0a1..c86a286f9 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -6,9 +6,12 @@ submodule(constitutive) constitutive_thermal type :: tDataContainer real(pReal), dimension(:), allocatable :: T end type tDataContainer + integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & + thermal_source type(tDataContainer), dimension(:), allocatable :: current + integer :: thermal_source_maxSizeDotState interface module function source_thermal_dissipation_init(source_length) result(mySources) @@ -60,30 +63,55 @@ module subroutine thermal_init(phases) class(tNode), pointer :: & phases + class(tNode), pointer :: & + phase, thermal, sources + integer :: & - ph, & + ph, so, & Nconstituents - print'(/,a)', ' <<<+- constitutive_mech init -+>>>' + print'(/,a)', ' <<<+- constitutive_thermal init -+>>>' allocate(current(phases%length)) + allocate(thermalState (phases%length)) + allocate(thermal_Nsources(phases%length),source = 0) do ph = 1, phases%length Nconstituents = count(material_phaseAt == ph) * discretization_nIPs allocate(current(ph)%T(Nconstituents)) + phase => phases%get(ph) + if(phase%contains('thermal')) then + thermal => phase%get('thermal') + sources => thermal%get('source',defaultVal=emptyList) + thermal_Nsources(ph) = sources%length + endif + allocate(thermalstate(ph)%p(thermal_Nsources(ph))) enddo -! initialize source mechanisms - if(maxval(phase_Nsources) /= 0) then - where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID - where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID + allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = SOURCE_undefined_ID) + + if(maxval(thermal_Nsources) /= 0) then + where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_dissipation_ID + where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_externalheat_ID endif + thermal_source_maxSizeDotState = 0 + PhaseLoop2:do ph = 1,phases%length + + do so = 1,thermal_Nsources(ph) + thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0 + thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%partitionedState0 + enddo + + thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & + maxval(thermalState(ph)%p%sizeDotState)) + enddo PhaseLoop2 + !-------------------------------------------------------------------------------------------------- !initialize kinematic mechanisms if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & @@ -123,8 +151,8 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, do co = 1, homogenization_Nconstituents(homog) ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - do so = 1, phase_Nsources(ph) - select case(phase_source(so,ph)) + do so = 1, thermal_Nsources(ph) + select case(thermal_source(so,ph)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & mech_S(ph,me),mech_L_p(ph,me), ph) @@ -145,6 +173,131 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, end subroutine constitutive_thermal_getRateAndItsTangents +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_thermal_collectDotState(ph,me) result(broken) + + integer, intent(in) :: ph, me + logical :: broken + + integer :: i + + + broken = .false. + + SourceLoop: do i = 1, thermal_Nsources(ph) + + if (thermal_source(i,ph) == SOURCE_thermal_externalheat_ID) & + call source_thermal_externalheat_dotState(ph,me) + + broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me))) + + enddo SourceLoop + +end function constitutive_thermal_collectDotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +module function integrateThermalState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + so + integer, dimension(maxval(thermal_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(thermal_source_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(thermal_source_maxSizeDotState,2,maxval(thermal_Nsources)) :: source_dotState + logical :: & + broken, converged_ + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_thermal_collectDotState(ph,me) + if(broken) return + + do so = 1, thermal_Nsources(ph) + size_so(so) = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) & + + thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do so = 1, thermal_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me) + enddo + + broken = constitutive_thermal_collectDotState(ph,me) + if(broken) exit iteration + + do so = 1, thermal_Nsources(ph) + zeta = damper(thermalState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) & + - thermalState(ph)%p(so)%subState0(1:size_so(so),me) & + - thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt + thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + thermalState(ph)%p(so)%state(1:size_so(so),me), & + thermalState(ph)%p(so)%atol(1:size_so(so))) + enddo + + if(converged_) exit iteration + + enddo iteration + + broken = broken .or. .not. converged_ + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end function integrateThermalState + + module subroutine thermal_initializeRestorationPoints(ph,me) @@ -153,7 +306,7 @@ module subroutine thermal_initializeRestorationPoints(ph,me) integer :: so - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) enddo @@ -168,7 +321,7 @@ module subroutine thermal_windForward(ph,me) integer :: so - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) enddo @@ -181,7 +334,7 @@ module subroutine thermal_forward() do ph = 1, size(thermalState) - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state enddo enddo @@ -200,7 +353,7 @@ module subroutine thermal_restore(ip,el) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) enddo @@ -237,4 +390,39 @@ module subroutine constitutive_thermal_setT(T,co,ip,el) end subroutine constitutive_thermal_setT + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a source mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function thermal_active(source_label,src_length) result(active_source) + + character(len=*), intent(in) :: source_label !< name of source mechanism + integer, intent(in) :: src_length !< max. number of sources in system + logical, dimension(:,:), allocatable :: active_source + + class(tNode), pointer :: & + phases, & + phase, & + sources, thermal, & + src + integer :: p,s + + phases => config_material%get('phase') + allocate(active_source(src_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + if (phase%contains('thermal')) then + thermal => phase%get('thermal',defaultVal=emptyList) + sources => thermal%get('source',defaultVal=emptyList) + do s = 1, sources%length + src => sources%get(s) + if(src%get_asString('type') == source_label) active_source(s,p) = .true. + enddo + endif + enddo + + +end function thermal_active + + end submodule constitutive_thermal diff --git a/src/constitutive_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 index 27653a9ef..44227536c 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -62,7 +62,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources) src => sources%get(sourceOffset) prm%kappa = src%get_asFloat('kappa') Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0) + call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0) end associate endif diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 index 3ef96790e..de1617efa 100644 --- a/src/constitutive_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -13,7 +13,7 @@ submodule(constitutive:constitutive_thermal) source_externalheat type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & - t_n, & + t_n, & f_T integer :: & nIntervals @@ -31,19 +31,20 @@ contains !-------------------------------------------------------------------------------------------------- module function source_thermal_externalheat_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & - sources, & - src + sources, thermal, & + src integer :: Ninstances,sourceOffset,Nconstituents,p print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' - mySources = source_active('thermal_externalheat',source_length) + mySources = thermal_active('externalheat',source_length) + Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -54,15 +55,16 @@ module function source_thermal_externalheat_init(source_length) result(mySources allocate(source_thermal_externalheat_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) + phase => phases%get(p) if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p)) if(count(mySources(:,p)) == 0) cycle - sources => phase%get('source') + thermal => phase%get('thermal') + sources => thermal%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_thermal_externalheat_offset(p) = sourceOffset associate(prm => param(source_thermal_externalheat_instance(p))) - src => sources%get(sourceOffset) + src => sources%get(sourceOffset) prm%t_n = src%get_asFloats('t_n') prm%nIntervals = size(prm%t_n) - 1 @@ -70,7 +72,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) + call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0) end associate endif @@ -95,7 +97,7 @@ module subroutine source_thermal_externalheat_dotState(phase, of) sourceOffset = source_thermal_externalheat_offset(phase) - sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time + thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time end subroutine source_thermal_externalheat_dotState @@ -121,7 +123,7 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d associate(prm => param(source_thermal_externalheat_instance(phase))) do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & + frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & diff --git a/src/homogenization.f90 b/src/homogenization.f90 index df7369096..8738ba6f1 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -69,6 +69,13 @@ module homogenization el !< element number end subroutine mech_partition + module subroutine thermal_partition(T,ip,el) + real(pReal), intent(in) :: T + integer, intent(in) :: & + ip, & !< integration point + el !< element number + end subroutine thermal_partition + module subroutine mech_homogenize(dt,ip,el) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -131,9 +138,10 @@ subroutine homogenization_init call mech_init(num_homog) + call thermal_init() - if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init - if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init + if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T) + if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T) if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 59e7357b6..f181d97fa 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -11,9 +11,11 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine thermal_init() + print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' - allocate(homogenization_T(discretization_nIPs*discretization_Nelems), source=0.0_pReal) + allocate(homogenization_T(discretization_nIPs*discretization_Nelems)) + end subroutine thermal_init diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index f98d36d3b..0cd1678e0 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -10,6 +10,7 @@ module thermal_conduction use results use constitutive use YAML_types + use discretization implicit none private @@ -38,25 +39,28 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init +subroutine thermal_conduction_init(T) - integer :: Ninstances,Nmaterialpoints,h + real(pReal), dimension(:), intent(inout) :: T + + integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce class(tNode), pointer :: & material_homogenization, & homog, & homogThermal - + + print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) Ninstances = count(thermal_type == THERMAL_conduction_ID) allocate(param(Ninstances)) material_homogenization => config_material%get('homogenization') - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_conduction_ID) cycle - homog => material_homogenization%get(h) + do ho = 1, size(material_name_homogenization) + if (thermal_type(ho) /= THERMAL_conduction_ID) cycle + homog => material_homogenization%get(ho) homogThermal => homog%get('thermal') - associate(prm => param(thermal_typeInstance(h))) + associate(prm => param(thermal_typeInstance(ho))) #if defined (__GFORTRAN__) prm%output = output_asStrings(homogThermal) @@ -64,14 +68,23 @@ subroutine thermal_conduction_init prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) #endif - Nmaterialpoints=count(material_homogenizationAt==h) + Nmaterialpoints=count(material_homogenizationAt==ho) - allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) - allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) + allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho)) + allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal) end associate enddo + ce = 0 + do el = 1, discretization_Nelems + do ip = 1, discretization_nIPs + ce = ce + 1 + ho = material_homogenizationAt(el) + if (thermal_type(ho) == THERMAL_conduction_ID) T(ce) = thermal_initialT(ho) + enddo + enddo + end subroutine thermal_conduction_init @@ -89,12 +102,12 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) Tdot, dTdot_dT integer :: & homog - + Tdot = 0.0_pReal dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) @@ -112,13 +125,13 @@ function thermal_conduction_getConductivity(ip,el) el !< element number real(pReal), dimension(3,3) :: & thermal_conduction_getConductivity - + integer :: & co thermal_conduction_getConductivity = 0.0_pReal - + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) thermal_conduction_getConductivity = thermal_conduction_getConductivity + & crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) @@ -168,7 +181,7 @@ function thermal_conduction_getMassDensity(ip,el) el !< element number real(pReal) :: & thermal_conduction_getMassDensity - + integer :: & co diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 2a41ada49..09e35931e 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -6,6 +6,7 @@ module thermal_isothermal use prec use config use material + use discretization implicit none public @@ -15,22 +16,33 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine thermal_isothermal_init +subroutine thermal_isothermal_init(T) - integer :: h,Nmaterialpoints + real(pReal), dimension(:), intent(inout) :: T + + integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_isothermal_ID) cycle + do ho = 1, size(thermal_type) + if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle - Nmaterialpoints = count(material_homogenizationAt == h) + Nmaterialpoints = count(material_homogenizationAt == ho) - allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) - allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) + allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho)) + allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal) enddo + ce = 0 + do el = 1, discretization_Nelems + do ip = 1, discretization_nIPs + ce = ce + 1 + ho = material_homogenizationAt(el) + if (thermal_type(ho) == THERMAL_isothermal_ID) T(ce) = thermal_initialT(ho) + enddo + enddo + end subroutine thermal_isothermal_init end module thermal_isothermal From 1df409376c884e1138dc3b2eeb0b7f9eee76481b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Jan 2021 23:32:54 +0100 Subject: [PATCH 11/21] sourceState is now damage state --- src/constitutive.f90 | 283 ++---------------------------------- src/constitutive_damage.f90 | 193 ++++++++++++++++++++++++ src/constitutive_mech.f90 | 82 ++++++++++- 3 files changed, 286 insertions(+), 272 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6a023022f..43d9b6b3f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -235,14 +235,22 @@ module constitutive ! == cleaned:end =================================================================================== module function integrateThermalState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt integer, intent(in) :: & el, & !< element index in element loop ip, & !< integration point index in ip loop co !< grain index in grain loop logical :: broken - end function + end function integrateThermalState + + module function integrateDamageState(dt,co,ip,el) result(broken) + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + end function integrateDamageState module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt @@ -395,7 +403,6 @@ module constitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_LiAndItsTangents, & constitutive_damage_getRateAndItsTangents, & constitutive_thermal_getRateAndItsTangents, & constitutive_results, & @@ -413,7 +420,8 @@ module constitutive crystallite_push33ToRef, & constitutive_restartWrite, & constitutive_restartRead, & - integrateSourceState, & + integrateThermalState, & + integrateDamageState, & constitutive_thermal_setT, & constitutive_mech_getP, & constitutive_mech_setF, & @@ -555,173 +563,6 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki end function kinematics_active -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: MD: S is Mi? -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Li !< intermediate velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dS, & !< derivative of Li with respect to S - dLi_dFi - - real(pReal), dimension(3,3) :: & - my_Li, & !< intermediate velocity gradient - FiInv, & - temp_33 - real(pReal), dimension(3,3,3,3) :: & - my_dLi_dS - real(pReal) :: & - detFi - integer :: & - k, i, j, & - instance, of - - Li = 0.0_pReal - dLi_dS = 0.0_pReal - dLi_dFi = 0.0_pReal - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - case (PLASTICITY_isotropic_ID) plasticityType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) - case default plasticityType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select plasticityType - - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_thermal_expansion_ID) kinematicsType - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - enddo KinematicsLoop - - FiInv = math_inv33(Fi) - detFi = math_det33(Fi) - Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration - temp_33 = matmul(FiInv,Li) - - do i = 1,3; do j = 1,3 - dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo - -end subroutine constitutive_LiAndItsTangents - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - integer :: & - so !< counter in source loop - logical :: broken - - - broken = .false. - - SourceLoop: do so = 1, phase_Nsources(ph) - - sourceType: select case (phase_source(so,ph)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& - co, ip, el) ! correct stress? - - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState(co, ip, el) - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState(co, ip, el) - - end select sourceType - - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) - - enddo SourceLoop - -end function constitutive_damage_collectDotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief for constitutive models having an instantaneous change of state -!> will return false if delta state is not needed/supported by the constitutive model -!-------------------------------------------------------------------------------------------------- -function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - real(pReal), intent(in), dimension(3,3) :: & - Fe !< elastic deformation gradient - integer :: & - so, & - myOffset, & - mySize - logical :: & - broken - - - broken = .false. - - sourceLoop: do so = 1, phase_Nsources(ph) - - sourceType: select case (phase_source(so,ph)) - - case (SOURCE_damage_isoBrittle_ID) sourceType - call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & - co, ip, el) - broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) - if(.not. broken) then - myOffset = sourceState(ph)%p(so)%offsetDeltaState - mySize = sourceState(ph)%p(so)%sizeDeltaState - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) - endif - - end select sourceType - - enddo SourceLoop - -end function constitutive_damage_deltaState - - !-------------------------------------------------------------------------------------------------- !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- @@ -1030,107 +871,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33) end function crystallite_push33ToRef -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -function integrateSourceState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - - integer :: & - NiterationState, & !< number of iterations in state loop - ph, & - me, & - so - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken, converged_ - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - converged_ = .true. - broken = constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) return - - do so = 1, phase_Nsources(ph) - size_so(so) = sourceState(ph)%p(so)%sizeDotState - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo - - iteration: do NiterationState = 1, num%nState - - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) - enddo - - broken = constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) exit iteration - - do so = 1, phase_Nsources(ph) - zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & - - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - sourceState(ph)%p(so)%state(1:size_so(so),me), & - sourceState(ph)%p(so)%atol(1:size_so(so))) - enddo - - if(converged_) then - broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) - exit iteration - endif - - enddo iteration - - broken = broken .or. .not. converged_ - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end function integrateSourceState !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 3ce614666..be47d92b6 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -214,6 +214,111 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +module function integrateDamageState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + so + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + converged_ + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) return + + do so = 1, phase_Nsources(ph) + size_so(so) = sourceState(ph)%p(so)%sizeDotState + sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & + + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do so = 1, phase_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) + enddo + + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) exit iteration + + do so = 1, phase_Nsources(ph) + zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & + - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & + - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + sourceState(ph)%p(so)%state(1:size_so(so),me), & + sourceState(ph)%p(so)%atol(1:size_so(so))) + enddo + + if(converged_) then + broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) + exit iteration + endif + + enddo iteration + + broken = broken .or. .not. converged_ + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end function integrateDamageState + + !---------------------------------------------------------------------------------------------- !< @brief writes damage sources results to HDF5 output file !---------------------------------------------------------------------------------------------- @@ -250,4 +355,92 @@ module subroutine damage_results(group,ph) end subroutine damage_results +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + of + integer :: & + so !< counter in source loop + logical :: broken + + + broken = .false. + + SourceLoop: do so = 1, phase_Nsources(ph) + + sourceType: select case (phase_source(so,ph)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& + co, ip, el) ! correct stress? + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_dotState(co, ip, el) + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_dotState(co, ip, el) + + end select sourceType + + broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) + + enddo SourceLoop + +end function constitutive_damage_collectDotState + + + +!-------------------------------------------------------------------------------------------------- +!> @brief for constitutive models having an instantaneous change of state +!> will return false if delta state is not needed/supported by the constitutive model +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + of + real(pReal), intent(in), dimension(3,3) :: & + Fe !< elastic deformation gradient + integer :: & + so, & + myOffset, & + mySize + logical :: & + broken + + + broken = .false. + + sourceLoop: do so = 1, phase_Nsources(ph) + + sourceType: select case (phase_source(so,ph)) + + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & + co, ip, el) + broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) + if(.not. broken) then + myOffset = sourceState(ph)%p(so)%offsetDeltaState + mySize = sourceState(ph)%p(so)%sizeDeltaState + sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & + sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) + endif + + end select sourceType + + enddo SourceLoop + +end function constitutive_damage_deltaState + + end submodule constitutive_damage diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 9a065a829..7c403eeea 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1653,7 +1653,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif @@ -1938,5 +1938,85 @@ module subroutine constitutive_mech_setF(F,co,ip,el) end subroutine constitutive_mech_setF + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: MD: S is Mi? +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + real(pReal), dimension(3,3) :: & + my_Li, & !< intermediate velocity gradient + FiInv, & + temp_33 + real(pReal), dimension(3,3,3,3) :: & + my_dLi_dS + real(pReal) :: & + detFi + integer :: & + k, i, j, & + instance, of + + Li = 0.0_pReal + dLi_dS = 0.0_pReal + dLi_dFi = 0.0_pReal + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_isotropic_ID) plasticityType + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + case default plasticityType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select plasticityType + + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) + kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) + case (KINEMATICS_cleavage_opening_ID) kinematicsType + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_slipplane_opening_ID) kinematicsType + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_thermal_expansion_ID) kinematicsType + call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) + case default kinematicsType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select kinematicsType + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + enddo KinematicsLoop + + FiInv = math_inv33(Fi) + detFi = math_det33(Fi) + Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration + temp_33 = matmul(FiInv,Li) + + do i = 1,3; do j = 1,3 + dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + enddo; enddo + +end subroutine constitutive_LiAndItsTangents + end submodule constitutive_mech From 5efa6c997a2f31d163a892800fc18d065104ea99 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Jan 2021 23:50:06 +0100 Subject: [PATCH 12/21] meaningful scope --- src/constitutive.f90 | 37 +---------------------------- src/constitutive_damage.f90 | 46 ++++++++++++++++++++++-------------- src/constitutive_mech.f90 | 20 +++++++++++++--- src/constitutive_thermal.f90 | 20 ++++++++++------ 4 files changed, 59 insertions(+), 64 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 43d9b6b3f..5f96e80e6 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -21,21 +21,6 @@ module constitutive private enum, bind(c); enumerator :: & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - SOURCE_UNDEFINED_ID ,& - SOURCE_THERMAL_DISSIPATION_ID, & - SOURCE_THERMAL_EXTERNALHEAT_ID, & - SOURCE_DAMAGE_ISOBRITTLE_ID, & - SOURCE_DAMAGE_ISODUCTILE_ID, & - SOURCE_DAMAGE_ANISOBRITTLE_ID, & - SOURCE_DAMAGE_ANISODUCTILE_ID, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & KINEMATICS_SLIPPLANE_OPENING_ID, & @@ -81,12 +66,7 @@ module constitutive type(tDebugOptions) :: debugCrystallite - - integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public :: & - phase_plasticity !< plasticity of each phase - - integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & - phase_source, & !< active sources mechanisms of each phase + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & phase_kinematics !< active kinematic mechanisms of each phase integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) @@ -428,21 +408,6 @@ module constitutive constitutive_mech_getF, & constitutive_initializeRestorationPoints, & constitutive_windForward, & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - SOURCE_UNDEFINED_ID ,& - SOURCE_THERMAL_DISSIPATION_ID, & - SOURCE_THERMAL_EXTERNALHEAT_ID, & - SOURCE_DAMAGE_ISOBRITTLE_ID, & - SOURCE_DAMAGE_ISODUCTILE_ID, & - SOURCE_DAMAGE_ANISOBRITTLE_ID, & - SOURCE_DAMAGE_ANISODUCTILE_ID, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & KINEMATICS_SLIPPLANE_OPENING_ID, & diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index be47d92b6..ea00b5c94 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -2,6 +2,16 @@ !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_damage + enum, bind(c); enumerator :: & + DAMAGE_UNDEFINED_ID, & + DAMAGE_ISOBRITTLE_ID, & + DAMAGE_ISODUCTILE_ID, & + DAMAGE_ANISOBRITTLE_ID, & + DAMAGE_ANISODUCTILE_ID + end enum + + integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & + phase_source !< active sources mechanisms of each phase interface @@ -129,14 +139,14 @@ module subroutine damage_init allocate(sourceState(ph)%p(phase_Nsources(ph))) enddo - allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) + allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then - where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoBrittle_ID - where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoDuctile_ID - where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoBrittle_ID - where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoDuctile_ID + where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID + where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID + where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID + where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID endif !-------------------------------------------------------------------------------------------------- @@ -189,16 +199,16 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi constituent = material_phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) + case (DAMAGE_ISOBRITTLE_ID) call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_isoDuctile_ID) + case (DAMAGE_ISODUCTILE_ID) call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_anisoBrittle_ID) + case (DAMAGE_ANISOBRITTLE_ID) call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_anisoDuctile_ID) + case (DAMAGE_ANISODUCTILE_ID) call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) case default @@ -331,21 +341,21 @@ module subroutine damage_results(group,ph) sourceLoop: do so = 1, phase_Nsources(ph) - if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) & + if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage' sourceType: select case (phase_source(so,ph)) - case (SOURCE_damage_anisoBrittle_ID) sourceType + case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_anisoBrittle_results(ph,group//'sources/') - case (SOURCE_damage_anisoDuctile_ID) sourceType + case (DAMAGE_ISODUCTILE_ID) sourceType call source_damage_anisoDuctile_results(ph,group//'sources/') - case (SOURCE_damage_isoBrittle_ID) sourceType + case (DAMAGE_ANISOBRITTLE_ID) sourceType call source_damage_isoBrittle_results(ph,group//'sources/') - case (SOURCE_damage_isoDuctile_ID) sourceType + case (DAMAGE_ANISODUCTILE_ID) sourceType call source_damage_isoDuctile_results(ph,group//'sources/') end select sourceType @@ -377,14 +387,14 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) sourceType: select case (phase_source(so,ph)) - case (SOURCE_damage_anisoBrittle_ID) sourceType + case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& co, ip, el) ! correct stress? - case (SOURCE_damage_isoDuctile_ID) sourceType + case (DAMAGE_ISODUCTILE_ID) sourceType call source_damage_isoDuctile_dotState(co, ip, el) - case (SOURCE_damage_anisoDuctile_ID) sourceType + case (DAMAGE_ANISODUCTILE_ID) sourceType call source_damage_anisoDuctile_dotState(co, ip, el) end select sourceType @@ -425,7 +435,7 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) sourceType: select case (phase_source(so,ph)) - case (SOURCE_damage_isoBrittle_ID) sourceType + case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & co, ip, el) broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 7c403eeea..9539d0b93 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -7,12 +7,20 @@ submodule(constitutive) constitutive_mech ELASTICITY_UNDEFINED_ID, & ELASTICITY_HOOKE_ID, & STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID + STIFFNESS_DEGRADATION_DAMAGE_ID, & + PLASTICITY_UNDEFINED_ID, & + PLASTICITY_NONE_ID, & + PLASTICITY_ISOTROPIC_ID, & + PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_KINEHARDENING_ID, & + PLASTICITY_DISLOTWIN_ID, & + PLASTICITY_DISLOTUNGSTEN_ID, & + PLASTICITY_NONLOCAL_ID end enum - integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & + integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase - integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & + integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase type(tTensorContainer), dimension(:), allocatable :: & @@ -41,6 +49,12 @@ submodule(constitutive) constitutive_mech constitutive_mech_partitionedS0 + + + integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & + phase_plasticity !< plasticity of each phase + + interface module function plastic_none_init() result(myPlasticity) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index c86a286f9..636d7e447 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -3,10 +3,16 @@ !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_thermal + enum, bind(c); enumerator :: & + THERMAL_UNDEFINED_ID ,& + THERMAL_DISSIPATION_ID, & + THERMAL_EXTERNALHEAT_ID + end enum + type :: tDataContainer real(pReal), dimension(:), allocatable :: T end type tDataContainer - integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & + integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & thermal_source type(tDataContainer), dimension(:), allocatable :: current @@ -93,11 +99,11 @@ module subroutine thermal_init(phases) allocate(thermalstate(ph)%p(thermal_Nsources(ph))) enddo - allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = SOURCE_undefined_ID) + allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) if(maxval(thermal_Nsources) /= 0) then - where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_dissipation_ID - where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_externalheat_ID + where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID + where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID endif thermal_source_maxSizeDotState = 0 @@ -153,11 +159,11 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, me = material_phasememberAt(co,ip,el) do so = 1, thermal_Nsources(ph) select case(thermal_source(so,ph)) - case (SOURCE_thermal_dissipation_ID) + case (THERMAL_DISSIPATION_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & mech_S(ph,me),mech_L_p(ph,me), ph) - case (SOURCE_thermal_externalheat_ID) + case (THERMAL_EXTERNALHEAT_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & ph, me) @@ -188,7 +194,7 @@ function constitutive_thermal_collectDotState(ph,me) result(broken) SourceLoop: do i = 1, thermal_Nsources(ph) - if (thermal_source(i,ph) == SOURCE_thermal_externalheat_ID) & + if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & call source_thermal_externalheat_dotState(ph,me) broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me))) From 7239c0b226188ac102c54d1a5167ef1a8c6076a9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Jan 2021 00:40:21 +0100 Subject: [PATCH 13/21] explicit Euler is ok (only state is current time) --- src/constitutive.f90 | 4 +- src/constitutive_thermal.f90 | 91 ++++++------------------------------ 2 files changed, 15 insertions(+), 80 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 5f96e80e6..9e9cfe423 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -214,8 +214,8 @@ module constitutive ! == cleaned:end =================================================================================== - module function integrateThermalState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt + module function integrateThermalState(Delta_t,co,ip,el) result(broken) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & el, & !< element index in element loop ip, & !< integration point index in ip loop diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 636d7e447..e6aa7c4cf 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -204,107 +204,42 @@ function constitutive_thermal_collectDotState(ph,me) result(broken) end function constitutive_thermal_collectDotState -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -module function integrateThermalState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt +!-------------------------------------------------------------------------------------------------- +!> @brief integrate state with 1st order explicit Euler method +!-------------------------------------------------------------------------------------------------- +function integrateThermalState(Delta_t,co,ip,el) result(broken) + + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & el, & !< element index in element loop ip, & !< integration point index in ip loop co !< grain index in grain loop + logical :: & + broken integer :: & - NiterationState, & !< number of iterations in state loop ph, & me, & - so - integer, dimension(maxval(thermal_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(thermal_source_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(thermal_source_maxSizeDotState,2,maxval(thermal_Nsources)) :: source_dotState - logical :: & - broken, converged_ + so, & + sizeDotState ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - converged_ = .true. broken = constitutive_thermal_collectDotState(ph,me) if(broken) return do so = 1, thermal_Nsources(ph) - size_so(so) = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) & - + thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal + sizeDotState = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) & + + thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t enddo - iteration: do NiterationState = 1, num%nState - - do so = 1, thermal_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me) - enddo - - broken = constitutive_thermal_collectDotState(ph,me) - if(broken) exit iteration - - do so = 1, thermal_Nsources(ph) - zeta = damper(thermalState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) & - - thermalState(ph)%p(so)%subState0(1:size_so(so),me) & - - thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt - thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - thermalState(ph)%p(so)%state(1:size_so(so),me), & - thermalState(ph)%p(so)%atol(1:size_so(so))) - enddo - - if(converged_) exit iteration - - enddo iteration - - broken = broken .or. .not. converged_ - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - end function integrateThermalState - module subroutine thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me From 88be08ae31f49fdd007843101c4491c04f96b69a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Jan 2021 00:44:16 +0100 Subject: [PATCH 14/21] modified structure for thermal tests, fixed damage branching --- PRIVATE | 2 +- src/constitutive_damage.f90 | 22 +++++++++++----------- src/constitutive_thermal.f90 | 2 +- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/PRIVATE b/PRIVATE index 591964dcf..7846c7112 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 591964dcf8521d95f6cccbfe840d462c430e63d9 +Subproject commit 7846c71126705cc5d41dd79f2d595f4864434068 diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index ea00b5c94..cc2b62002 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -347,17 +347,17 @@ module subroutine damage_results(group,ph) sourceType: select case (phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_anisoBrittle_results(ph,group//'sources/') - - case (DAMAGE_ISODUCTILE_ID) sourceType - call source_damage_anisoDuctile_results(ph,group//'sources/') - - case (DAMAGE_ANISOBRITTLE_ID) sourceType call source_damage_isoBrittle_results(ph,group//'sources/') - case (DAMAGE_ANISODUCTILE_ID) sourceType + case (DAMAGE_ISODUCTILE_ID) sourceType call source_damage_isoDuctile_results(ph,group//'sources/') + case (DAMAGE_ANISOBRITTLE_ID) sourceType + call source_damage_anisoBrittle_results(ph,group//'sources/') + + case (DAMAGE_ANISODUCTILE_ID) sourceType + call source_damage_anisoDuctile_results(ph,group//'sources/') + end select sourceType enddo SourceLoop @@ -387,16 +387,16 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) sourceType: select case (phase_source(so,ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& - co, ip, el) ! correct stress? - case (DAMAGE_ISODUCTILE_ID) sourceType call source_damage_isoDuctile_dotState(co, ip, el) case (DAMAGE_ANISODUCTILE_ID) sourceType call source_damage_anisoDuctile_dotState(co, ip, el) + case (DAMAGE_ANISOBRITTLE_ID) sourceType + call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& + co, ip, el) ! correct stress? + end select sourceType broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index e6aa7c4cf..5017904df 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -208,7 +208,7 @@ end function constitutive_thermal_collectDotState !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateThermalState(Delta_t,co,ip,el) result(broken) +module function integrateThermalState(Delta_t,co,ip,el) result(broken) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & From 65bd880fdf79394e2693a49e3c681c34574fc9ba Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Jan 2021 07:10:38 +0100 Subject: [PATCH 15/21] clearerer names --- src/constitutive.f90 | 24 +++++++++--------- src/constitutive_damage.f90 | 40 +++++++++++++++--------------- src/constitutive_mech.f90 | 6 ++--- src/damage_none.f90 | 8 +++--- src/damage_nonlocal.f90 | 10 ++++---- src/homogenization.f90 | 8 +++--- src/material.f90 | 4 +-- src/source_damage_anisoBrittle.f90 | 16 ++++++------ src/source_damage_anisoDuctile.f90 | 12 ++++----- src/source_damage_isoBrittle.f90 | 24 +++++++++--------- src/source_damage_isoDuctile.f90 | 12 ++++----- 11 files changed, 82 insertions(+), 82 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 9e9cfe423..0d8e35ba3 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -83,7 +83,7 @@ module constitutive type(tPlasticState), allocatable, dimension(:), public :: & plasticState type(tSourceState), allocatable, dimension(:), public :: & - sourceState, thermalState + damageState, thermalState integer, public, protected :: & @@ -454,12 +454,12 @@ subroutine constitutive_init plasticState(ph)%partitionedState0 = plasticState(ph)%state0 plasticState(ph)%state = plasticState(ph)%partitionedState0 forall(so = 1:phase_Nsources(ph)) - sourceState(ph)%p(so)%partitionedState0 = sourceState(ph)%p(so)%state0 - sourceState(ph)%p(so)%state = sourceState(ph)%p(so)%partitionedState0 + damageState(ph)%p(so)%partitionedState0 = damageState(ph)%p(so)%state0 + damageState(ph)%p(so)%state = damageState(ph)%p(so)%partitionedState0 end forall constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & - maxval(sourceState(ph)%p%sizeDotState)) + maxval(damageState(ph)%p%sizeDotState)) enddo PhaseLoop2 constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) @@ -578,8 +578,8 @@ subroutine constitutive_restore(ip,el,includeL) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do so = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = & - sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) + damageState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = & + damageState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) enddo enddo @@ -601,9 +601,9 @@ subroutine constitutive_forward() call mech_forward() call thermal_forward() - do ph = 1, size(sourceState) + do ph = 1, size(damageState) do so = 1,phase_Nsources(ph) - sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state + damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state enddo; enddo end subroutine constitutive_forward @@ -704,7 +704,7 @@ subroutine crystallite_init() do ph = 1, phases%length do so = 1, phase_Nsources(ph) - allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack + allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack enddo do so = 1, thermal_Nsources(ph) allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack @@ -753,8 +753,8 @@ subroutine constitutive_initializeRestorationPoints(ip,el) call mech_initializeRestorationPoints(ph,me) call thermal_initializeRestorationPoints(ph,me) - do so = 1, size(sourceState(ph)%p) - sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me) + do so = 1, size(damageState(ph)%p) + damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) enddo enddo @@ -784,7 +784,7 @@ subroutine constitutive_windForward(ip,el) call thermal_windForward(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state(:,me) + damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me) enddo enddo diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index cc2b62002..85500e260 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -129,14 +129,14 @@ module subroutine damage_init phases => config_material%get('phase') - allocate(sourceState (phases%length)) + allocate(damageState (phases%length)) allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics do ph = 1,phases%length phase => phases%get(ph) sources => phase%get('source',defaultVal=emptyList) phase_Nsources(ph) = sources%length - allocate(sourceState(ph)%p(phase_Nsources(ph))) + allocate(damageState(ph)%p(phase_Nsources(ph))) enddo allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) @@ -262,9 +262,9 @@ module function integrateDamageState(dt,co,ip,el) result(broken) if(broken) return do so = 1, phase_Nsources(ph) - size_so(so) = sourceState(ph)%p(so)%sizeDotState - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + size_so(so) = damageState(ph)%p(so)%sizeDotState + damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) & + + damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt source_dotState(1:size_so(so),2,so) = 0.0_pReal enddo @@ -272,26 +272,26 @@ module function integrateDamageState(dt,co,ip,el) result(broken) do so = 1, phase_Nsources(ph) if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) + source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) enddo broken = constitutive_damage_collectDotState(co,ip,el,ph,me) if(broken) exit iteration do so = 1, phase_Nsources(ph) - zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & + zeta = damper(damageState(ph)%p(so)%dotState(:,me), & source_dotState(1:size_so(so),1,so),& source_dotState(1:size_so(so),2,so)) - sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & + damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta & + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & - - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & + r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) & + - damageState(ph)%p(so)%subState0(1:size_so(so),me) & + - damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt + damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) & - r(1:size_so(so)) converged_ = converged_ .and. converged(r(1:size_so(so)), & - sourceState(ph)%p(so)%state(1:size_so(so),me), & - sourceState(ph)%p(so)%atol(1:size_so(so))) + damageState(ph)%p(so)%state(1:size_so(so),me), & + damageState(ph)%p(so)%atol(1:size_so(so))) enddo if(converged_) then @@ -399,7 +399,7 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) end select sourceType - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,of))) enddo SourceLoop @@ -438,12 +438,12 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & co, ip, el) - broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) + broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of))) if(.not. broken) then - myOffset = sourceState(ph)%p(so)%offsetDeltaState - mySize = sourceState(ph)%p(so)%sizeDeltaState - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) + myOffset = damageState(ph)%p(so)%offsetDeltaState + mySize = damageState(ph)%p(so)%sizeDeltaState + damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & + damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + damageState(ph)%p(so)%deltaState(1:mySize,of) endif end select sourceType diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 9539d0b93..8bc85354f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1600,7 +1600,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) - sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me) + damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me) enddo do so = 1, thermal_Nsources(ph) thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) @@ -1631,7 +1631,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) do so = 1, phase_Nsources(ph) - sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) + damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) enddo do so = 1, thermal_Nsources(ph) thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me) @@ -1650,7 +1650,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) endif plasticState(ph)%state(:,me) = subState0 do so = 1, phase_Nsources(ph) - sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) + damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) enddo do so = 1, thermal_Nsources(ph) thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 3f1144833..078d42af7 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -25,10 +25,10 @@ subroutine damage_none_init if (damage_type(h) /= DAMAGE_NONE_ID) cycle Nmaterialpoints = count(material_homogenizationAt == h) - damageState(h)%sizeState = 0 - allocate(damageState(h)%state0 (0,Nmaterialpoints)) - allocate(damageState(h)%subState0(0,Nmaterialpoints)) - allocate(damageState(h)%state (0,Nmaterialpoints)) + damageState_h(h)%sizeState = 0 + allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) + allocate(damageState_h(h)%subState0(0,Nmaterialpoints)) + allocate(damageState_h(h)%state (0,Nmaterialpoints)) allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 3db63cab2..4423c1e3a 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -76,12 +76,12 @@ subroutine damage_nonlocal_init #endif Nmaterialpoints = count(material_homogenizationAt == h) - damageState(h)%sizeState = 1 - allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) + damageState_h(h)%sizeState = 1 + allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) - damage(h)%p => damageState(h)%state(1,:) + damage(h)%p => damageState_h(h)%state(1,:) end associate enddo diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 8738ba6f1..9112562b9 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -186,7 +186,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) - if (damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me) + if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me) cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) @@ -200,7 +200,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call constitutive_windForward(ip,el) if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) - if(damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me) endif steppingNeeded elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite @@ -215,7 +215,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call constitutive_restore(ip,el,subStep < 1.0_pReal) if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) - if(damageState(ho)%sizeState > 0) damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me) endif if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] @@ -326,7 +326,7 @@ subroutine homogenization_forward do ho = 1, size(material_name_homogenization) homogState (ho)%state0 = homogState (ho)%state - damageState(ho)%state0 = damageState(ho)%state + damageState_h(ho)%state0 = damageState_h(ho)%state enddo end subroutine homogenization_forward diff --git a/src/material.f90 b/src/material.f90 index 581182d22..16116ca91 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -61,7 +61,7 @@ module material type(tState), allocatable, dimension(:), public :: & homogState, & - damageState + damageState_h type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element @@ -101,7 +101,7 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) - allocate(damageState (size(material_name_homogenization))) + allocate(damageState_h (size(material_name_homogenization))) allocate(temperature (size(material_name_homogenization))) allocate(damage (size(material_name_homogenization))) diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 0f923ceba..7c00c6580 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -101,9 +101,9 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' end associate @@ -146,7 +146,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_anisoBrittle_instance(phase))) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal + damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal do i = 1, prm%sum_N_cl traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) @@ -154,8 +154,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & - = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + damageState(phase)%p(sourceOffset)%dotState(1,constituent) & + = damageState(phase)%p(sourceOffset)%dotState(1,constituent) & + prm%dot_o / prm%s_crit(i) & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & @@ -185,7 +185,7 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d sourceOffset = source_damage_anisoBrittle_offset(phase) - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -204,7 +204,7 @@ module subroutine source_damage_anisoBrittle_results(phase,group) integer :: o associate(prm => param(source_damage_anisoBrittle_instance(phase)), & - stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) + stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 6f71fc145..7ec06cb62 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -87,9 +87,9 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -128,7 +128,7 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el) damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_anisoDuctile_instance(phase))) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + damageState(phase)%p(sourceOffset)%dotState(1,constituent) & = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) end associate @@ -154,7 +154,7 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d sourceOffset = source_damage_anisoDuctile_offset(phase) - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -173,7 +173,7 @@ module subroutine source_damage_anisoDuctile_results(phase,group) integer :: o associate(prm => param(source_damage_anisoDuctile_instance(phase)), & - stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) + stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 8c768b08d..1721b0201 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -74,9 +74,9 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,1) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -124,13 +124,13 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit - if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent) + if (strainenergy > damageState(phase)%p(sourceOffset)%subState0(1,constituent)) then + damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & + strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent) else - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - & - sourceState(phase)%p(sourceOffset)%state(1,constituent) + damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & + damageState(phase)%p(sourceOffset)%subState0(1,constituent) - & + damageState(phase)%p(sourceOffset)%state(1,constituent) endif end associate @@ -158,8 +158,8 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo associate(prm => param(source_damage_isoBrittle_instance(phase))) localphiDot = 1.0_pReal & - - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - sourceState(phase)%p(sourceOffset)%state(1,constituent) + - phi*damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) end associate end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -176,7 +176,7 @@ module subroutine source_damage_isoBrittle_results(phase,group) integer :: o associate(prm => param(source_damage_isoBrittle_instance(phase)), & - stt => sourceState(phase)%p(source_damage_isoBrittle_offset(phase))%state) + stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 86222bbf9..dd2910182 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -78,9 +78,9 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -119,7 +119,7 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el) damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_isoDuctile_instance(phase))) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & + damageState(phase)%p(sourceOffset)%dotState(1,constituent) = & sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit end associate @@ -145,7 +145,7 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo sourceOffset = source_damage_isoDuctile_offset(phase) - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -164,7 +164,7 @@ module subroutine source_damage_isoDuctile_results(phase,group) integer :: o associate(prm => param(source_damage_isoDuctile_instance(phase)), & - stt => sourceState(phase)%p(source_damage_isoDuctile_offset(phase))%state) + stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') From 6c62e186deccb414cec8681b298816e35e33df37 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Jan 2021 07:37:51 +0100 Subject: [PATCH 16/21] separte functionality --- src/constitutive.f90 | 31 +---------------------- src/constitutive_damage.f90 | 31 +++++++++++++++++++++++ src/constitutive_thermal_dissipation.f90 | 20 ++++++++------- src/constitutive_thermal_externalheat.f90 | 3 +-- 4 files changed, 44 insertions(+), 41 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 0d8e35ba3..111e68fdf 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -369,7 +369,7 @@ module constitutive module subroutine constitutive_plastic_dependentState(co,ip,el) integer, intent(in) :: & - co, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element end subroutine constitutive_plastic_dependentState @@ -390,7 +390,6 @@ module constitutive constitutive_forward, & constitutive_restore, & plastic_nonlocal_updateCompatibility, & - source_active, & kinematics_active, & converged, & crystallite_init, & @@ -466,35 +465,7 @@ subroutine constitutive_init end subroutine constitutive_init -!-------------------------------------------------------------------------------------------------- -!> @brief checks if a source mechanism is active or not -!-------------------------------------------------------------------------------------------------- -function source_active(source_label,src_length) result(active_source) - character(len=*), intent(in) :: source_label !< name of source mechanism - integer, intent(in) :: src_length !< max. number of sources in system - logical, dimension(:,:), allocatable :: active_source - - class(tNode), pointer :: & - phases, & - phase, & - sources, & - src - integer :: p,s - - phases => config_material%get('phase') - allocate(active_source(src_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) - sources => phase%get('source',defaultVal=emptyList) - do s = 1, sources%length - src => sources%get(s) - if(src%get_asString('type') == source_label) active_source(s,p) = .true. - enddo - enddo - - -end function source_active !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 85500e260..8c9104946 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -453,4 +453,35 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) end function constitutive_damage_deltaState +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a source mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function source_active(source_label,src_length) result(active_source) + + character(len=*), intent(in) :: source_label !< name of source mechanism + integer, intent(in) :: src_length !< max. number of sources in system + logical, dimension(:,:), allocatable :: active_source + + class(tNode), pointer :: & + phases, & + phase, & + sources, & + src + integer :: p,s + + phases => config_material%get('phase') + allocate(active_source(src_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + sources => phase%get('source',defaultVal=emptyList) + do s = 1, sources%length + src => sources%get(s) + if(src%get_asString('type') == source_label) active_source(s,p) = .true. + enddo + enddo + + +end function source_active + + end submodule constitutive_damage diff --git a/src/constitutive_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 index 44227536c..f15d1cfe9 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -27,19 +27,20 @@ contains !-------------------------------------------------------------------------------------------------- module function source_thermal_dissipation_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & - sources, & - src + sources, thermal, & + src integer :: Ninstances,sourceOffset,Nconstituents,p - print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>' + print'(/,a)', ' <<<+- thermal_externalheat init -+>>>' + + mySources = thermal_active('dissipation',source_length) - mySources = source_active('thermal_dissipation',source_length) Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -50,16 +51,17 @@ module function source_thermal_dissipation_init(source_length) result(mySources) allocate(source_thermal_dissipation_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle + phase => phases%get(p) if(any(mySources(:,p))) source_thermal_dissipation_instance(p) = count(mySources(:,1:p)) - sources => phase%get('source') + if(count(mySources(:,p)) == 0) cycle + thermal => phase%get('thermal') + sources => thermal%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_thermal_dissipation_offset(p) = sourceOffset associate(prm => param(source_thermal_dissipation_instance(p))) + src => sources%get(sourceOffset) - src => sources%get(sourceOffset) prm%kappa = src%get_asFloat('kappa') Nconstituents = count(material_phaseAt==p) * discretization_nIPs call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0) diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 index de1617efa..2a3ec7362 100644 --- a/src/constitutive_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -41,7 +41,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources src integer :: Ninstances,sourceOffset,Nconstituents,p - print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' + print'(/,a)', ' <<<+- thermal_externalheat init -+>>>' mySources = thermal_active('externalheat',source_length) @@ -74,7 +74,6 @@ module function source_thermal_externalheat_init(source_length) result(mySources Nconstituents = count(material_phaseAt==p) * discretization_nIPs call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0) end associate - endif enddo enddo From d494c2d81e66ceb770d7e902f262ba337c3fd6d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Jan 2021 08:26:17 +0100 Subject: [PATCH 17/21] better to read --- src/constitutive.f90 | 3 +-- src/constitutive_thermal.f90 | 2 +- src/thermal_conduction.f90 | 5 +---- 3 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 111e68fdf..250f7f99e 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -305,6 +305,7 @@ module constitutive orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -466,8 +467,6 @@ end subroutine constitutive_init - - !-------------------------------------------------------------------------------------------------- !> @brief checks if a kinematic mechanism is active or not !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 5017904df..a716a0c55 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -136,7 +136,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, el !< element number real(pReal), intent(in) :: & T !< plastic velocity gradient - real(pReal), intent(inout) :: & + real(pReal), intent(out) :: & TDot, & dTDot_dT diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 0cd1678e0..02649b1ad 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -103,10 +103,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) integer :: & homog - Tdot = 0.0_pReal - dTdot_dT = 0.0_pReal - - homog = material_homogenizationAt(el) + homog = material_homogenizationAt(el) call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) From 350466dd5f18aad38d11bb1cf49dd71591524d6f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Jan 2021 08:57:30 +0100 Subject: [PATCH 18/21] not needed --- src/DAMASK_marc.f90 | 3 ++- src/constitutive.f90 | 7 +++--- src/constitutive_thermal.f90 | 27 +++++++++-------------- src/constitutive_thermal_dissipation.f90 | 6 ++--- src/constitutive_thermal_externalheat.f90 | 6 ++--- src/grid/grid_thermal_spectral.f90 | 4 ++-- src/thermal_conduction.f90 | 10 ++++----- 7 files changed, 26 insertions(+), 37 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 0ad68445c..7c002e63c 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -364,7 +364,8 @@ subroutine flux(f,ts,n,time) real(pReal), dimension(2), intent(out) :: & f - call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) + f(2) = 0.0_pReal + call thermal_conduction_getSourceAndItsTangent(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) end subroutine flux diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 250f7f99e..0b58e524f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -283,15 +283,14 @@ module constitutive dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents - module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + module subroutine constitutive_thermal_getRateAndItsTangents(TDot, T,ip,el) integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), intent(in) :: & T - real(pReal), intent(inout) :: & - TDot, & - dTDot_dT + real(pReal), intent(out) :: & + TDot end subroutine constitutive_thermal_getRateAndItsTangents diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index a716a0c55..721d6925c 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -36,7 +36,7 @@ submodule(constitutive) constitutive_thermal end function kinematics_thermal_expansion_init - module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) + module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, Tstar,Lp,phase) integer, intent(in) :: & phase !< phase ID of element real(pReal), intent(in), dimension(3,3) :: & @@ -44,17 +44,15 @@ submodule(constitutive) constitutive_thermal real(pReal), intent(in), dimension(3,3) :: & Lp !< plastic velocuty gradient for a given element real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot end subroutine source_thermal_dissipation_getRateAndItsTangent - module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) + module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, phase,of) integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot end subroutine source_thermal_externalheat_getRateAndItsTangent end interface @@ -129,7 +127,7 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) +module subroutine constitutive_thermal_getRateAndItsTangents(TDot, T, ip, el) integer, intent(in) :: & ip, & !< integration point number @@ -137,12 +135,10 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, real(pReal), intent(in) :: & T !< plastic velocity gradient real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot real(pReal) :: & - my_Tdot, & - my_dTdot_dT + my_Tdot integer :: & ph, & homog, & @@ -154,25 +150,22 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, homog = material_homogenizationAt(el) instance = thermal_typeInstance(homog) + TDot = 0.0_pReal do co = 1, homogenization_Nconstituents(homog) ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) do so = 1, thermal_Nsources(ph) select case(thermal_source(so,ph)) case (THERMAL_DISSIPATION_ID) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - mech_S(ph,me),mech_L_p(ph,me), ph) + call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph) case (THERMAL_EXTERNALHEAT_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - ph, me) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, ph,me) case default my_Tdot = 0.0_pReal - my_dTdot_dT = 0.0_pReal end select Tdot = Tdot + my_Tdot - dTdot_dT = dTdot_dT + my_dTdot_dT enddo enddo diff --git a/src/constitutive_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 index f15d1cfe9..88b170f27 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -78,7 +78,7 @@ end function source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) +module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -88,12 +88,10 @@ module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT Lp real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot associate(prm => param(source_thermal_dissipation_instance(phase))) TDot = prm%kappa*sum(abs(Tstar*Lp)) - dTDot_dT = 0.0_pReal end associate end subroutine source_thermal_dissipation_getRateAndItsTangent diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 index 2a3ec7362..853f1e7dd 100644 --- a/src/constitutive_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -104,14 +104,13 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) +module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, phase, of) integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot integer :: & sourceOffset, interval @@ -131,7 +130,6 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside of bounds enddo - dTDot_dT = 0.0 end associate end subroutine source_thermal_externalheat_getRateAndItsTangent diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 259b45f33..aa3c38e4c 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -256,7 +256,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr integer :: i, j, k, cell - real(pReal) :: Tdot, dTdot_dT + real(pReal) :: Tdot T_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -278,7 +278,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T_current(i,j,k), 1, cell) + call thermal_conduction_getSourceAndItsTangent(Tdot, T_current(i,j,k), 1, cell) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & + thermal_conduction_getMassDensity (1,cell)* & thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - & diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 02649b1ad..4d6869d04 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -91,7 +91,7 @@ end subroutine thermal_conduction_init !-------------------------------------------------------------------------------------------------- !> @brief return heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) +subroutine thermal_conduction_getSourceAndItsTangent(Tdot, T,ip,el) integer, intent(in) :: & ip, & !< integration point number @@ -99,15 +99,15 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) real(pReal), intent(in) :: & T real(pReal), intent(out) :: & - Tdot, dTdot_dT - integer :: & + Tdot + + integer :: & homog homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + call constitutive_thermal_getRateAndItsTangents(TDot, T,ip,el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) - dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) end subroutine thermal_conduction_getSourceAndItsTangent From fc7f919c231bb0b6e0973efcdecc8836ab496259 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 8 Jan 2021 21:37:15 +0100 Subject: [PATCH 19/21] [skip ci] updated version information after successful test of v3.0.0-alpha2-230-g0fc670d01 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 62c706093..f123674b3 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-173-g584c7cc3a +v3.0.0-alpha2-230-g0fc670d01 From 209d59534aa865ab474d3544238e5d3966538ad0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 9 Jan 2021 17:19:48 +0100 Subject: [PATCH 20/21] copy and paste error --- src/constitutive_thermal_dissipation.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 index 88b170f27..9153915ca 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -37,7 +37,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources) src integer :: Ninstances,sourceOffset,Nconstituents,p - print'(/,a)', ' <<<+- thermal_externalheat init -+>>>' + print'(/,a)', ' <<<+- thermal_dissipation init -+>>>' mySources = thermal_active('dissipation',source_length) From b5bfb1dba9c88cb26c5d6040056d798d51556437 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 11 Jan 2021 16:13:59 +0100 Subject: [PATCH 21/21] tangent is not included anymore --- src/DAMASK_marc.f90 | 2 +- src/constitutive.f90 | 6 +++--- src/constitutive_thermal.f90 | 16 ++++++++-------- src/constitutive_thermal_dissipation.f90 | 4 ++-- src/constitutive_thermal_externalheat.f90 | 4 ++-- src/grid/grid_thermal_spectral.f90 | 2 +- src/thermal_conduction.f90 | 8 ++++---- 7 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 7c002e63c..0f9d37ddb 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -365,7 +365,7 @@ subroutine flux(f,ts,n,time) f f(2) = 0.0_pReal - call thermal_conduction_getSourceAndItsTangent(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) + call thermal_conduction_getSource(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) end subroutine flux diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 0b58e524f..c6415a883 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -283,7 +283,7 @@ module constitutive dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents - module subroutine constitutive_thermal_getRateAndItsTangents(TDot, T,ip,el) + module subroutine constitutive_thermal_getRate(TDot, T,ip,el) integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -291,7 +291,7 @@ module constitutive T real(pReal), intent(out) :: & TDot - end subroutine constitutive_thermal_getRateAndItsTangents + end subroutine constitutive_thermal_getRate @@ -384,7 +384,7 @@ module constitutive constitutive_init, & constitutive_homogenizedC, & constitutive_damage_getRateAndItsTangents, & - constitutive_thermal_getRateAndItsTangents, & + constitutive_thermal_getRate, & constitutive_results, & constitutive_allocateState, & constitutive_forward, & diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 721d6925c..9787cb0e4 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -36,7 +36,7 @@ submodule(constitutive) constitutive_thermal end function kinematics_thermal_expansion_init - module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, Tstar,Lp,phase) + module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase) integer, intent(in) :: & phase !< phase ID of element real(pReal), intent(in), dimension(3,3) :: & @@ -45,15 +45,15 @@ submodule(constitutive) constitutive_thermal Lp !< plastic velocuty gradient for a given element real(pReal), intent(out) :: & TDot - end subroutine source_thermal_dissipation_getRateAndItsTangent + end subroutine thermal_dissipation_getRate - module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, phase,of) + module subroutine thermal_externalheat_getRate(TDot, phase,of) integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & TDot - end subroutine source_thermal_externalheat_getRateAndItsTangent + end subroutine thermal_externalheat_getRate end interface @@ -127,7 +127,7 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRateAndItsTangents(TDot, T, ip, el) +module subroutine constitutive_thermal_getRate(TDot, T, ip, el) integer, intent(in) :: & ip, & !< integration point number @@ -157,10 +157,10 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, T, ip, el) do so = 1, thermal_Nsources(ph) select case(thermal_source(so,ph)) case (THERMAL_DISSIPATION_ID) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph) + call thermal_dissipation_getRate(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph) case (THERMAL_EXTERNALHEAT_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, ph,me) + call thermal_externalheat_getRate(my_Tdot, ph,me) case default my_Tdot = 0.0_pReal @@ -169,7 +169,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, T, ip, el) enddo enddo -end subroutine constitutive_thermal_getRateAndItsTangents +end subroutine constitutive_thermal_getRate !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 index 9153915ca..ae2d5735e 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -78,7 +78,7 @@ end function source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, Tstar, Lp, phase) +module subroutine thermal_dissipation_getRate(TDot, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -94,6 +94,6 @@ module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, Tstar, L TDot = prm%kappa*sum(abs(Tstar*Lp)) end associate -end subroutine source_thermal_dissipation_getRateAndItsTangent +end subroutine thermal_dissipation_getRate end submodule source_dissipation diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 index 853f1e7dd..2e8c02f8c 100644 --- a/src/constitutive_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -104,7 +104,7 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, phase, of) +module subroutine thermal_externalheat_getRate(TDot, phase, of) integer, intent(in) :: & phase, & @@ -132,6 +132,6 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, phase, enddo end associate -end subroutine source_thermal_externalheat_getRateAndItsTangent +end subroutine thermal_externalheat_getRate end submodule source_externalheat diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index aa3c38e4c..9d804ec56 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -278,7 +278,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - call thermal_conduction_getSourceAndItsTangent(Tdot, T_current(i,j,k), 1, cell) + call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & + thermal_conduction_getMassDensity (1,cell)* & thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - & diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 4d6869d04..79fe0d6cd 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -25,7 +25,7 @@ module thermal_conduction public :: & thermal_conduction_init, & - thermal_conduction_getSourceAndItsTangent, & + thermal_conduction_getSource, & thermal_conduction_getConductivity, & thermal_conduction_getSpecificHeat, & thermal_conduction_getMassDensity, & @@ -91,7 +91,7 @@ end subroutine thermal_conduction_init !-------------------------------------------------------------------------------------------------- !> @brief return heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_getSourceAndItsTangent(Tdot, T,ip,el) +subroutine thermal_conduction_getSource(Tdot, T,ip,el) integer, intent(in) :: & ip, & !< integration point number @@ -105,11 +105,11 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, T,ip,el) homog homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, T,ip,el) + call constitutive_thermal_getRate(TDot, T,ip,el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) -end subroutine thermal_conduction_getSourceAndItsTangent +end subroutine thermal_conduction_getSource !--------------------------------------------------------------------------------------------------