From b3dde6d7227b60a830c2486d068ca8b6b5b0343c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 08:32:43 +0100 Subject: [PATCH] only one damage mechanism per phase material.yaml specification is designed to allow more than one, but that requires to have two phase fields etc. For the moment, keep it as simple as possible. --- PRIVATE | 2 +- src/phase.f90 | 43 +++++++++---------- src/phase_damage.f90 | 71 +++++++++++++------------------ src/phase_damage_anisobrittle.f90 | 30 +++++-------- src/phase_damage_anisoductile.f90 | 29 ++++--------- src/phase_damage_isobrittle.f90 | 37 ++++++---------- src/phase_damage_isoductile.f90 | 29 ++++--------- src/phase_mechanical.f90 | 18 ++++---- 8 files changed, 97 insertions(+), 162 deletions(-) diff --git a/PRIVATE b/PRIVATE index 2a0f30ec4..571b06ac4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2a0f30ec47473f42e0da922055b72b527730378d +Subproject commit 571b06ac4b42ed217eacfa5c71281e571bdab1df diff --git a/src/phase.f90 b/src/phase.f90 index 1d08f680d..8ea1fa5f7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -70,7 +70,7 @@ module phase type(tPlasticState), allocatable, dimension(:), public :: & plasticState - type(tSourceState), allocatable, dimension(:), public :: & + type(tState), allocatable, dimension(:), public :: & damageState @@ -350,14 +350,12 @@ subroutine phase_init PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(ph)%state = plasticState(ph)%state0 - forall(so = 1:phase_Nsources(ph)) - damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 - end forall - - phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, & - maxval(damageState(ph)%p%sizeDotState)) + plasticState(ph)%state = plasticState(ph)%state0 + if(phase_Nsources(ph) > 0) & + damageState(ph)%state = damageState(ph)%state0 enddo PhaseLoop2 + + phase_source_maxSizeDotState = maxval(damageState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) end subroutine phase_init @@ -404,15 +402,13 @@ subroutine phase_restore(ce,includeL) integer, intent(in) :: ce integer :: & - co, & !< constituent number - so + co do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - do so = 1, phase_Nsources(material_phaseAt2(co,ce)) - damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) - enddo + if (phase_Nsources(material_phaseAt2(co,ce)) > 0) & + damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) enddo call mechanical_restore(ce,includeL) @@ -426,16 +422,16 @@ end subroutine phase_restore !-------------------------------------------------------------------------------------------------- subroutine phase_forward() - integer :: ph, so + integer :: ph call mechanical_forward() call thermal_forward() do ph = 1, size(damageState) - do so = 1,phase_Nsources(ph) - damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state - enddo; enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state0 = damageState(ph)%state + enddo end subroutine phase_forward @@ -531,9 +527,8 @@ subroutine crystallite_init() phases => config_material%get('phase') do ph = 1, phases%length - do so = 1, phase_Nsources(ph) - allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack - enddo + if (phase_Nsources(ph) > 0) & + allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -579,9 +574,9 @@ subroutine phase_windForward(ip,el) call mechanical_windForward(ph,me) - do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) + enddo diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 8defcf3a1..5f323c3a2 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -176,8 +176,8 @@ module subroutine damage_init phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) - phase_Nsources(ph) = sources%length - allocate(damageState(ph)%p(phase_Nsources(ph))) + if (sources%length > 1) error stop + enddo allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) @@ -267,57 +267,48 @@ module function integrateDamageState(dt,co,ip,el) result(broken) NiterationState, & !< number of iterations in state loop ph, & me, & - so - integer, dimension(maxval(phase_Nsources)) :: & size_so real(pReal) :: & zeta real(pReal), dimension(phase_source_maxSizeDotState) :: & r ! state residuum - real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState logical :: & converged_ - ph = material_phaseAt(co,el) + if (phase_Nsources(ph) == 0) return me = material_phaseMemberAt(co,ip,el) converged_ = .true. broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) return - do so = 1, phase_Nsources(ph) - 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 + size_so = damageState(ph)%sizeDotState + damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & + + damageState(ph)%dotState (1:size_so,me) * dt + source_dotState(1:size_so,2) = 0.0_pReal 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) = damageState(ph)%p(so)%dotState(:,me) - enddo + 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) if(broken) exit iteration - do so = 1, phase_Nsources(ph) - zeta = damper(damageState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - 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)) = 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)), & - damageState(ph)%p(so)%state(1:size_so(so),me), & - damageState(ph)%p(so)%atol(1:size_so(so))) - enddo + + zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) + damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & + + source_dotState(1:size_so,1)* (1.0_pReal - zeta) + r(1:size_so) = damageState(ph)%state (1:size_so,me) & + - damageState(ph)%subState0(1:size_so,me) & + - damageState(ph)%dotState (1:size_so,me) * dt + damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) + converged_ = converged_ .and. converged(r(1:size_so), & + damageState(ph)%state(1:size_so,me), & + damageState(ph)%atol(1:size_so)) + if(converged_) then broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me) @@ -423,7 +414,7 @@ function phase_damage_collectDotState(co,ip,el,ph,me) result(broken) end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) enddo SourceLoop @@ -443,7 +434,6 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & - so, & myOffset, & mySize logical :: & @@ -452,23 +442,22 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) broken = .false. - sourceLoop: do so = 1, phase_Nsources(ph) + if (phase_Nsources(ph) == 0) return - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(1,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) - broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) + broken = any(IEEE_is_NaN(damageState(ph)%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,me) = & - damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me) + myOffset = damageState(ph)%offsetDeltaState + mySize = damageState(ph)%sizeDeltaState + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) endif end select sourceType - enddo SourceLoop end function phase_damage_deltaState diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index ad7f80604..422b2a91a 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -6,8 +6,7 @@ !-------------------------------------------------------------------------------------------------- submodule (phase:damagee) anisobrittle - integer, dimension(:), allocatable :: & - source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? + integer, dimension(:), allocatable :: & source_damage_anisoBrittle_instance !< instance of source mechanism type :: tParameters !< container type for internal constitutive parameters @@ -58,7 +57,6 @@ module function anisobrittle_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_anisoBrittle_offset (phases%length), source=0) allocate(source_damage_anisoBrittle_instance(phases%length), source=0) do p = 1, phases%length @@ -68,7 +66,6 @@ module function anisobrittle_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_anisoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_anisoBrittle_instance(p))) src => sources%get(sourceOffset) @@ -101,9 +98,9 @@ module function 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 phase_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' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' end associate @@ -139,14 +136,13 @@ module subroutine anisobrittle_dotState(S, co, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit + 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(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal + damageState(ph)%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 +150,8 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = damageState(ph)%p(sourceOffset)%dotState(1,me) & + damageState(ph)%dotState(1,me) & + = damageState(ph)%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 + & @@ -180,12 +176,8 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoBrittle_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -204,7 +196,7 @@ module subroutine anisobrittle_results(phase,group) integer :: o associate(prm => param(source_damage_anisoBrittle_instance(phase)), & - stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) + 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_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index b22f5b3d6..bafcf80a2 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) anisoductile integer, dimension(:), allocatable :: & - source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism type :: tParameters !< container type for internal constitutive parameters @@ -53,7 +52,6 @@ module function anisoductile_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_anisoDuctile_offset (phases%length), source=0) allocate(source_damage_anisoDuctile_instance(phases%length), source=0) do p = 1, phases%length @@ -65,7 +63,6 @@ module function anisoductile_init(source_length) result(mySources) sources => phase%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_anisoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_anisoDuctile_instance(p))) src => sources%get(sourceOffset) @@ -87,9 +84,9 @@ module function 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 phase_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' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -116,20 +113,14 @@ module subroutine anisoductile_dotState(co, ip, el) integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me 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(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) end associate end subroutine anisoductile_dotState @@ -149,12 +140,8 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -173,7 +160,7 @@ module subroutine anisoductile_results(phase,group) integer :: o associate(prm => param(source_damage_anisoDuctile_instance(phase)), & - stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) + 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 a2ae7dfb7..0a87493a6 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) isobrittle integer, dimension(:), allocatable :: & - source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance type :: tParameters !< container type for internal constitutive parameters @@ -48,7 +47,6 @@ module function isobrittle_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_isoBrittle_offset (phases%length), source=0) allocate(source_damage_isoBrittle_instance(phases%length), source=0) do p = 1, phases%length @@ -58,7 +56,6 @@ module function isobrittle_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_isoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_isoBrittle_instance(p))) src => sources%get(sourceOffset) @@ -74,9 +71,9 @@ module function 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 phase_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' + call phase_allocateState(damageState(p),Nconstituents,1,1,1) + damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -102,14 +99,11 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) real(pReal), intent(in), dimension(6,6) :: & C - integer :: & - sourceOffset real(pReal), dimension(6) :: & strain real(pReal) :: & strainenergy - sourceOffset = source_damage_isoBrittle_offset(ph) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) @@ -117,14 +111,11 @@ 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)%p(sourceOffset)%subState0(1,me)) then - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - strainenergy - damageState(ph)%p(sourceOffset)%state(1,me) - else - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - damageState(ph)%p(sourceOffset)%subState0(1,me) - & - damageState(ph)%p(sourceOffset)%state(1,me) - endif + 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 end associate end subroutine source_damage_isoBrittle_deltaState @@ -144,15 +135,11 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - - sourceOffset = source_damage_isoBrittle_offset(phase) associate(prm => param(source_damage_isoBrittle_instance(phase))) - localphiDot = 1.0_pReal & - - phi*damageState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) + localphiDot = 1.0_pReal & + - phi*damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent) end associate end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -169,7 +156,7 @@ module subroutine isobrittle_results(phase,group) integer :: o associate(prm => param(source_damage_isoBrittle_instance(phase)), & - stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state) + 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_isoductile.f90 b/src/phase_damage_isoductile.f90 index f378b7ad2..657f9da39 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) isoductile integer, dimension(:), allocatable :: & - source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism type:: tParameters !< container type for internal constitutive parameters @@ -50,7 +49,6 @@ module function isoductile_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_isoDuctile_offset (phases%length), source=0) allocate(source_damage_isoDuctile_instance(phases%length), source=0) do p = 1, phases%length @@ -60,7 +58,6 @@ module function isoductile_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_isoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_isoDuctile_instance(p))) src => sources%get(sourceOffset) @@ -78,9 +75,9 @@ module function 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 phase_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' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -107,20 +104,14 @@ module subroutine isoductile_dotState(co, ip, el) integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me 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(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = & - sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit end associate end subroutine isoductile_dotState @@ -140,12 +131,8 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_isoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -164,7 +151,7 @@ module subroutine isoductile_results(phase,group) integer :: o associate(prm => param(source_damage_isoDuctile_instance(phase)), & - stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state) + 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_mechanical.f90 b/src/phase_mechanical.f90 index f7631a38d..bc4b2cd47 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1138,7 +1138,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal) :: & formerSubStep integer :: & - so, ph, me, sizeDotState + ph, me, sizeDotState logical :: todo real(pReal) :: subFrac,subStep real(pReal), dimension(3,3) :: & @@ -1159,10 +1159,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%State0(:,me) + if (phase_Nsources(ph) > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) - enddo subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me) @@ -1188,9 +1187,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) + endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1204,9 +1203,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 endif plasticState(ph)%state(:,me) = subState0 - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif