From 494ed244a08c510e315856afd00078c5df68487e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Jan 2021 10:14:34 +0100 Subject: [PATCH 1/8] consistent access via (ph)ase and phase(me)mber --- src/constitutive.f90 | 4 ++-- src/constitutive_damage.f90 | 3 ++- src/constitutive_mech.f90 | 31 ++++++++++---------------- src/constitutive_plastic_dislotwin.f90 | 27 ++++++++++------------ src/homogenization_mech_RGC.f90 | 5 ++++- 5 files changed, 32 insertions(+), 38 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 29796cfbc..1564e610b 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -229,8 +229,8 @@ module constitutive logical :: converged_ end function crystallite_stress - module function constitutive_homogenizedC(co,ip,el) result(C) - integer, intent(in) :: co, ip, el + module function constitutive_homogenizedC(ph,me) result(C) + integer, intent(in) :: ph, me real(pReal), dimension(6,6) :: C end function constitutive_homogenizedC diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 8c9104946..4fa7a9627 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -436,7 +436,8 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) sourceType: select case (phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & + call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(material_phaseAt(co,el), & + material_phaseMemberAt(co,ip,el)), Fe, & co, ip, el) broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of))) if(.not. broken) then diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 50afe2e87..cd1dedb61 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -307,13 +307,9 @@ submodule(constitutive) constitutive_mech character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results - module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC) - real(pReal), dimension(6,6) :: & - homogenizedC - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) + real(pReal), dimension(6,6) :: homogenizedC + integer, intent(in) :: ph,me end function plastic_dislotwin_homogenizedC @@ -559,16 +555,16 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & 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 + real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3,3,3) :: C integer :: & ho, & !< homogenization - d !< counter in degradation loop - integer :: & - i, j + d, & !< counter in degradation loop + i, j, ph, me ho = material_homogenizationAt(el) - C = math_66toSym3333(constitutive_homogenizedC(co,ip,el)) + C = math_66toSym3333(constitutive_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))) DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) @@ -1535,19 +1531,16 @@ end subroutine mech_forward !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- -module function constitutive_homogenizedC(co,ip,el) result(C) +module function constitutive_homogenizedC(ph,me) result(C) real(pReal), dimension(6,6) :: C - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer, intent(in) :: ph, me - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + plasticityType: select case (phase_plasticity(ph)) case (PLASTICITY_DISLOTWIN_ID) plasticityType - C = plastic_dislotwin_homogenizedC(co,ip,el) + C = plastic_dislotwin_homogenizedC(ph,me) case default plasticityType - C = lattice_C66(1:6,1:6,material_phaseAt(co,el)) + C = lattice_C66(1:6,1:6,ph) end select plasticityType end function constitutive_homogenizedC diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 0474427fe..dea84f751 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -150,7 +150,7 @@ module function plastic_dislotwin_init() result(myPlasticity) Ninstances = count(myPlasticity) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return - + print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL @@ -485,35 +485,32 @@ end function plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- !> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- -module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC) +module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) + integer, intent(in) :: & + ph, me real(pReal), dimension(6,6) :: & homogenizedC - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - integer :: i, & - of + integer :: i real(pReal) :: f_unrotated - of = material_phasememberAt(co,ip,el) - associate(prm => param(phase_plasticityInstance(material_phaseAt(co,el))),& - stt => state(phase_plasticityInstance(material_phaseAT(co,el)))) + + associate(prm => param(phase_plasticityInstance(ph)),& + stt => state(phase_plasticityInstance(ph))) f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - - sum(stt%f_tr(1:prm%sum_N_tr,of)) + - sum(stt%f_tw(1:prm%sum_N_tw,me)) & + - sum(stt%f_tr(1:prm%sum_N_tr,me)) homogenizedC = f_unrotated * prm%C66 do i=1,prm%sum_N_tw homogenizedC = homogenizedC & - + stt%f_tw(i,of)*prm%C66_tw(1:6,1:6,i) + + stt%f_tw(i,me)*prm%C66_tw(1:6,1:6,i) enddo do i=1,prm%sum_N_tr homogenizedC = homogenizedC & - + stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) + + stt%f_tr(i,me)*prm%C66_tr(1:6,1:6,i) enddo end associate diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 580bb0268..24d14e059 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -656,8 +656,11 @@ module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy ip, & !< integration point number el !< element number + real(pReal), dimension(6,6) :: C - equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt') + + C = constitutive_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el)) + equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu From 261f32d7c9d2a912fae4d1add0b52e61389bbfd7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Jan 2021 10:21:51 +0100 Subject: [PATCH 2/8] consistent names --- src/constitutive_damage.f90 | 16 +++++++--------- src/source_damage_isoBrittle.f90 | 32 ++++++++++++++++---------------- 2 files changed, 23 insertions(+), 25 deletions(-) diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 4fa7a9627..46f2a3057 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -411,14 +411,14 @@ 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) +function constitutive_damage_deltaState(Fe, co, ip, el, ph, me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point + co, & !< component-ID me integration point ip, & !< integration point el, & !< element ph, & - of + me real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & @@ -436,15 +436,13 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) sourceType: select case (phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(material_phaseAt(co,el), & - material_phaseMemberAt(co,ip,el)), Fe, & - co, ip, el) - broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of))) + call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, co, ip, el) + broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) if(.not. broken) then 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) + damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = & + damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me) endif end select sourceType diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 1721b0201..b2da8411f 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -28,7 +28,7 @@ contains !-------------------------------------------------------------------------------------------------- module function source_damage_isoBrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & @@ -52,7 +52,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) allocate(source_damage_isoBrittle_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) + phase => phases%get(p) if(any(mySources(:,p))) source_damage_isoBrittle_instance(p) = count(mySources(:,1:p)) if(count(mySources(:,p)) == 0) cycle sources => phase%get('source') @@ -60,7 +60,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) if(mySources(sourceOffset,p)) then source_damage_isoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_isoBrittle_instance(p))) - src => sources%get(sourceOffset) + src => sources%get(sourceOffset) prm%W_crit = src%get_asFloat('W_crit') @@ -69,7 +69,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) #else prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) #endif - + ! sanity checks if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' @@ -106,31 +106,31 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) C integer :: & - phase, & - constituent, & + ph, & + me, & sourceOffset real(pReal), dimension(6) :: & strain real(pReal) :: & strainenergy - phase = material_phaseAt(co,el) !< phase ID at co,ip,el - constituent = material_phasememberAt(co,ip,el) !< state array offset for phase ID at co,ip,el - sourceOffset = source_damage_isoBrittle_offset(phase) + ph = material_phaseAt(co,el) !< ph ID at co,ip,el + me = material_phasememberAt(co,ip,el) !< state array offset for ph ID at co,ip,el + sourceOffset = source_damage_isoBrittle_offset(ph) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - associate(prm => param(source_damage_isoBrittle_instance(phase))) + associate(prm => param(source_damage_isoBrittle_instance(ph))) 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 > damageState(phase)%p(sourceOffset)%subState0(1,constituent)) then - damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent) + if (strainenergy > damageState(ph)%p(sourceOffset)%subState0(1,me)) then + damageState(ph)%p(sourceOffset)%deltaState(1,me) = & + strainenergy - damageState(ph)%p(sourceOffset)%state(1,me) else - damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - damageState(phase)%p(sourceOffset)%subState0(1,constituent) - & - damageState(phase)%p(sourceOffset)%state(1,constituent) + damageState(ph)%p(sourceOffset)%deltaState(1,me) = & + damageState(ph)%p(sourceOffset)%subState0(1,me) - & + damageState(ph)%p(sourceOffset)%state(1,me) endif end associate From 47724be32b4e85ea4c564da8e4b6435fe3534699 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Jan 2021 10:25:52 +0100 Subject: [PATCH 3/8] simplified --- src/constitutive.f90 | 13 ------------- src/constitutive_damage.f90 | 17 +++++++++++------ src/source_damage_isoBrittle.f90 | 11 ++--------- 3 files changed, 13 insertions(+), 28 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 1564e610b..d54f55158 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -343,19 +343,6 @@ module constitutive dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) end subroutine kinematics_thermal_expansion_LiAndItsTangent - - module subroutine source_damage_isoBrittle_deltaState(C, Fe, 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 - real(pReal), intent(in), dimension(6,6) :: & - C - end subroutine source_damage_isoBrittle_deltaState - - module subroutine constitutive_plastic_dependentState(co,ip,el) integer, intent(in) :: & co, & !< component-ID of integration point diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 46f2a3057..235fb8967 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -45,6 +45,14 @@ submodule(constitutive) constitutive_damage logical, dimension(:,:), allocatable :: myKinematics end function kinematics_slipplane_opening_init + module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me) + integer, intent(in) :: ph,me + real(pReal), intent(in), dimension(3,3) :: & + Fe + real(pReal), intent(in), dimension(6,6) :: & + C + end subroutine source_damage_isoBrittle_deltaState + module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & @@ -295,7 +303,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) enddo if(converged_) then - broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) + broken = constitutive_damage_deltaState(mech_F_e(ph,me),ph,me) exit iteration endif @@ -411,12 +419,9 @@ 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, me) result(broken) +function constitutive_damage_deltaState(Fe, ph, me) result(broken) integer, intent(in) :: & - co, & !< component-ID me integration point - ip, & !< integration point - el, & !< element ph, & me real(pReal), intent(in), dimension(3,3) :: & @@ -436,7 +441,7 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, me) result(broken) sourceType: select case (phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, co, ip, el) + call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, ph,me) broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) if(.not. broken) then myOffset = damageState(ph)%p(so)%offsetDeltaState diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index b2da8411f..cff993e7f 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -94,28 +94,21 @@ end function source_damage_isoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) +module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & C integer :: & - ph, & - me, & sourceOffset real(pReal), dimension(6) :: & strain real(pReal) :: & strainenergy - ph = material_phaseAt(co,el) !< ph ID at co,ip,el - me = material_phasememberAt(co,ip,el) !< state array offset for ph ID at co,ip,el sourceOffset = source_damage_isoBrittle_offset(ph) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) From ff1dbfbb952a819c19a885689f52183cfb737c65 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Jan 2021 10:30:10 +0100 Subject: [PATCH 4/8] only needed at lower level --- src/constitutive.f90 | 25 ------------------------- src/constitutive_damage.f90 | 24 ++++++++++++++++++++++++ src/constitutive_thermal.f90 | 7 +++++++ 3 files changed, 31 insertions(+), 25 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index d54f55158..7ecff0eef 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -234,34 +234,9 @@ module constitutive real(pReal), dimension(6,6) :: C end function constitutive_homogenizedC - module subroutine source_damage_anisoBrittle_dotState(S, 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 - end subroutine source_damage_anisoBrittle_dotState - module subroutine source_damage_anisoDuctile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_anisoDuctile_dotState - module subroutine source_damage_isoDuctile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_isoDuctile_dotState - module subroutine source_thermal_externalheat_dotState(phase, of) - integer, intent(in) :: & - phase, & - of - end subroutine source_thermal_externalheat_dotState module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) integer, intent(in) :: & diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 235fb8967..6565e1813 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -54,6 +54,30 @@ submodule(constitutive) constitutive_damage end subroutine source_damage_isoBrittle_deltaState + module subroutine source_damage_anisoBrittle_dotState(S, 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 + end subroutine source_damage_anisoBrittle_dotState + + module subroutine source_damage_anisoDuctile_dotState(co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_anisoDuctile_dotState + + module subroutine source_damage_isoDuctile_dotState(co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_isoDuctile_dotState + + module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & !< phase ID of element diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 9007ef729..511d44f8e 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -36,6 +36,13 @@ submodule(constitutive) constitutive_thermal end function kinematics_thermal_expansion_init + module subroutine source_thermal_externalheat_dotState(phase, of) + integer, intent(in) :: & + phase, & + of + end subroutine source_thermal_externalheat_dotState + + module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase) integer, intent(in) :: & phase !< phase ID of element From b1674b6835067341180c21407bec86d996ad0141 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Jan 2021 10:32:56 +0100 Subject: [PATCH 5/8] consistent names --- src/constitutive_thermal.f90 | 12 ++++++------ src/constitutive_thermal_externalheat.f90 | 24 +++++++++++------------ 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 511d44f8e..2746ff889 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -36,10 +36,10 @@ submodule(constitutive) constitutive_thermal end function kinematics_thermal_expansion_init - module subroutine source_thermal_externalheat_dotState(phase, of) + module subroutine source_thermal_externalheat_dotState(ph, me) integer, intent(in) :: & - phase, & - of + ph, & + me end subroutine source_thermal_externalheat_dotState @@ -54,10 +54,10 @@ submodule(constitutive) constitutive_thermal TDot end subroutine thermal_dissipation_getRate - module subroutine thermal_externalheat_getRate(TDot, phase,of) + module subroutine thermal_externalheat_getRate(TDot, ph,me) integer, intent(in) :: & - phase, & - of + ph, & + me real(pReal), intent(out) :: & TDot end subroutine thermal_externalheat_getRate diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 index 2e8c02f8c..16b2bdb65 100644 --- a/src/constitutive_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -85,18 +85,18 @@ end function source_thermal_externalheat_init !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_externalheat_dotState(phase, of) +module subroutine source_thermal_externalheat_dotState(ph, me) integer, intent(in) :: & - phase, & - of + ph, & + me integer :: & sourceOffset - sourceOffset = source_thermal_externalheat_offset(phase) + sourceOffset = source_thermal_externalheat_offset(ph) - thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time + thermalState(ph)%p(sourceOffset)%dotState(1,me) = 1.0_pReal ! state is current time end subroutine source_thermal_externalheat_dotState @@ -104,11 +104,11 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module subroutine thermal_externalheat_getRate(TDot, phase, of) +module subroutine thermal_externalheat_getRate(TDot, ph, me) integer, intent(in) :: & - phase, & - of + ph, & + me real(pReal), intent(out) :: & TDot @@ -117,18 +117,18 @@ module subroutine thermal_externalheat_getRate(TDot, phase, of) real(pReal) :: & frac_time - sourceOffset = source_thermal_externalheat_offset(phase) + sourceOffset = source_thermal_externalheat_offset(ph) - associate(prm => param(source_thermal_externalheat_instance(phase))) + associate(prm => param(source_thermal_externalheat_instance(ph))) do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & + frac_time = (thermalState(ph)%p(sourceOffset)%state(1,me) - 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) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... - ! ...or extrapolate if outside of bounds + ! ...or extrapolate if outside me bounds enddo end associate From 440790ca01a38db9029f2f966f799d115b6e8f37 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Jan 2021 10:39:16 +0100 Subject: [PATCH 6/8] consistent names --- src/constitutive_damage.f90 | 11 ++++----- src/source_damage_anisoBrittle.f90 | 38 +++++++++++++++--------------- src/source_damage_anisoDuctile.f90 | 22 ++++++++--------- src/source_damage_isoDuctile.f90 | 24 +++++++++---------- 4 files changed, 47 insertions(+), 48 deletions(-) diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 6565e1813..a8792d5e2 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -400,14 +400,14 @@ 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) +function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point + co, & !< component-ID me integration point ip, & !< integration point el, & !< element ph, & - of + me integer :: & so !< counter in source loop logical :: broken @@ -426,12 +426,11 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) 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? + call source_damage_anisoBrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress? end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,of))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me))) enddo SourceLoop diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 7c00c6580..8da37c5d5 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -15,7 +15,7 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle dot_o, & !< opening rate of cleavage planes q !< damage rate sensitivity real(pReal), dimension(:), allocatable :: & - s_crit, & !< critical displacement + s_crit, & !< critical displacement g_crit !< critical load real(pReal), dimension(:,:,:,:), allocatable :: & cleavage_systems @@ -37,7 +37,7 @@ contains !-------------------------------------------------------------------------------------------------- module function source_damage_anisoBrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & @@ -62,7 +62,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) allocate(source_damage_anisoBrittle_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) + phase => phases%get(p) if(any(mySources(:,p))) source_damage_anisoBrittle_instance(p) = count(mySources(:,1:p)) if(count(mySources(:,p)) == 0) cycle sources => phase%get('source') @@ -70,20 +70,20 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) if(mySources(sourceOffset,p)) then source_damage_anisoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_anisoBrittle_instance(p))) - src => sources%get(sourceOffset) - + src => sources%get(sourceOffset) + N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) prm%sum_N_cl = sum(abs(N_cl)) - + prm%q = src%get_asFloat('q') prm%dot_o = src%get_asFloat('dot_o') - + prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl)) prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl)) - + prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - + ! expand: family => system prm%s_crit = math_expand(prm%s_crit,N_cl) prm%g_crit = math_expand(prm%g_crit,N_cl) @@ -93,7 +93,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) #else prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) #endif - + ! sanity checks if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' @@ -130,8 +130,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) S integer :: & - phase, & - constituent, & + ph, & + me, & sourceOffset, & damageOffset, & homog, & @@ -139,14 +139,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit - phase = material_phaseAt(co,el) - constituent = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_anisoBrittle_offset(phase) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + sourceOffset = source_damage_anisoBrittle_offset(ph) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(source_damage_anisoBrittle_instance(phase))) - damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal + associate(prm => param(source_damage_anisoBrittle_instance(ph))) + damageState(ph)%p(sourceOffset)%dotState(1,me) = 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 - damageState(phase)%p(sourceOffset)%dotState(1,constituent) & - = damageState(phase)%p(sourceOffset)%dotState(1,constituent) & + damageState(ph)%p(sourceOffset)%dotState(1,me) & + = damageState(ph)%p(sourceOffset)%dotState(1,me) & + 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 + & diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 7ec06cb62..b5e6fd17c 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -30,7 +30,7 @@ contains !-------------------------------------------------------------------------------------------------- module function source_damage_anisoDuctile_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & @@ -67,7 +67,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) if(mySources(sourceOffset,p)) then source_damage_anisoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_anisoDuctile_instance(p))) - src => sources%get(sourceOffset) + src => sources%get(sourceOffset) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) prm%q = src%get_asFloat('q') @@ -81,7 +81,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) #else prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) #endif - + ! sanity checks if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' @@ -115,21 +115,21 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el) el !< element integer :: & - phase, & - constituent, & + ph, & + me, & sourceOffset, & damageOffset, & homog - phase = material_phaseAt(co,el) - constituent = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_anisoDuctile_offset(phase) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + sourceOffset = source_damage_anisoDuctile_offset(ph) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(source_damage_anisoDuctile_instance(phase))) - damageState(phase)%p(sourceOffset)%dotState(1,constituent) & - = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) + associate(prm => param(source_damage_anisoDuctile_instance(ph))) + damageState(ph)%p(sourceOffset)%dotState(1,me) & + = sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) end associate end subroutine source_damage_anisoDuctile_dotState diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index dd2910182..3cebc474a 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -30,7 +30,7 @@ contains !-------------------------------------------------------------------------------------------------- module function source_damage_isoDuctile_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & @@ -54,7 +54,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) allocate(source_damage_isoDuctile_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) + phase => phases%get(p) if(count(mySources(:,p)) == 0) cycle if(any(mySources(:,p))) source_damage_isoDuctile_instance(p) = count(mySources(:,1:p)) sources => phase%get('source') @@ -62,7 +62,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) if(mySources(sourceOffset,p)) then source_damage_isoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_isoDuctile_instance(p))) - src => sources%get(sourceOffset) + src => sources%get(sourceOffset) prm%q = src%get_asFloat('q') prm%gamma_crit = src%get_asFloat('gamma_crit') @@ -72,7 +72,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) #else prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) #endif - + ! sanity checks if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' @@ -106,21 +106,21 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el) el !< element integer :: & - phase, & - constituent, & + ph, & + me, & sourceOffset, & damageOffset, & homog - phase = material_phaseAt(co,el) - constituent = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_isoDuctile_offset(phase) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + sourceOffset = source_damage_isoDuctile_offset(ph) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(source_damage_isoDuctile_instance(phase))) - damageState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit + associate(prm => param(source_damage_isoDuctile_instance(ph))) + damageState(ph)%p(sourceOffset)%dotState(1,me) = & + sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit end associate end subroutine source_damage_isoDuctile_dotState From 219529a4029ab92e82a8f43a02812128dd13a5fb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Jan 2021 20:40:14 +0100 Subject: [PATCH 7/8] field variables --- src/homogenization.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index a63487b6f..d576c5672 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -29,7 +29,9 @@ module homogenization ! General variables for the homogenization at a material point real(pReal), dimension(:), allocatable, public :: & homogenization_T, & - homogenization_dot_T + homogenization_dot_T, & + homogenization_phi, & + homogenization_dot_phi 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 From 32bb0d8c6e5b4b223a8591133b0d2283b755dcd6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Jan 2021 20:54:31 +0100 Subject: [PATCH 8/8] new variables for damage --- PRIVATE | 2 +- src/commercialFEM_fileList.f90 | 1 + src/constitutive.f90 | 12 ++++++ src/constitutive_damage.f90 | 42 ++++++++++++++++++- src/grid/grid_damage_spectral.f90 | 69 ++++++++++++++++--------------- src/homogenization.f90 | 11 ++++- src/homogenization_damage.f90 | 40 ++++++++++++++++++ 7 files changed, 141 insertions(+), 36 deletions(-) create mode 100644 src/homogenization_damage.f90 diff --git a/PRIVATE b/PRIVATE index 9e8011876..5fd23c60c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 9e8011876b16896a18902737790eba8018e94f85 +Subproject commit 5fd23c60c721aada27bcd3a4ba96301753b7bd0d diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index e04936a3b..963871c7a 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -52,4 +52,5 @@ #include "homogenization_mech_isostrain.f90" #include "homogenization_mech_RGC.f90" #include "homogenization_thermal.f90" +#include "homogenization_damage.f90" #include "CPFEM.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7ecff0eef..16f3d3b59 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -188,6 +188,11 @@ module constitutive real(pReal), dimension(3,3) :: P end function constitutive_mech_getP + module function constitutive_damage_get_phi(co,ip,el) result(phi) + integer, intent(in) :: co, ip, el + real(pReal) :: phi + end function constitutive_damage_get_phi + module function thermal_T(ph,me) result(T) integer, intent(in) :: ph,me real(pReal) :: T @@ -203,6 +208,11 @@ module constitutive real(pReal), intent(in) :: T integer, intent(in) :: co, ce end subroutine constitutive_thermal_setT + + module subroutine constitutive_damage_set_phi(phi,co,ce) + real(pReal), intent(in) :: phi + integer, intent(in) :: co, ce + end subroutine constitutive_damage_set_phi ! == cleaned:end =================================================================================== @@ -353,6 +363,8 @@ module constitutive constitutive_restartRead, & integrateDamageState, & constitutive_thermal_setT, & + constitutive_damage_set_phi, & + constitutive_damage_get_phi, & constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index a8792d5e2..07bfbeeca 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -10,9 +10,16 @@ submodule(constitutive) constitutive_damage DAMAGE_ANISODUCTILE_ID end enum + + type :: tDataContainer + real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi + end type tDataContainer + integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & phase_source !< active sources mechanisms of each phase + type(tDataContainer), dimension(:), allocatable :: current + interface module function source_damage_anisoBrittle_init(source_length) result(mySources) @@ -152,7 +159,8 @@ contains module subroutine damage_init integer :: & - ph !< counter in phase loop + ph, & !< counter in phase loop + Nconstituents class(tNode), pointer :: & phases, & phase, & @@ -161,10 +169,18 @@ module subroutine damage_init phases => config_material%get('phase') + allocate(current(phases%length)) + allocate(damageState (phases%length)) allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics do ph = 1,phases%length + + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + + allocate(current(ph)%phi(Nconstituents),source=1.0_pReal) + allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) + phase => phases%get(ph) sources => phase%get('source',defaultVal=emptyList) phase_Nsources(ph) = sources%length @@ -511,4 +527,28 @@ function source_active(source_label,src_length) result(active_source) end function source_active +!---------------------------------------------------------------------------------------------- +!< @brief Set damage parameter +!---------------------------------------------------------------------------------------------- +module subroutine constitutive_damage_set_phi(phi,co,ce) + + real(pReal), intent(in) :: phi + integer, intent(in) :: ce, co + + + current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi + +end subroutine constitutive_damage_set_phi + + +module function constitutive_damage_get_phi(co,ip,el) result(phi) + + integer, intent(in) :: co, ip, el + real(pReal) :: phi + + phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el)) + +end function constitutive_damage_get_phi + + end submodule constitutive_damage diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 79437945b..a891f070d 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -17,8 +17,9 @@ module grid_damage_spectral use discretization_grid use damage_nonlocal use YAML_types + use homogenization use config - + implicit none private @@ -43,13 +44,13 @@ module grid_damage_spectral phi_current, & !< field of current damage phi_lastInc, & !< field of previous damage phi_stagInc !< field of staggered damage - + !-------------------------------------------------------------------------------------------------- -! reference diffusion tensor, mobility etc. +! reference diffusion tensor, mobility etc. integer :: totalIter = 0 !< total iteration in current increment real(pReal), dimension(3,3) :: K_ref real(pReal) :: mu_ref - + public :: & grid_damage_spectral_init, & grid_damage_spectral_solution, & @@ -62,8 +63,8 @@ contains ! ToDo: Restart not implemented !-------------------------------------------------------------------------------------------------- subroutine grid_damage_spectral_init - - PetscInt, dimension(0:worldsize-1) :: localK + + PetscInt, dimension(0:worldsize-1) :: localK DM :: damage_grid Vec :: uBound, lBound PetscErrorCode :: ierr @@ -72,22 +73,22 @@ subroutine grid_damage_spectral_init num_generic character(len=pStringLen) :: & snes_type - + print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal) num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal) - + num_generic => config_numerics%get('generic',defaultVal=emptyDict) num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) - + if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol') @@ -104,7 +105,7 @@ subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -122,7 +123,7 @@ subroutine grid_damage_spectral_init call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) + CHKERRQ(ierr) call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) if (trim(snes_type) == 'vinewtonrsls' .or. & @@ -137,12 +138,12 @@ subroutine grid_damage_spectral_init endif !-------------------------------------------------------------------------------------------------- -! init fields +! init fields call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) CHKERRQ(ierr) xend = xstart + xend - 1 yend = ystart + yend - 1 - zend = zstart + zend - 1 + zend = zstart + zend - 1 allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) @@ -157,26 +158,26 @@ end subroutine grid_damage_spectral_init !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_damage_spectral_solution(timeinc) result(solution) - + real(pReal), intent(in) :: & timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm, solnNorm - + PetscErrorCode :: ierr SNESConvergedReason :: reason solution%converged =.false. - + !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data params%timeinc = timeinc - + call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) - + if (reason < 1) then solution%converged = .false. solution%iterationsNeeded = num%itmax @@ -192,20 +193,21 @@ function grid_damage_spectral_solution(timeinc) result(solution) solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm) !-------------------------------------------------------------------------------------------------- -! updating damage state +! updating damage state cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo - + call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) if (solution%converged) & print'(/,a)', ' ... nonlocal damage converged .....................................' print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,a)', ' ===========================================================================' - flush(IO_STDOUT) + flush(IO_STDOUT) end function grid_damage_spectral_solution @@ -214,31 +216,32 @@ end function grid_damage_spectral_solution !> @brief spectral damage forwarding routine !-------------------------------------------------------------------------------------------------- subroutine grid_damage_spectral_forward(cutBack) - + logical, intent(in) :: cutBack integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr - if (cutBack) then + if (cutBack) then phi_current = phi_lastInc phi_stagInc = phi_lastInc !-------------------------------------------------------------------------------------------------- -! reverting damage field state +! reverting damage field state cell = 0 call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 + cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo else phi_lastInc = phi_current call updateReference - endif + endif end subroutine grid_damage_spectral_forward @@ -247,7 +250,7 @@ end subroutine grid_damage_spectral_forward !> @brief forms the spectral damage residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,f_scal,dummy,ierr) - + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & @@ -261,11 +264,11 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) integer :: i, j, k, cell real(pReal) :: phiDot, dPhiDot_dPhi, mobility - phi_current = x_scal + phi_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current + scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_FFTvectorBackward @@ -297,7 +300,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness - + !-------------------------------------------------------------------------------------------------- ! constructing residual f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current @@ -311,7 +314,7 @@ end subroutine formResidual subroutine updateReference integer :: i,j,k,cell,ierr - + cell = 0 K_ref = 0.0_pReal mu_ref = 0.0_pReal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index d576c5672..ed9b7fead 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -64,6 +64,9 @@ module homogenization module subroutine thermal_init end subroutine thermal_init + module subroutine damage_init + end subroutine damage_init + module subroutine mech_partition(subF,ip,el) real(pReal), intent(in), dimension(3,3) :: & subF @@ -77,6 +80,11 @@ module homogenization integer, intent(in) :: ce end subroutine thermal_partition + module subroutine damage_partition(phi,ce) + real(pReal), intent(in) :: phi + integer, intent(in) :: ce + end subroutine damage_partition + module subroutine thermal_homogenize(ip,el) integer, intent(in) :: ip,el end subroutine thermal_homogenize @@ -120,7 +128,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !-------------------------------------------------------------------------------------------------- -subroutine homogenization_init +subroutine homogenization_init() class (tNode) , pointer :: & num_homog, & @@ -144,6 +152,7 @@ subroutine homogenization_init call mech_init(num_homog) call thermal_init() + call damage_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) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 new file mode 100644 index 000000000..479318340 --- /dev/null +++ b/src/homogenization_damage.f90 @@ -0,0 +1,40 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!-------------------------------------------------------------------------------------------------- +submodule(homogenization) homogenization_damage + + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief Allocate variables and set parameters. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_init() + + + print'(/,a)', ' <<<+- homogenization_damage init -+>>>' + + allocate(homogenization_phi(discretization_nIPs*discretization_Nelems)) + allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems)) + +end subroutine damage_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Partition temperature onto the individual constituents. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_partition(phi,ce) + + real(pReal), intent(in) :: phi + integer, intent(in) :: ce + + integer :: co + + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + call constitutive_damage_set_phi(phi,co,ce) + enddo + +end subroutine damage_partition + + +end submodule homogenization_damage