diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5f323c3a2..803a9251d 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -43,82 +43,65 @@ submodule(phase) damagee end function isoductile_init - module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me) + module subroutine 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 + end subroutine isobrittle_deltaState - module subroutine anisobrittle_dotState(S, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine anisobrittle_dotState(S, ph, me) + integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & S end subroutine anisobrittle_dotState - module subroutine anisoductile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine anisoductile_dotState(ph,me) + integer, intent(in) :: ph,me end subroutine anisoductile_dotState - module subroutine isoductile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine isoductile_dotState(ph,me) + integer, intent(in) :: ph,me end subroutine isoductile_dotState - module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoBrittle_getRateAndItsTangent + end subroutine anisobrittle_getRateAndItsTangent - module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoDuctile_getRateAndItsTangent + end subroutine anisoductile_getRateAndItsTangent - module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoBrittle_getRateAndItsTangent + end subroutine isobrittle_getRateAndItsTangent - module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoDuctile_getRateAndItsTangent + end subroutine isoductile_getRateAndItsTangent module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase @@ -225,16 +208,16 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, do so = 1, phase_Nsources(ph) select case(phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ISODUCTILE_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISOBRITTLE_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISODUCTILE_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case default localphiDot = 0.0_pReal @@ -281,7 +264,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) me = material_phaseMemberAt(co,ip,el) converged_ = .true. - broken = phase_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(ph,me) if(broken) return size_so = damageState(ph)%sizeDotState @@ -294,7 +277,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) - broken = phase_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(ph,me) if(broken) exit iteration @@ -384,39 +367,34 @@ end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(co,ip,el,ph,me) result(broken) +function phase_damage_collectDotState(ph,me) result(broken) integer, intent(in) :: & - co, & !< component-ID me integration point - ip, & !< integration point - el, & !< element ph, & - me - integer :: & - so !< counter in source loop + me !< counter in source loop logical :: broken broken = .false. - SourceLoop: do so = 1, phase_Nsources(ph) + if (phase_Nsources(ph)==1) then - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(1,ph)) case (DAMAGE_ISODUCTILE_ID) sourceType - call isoductile_dotState(co, ip, el) + call isoductile_dotState(ph,me) case (DAMAGE_ANISODUCTILE_ID) sourceType - call anisoductile_dotState(co, ip, el) + call anisoductile_dotState(ph,me) case (DAMAGE_ANISOBRITTLE_ID) sourceType - call anisobrittle_dotState(mechanical_S(ph,me),co, ip, el) ! correct stress? + call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? end select sourceType broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) - enddo SourceLoop + endif end function phase_damage_collectDotState @@ -447,7 +425,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) sourceType: select case (phase_source(1,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) + call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) if(.not. broken) then myOffset = damageState(ph)%offsetDeltaState diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 81072e16d..a121ee932 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -113,18 +113,14 @@ end function anisobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_dotState(S, co, ip, el) +module subroutine anisobrittle_dotState(S, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S integer :: & - ph, & - me, & sourceOffset, & damageOffset, & homog, & @@ -133,10 +129,6 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_d, traction_t, traction_n, traction_crit - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - - associate(prm => param(ph)) damageState(ph)%dotState(1,me) = 0.0_pReal do i = 1, prm%sum_N_cl @@ -144,7 +136,7 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) & + prm%dot_o / prm%s_crit(i) & @@ -160,11 +152,11 @@ end subroutine anisobrittle_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -172,12 +164,12 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d dLocalphiDot_dPhi - dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_anisoBrittle_getRateAndItsTangent +end subroutine anisobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -190,6 +182,7 @@ module subroutine anisobrittle_results(phase,group) integer :: o + associate(prm => param(phase), stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index 4d2392589..cf55e398b 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -99,20 +99,12 @@ end function anisoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine anisoductile_dotState(co, ip, el) +module subroutine anisoductile_dotState(ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & me - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - associate(prm => param(ph)) damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(damage_phi(ph,me)**prm%q)/prm%gamma_crit) @@ -124,11 +116,11 @@ end subroutine anisoductile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -136,12 +128,12 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d dLocalphiDot_dPhi - dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_anisoDuctile_getRateAndItsTangent +end subroutine anisoductile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -154,8 +146,8 @@ module subroutine anisoductile_results(phase,group) integer :: o - associate(prm => param(phase), & - stt => damageState(phase)%state) + + associate(prm => param(phase), stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 4e37f3ee9..e267ffcc9 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -86,7 +86,7 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) +module subroutine isobrittle_deltaState(C, Fe, ph,me) integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & @@ -106,24 +106,21 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) 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(ph)%subState0(1,me)) then - damageState(ph)%deltaState(1,me) = strainenergy - damageState(ph)%state(1,me) - else - damageState(ph)%deltaState(1,me) = damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me) - endif + damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), & + damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), & + strainenergy > damageState(ph)%subState0(1,me)) end associate -end subroutine source_damage_isoBrittle_deltaState +end subroutine isobrittle_deltaState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -131,13 +128,13 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo dLocalphiDot_dPhi - associate(prm => param(phase)) + associate(prm => param(ph)) localphiDot = 1.0_pReal & - - phi*damageState(phase)%state(1,constituent) - dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent) + - phi*damageState(ph)%state(1,me) + dLocalphiDot_dPhi = - damageState(ph)%state(1,me) end associate -end subroutine source_damage_isoBrittle_getRateAndItsTangent +end subroutine isobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 89211b746..f59952222 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -90,23 +90,16 @@ end function isoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine isoductile_dotState(co, ip, el) +module subroutine isoductile_dotState(ph, me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & me - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - associate(prm => param(ph)) - damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(damage_phi(ph,me)**prm%q)/prm%gamma_crit + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)) & + / (prm%gamma_crit*damage_phi(ph,me)**prm%q) end associate end subroutine isoductile_dotState @@ -115,11 +108,11 @@ end subroutine isoductile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -127,12 +120,12 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo dLocalphiDot_dPhi - dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_isoDuctile_getRateAndItsTangent +end subroutine isoductile_getRateAndItsTangent !--------------------------------------------------------------------------------------------------