From c09c2a6c8e123d42ee33527675a5096275cdf236 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 00:50:42 +0100 Subject: [PATCH 1/7] easier to read without instance --- src/phase_mechanical.f90 | 12 +-- src/phase_mechanical_plastic.f90 | 6 +- ...phase_mechanical_plastic_kinehardening.f90 | 73 ++++++++----------- ...phase_mechanical_plastic_phenopowerlaw.f90 | 69 +++++++++--------- 4 files changed, 74 insertions(+), 86 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index c10a916da..ae3967b26 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -135,13 +135,13 @@ submodule(phase) mechanical character(len=*), intent(in) :: group end subroutine plastic_isotropic_results - module subroutine plastic_phenopowerlaw_results(instance,group) - integer, intent(in) :: instance + module subroutine plastic_phenopowerlaw_results(ph,group) + integer, intent(in) :: ph character(len=*), intent(in) :: group end subroutine plastic_phenopowerlaw_results - module subroutine plastic_kinehardening_results(instance,group) - integer, intent(in) :: instance + module subroutine plastic_kinehardening_results(ph,group) + integer, intent(in) :: ph character(len=*), intent(in) :: group end subroutine plastic_kinehardening_results @@ -406,10 +406,10 @@ module subroutine mechanical_results(group,ph) call plastic_isotropic_results(phase_plasticInstance(ph),group//'plastic/') case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticInstance(ph),group//'plastic/') + call plastic_phenopowerlaw_results(ph,group//'plastic/') case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticInstance(ph),group//'plastic/') + call plastic_kinehardening_results(ph,group//'plastic/') case(PLASTICITY_DISLOTWIN_ID) call plastic_dislotwin_results(phase_plasticInstance(ph),group//'plastic/') diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 843273c72..b08e27276 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -201,11 +201,11 @@ submodule(phase:mechanical) plastic el !< current element number end subroutine nonlocal_dependentState - module subroutine plastic_kinehardening_deltaState(Mp,instance,me) + module subroutine plastic_kinehardening_deltaState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me end subroutine plastic_kinehardening_deltaState @@ -417,7 +417,7 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken) plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_KINEHARDENING_ID) plasticType - call plastic_kinehardening_deltaState(Mp,instance,me) + call plastic_kinehardening_deltaState(Mp,ph,me) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case (PLASTICITY_NONLOCAL_ID) plasticType diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 1244988de..14aa58fb9 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -62,7 +62,6 @@ module function plastic_kinehardening_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, i, o, & Nconstituents, & sizeState, sizeDeltaState, sizeDotState, & @@ -81,25 +80,24 @@ module function plastic_kinehardening_init() result(myPlasticity) pl myPlasticity = plastic_active('kinehardening') - Ninstances = count(myPlasticity) - if(Ninstances == 0) return + if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanical:plastic:kinehardening init -+>>>' - print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(dotState(Ninstances)) - allocate(deltaState(Ninstances)) - phases => config_material%get('phase') - i = 0 + allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(dotState(phases%length)) + allocate(deltaState(phases%length)) + + do p = 1, phases%length + if(.not. myPlasticity(p)) cycle phase => phases%get(p) mech => phase%get('mechanics') - if(.not. myPlasticity(p)) cycle - i = i + 1 + i = p associate(prm => param(i), & dot => dotState(i), & dlt => deltaState(i), & @@ -256,16 +254,16 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) integer :: & i,k,l,m,n - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(phase_plasticInstance(ph))) + associate(prm => param(ph)) - call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,ph,me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) @@ -293,14 +291,14 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,me) real(pReal) :: & sumGamma - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & gdot_pos,gdot_neg - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),& - dot => dotState(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph),& + dot => dotState(ph)) - call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg) + call kinetics(Mp,ph,me,gdot_pos,gdot_neg) dot%accshear(:,me) = abs(gdot_pos+gdot_neg) sumGamma = sum(stt%accshear(:,me)) @@ -326,32 +324,25 @@ end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate (instantaneous) incremental change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_deltaState(Mp,instance,me) +module subroutine plastic_kinehardening_deltaState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & gdot_pos,gdot_neg, & sense - associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) + associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) - call kinetics(Mp,instance,me,gdot_pos,gdot_neg) - sense = merge(state(instance)%sense(:,me), & ! keep existing... + call kinetics(Mp,ph,me,gdot_pos,gdot_neg) + sense = merge(state(ph)%sense(:,me), & ! keep existing... sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction -#ifdef DEBUG - if (debugConstitutive%extensive & - .and. (me == prm%of_debug .or. .not. debugConstitutive%selective)) then - print*, '======= kinehardening delta state =======' - print*, sense,state(instance)%sense(:,me) - endif -#endif !-------------------------------------------------------------------------------------------------- ! switch in sense me shear? @@ -373,14 +364,14 @@ end subroutine plastic_kinehardening_deltaState !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_results(instance,group) +module subroutine plastic_kinehardening_results(ph,group) - integer, intent(in) :: instance + integer, intent(in) :: ph character(len=*), intent(in) :: group integer :: o - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case('xi') @@ -415,28 +406,28 @@ end subroutine plastic_kinehardening_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,instance,me, & +pure subroutine kinetics(Mp,ph,me, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me - real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & + real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & gdot_pos, & gdot_neg - real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & + real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & dgdot_dtau_pos, & dgdot_dtau_neg - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & tau_pos, & tau_neg integer :: i - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) do i = 1, prm%sum_N_sl tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,me) diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index bea7ef06b..94214f14a 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -70,7 +70,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, i, & Nconstituents, & sizeState, sizeDotState, & @@ -91,24 +90,22 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) myPlasticity = plastic_active('phenopowerlaw') - Ninstances = count(myPlasticity) - if(Ninstances == 0) return + if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' - print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(dotState(Ninstances)) - phases => config_material%get('phase') - i = 0 + allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(dotState(phases%length)) + do p = 1, phases%length + if(.not. myPlasticity(p)) cycle phase => phases%get(p) mech => phase%get('mechanics') - if(.not. myPlasticity(p)) cycle - i = i + 1 + i = p associate(prm => param(i), & dot => dotState(i), & stt => state(i)) @@ -302,18 +299,18 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) integer :: & i,k,l,m,n - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & gdot_slip_pos,gdot_slip_neg, & dgdot_dtauslip_pos,dgdot_dtauslip_neg - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: & + real(pReal), dimension(param(ph)%sum_N_tw) :: & gdot_twin,dgdot_dtautwin Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(phase_plasticInstance(ph))) + associate(prm => param(ph)) - call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%sum_N_sl Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -322,7 +319,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - call kinetics_twin(Mp,phase_plasticInstance(ph),me,gdot_twin,dgdot_dtautwin) + call kinetics_twin(Mp,ph,me,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1, prm%sum_N_tw Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -350,12 +347,12 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me) c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& sumGamma,sumF - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), & - dot => dotState(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph), & + dot => dotState(ph)) sumGamma = sum(stt%gamma_slip(:,me)) sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char) @@ -375,9 +372,9 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me) !-------------------------------------------------------------------------------------------------- ! shear rates - call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg) + call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,phase_plasticInstance(ph),me,dot%gamma_twin(:,me)) + call kinetics_twin(Mp,ph,me,dot%gamma_twin(:,me)) !-------------------------------------------------------------------------------------------------- ! hardening @@ -395,14 +392,14 @@ end subroutine phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_phenopowerlaw_results(instance,group) +module subroutine plastic_phenopowerlaw_results(ph,group) - integer, intent(in) :: instance + integer, intent(in) :: ph character(len=*), intent(in) :: group integer :: o - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) @@ -434,28 +431,28 @@ end subroutine plastic_phenopowerlaw_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,instance,me, & +pure subroutine kinetics_slip(Mp,ph,me, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me - real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & + real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & gdot_slip_pos, & gdot_slip_neg - real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & + real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & tau_slip_pos, & tau_slip_neg integer :: i - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) do i = 1, prm%sum_N_sl tau_slip_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) @@ -503,25 +500,25 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,instance,me,& +pure subroutine kinetics_twin(Mp,ph,me,& gdot_twin,dgdot_dtau_twin) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me - real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & gdot_twin - real(pReal), dimension(param(instance)%sum_N_tw), intent(out), optional :: & + real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & dgdot_dtau_twin - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(ph)%sum_N_tw) :: & tau_twin integer :: i - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) do i = 1, prm%sum_N_tw tau_twin(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) From 18971d7d8b54900fcca829cfa2253be4c2da27ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 07:06:14 +0100 Subject: [PATCH 2/7] separation by instance does not add any value --- src/phase_mechanical.f90 | 22 ++-- src/phase_mechanical_eigen.f90 | 2 +- src/phase_mechanical_plastic.f90 | 12 +- ...phase_mechanical_plastic_dislotungsten.f90 | 63 +++++---- src/phase_mechanical_plastic_dislotwin.f90 | 122 +++++++++--------- src/phase_mechanical_plastic_isotropic.f90 | 39 +++--- src/phase_mechanical_plastic_nonlocal.f90 | 2 +- 7 files changed, 129 insertions(+), 133 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index ae3967b26..394c3d3ba 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -60,7 +60,7 @@ submodule(phase) mechanical module subroutine plastic_init end subroutine plastic_init - module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -68,7 +68,7 @@ submodule(phase) mechanical real(pReal), dimension(3,3), intent(in) :: & Mi !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me end subroutine plastic_isotropic_LiAndItsTangent @@ -130,8 +130,8 @@ submodule(phase) mechanical end subroutine plastic_LpAndItsTangents - module subroutine plastic_isotropic_results(instance,group) - integer, intent(in) :: instance + module subroutine plastic_isotropic_results(ph,group) + integer, intent(in) :: ph character(len=*), intent(in) :: group end subroutine plastic_isotropic_results @@ -145,13 +145,13 @@ submodule(phase) mechanical character(len=*), intent(in) :: group end subroutine plastic_kinehardening_results - module subroutine plastic_dislotwin_results(instance,group) - integer, intent(in) :: instance + module subroutine plastic_dislotwin_results(ph,group) + integer, intent(in) :: ph character(len=*), intent(in) :: group end subroutine plastic_dislotwin_results - module subroutine plastic_dislotungsten_results(instance,group) - integer, intent(in) :: instance + module subroutine plastic_dislotungsten_results(ph,group) + integer, intent(in) :: ph character(len=*), intent(in) :: group end subroutine plastic_dislotungsten_results @@ -403,7 +403,7 @@ module subroutine mechanical_results(group,ph) select case(phase_plasticity(ph)) case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticInstance(ph),group//'plastic/') + call plastic_isotropic_results(ph,group//'plastic/') case(PLASTICITY_PHENOPOWERLAW_ID) call plastic_phenopowerlaw_results(ph,group//'plastic/') @@ -412,10 +412,10 @@ module subroutine mechanical_results(group,ph) call plastic_kinehardening_results(ph,group//'plastic/') case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticInstance(ph),group//'plastic/') + call plastic_dislotwin_results(ph,group//'plastic/') case(PLASTICITY_DISLOTUNGSTEN_ID) - call plastic_dislotungsten_results(phase_plasticInstance(ph),group//'plastic/') + call plastic_dislotungsten_results(ph,group//'plastic/') case(PLASTICITY_NONLOCAL_ID) call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/') diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index bee7dc2d2..ce974a6a6 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -180,7 +180,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_isotropic_ID) plasticType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,me) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index b08e27276..a13b5c73f 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -179,17 +179,17 @@ submodule(phase:mechanical) plastic el !< current element number end subroutine nonlocal_dotState - module subroutine dislotwin_dependentState(T,instance,me) + module subroutine dislotwin_dependentState(T,ph,me) integer, intent(in) :: & - instance, & + ph, & me real(pReal), intent(in) :: & T end subroutine dislotwin_dependentState - module subroutine dislotungsten_dependentState(instance,me) + module subroutine dislotungsten_dependentState(ph,me) integer, intent(in) :: & - instance, & + ph, & me end subroutine dislotungsten_dependentState @@ -374,10 +374,10 @@ module subroutine plastic_dependentState(co, ip, el) plasticType: select case (phase_plasticity(material_phaseAt(co,el))) case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_dependentState(thermal_T(ph,me),instance,me) + call dislotwin_dependentState(thermal_T(ph,me),ph,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_dependentState(instance,me) + call dislotungsten_dependentState(ph,me) case (PLASTICITY_NONLOCAL_ID) plasticType call nonlocal_dependentState(instance,me,ip,el) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 6eb1e4abe..81c4605ac 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -78,7 +78,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, i, & Nconstituents, & sizeState, sizeDotState, & @@ -97,29 +96,29 @@ module function plastic_dislotungsten_init() result(myPlasticity) mech, & pl + myPlasticity = plastic_active('dislotungsten') - Ninstances = count(myPlasticity) - if(Ninstances == 0) return + if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotungsten init -+>>>' - print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) - + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002' - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(dotState(Ninstances)) - allocate(dependentState(Ninstances)) phases => config_material%get('phase') - i = 0 + allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(dotState(phases%length)) + allocate(dependentState(phases%length)) + + do p = 1, phases%length + if(.not. myPlasticity(p)) cycle phase => phases%get(p) mech => phase%get('mechanics') - if(.not. myPlasticity(p)) cycle - i = i + 1 + i = p associate(prm => param(i), & dot => dotState(i), & stt => state(i), & @@ -290,16 +289,16 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & integer :: & i,k,l,m,n - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(phase_plasticInstance(ph))) + associate(prm => param(ph)) - call kinetics(Mp,T,phase_plasticInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + call kinetics(Mp,T,ph,me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -328,7 +327,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) real(pReal) :: & VacancyDiffusion - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & gdot_pos, gdot_neg,& tau_pos,& tau_neg, & @@ -337,10 +336,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) dot_rho_dip_climb, & dip_distance - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),& - dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph),& + dot => dotState(ph), dst => dependentState(ph)) - call kinetics(Mp,T,phase_plasticInstance(ph),me,& + call kinetics(Mp,T,ph,me,& gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) @@ -377,16 +376,16 @@ end subroutine dislotungsten_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine dislotungsten_dependentState(instance,me) +module subroutine dislotungsten_dependentState(ph,me) integer, intent(in) :: & - instance, & + ph, & me - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & dislocationSpacing - associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me))) dst%threshold_stress(:,me) = prm%mu*prm%b_sl & @@ -402,14 +401,14 @@ end subroutine dislotungsten_dependentState !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotungsten_results(instance,group) +module subroutine plastic_dislotungsten_results(ph,group) - integer, intent(in) :: instance + integer, intent(in) :: ph character(len=*), intent(in) :: group integer :: o - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case('rho_mob') @@ -441,7 +440,7 @@ end subroutine plastic_dislotungsten_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,T,instance,me, & +pure subroutine kinetics(Mp,T,ph,me, & dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) real(pReal), dimension(3,3), intent(in) :: & @@ -449,18 +448,18 @@ pure subroutine kinetics(Mp,T,instance,me, & real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & + ph, & me - real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & + real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos, & dot_gamma_neg - real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & + real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & ddot_gamma_dtau_pos, & ddot_gamma_dtau_neg, & tau_pos_out, & tau_neg_out - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & StressRatio, & StressRatio_p,StressRatio_pminus1, & dvel, vel, & @@ -469,7 +468,7 @@ pure subroutine kinetics(Mp,T,instance,me, & needsGoodName ! ToDo: @Karo: any idea? integer :: j - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) do j = 1, prm%sum_N_sl tau_pos(j) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,j)) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index ca385f417..3b33196d9 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -48,7 +48,7 @@ submodule(phase:plastic) dislotwin dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system t_tw, & !< twin thickness [m] for each twin system i_sl, & !< Adj. parameter for distance between 2 forest dislocations for each slip system - t_tr, & !< martensite lamellar thickness [m] for each trans system and instance + t_tr, & !< martensite lamellar thickness [m] for each trans system p, & !< p-exponent in glide velocity q, & !< q-exponent in glide velocity r, & !< r-exponent in twin nucleation rate @@ -126,7 +126,6 @@ module function plastic_dislotwin_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, i, & Nconstituents, & sizeState, sizeDotState, & @@ -144,12 +143,12 @@ module function plastic_dislotwin_init() result(myPlasticity) mech, & pl + myPlasticity = plastic_active('dislotwin') - Ninstances = count(myPlasticity) - if(Ninstances == 0) return + if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotwin init -+>>>' - print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL @@ -160,18 +159,19 @@ module function plastic_dislotwin_init() result(myPlasticity) print*, 'Wong et al., Acta Materialia 118:140–151, 2016' print*, 'https://doi.org/10.1016/j.actamat.2016.07.032' - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(dotState(Ninstances)) - allocate(dependentState(Ninstances)) phases => config_material%get('phase') - i = 0 + allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(dotState(phases%length)) + allocate(dependentState(phases%length)) + + do p = 1, phases%length + if(.not. myPlasticity(p)) cycle phase => phases%get(p) mech => phase%get('mechanics') - if(.not. myPlasticity(p)) cycle - i = i + 1 + i = p associate(prm => param(i), & dot => dotState(i), & stt => state(i), & @@ -496,8 +496,8 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) real(pReal) :: f_unrotated - associate(prm => param(phase_plasticInstance(ph)),& - stt => state(phase_plasticInstance(ph))) + associate(prm => param(ph),& + stt => state(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,me)) & @@ -535,11 +535,11 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) BoltzmannRatio, & ddot_gamma_dtau, & tau - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_slip - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: & + real(pReal), dimension(param(ph)%sum_N_tw) :: & dot_gamma_twin,ddot_gamma_dtau_twin - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: & + real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr,ddot_gamma_dtau_trans real(pReal):: dot_gamma_sb real(pReal), dimension(3,3) :: eigVectors, P_sb @@ -564,7 +564,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) 0, 1, 1 & ],pReal),[ 3,6]) - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,me)) & @@ -573,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip) + call kinetics_slip(Mp,T,ph,me,dot_gamma_sl,ddot_gamma_dtau_slip) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -581,7 +581,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin) + call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin,ddot_gamma_dtau_twin) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -589,7 +589,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinContibution - call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans) + call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_trans) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -653,24 +653,24 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) tau, & sigma_cl, & !< climb stress b_d !< ratio of Burgers vector to stacking fault width - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_rho_dip_formation, & dot_rho_dip_climb, & rho_dip_distance_min, & dot_gamma_sl - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: & + real(pReal), dimension(param(ph)%sum_N_tw) :: & dot_gamma_twin - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: & + real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), & - dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph), & + dot => dotState(ph), dst => dependentState(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tr(1:prm%sum_N_tr,me)) - call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl) + call kinetics_slip(Mp,T,ph,me,dot_gamma_sl) dot%gamma_sl(:,me) = abs(dot_gamma_sl) rho_dip_distance_min = prm%D_a*prm%b_sl @@ -721,10 +721,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin) + call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin) dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char - call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr) + call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr) dot%f_tr(:,me) = f_unrotated*dot_gamma_tr end associate @@ -735,33 +735,33 @@ end subroutine dislotwin_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_dependentState(T,instance,me) +module subroutine dislotwin_dependentState(T,ph,me) integer, intent(in) :: & - instance, & + ph, & me real(pReal), intent(in) :: & T real(pReal) :: & sumf_twin,Gamma,sumf_trans - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(ph)%sum_N_tw) :: & inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin f_over_t_tw - real(pReal), dimension(param(instance)%sum_N_tr) :: & + real(pReal), dimension(param(ph)%sum_N_tr) :: & inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite f_over_t_tr real(pReal), dimension(:), allocatable :: & x0 - associate(prm => param(instance),& - stt => state(instance),& - dst => dependentState(instance)) + associate(prm => param(ph),& + stt => state(ph),& + dst => dependentState(ph)) sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me)) @@ -827,14 +827,14 @@ end subroutine dislotwin_dependentState !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotwin_results(instance,group) +module subroutine plastic_dislotwin_results(ph,group) - integer, intent(in) :: instance + integer, intent(in) :: ph character(len=*), intent(in) :: group integer :: o - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) @@ -882,7 +882,7 @@ end subroutine plastic_dislotwin_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,T,instance,me, & +pure subroutine kinetics_slip(Mp,T,ph,me, & dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) real(pReal), dimension(3,3), intent(in) :: & @@ -890,18 +890,18 @@ pure subroutine kinetics_slip(Mp,T,instance,me, & real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & + ph, & me - real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: & ddot_gamma_dtau_slip, & tau_slip - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & ddot_gamma_dtau - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & stressRatio, & StressRatio_p, & @@ -914,7 +914,7 @@ pure subroutine kinetics_slip(Mp,T,instance,me, & tau_eff !< effective resolved stress integer :: i - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) do i = 1, prm%sum_N_sl tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) @@ -959,7 +959,7 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,& +pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& dot_gamma_twin,ddot_gamma_dtau_twin) real(pReal), dimension(3,3), intent(in) :: & @@ -967,17 +967,17 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,& real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & + ph, & me - real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & + real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & dot_gamma_twin - real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & ddot_gamma_dtau_twin - real, dimension(param(instance)%sum_N_tw) :: & + real, dimension(param(ph)%sum_N_tw) :: & tau, & Ndot0, & stressRatio_r, & @@ -985,7 +985,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,& integer :: i,s1,s2 - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) do i = 1, prm%sum_N_tw tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) @@ -1028,7 +1028,7 @@ end subroutine kinetics_twin ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,me,& +pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& dot_gamma_tr,ddot_gamma_dtau_trans) real(pReal), dimension(3,3), intent(in) :: & @@ -1036,24 +1036,24 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,me,& real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & + ph, & me - real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & + real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: & dot_gamma_tr - real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & ddot_gamma_dtau_trans - real, dimension(param(instance)%sum_N_tr) :: & + real, dimension(param(ph)%sum_N_tr) :: & tau, & Ndot0, & stressRatio_s, & ddot_gamma_dtau integer :: i,s1,s2 - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) do i = 1, prm%sum_N_tr tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index e315a38fb..39a4b3a7c 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -53,7 +53,6 @@ module function plastic_isotropic_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, & i, & Nconstituents, & @@ -68,28 +67,26 @@ module function plastic_isotropic_init() result(myPlasticity) mech, & pl + myPlasticity = plastic_active('isotropic') - Ninstances = count(myPlasticity) - if(Ninstances == 0) return + if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanical:plastic:isotropic init -+>>>' - print'(a,i0)', ' # phses: ',Ninstances; flush(IO_STDOUT) - + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(dotState(Ninstances)) - phases => config_material%get('phase') - i = 0 + allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(dotState(phases%length)) + do p = 1, phases%length + if(.not. myPlasticity(p)) cycle phase => phases%get(p) mech => phase%get('mechanics') - if(.not. myPlasticity(p)) cycle - i = i + 1 + i = p associate(prm => param(i), & dot => dotState(i), & stt => state(i)) @@ -191,7 +188,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) integer :: & k, l, m, n - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph)) Mp_dev = math_deviatoric33(Mp) squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) @@ -221,7 +218,7 @@ end subroutine isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate inelastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) +module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -231,7 +228,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) real(pReal), dimension(3,3), intent(in) :: & Mi !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me real(pReal) :: & @@ -239,7 +236,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) integer :: & k, l, m, n - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) tr=math_trace33(math_spherical33(Mi)) @@ -276,8 +273,8 @@ module subroutine isotropic_dotState(Mp,ph,me) xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), & - dot => dotState(phase_plasticInstance(ph))) + associate(prm => param(ph), stt => state(ph), & + dot => dotState(ph)) if (prm%dilatation) then norm_Mp = sqrt(math_tensordot(Mp,Mp)) @@ -313,14 +310,14 @@ end subroutine isotropic_dotState !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_isotropic_results(instance,group) +module subroutine plastic_isotropic_results(ph,group) - integer, intent(in) :: instance + integer, intent(in) :: ph character(len=*), intent(in) :: group integer :: o - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(ph), stt => state(ph)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('xi') diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 4e94f83ae..285217e1b 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -521,7 +521,7 @@ module function plastic_nonlocal_init() result(myPlasticity) if(.not. myPlasticity(p)) cycle i = i + 1 - Nconstituents = count(material_phaseAt==p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) l = 0 do t = 1,4 do s = 1,param(i)%sum_N_sl From 5a1ca012f8f0abb1d9ca1ff6a4bbdecde24145b3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 14:36:56 +0100 Subject: [PATCH 3/7] more suitable data structure no need to know (ip,el) at the constitutive level --- src/phase_mechanical_eigen.f90 | 1 - src/phase_mechanical_plastic_nonlocal.f90 | 38 +++++++++++++++++++++-- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index ce974a6a6..eb6d4f219 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -186,7 +186,6 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & active = .true. end select plasticType - KinematicsLoop: do k = 1, Nmodels(ph) kinematicsType: select case (model(k,ph)) case (KINEMATICS_thermal_expansion_ID) kinematicsType diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 285217e1b..8f8e6cf3d 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -13,6 +13,12 @@ submodule(phase:plastic) nonlocal IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & geometry_plastic_nonlocal_disable + type :: tGeometry + real(pReal), dimension(:), allocatable :: V_0 + end type tGeometry + + type(tGeometry), dimension(:), allocatable :: geom + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -203,6 +209,10 @@ module function plastic_nonlocal_init() result(myPlasticity) print*, 'Kords, Dissertation RWTH Aachen, 2014' print*, 'http://publications.rwth-aachen.de/record/229993' + + phases => config_material%get('phase') + allocate(geom(phases%length)) + allocate(param(Ninstances)) allocate(state(Ninstances)) allocate(state0(Ninstances)) @@ -210,7 +220,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(deltaState(Ninstances)) allocate(microstructure(Ninstances)) - phases => config_material%get('phase') + i = 0 do p = 1, phases%length phase => phases%get(p) @@ -225,7 +235,7 @@ module function plastic_nonlocal_init() result(myPlasticity) dst => microstructure(i)) pl => mech%get('plasticity') - phase_localPlasticity(p) = .not. pl%contains('nonlocal') + phase_localPlasticity(p) = .not. pl%contains('nonlocal') #if defined (__GFORTRAN__) prm%output = output_asStrings(pl) @@ -409,6 +419,9 @@ module function plastic_nonlocal_init() result(myPlasticity) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) + allocate(geom(p)%V_0(Nconstituents)) + call storeGeometry(p) + plasticState(p)%nonlocal = pl%get_asBool('nonlocal') if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) & call IO_error(212,ext_msg='IPneighborhood does not exist') @@ -1831,4 +1844,25 @@ pure function getRho0(instance,me,ip,el) end function getRho0 + +subroutine storeGeometry(ph) + + integer, intent(in) :: ph + + integer :: ip, el, ce, co + + ce = 0 + do el = 1, size(material_homogenizationMemberAt,2) + do ip = 1, size(material_homogenizationMemberAt,1) + ce = ce + 1 + do co = 1, homogenization_maxNconstituents + if(material_phaseAt2(co,ce) == ph) then + geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = IPvolume(ip,el) + endif + enddo + enddo + enddo + +end subroutine + end submodule nonlocal From 4026881e5a08204c2e117cd1bcb8e9eaf62f2954 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 15:29:10 +0100 Subject: [PATCH 4/7] clean interface still need to get rid of internal converstion to instance and el,ip arguments --- src/phase.f90 | 6 +- src/phase_mechanical.f90 | 6 +- src/phase_mechanical_plastic.f90 | 12 +- src/phase_mechanical_plastic_nonlocal.f90 | 195 +++++++++++----------- 4 files changed, 107 insertions(+), 112 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index f129ba92f..39d5cc340 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -247,9 +247,9 @@ module phase TDot end subroutine phase_thermal_getRate - module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) + module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) integer, intent(in) :: & - instance, & + ph, & i, & e type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & @@ -616,7 +616,7 @@ subroutine crystallite_orientations(co,ip,el) if (plasticState(material_phaseAt(1,el))%nonlocal) & call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticInstance(material_phaseAt(1,el)),ip,el) + material_phaseAt(1,el),ip,el) end subroutine crystallite_orientations diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 394c3d3ba..824c1e257 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -155,8 +155,8 @@ submodule(phase) mechanical character(len=*), intent(in) :: group end subroutine plastic_dislotungsten_results - module subroutine plastic_nonlocal_results(instance,group) - integer, intent(in) :: instance + module subroutine plastic_nonlocal_results(ph,group) + integer, intent(in) :: ph character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results @@ -418,7 +418,7 @@ module subroutine mechanical_results(group,ph) call plastic_dislotungsten_results(ph,group//'plastic/') case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/') + call plastic_nonlocal_results(ph,group//'plastic/') end select diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index a13b5c73f..cf9faf6c9 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -193,9 +193,9 @@ submodule(phase:mechanical) plastic me end subroutine dislotungsten_dependentState - module subroutine nonlocal_dependentState(instance, me, ip, el) + module subroutine nonlocal_dependentState(ph, me, ip, el) integer, intent(in) :: & - instance, & + ph, & me, & ip, & !< current integration point el !< current element number @@ -209,11 +209,11 @@ submodule(phase:mechanical) plastic me end subroutine plastic_kinehardening_deltaState - module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) + module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & - instance, & + ph, & me, & ip, & el @@ -380,7 +380,7 @@ module subroutine plastic_dependentState(co, ip, el) call dislotungsten_dependentState(ph,me) case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_dependentState(instance,me,ip,el) + call nonlocal_dependentState(ph,me,ip,el) end select plasticType @@ -421,7 +421,7 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case (PLASTICITY_NONLOCAL_ID) plasticType - call plastic_nonlocal_deltaState(Mp,instance,me,ip,el) + call plastic_nonlocal_deltaState(Mp,ph,me,ip,el) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case default diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 8f8e6cf3d..1f6fce029 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -12,6 +12,8 @@ submodule(phase:plastic) nonlocal IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & geometry_plastic_nonlocal_disable + use phase, & + ins => phase_plasticInstance type :: tGeometry real(pReal), dimension(:), allocatable :: V_0 @@ -160,7 +162,7 @@ submodule(phase:plastic) nonlocal state, & state0 - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure @@ -510,7 +512,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) end associate - if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i) + if (Nconstituents > 0) call stateInit(ini,p,Nconstituents) plasticState(p)%state0 = plasticState(p)%state !-------------------------------------------------------------------------------------------------- @@ -565,10 +567,10 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dependentState(instance, me, ip, el) +module subroutine nonlocal_dependentState(ph, me, ip, el) integer, intent(in) :: & - instance, & + ph, & me, & ip, & el @@ -602,29 +604,29 @@ module subroutine nonlocal_dependentState(instance, me, ip, el) invConnections real(pReal), dimension(3,nIPneighbors) :: & connection_latticeConf - real(pReal), dimension(2,param(instance)%sum_N_sl) :: & + real(pReal), dimension(2,param(ins(ph))%sum_N_sl) :: & rhoExcess - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & rho_edg_delta, & rho_scr_delta - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & rho, & rho0, & rho_neighbor0 - real(pReal), dimension(param(instance)%sum_N_sl,param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,param(ins(ph))%sum_N_sl) :: & myInteractionMatrix ! corrected slip interaction matrix - real(pReal), dimension(param(instance)%sum_N_sl,nIPneighbors) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,nIPneighbors) :: & rho_edg_delta_neighbor, & rho_scr_delta_neighbor real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: & neighbor_rhoExcess, & ! excess density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point - real(pReal), dimension(3,param(instance)%sum_N_sl,2) :: & + real(pReal), dimension(3,param(ins(ph))%sum_N_sl,2) :: & m ! direction of dislocation motion - associate(prm => param(instance),dst => microstructure(instance), stt => state(instance)) + associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph))) - rho = getRho(instance,me,ip,el) + rho = getRho(ph,me,ip,el) stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) @@ -652,9 +654,8 @@ module subroutine nonlocal_dependentState(instance, me, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - rho0 = getRho0(instance,me,ip,el) + rho0 = getRho0(ph,me,ip,el) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then - ph = material_phaseAt(1,el) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me)) @@ -675,11 +676,11 @@ module subroutine nonlocal_dependentState(instance, me, ip, el) neighbor_ip = IPneighborhood(2,n,ip,el) no = material_phasememberAt(1,neighbor_ip,neighbor_el) if (neighbor_el > 0 .and. neighbor_ip > 0) then - neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el)) - if (neighbor_instance == instance) then + neighbor_instance = ins(material_phaseAt(1,neighbor_el)) + if (neighbor_instance == ins(ph)) then nRealNeighbors = nRealNeighbors + 1.0_pReal - rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el) + rho_neighbor0 = getRho0(ph,no,neighbor_ip,neighbor_el) rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) @@ -795,25 +796,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & l, & t, & !< dislocation type s !< index of my current slip system - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: & rhoSgl !< single dislocation densities (including blocked) - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & rho - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - associate(prm => param(phase_plasticInstance(ph)),dst=>microstructure(phase_plasticInstance(ph)),& - stt=>state(phase_plasticInstance(ph))) + associate(prm => param(ins(ph)),dst=>microstructure(ins(ph)),& + stt=>state(ins(ph))) ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(phase_plasticInstance(ph),me,ip,el) + rho = getRho(ph,me,ip,el) rhoSgl = rho(:,sgl) do s = 1,ns @@ -833,7 +834,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & ! edges call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & - tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticInstance(ph)) + tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, ph) v(:,2) = v(:,1) dv_dtau(:,2) = dv_dtau(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1) @@ -846,7 +847,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & else do t = 3,4 call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & - tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticInstance(ph)) + tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, ph) enddo endif @@ -879,12 +880,12 @@ end subroutine nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) +module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress integer, intent(in) :: & - instance, & ! current instance of this plasticity + ph, & me, & !< offset ip, & el @@ -895,31 +896,29 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) c, & ! character of dislocation t, & ! type of dislocation s ! index of my current slip system - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & deltaRhoRemobilization, & ! density increment by remobilization deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & rho ! current dislocation densities - real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & v ! dislocation glide velocity - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & tau ! current resolved shear stress - real(pReal), dimension(param(instance)%sum_N_sl,2) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,2) :: & rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - ph = material_phaseAt(1,el) - - associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance)) - ns = prm%sum_N_sl + associate(prm => param(ins(ph)),dst => microstructure(ins(ph)),del => deltaState(ins(ph))) + ns = prm%sum_N_sl !*** shortcut to state variables - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) - forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) + forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ins(ph)),me) - rho = getRho(instance,me,ip,el) + rho = getRho(ph,me,ip,el) rhoDip = rho(:,dip) !**************************************************************************** @@ -964,7 +963,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9) - forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),me) = dUpper(s,c) + forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ins(ph)),me) = dUpper(s,c) plasticState(ph)%deltaState(:,me) = 0.0_pReal del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) @@ -1005,7 +1004,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & c, & !< character of dislocation t, & !< type of dislocation s !< index of my current slip system - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution @@ -1013,17 +1012,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & gdot !< shear rates - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,2) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws @@ -1035,22 +1034,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & return endif - associate(prm => param(phase_plasticInstance(ph)), & - dst => microstructure(phase_plasticInstance(ph)), & - dot => dotState(phase_plasticInstance(ph)), & - stt => state(phase_plasticInstance(ph))) + associate(prm => param(ins(ph)), & + dst => microstructure(ins(ph)), & + dot => dotState(ins(ph)), & + stt => state(ins(ph))) ns = prm%sum_N_sl tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(phase_plasticInstance(ph),me,ip,el) + rho = getRho(ph,me,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - rho0 = getRho0(phase_plasticInstance(ph),me,ip,el) + rho0 = getRho0(ph,me,ip,el) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticInstance(ph)),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) #ifdef DEBUG @@ -1099,7 +1098,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) endif isBCC - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticInstance(ph)),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ins(ph)),me) !**************************************************************************** @@ -1155,7 +1154,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(timestep, phase_plasticInstance(ph),me,ip,el) & + rhoDot = rhoDotFlux(timestep, ph,me,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1184,12 +1183,12 @@ end subroutine nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -function rhoDotFlux(timestep,instance,me,ip,el) +function rhoDotFlux(timestep,ph,me,ip,el) real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & - instance, & + ph, & me, & ip, & !< current integration point el !< current element number @@ -1212,20 +1211,20 @@ function rhoDotFlux(timestep,instance,me,ip,el) np,& !< neighbor phase shortcut topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDotFlux !< density evolution by flux - real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) - real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates - real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(3,param(ins(ph))%sum_N_sl,4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: & my_F, & !< my total deformation gradient @@ -1243,26 +1242,25 @@ function rhoDotFlux(timestep,instance,me,ip,el) transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength !< dislocation line length leaving the current interface - ph = material_phaseAt(1,el) - associate(prm => param(instance), & - dst => microstructure(instance), & - dot => dotState(instance), & - stt => state(instance)) + associate(prm => param(ins(ph)), & + dst => microstructure(ins(ph)), & + dot => dotState(ins(ph)), & + stt => state(ins(ph))) ns = prm%sum_N_sl gdot = 0.0_pReal - rho = getRho(instance,me,ip,el) + rho = getRho(ph,me,ip,el) rhoSgl = rho(:,sgl) - rho0 = getRho0(instance,me,ip,el) + rho0 = getRho0(ph,me,ip,el) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) !ToDo: MD: I think we should use state0 here gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ins(ph)),me) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) @@ -1314,7 +1312,7 @@ function rhoDotFlux(timestep,instance,me,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el)) + neighbor_instance = ins(material_phaseAt(1,neighbor_el)) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) Favg = 0.5_pReal * (my_F + neighbor_F) @@ -1421,12 +1419,12 @@ end function rhoDotFlux ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) +module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & - instance, & + ph, & i, & e @@ -1439,19 +1437,17 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) ns, & ! number of active slip systems s1, & ! slip system index (me) s2 ! slip system index (my neighbor) - real(pReal), dimension(2,param(instance)%sum_N_sl,param(instance)%sum_N_sl,nIPneighbors) :: & + real(pReal), dimension(2,param(ins(ph))%sum_N_sl,param(ins(ph))%sum_N_sl,nIPneighbors) :: & my_compatibility ! my_compatibility for current element and ip real(pReal) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(param(instance)%sum_N_sl) :: & + logical, dimension(param(ins(ph))%sum_N_sl) :: & belowThreshold type(rotation) :: mis - ph = material_phaseAt(1,e) - - associate(prm => param(instance)) + associate(prm => param(ins(ph))) ns = prm%sum_N_sl !*** start out fully compatible @@ -1537,14 +1533,14 @@ end subroutine plastic_nonlocal_updateCompatibility !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_results(instance,group) +module subroutine plastic_nonlocal_results(ph,group) - integer, intent(in) :: instance + integer, intent(in) :: ph character(len=*),intent(in) :: group integer :: o - associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) + associate(prm => param(ins(ph)),dst => microstructure(ins(ph)),stt=>state(ins(ph))) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case('rho_u_ed_pos') @@ -1608,14 +1604,13 @@ end subroutine plastic_nonlocal_results !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- -subroutine stateInit(ini,phase,Nconstituents,instance) +subroutine stateInit(ini,phase,Nconstituents) type(tInitialParameters) :: & ini integer,intent(in) :: & phase, & - Nconstituents, & - instance + Nconstituents integer :: & e, & i, & @@ -1636,7 +1631,7 @@ subroutine stateInit(ini,phase,Nconstituents,instance) volume - associate(stt => state(instance)) + associate(stt => state(ins(phase))) if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume do e = 1,discretization_Nelems @@ -1684,18 +1679,18 @@ end subroutine stateInit !-------------------------------------------------------------------------------------------------- !> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance) +pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, ph) integer, intent(in) :: & c, & !< dislocation character (1:edge, 2:screw) - instance + ph real(pReal), intent(in) :: & Temperature !< temperature - real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl), intent(in) :: & tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) tauThreshold !< threshold shear stress - real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & + real(pReal), dimension(param(ins(ph))%sum_N_sl), intent(out) :: & v, & !< velocity dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) @@ -1726,7 +1721,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem criticalStress_S, & !< maximum obstacle strength mobility !< dislocation mobility - associate(prm => param(instance)) + associate(prm => param(ins(ph))) ns = prm%sum_N_sl v = 0.0_pReal dv_dtau = 0.0_pReal @@ -1799,14 +1794,14 @@ end subroutine kinetics !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho(instance,me,ip,el) +pure function getRho(ph,me,ip,el) - integer, intent(in) :: instance, me,ip,el - real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho + integer, intent(in) :: ph, me,ip,el + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho - associate(prm => param(instance)) + associate(prm => param(ins(ph))) - getRho = reshape(state(instance)%rho(:,me),[prm%sum_N_sl,10]) + getRho = reshape(state(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) @@ -1824,14 +1819,14 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho0(instance,me,ip,el) +pure function getRho0(ph,me,ip,el) - integer, intent(in) :: instance, me,ip,el - real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0 + integer, intent(in) :: ph, me,ip,el + real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0 - associate(prm => param(instance)) + associate(prm => param(ins(ph))) - getRho0 = reshape(state0(instance)%rho(:,me),[prm%sum_N_sl,10]) + getRho0 = reshape(state0(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) From f46d212e47bdc3e675092dc236ee99a0dd3e3b66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 17:29:23 +0100 Subject: [PATCH 5/7] simplified access pattern --- src/phase_mechanical.f90 | 23 +++--- src/phase_mechanical_plastic.f90 | 34 ++++----- src/phase_mechanical_plastic_nonlocal.f90 | 85 ++++++++++------------- 3 files changed, 56 insertions(+), 86 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 824c1e257..381820dcc 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -85,11 +85,8 @@ submodule(phase) mechanical logical :: broken end function plastic_dotState - module function plastic_deltaState(co, ip, el, ph, me) result(broken) + module function plastic_deltaState(ph, me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & me logical :: & @@ -114,11 +111,9 @@ submodule(phase) mechanical module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) + S, Fi, 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, & !< 2nd Piola-Kirchhoff stress Fi !< intermediate deformation gradient @@ -538,7 +533,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) Fe, Fi_new, co, ip, el) call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, co, ip, el) + S, Fi_new, ph,me) !* update current residuum and check for convergence of loop atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -702,7 +697,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & - r(1:sizeDotState) if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) exit iteration endif @@ -765,7 +760,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -807,7 +802,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -956,7 +951,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D if(broken) return - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -1329,7 +1324,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & phase_mechanical_S(ph)%data(1:3,1:3,me), & - phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el) + phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index cf9faf6c9..c8e9b6271 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -104,7 +104,7 @@ submodule(phase:mechanical) plastic end subroutine dislotungsten_LpAndItsTangent module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,me,ip,el) + Mp,Temperature,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -116,9 +116,7 @@ submodule(phase:mechanical) plastic Temperature integer, intent(in) :: & ph, & - me, & - ip, & !< current integration point - el !< current element number + me end subroutine nonlocal_LpAndItsTangent @@ -209,14 +207,12 @@ submodule(phase:mechanical) plastic me end subroutine plastic_kinehardening_deltaState - module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) + module subroutine plastic_nonlocal_deltaState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & ph, & - me, & - ip, & - el + me end subroutine plastic_nonlocal_deltaState end interface @@ -244,11 +240,9 @@ end subroutine plastic_init ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) + S, Fi, 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, & !< 2nd Piola-Kirchhoff stress Fi !< intermediate deformation gradient @@ -263,14 +257,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & real(pReal), dimension(3,3) :: & Mp !< Mandel stress work conjugate with Lp integer :: & - i, j, me, ph + i, j Mp = matmul(matmul(transpose(Fi),Fi),S) - me = material_phasememberAt(co,ip,el) - ph = material_phaseAt(co,el) - plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + + plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_NONE_ID) plasticType Lp = 0.0_pReal @@ -286,7 +279,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) case (PLASTICITY_DISLOTWIN_ID) plasticType call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) @@ -391,12 +384,9 @@ end subroutine plastic_dependentState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -module function plastic_deltaState(co, ip, el, ph, me) result(broken) +module function plastic_deltaState(ph, me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & me logical :: & @@ -421,7 +411,7 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case (PLASTICITY_NONLOCAL_ID) plasticType - call plastic_nonlocal_deltaState(Mp,ph,me,ip,el) + call plastic_nonlocal_deltaState(Mp,ph,me) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case default diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 1f6fce029..b579a3164 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -512,7 +512,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) end associate - if (Nconstituents > 0) call stateInit(ini,p,Nconstituents) + if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i) plasticState(p)%state0 = plasticState(p)%state !-------------------------------------------------------------------------------------------------- @@ -576,7 +576,6 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) el integer :: & - ph, & no, & !< neighbor offset neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point @@ -626,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph))) - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) @@ -654,7 +653,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - rho0 = getRho0(ph,me,ip,el) + rho0 = getRho0(ph,me) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me)) @@ -665,7 +664,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) rhoExcess(1,:) = rho_edg_delta rhoExcess(2,:) = rho_scr_delta - FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal) + FVsize = geom(ph)%V_0(me) ** (1.0_pReal/3.0_pReal) !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities @@ -680,7 +679,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) if (neighbor_instance == ins(ph)) then nRealNeighbors = nRealNeighbors + 1.0_pReal - rho_neighbor0 = getRho0(ph,no,neighbor_ip,neighbor_el) + rho_neighbor0 = getRho0(ph,no) rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) @@ -772,16 +771,14 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,me,ip,el) + Mp,Temperature,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp integer, intent(in) :: & ph, & - me, & - ip, & !< current integration point - el !< current element number + me real(pReal), intent(in) :: & Temperature !< temperature @@ -814,7 +811,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoSgl = rho(:,sgl) do s = 1,ns @@ -880,18 +877,15 @@ end subroutine nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) +module subroutine plastic_nonlocal_deltaState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress integer, intent(in) :: & ph, & - me, & !< offset - ip, & - el + me integer :: & - ph, & !< phase ns, & ! short notation for the total number of active slip systems c, & ! character of dislocation t, & ! type of dislocation @@ -918,7 +912,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ins(ph)),me) - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoDip = rho(:,dip) !**************************************************************************** @@ -968,15 +962,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) plasticState(ph)%deltaState(:,me) = 0.0_pReal del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) -#ifdef DEBUG - if (debugConstitutive%extensive & - .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& - .or. .not. debugConstitutive%selective)) then - print'(a,/,8(12x,12(e12.5,1x),/))', '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8) - print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress - endif -#endif - end associate end subroutine plastic_nonlocal_deltaState @@ -1043,10 +1028,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - rho0 = getRho0(ph,me,ip,el) + rho0 = getRho0(ph,me) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) @@ -1194,7 +1179,6 @@ function rhoDotFlux(timestep,ph,me,ip,el) el !< current element number integer :: & - ph, & neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation @@ -1251,9 +1235,9 @@ function rhoDotFlux(timestep,ph,me,ip,el) gdot = 0.0_pReal - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoSgl = rho(:,sgl) - rho0 = getRho0(ph,me,ip,el) + rho0 = getRho0(ph,me) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) !ToDo: MD: I think we should use state0 here @@ -1432,7 +1416,6 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor - ph, & neighbor_phase, & ns, & ! number of active slip systems s1, & ! slip system index (me) @@ -1604,7 +1587,7 @@ end subroutine plastic_nonlocal_results !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- -subroutine stateInit(ini,phase,Nconstituents) +subroutine stateInit(ini,phase,Nconstituents,i) type(tInitialParameters) :: & ini @@ -1631,7 +1614,7 @@ subroutine stateInit(ini,phase,Nconstituents) volume - associate(stt => state(ins(phase))) + associate(stt => state(i)) if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume do e = 1,discretization_Nelems @@ -1794,21 +1777,22 @@ end subroutine kinetics !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho(ph,me,ip,el) +pure function getRho(ph,me) - integer, intent(in) :: ph, me,ip,el + integer, intent(in) :: ph, me real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho + associate(prm => param(ins(ph))) - getRho = reshape(state(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) + getRho = reshape(state(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) - ! ensure positive densities (not for imm, they have a sign) - getRho(:,mob) = max(getRho(:,mob),0.0_pReal) - getRho(:,dip) = max(getRho(:,dip),0.0_pReal) + ! ensure positive densities (not for imm, they have a sign) + getRho(:,mob) = max(getRho(:,mob),0.0_pReal) + getRho(:,dip) = max(getRho(:,dip),0.0_pReal) - where(abs(getRho) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & - getRho = 0.0_pReal + where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + getRho = 0.0_pReal end associate @@ -1819,21 +1803,22 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho0(ph,me,ip,el) +pure function getRho0(ph,me) - integer, intent(in) :: ph, me,ip,el + integer, intent(in) :: ph, me real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0 + associate(prm => param(ins(ph))) - getRho0 = reshape(state0(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) + getRho0 = reshape(state0(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) - ! ensure positive densities (not for imm, they have a sign) - getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) - getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) + ! ensure positive densities (not for imm, they have a sign) + getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) + getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) - where(abs(getRho0) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & - getRho0 = 0.0_pReal + where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + getRho0 = 0.0_pReal end associate From 341e8ddd6a7d60631d0a212d3c30009fa75e7be6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 18:30:57 +0100 Subject: [PATCH 6/7] storing per instance does not add any value --- src/phase.f90 | 3 +- src/phase_mechanical.f90 | 2 - src/phase_mechanical_plastic.f90 | 5 +- src/phase_mechanical_plastic_nonlocal.f90 | 166 +++++++++++----------- 4 files changed, 82 insertions(+), 94 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 39d5cc340..a36ddc157 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -59,8 +59,7 @@ module phase integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) phase_elasticityInstance, & - phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase - phase_plasticInstance + phase_NstiffnessDegradations logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) phase_localPlasticity !< flags phases with local constitutive law diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 381820dcc..57d4ea5ec 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -298,14 +298,12 @@ module subroutine mechanical_init(phases) ! initialize plasticity allocate(plasticState(phases%length)) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) - allocate(phase_plasticInstance(phases%length),source = 0) allocate(phase_localPlasticity(phases%length), source=.true.) call plastic_init() do ph = 1, phases%length phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) - phase_plasticInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph)) enddo num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index c8e9b6271..136f0884c 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -357,12 +357,11 @@ module subroutine plastic_dependentState(co, ip, el) integer :: & ph, & - instance, me + me ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - instance = phase_plasticInstance(ph) plasticType: select case (phase_plasticity(material_phaseAt(co,el))) @@ -395,14 +394,12 @@ module function plastic_deltaState(ph, me) result(broken) real(pReal), dimension(3,3) :: & Mp integer :: & - instance, & myOffset, & mySize Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),& phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me)) - instance = phase_plasticInstance(ph) plasticType: select case (phase_plasticity(ph)) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index b579a3164..0212baee8 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -12,8 +12,6 @@ submodule(phase:plastic) nonlocal IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & geometry_plastic_nonlocal_disable - use phase, & - ins => phase_plasticInstance type :: tGeometry real(pReal), dimension(:), allocatable :: V_0 @@ -215,20 +213,19 @@ module function plastic_nonlocal_init() result(myPlasticity) phases => config_material%get('phase') allocate(geom(phases%length)) - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(state0(Ninstances)) - allocate(dotState(Ninstances)) - allocate(deltaState(Ninstances)) - allocate(microstructure(Ninstances)) + allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(state0(phases%length)) + allocate(dotState(phases%length)) + allocate(deltaState(phases%length)) + allocate(microstructure(phases%length)) - - i = 0 do p = 1, phases%length + if(.not. myPlasticity(p)) cycle phase => phases%get(p) mech => phase%get('mechanics') - if(.not. myPlasticity(p)) cycle - i = i + 1 + + i = p associate(prm => param(i), & dot => dotState(i), & stt => state(i), & @@ -512,7 +509,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) end associate - if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i) + if (Nconstituents > 0) call stateInit(ini,p,Nconstituents) plasticState(p)%state0 = plasticState(p)%state !-------------------------------------------------------------------------------------------------- @@ -525,16 +522,15 @@ module function plastic_nonlocal_init() result(myPlasticity) discretization_nIPs,discretization_Nelems), source=0.0_pReal) ! BEGIN DEPRECATED---------------------------------------------------------------------------------- - allocate(iRhoU(maxval(param%sum_N_sl),4,Ninstances), source=0) - allocate(iV(maxval(param%sum_N_sl),4,Ninstances), source=0) - allocate(iD(maxval(param%sum_N_sl),2,Ninstances), source=0) + allocate(iRhoU(maxval(param%sum_N_sl),4,phases%length), source=0) + allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0) + allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0) - i = 0 do p = 1, phases%length phase => phases%get(p) if(.not. myPlasticity(p)) cycle - i = i + 1 + i = p Nconstituents = count(material_phaseAt2 == p) l = 0 @@ -579,7 +575,6 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) no, & !< neighbor offset neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point - neighbor_instance, & ! instance of this plasticity of neighboring material point c, & ! index of dilsocation character (edge, screw) s, & ! slip system index dir, & @@ -603,27 +598,27 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) invConnections real(pReal), dimension(3,nIPneighbors) :: & connection_latticeConf - real(pReal), dimension(2,param(ins(ph))%sum_N_sl) :: & + real(pReal), dimension(2,param(ph)%sum_N_sl) :: & rhoExcess - real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & rho_edg_delta, & rho_scr_delta - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ph)%sum_N_sl,10) :: & rho, & rho0, & rho_neighbor0 - real(pReal), dimension(param(ins(ph))%sum_N_sl,param(ins(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: & myInteractionMatrix ! corrected slip interaction matrix - real(pReal), dimension(param(ins(ph))%sum_N_sl,nIPneighbors) :: & + real(pReal), dimension(param(ph)%sum_N_sl,nIPneighbors) :: & rho_edg_delta_neighbor, & rho_scr_delta_neighbor real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: & neighbor_rhoExcess, & ! excess density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point - real(pReal), dimension(3,param(ins(ph))%sum_N_sl,2) :: & + real(pReal), dimension(3,param(ph)%sum_N_sl,2) :: & m ! direction of dislocation motion - associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph))) + associate(prm => param(ph),dst => microstructure(ph), stt => state(ph)) rho = getRho(ph,me) @@ -675,8 +670,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) neighbor_ip = IPneighborhood(2,n,ip,el) no = material_phasememberAt(1,neighbor_ip,neighbor_el) if (neighbor_el > 0 .and. neighbor_ip > 0) then - neighbor_instance = ins(material_phaseAt(1,neighbor_el)) - if (neighbor_instance == ins(ph)) then + if (material_phaseAt(1,neighbor_el) == ph) then nRealNeighbors = nRealNeighbors + 1.0_pReal rho_neighbor0 = getRho0(ph,no) @@ -793,21 +787,21 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & l, & t, & !< dislocation type s !< index of my current slip system - real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: & + real(pReal), dimension(param(ph)%sum_N_sl,8) :: & rhoSgl !< single dislocation densities (including blocked) - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ph)%sum_N_sl,10) :: & rho - real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & + real(pReal), dimension(param(ph)%sum_N_sl,4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress - real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - associate(prm => param(ins(ph)),dst=>microstructure(ins(ph)),& - stt=>state(ins(ph))) + associate(prm => param(ph),dst=>microstructure(ph),& + stt=>state(ph)) ns = prm%sum_N_sl !*** shortcut to state variables @@ -890,27 +884,27 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me) c, & ! character of dislocation t, & ! type of dislocation s ! index of my current slip system - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ph)%sum_N_sl,10) :: & deltaRhoRemobilization, & ! density increment by remobilization deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ph)%sum_N_sl,10) :: & rho ! current dislocation densities - real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & + real(pReal), dimension(param(ph)%sum_N_sl,4) :: & v ! dislocation glide velocity - real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & tau ! current resolved shear stress - real(pReal), dimension(param(ins(ph))%sum_N_sl,2) :: & + real(pReal), dimension(param(ph)%sum_N_sl,2) :: & rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - associate(prm => param(ins(ph)),dst => microstructure(ins(ph)),del => deltaState(ins(ph))) + associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph)) ns = prm%sum_N_sl !*** shortcut to state variables - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) - forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ins(ph)),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) + forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),me) rho = getRho(ph,me) rhoDip = rho(:,dip) @@ -957,7 +951,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9) - forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ins(ph)),me) = dUpper(s,c) + forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),me) = dUpper(s,c) plasticState(ph)%deltaState(:,me) = 0.0_pReal del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) @@ -989,7 +983,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & c, & !< character of dislocation t, & !< type of dislocation s !< index of my current slip system - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ph)%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution @@ -997,17 +991,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation - real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: & + real(pReal), dimension(param(ph)%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) - real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & + real(pReal), dimension(param(ph)%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & gdot !< shear rates - real(pReal), dimension(param(ins(ph))%sum_N_sl) :: & + real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - real(pReal), dimension(param(ins(ph))%sum_N_sl,2) :: & + real(pReal), dimension(param(ph)%sum_N_sl,2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws @@ -1019,10 +1013,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & return endif - associate(prm => param(ins(ph)), & - dst => microstructure(ins(ph)), & - dot => dotState(ins(ph)), & - stt => state(ins(ph))) + associate(prm => param(ph), & + dst => microstructure(ph), & + dot => dotState(ph), & + stt => state(ph)) ns = prm%sum_N_sl tau = 0.0_pReal @@ -1034,7 +1028,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rho0 = getRho0(ph,me) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) #ifdef DEBUG @@ -1083,7 +1077,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) endif isBCC - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ins(ph)),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me) !**************************************************************************** @@ -1179,7 +1173,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) el !< current element number integer :: & - neighbor_instance, & !< instance of my neighbor's plasticity + neighbor_ph, & !< phase of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation n, & !< index of my current neighbor @@ -1195,20 +1189,20 @@ function rhoDotFlux(timestep,ph,me,ip,el) np,& !< neighbor phase shortcut topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(ph)%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDotFlux !< density evolution by flux - real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: & + real(pReal), dimension(param(ph)%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) - real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: & + real(pReal), dimension(param(ph)%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates - real(pReal), dimension(3,param(ins(ph))%sum_N_sl,4) :: & + real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: & my_F, & !< my total deformation gradient @@ -1227,10 +1221,10 @@ function rhoDotFlux(timestep,ph,me,ip,el) lineLength !< dislocation line length leaving the current interface - associate(prm => param(ins(ph)), & - dst => microstructure(ins(ph)), & - dot => dotState(ins(ph)), & - stt => state(ins(ph))) + associate(prm => param(ph), & + dst => microstructure(ph), & + dot => dotState(ph), & + stt => state(ph)) ns = prm%sum_N_sl gdot = 0.0_pReal @@ -1240,11 +1234,11 @@ function rhoDotFlux(timestep,ph,me,ip,el) rho0 = getRho0(ph,me) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) !ToDo: MD: I think we should use state0 here gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ins(ph)),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) @@ -1296,7 +1290,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_instance = ins(material_phaseAt(1,neighbor_el)) + neighbor_ph = material_phaseAt(1,neighbor_el) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) Favg = 0.5_pReal * (my_F + neighbor_F) @@ -1319,8 +1313,8 @@ function rhoDotFlux(timestep,ph,me,ip,el) any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then forall (s = 1:ns, t = 1:4) - neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) - neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal) + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no) + neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal) endforall where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min & @@ -1420,17 +1414,17 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) ns, & ! number of active slip systems s1, & ! slip system index (me) s2 ! slip system index (my neighbor) - real(pReal), dimension(2,param(ins(ph))%sum_N_sl,param(ins(ph))%sum_N_sl,nIPneighbors) :: & + real(pReal), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: & my_compatibility ! my_compatibility for current element and ip real(pReal) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(param(ins(ph))%sum_N_sl) :: & + logical, dimension(param(ph)%sum_N_sl) :: & belowThreshold type(rotation) :: mis - associate(prm => param(ins(ph))) + associate(prm => param(ph)) ns = prm%sum_N_sl !*** start out fully compatible @@ -1523,7 +1517,7 @@ module subroutine plastic_nonlocal_results(ph,group) integer :: o - associate(prm => param(ins(ph)),dst => microstructure(ins(ph)),stt=>state(ins(ph))) + associate(prm => param(ph),dst => microstructure(ph),stt=>state(ph)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case('rho_u_ed_pos') @@ -1587,7 +1581,7 @@ end subroutine plastic_nonlocal_results !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- -subroutine stateInit(ini,phase,Nconstituents,i) +subroutine stateInit(ini,phase,Nconstituents) type(tInitialParameters) :: & ini @@ -1595,8 +1589,8 @@ subroutine stateInit(ini,phase,Nconstituents,i) phase, & Nconstituents integer :: & - e, & i, & + e, & f, & from, & upto, & @@ -1614,7 +1608,7 @@ subroutine stateInit(ini,phase,Nconstituents,i) volume - associate(stt => state(i)) + associate(stt => state(phase)) if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume do e = 1,discretization_Nelems @@ -1643,8 +1637,8 @@ subroutine stateInit(ini,phase,Nconstituents,i) do s = from,upto noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] - stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) - stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) + stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) + stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) enddo @@ -1669,11 +1663,11 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem ph real(pReal), intent(in) :: & Temperature !< temperature - real(pReal), dimension(param(ins(ph))%sum_N_sl), intent(in) :: & + real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) tauThreshold !< threshold shear stress - real(pReal), dimension(param(ins(ph))%sum_N_sl), intent(out) :: & + real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & v, & !< velocity dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) @@ -1704,7 +1698,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem criticalStress_S, & !< maximum obstacle strength mobility !< dislocation mobility - associate(prm => param(ins(ph))) + associate(prm => param(ph)) ns = prm%sum_N_sl v = 0.0_pReal dv_dtau = 0.0_pReal @@ -1780,12 +1774,12 @@ end subroutine kinetics pure function getRho(ph,me) integer, intent(in) :: ph, me - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho + real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho - associate(prm => param(ins(ph))) + associate(prm => param(ph)) - getRho = reshape(state(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) + getRho = reshape(state(ph)%rho(:,me),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) @@ -1806,12 +1800,12 @@ end function getRho pure function getRho0(ph,me) integer, intent(in) :: ph, me - real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0 + real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho0 - associate(prm => param(ins(ph))) + associate(prm => param(ph)) - getRho0 = reshape(state0(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) + getRho0 = reshape(state0(ph)%rho(:,me),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) From 53bab41b4762eb3c8da911d3f409da800eb9ed16 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Feb 2021 16:06:09 +0100 Subject: [PATCH 7/7] consistent name 'ph' and cleaning --- ...phase_mechanical_plastic_dislotungsten.f90 | 46 ++--- src/phase_mechanical_plastic_dislotwin.f90 | 75 ++++--- src/phase_mechanical_plastic_isotropic.f90 | 54 ++--- ...phase_mechanical_plastic_kinehardening.f90 | 73 +++---- src/phase_mechanical_plastic_none.f90 | 10 +- src/phase_mechanical_plastic_nonlocal.f90 | 189 +++++++++--------- ...phase_mechanical_plastic_phenopowerlaw.f90 | 57 +++--- 7 files changed, 235 insertions(+), 269 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 81c4605ac..9c8baf0a5 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -78,7 +78,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - p, i, & + ph, i, & Nconstituents, & sizeState, sizeDotState, & startIndex, endIndex @@ -114,15 +114,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) allocate(dependentState(phases%length)) - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - phase => phases%get(p) + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle + + associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) + + phase => phases%get(ph) mech => phase%get('mechanics') - i = p - associate(prm => param(i), & - dot => dotState(i), & - stt => state(i), & - dst => dependentState(i)) pl => mech%get('plasticity') #if defined (__GFORTRAN__) @@ -132,7 +130,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) #endif ! This data is read in already in lattice - prm%mu = lattice_mu(p) + prm%mu = lattice_mu(ph) !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -222,41 +220,41 @@ module function plastic_dislotungsten_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == p) + Nconstituents = count(material_phaseAt2 == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:) + stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob = spread(rho_mob_0,2,Nconstituents) - dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) - if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' + dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:) + stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip = spread(rho_dip_0,2,Nconstituents) - dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) + dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_sl => plasticState(p)%state(startIndex:endIndex,:) - dot%gamma_sl => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = 1.0e-2_pReal + stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) + dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal) - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally end associate diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 3b33196d9..0c52176ce 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -126,7 +126,7 @@ module function plastic_dislotwin_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - p, i, & + ph, i, & Nconstituents, & sizeState, sizeDotState, & startIndex, endIndex @@ -167,15 +167,13 @@ module function plastic_dislotwin_init() result(myPlasticity) allocate(dependentState(phases%length)) - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - phase => phases%get(p) + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle + + associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) + + phase => phases%get(ph) mech => phase%get('mechanics') - i = p - associate(prm => param(i), & - dot => dotState(i), & - stt => state(i), & - dst => dependentState(i)) pl => mech%get('plasticity') #if defined (__GFORTRAN__) @@ -185,9 +183,9 @@ module function plastic_dislotwin_init() result(myPlasticity) #endif ! This data is read in already in lattice - prm%mu = lattice_mu(p) - prm%nu = lattice_nu(p) - prm%C66 = lattice_C66(1:6,1:6,p) + prm%mu = lattice_mu(ph) + prm%nu = lattice_nu(ph) + prm%C66 = lattice_C66(1:6,1:6,ph) !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -204,8 +202,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),& phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == lattice_FCC_ID) & - .and. (N_sl(1) == 12) + prm%fccTwinTransNucleation = lattice_structure(ph) == lattice_FCC_ID .and. (N_sl(1) == 12) if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl)) @@ -234,7 +231,7 @@ module function plastic_dislotwin_init() result(myPlasticity) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 prm%omega = pl%get_asFloat('omega', defaultVal = 1000.0_pReal) & - * merge(12.0_pReal,8.0_pReal,any(lattice_structure(p) == [lattice_FCC_ID,lattice_HEX_ID])) + * merge(12.0_pReal,8.0_pReal,any(lattice_structure(ph) == [lattice_FCC_ID,lattice_HEX_ID])) ! expand: family => system rho_mob_0 = math_expand(rho_mob_0, N_sl) @@ -342,7 +339,7 @@ module function plastic_dislotwin_init() result(myPlasticity) pl%get_asFloat('a_cI', defaultVal=0.0_pReal), & pl%get_asFloat('a_cF', defaultVal=0.0_pReal)) - if (lattice_structure(p) /= lattice_FCC_ID) then + if (lattice_structure(ph) /= lattice_FCC_ID) then prm%dot_N_0_tr = pl%get_asFloats('dot_N_0_tr') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr) endif @@ -357,7 +354,7 @@ module function plastic_dislotwin_init() result(myPlasticity) if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr' - if (lattice_structure(p) /= lattice_FCC_ID) then + if (lattice_structure(ph) /= lattice_FCC_ID) then if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr' endif else transActive @@ -408,53 +405,53 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == p) + Nconstituents = count(material_phaseAt2 == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and atol startIndex = 1 endIndex = prm%sum_N_sl - stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob= spread(rho_mob_0,2,Nconstituents) - dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) - if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' + dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip= spread(rho_dip_0,2,Nconstituents) - dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) + dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) - dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = 1.0e-2_pReal + stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) + dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:) - dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal) - if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin' + stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) + dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal) + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tr - stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:) - dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) - if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' + stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:) + dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal) @@ -469,7 +466,7 @@ module function plastic_dislotwin_init() result(myPlasticity) allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally end associate diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 39a4b3a7c..245f4293a 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -22,8 +22,6 @@ submodule(phase:plastic) isotropic c_4, & c_3, & c_2 - integer :: & - of_debug = 0 logical :: & dilatation character(len=pStringLen), allocatable, dimension(:) :: & @@ -53,8 +51,7 @@ module function plastic_isotropic_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - p, & - i, & + ph, & Nconstituents, & sizeState, sizeDotState real(pReal) :: & @@ -82,16 +79,14 @@ module function plastic_isotropic_init() result(myPlasticity) allocate(state(phases%length)) allocate(dotState(phases%length)) - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - phase => phases%get(p) - mech => phase%get('mechanics') - i = p - associate(prm => param(i), & - dot => dotState(i), & - stt => state(i)) - pl => mech%get('plasticity') + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle + associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) + + phase => phases%get(ph) + mech => phase%get('mechanics') + pl => mech%get('plasticity') #if defined (__GFORTRAN__) prm%output = output_asStrings(pl) @@ -99,11 +94,6 @@ module function plastic_isotropic_init() result(myPlasticity) prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray) #endif -#ifdef DEBUG - if (p==material_phaseAt(debugConstitutive%grain,debugConstitutive%element)) & - prm%of_debug = material_phasememberAt(debugConstitutive%grain,debugConstitutive%ip,debugConstitutive%element) -#endif - xi_0 = pl%get_asFloat('xi_0') prm%xi_inf = pl%get_asFloat('xi_inf') prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0') @@ -129,28 +119,28 @@ module function plastic_isotropic_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == p) + Nconstituents = count(material_phaseAt2 == ph) sizeDotState = size(['xi ','gamma']) sizeState = sizeDotState - call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization - stt%xi => plasticState(p)%state (1,:) + stt%xi => plasticState(ph)%state (1,:) stt%xi = xi_0 - dot%xi => plasticState(p)%dotState(1,:) - plasticState(p)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) - if (plasticState(p)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi' + dot%xi => plasticState(ph)%dotState(1,:) + plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) + if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi' - stt%gamma => plasticState(p)%state (2,:) - dot%gamma => plasticState(p)%dotState(2,:) - plasticState(p)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if (plasticState(p)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma' + stt%gamma => plasticState(ph)%state (2,:) + dot%gamma => plasticState(ph)%dotState(2,:) + plasticState(ph)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + if (plasticState(ph)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma' ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:) + plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:) - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally end associate @@ -240,13 +230,11 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) tr=math_trace33(math_spherical33(Mi)) - if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero + if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero Li = math_I3 & * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) & * tr * abs(tr)**(prm%n-1.0_pReal) - forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) - else Li = 0.0_pReal dLi_dMi = 0.0_pReal diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 14aa58fb9..8acf42710 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -25,8 +25,7 @@ submodule(phase:plastic) kinehardening nonSchmid_pos, & nonSchmid_neg integer :: & - sum_N_sl, & !< total number of active slip system - of_debug = 0 + sum_N_sl logical :: & nonSchmidActive = .false. character(len=pStringLen), allocatable, dimension(:) :: & @@ -62,7 +61,7 @@ module function plastic_kinehardening_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - p, i, o, & + ph, o, & Nconstituents, & sizeState, sizeDeltaState, sizeDotState, & startIndex, endIndex @@ -93,15 +92,13 @@ module function plastic_kinehardening_init() result(myPlasticity) allocate(deltaState(phases%length)) - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - phase => phases%get(p) + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle + + associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph)) + + phase => phases%get(ph) mech => phase%get('mechanics') - i = p - associate(prm => param(i), & - dot => dotState(i), & - dlt => deltaState(i), & - stt => state(i)) pl => mech%get('plasticity') #if defined (__GFORTRAN__) @@ -110,12 +107,6 @@ module function plastic_kinehardening_init() result(myPlasticity) prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray) #endif -#ifdef DEBUG - if (p==material_phaseAt(debugConstitutive%grain,debugConstitutive%element)) then - prm%of_debug = material_phasememberAt(debugConstitutive%grain,debugConstitutive%ip,debugConstitutive%element) - endif -#endif - !-------------------------------------------------------------------------------------------------- ! slip related parameters N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) @@ -174,55 +165,55 @@ module function plastic_kinehardening_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == p) - sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip - sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names + Nconstituents = count(material_phaseAt2 == ph) + sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml + sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeState = sizeDotState + sizeDeltaState - call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%crss => plasticState(p)%state (startIndex:endIndex,:) + stt%crss => plasticState(ph)%state (startIndex:endIndex,:) stt%crss = spread(xi_0, 2, Nconstituents) - dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) - if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' + dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) - dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) + stt%crss_back => plasticState(ph)%state (startIndex:endIndex,:) + dot%crss_back => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%accshear => plasticState(p)%state (startIndex:endIndex,:) - dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + stt%accshear => plasticState(ph)%state (startIndex:endIndex,:) + dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) - o = plasticState(p)%offsetDeltaState + o = plasticState(ph)%offsetDeltaState startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%sense => plasticState(p)%state (startIndex :endIndex ,:) - dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + stt%sense => plasticState(ph)%state (startIndex :endIndex ,:) + dlt%sense => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) - dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + stt%chi0 => plasticState(ph)%state (startIndex :endIndex ,:) + dlt%chi0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) - dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:) + dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally end associate diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index 1510262cc..28e9fbc7c 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -16,8 +16,7 @@ module function plastic_none_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - p, & - Nconstituents + ph class(tNode), pointer :: & phases @@ -29,10 +28,9 @@ module function plastic_none_init() result(myPlasticity) print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - Nconstituents = count(material_phaseAt2 == p) - call phase_allocateState(plasticState(p),Nconstituents,0,0,0) + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle + call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0) enddo end function plastic_none_init diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 0212baee8..fbfcfa5af 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -176,7 +176,7 @@ module function plastic_nonlocal_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & Ninstances, & - p, i, & + ph, & Nconstituents, & sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & s1, s2, & @@ -220,21 +220,17 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(deltaState(phases%length)) allocate(microstructure(phases%length)) - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - phase => phases%get(p) - mech => phase%get('mechanics') + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle - i = p - associate(prm => param(i), & - dot => dotState(i), & - stt => state(i), & - st0 => state0(i), & - del => deltaState(i), & - dst => microstructure(i)) + associate(prm => param(ph), dot => dotState(ph), stt => state(ph), & + st0 => state0(ph), del => deltaState(ph), dst => microstructure(ph)) + + phase => phases%get(ph) + mech => phase%get('mechanics') pl => mech%get('plasticity') - phase_localPlasticity(p) = .not. pl%contains('nonlocal') + phase_localPlasticity(ph) = .not. pl%contains('nonlocal') #if defined (__GFORTRAN__) prm%output = output_asStrings(pl) @@ -245,8 +241,8 @@ module function plastic_nonlocal_init() result(myPlasticity) prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0e4_pReal) ! This data is read in already in lattice - prm%mu = lattice_mu(p) - prm%nu = lattice_nu(p) + prm%mu = lattice_mu(ph) + prm%nu = lattice_nu(ph) ini%N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(ini%N_sl)) @@ -402,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == p) + Nconstituents = count(material_phaseAt2 == ph) sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & @@ -416,101 +412,101 @@ module function plastic_nonlocal_init() result(myPlasticity) 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState) - allocate(geom(p)%V_0(Nconstituents)) - call storeGeometry(p) + allocate(geom(ph)%V_0(Nconstituents)) + call storeGeometry(ph) - plasticState(p)%nonlocal = pl%get_asBool('nonlocal') - if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) & + plasticState(ph)%nonlocal = pl%get_asBool('nonlocal') + if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & call IO_error(212,ext_msg='IPneighborhood does not exist') - plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention + plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - stt%rho => plasticState(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - dot%rho => plasticState(p)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - del%rho => plasticState(p)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - plasticState(p)%atol(1:10*prm%sum_N_sl) = prm%atol_rho + st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho - stt%rhoSgl => plasticState(p)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - dot%rhoSgl => plasticState(p)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - del%rhoSgl => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + stt%rhoSgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + dot%rhoSgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + del%rhoSgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - stt%rhoSglMobile => plasticState(p)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - dot%rhoSglMobile => plasticState(p)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - del%rhoSglMobile => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + stt%rhoSglMobile => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + dot%rhoSglMobile => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + del%rhoSglMobile => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) - dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) - del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) + stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) + dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) + del%rho_sgl_mob_edg_pos => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) - stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:) - dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:) - del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:) + stt%rho_sgl_mob_edg_neg => plasticState(ph)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:) + dot%rho_sgl_mob_edg_neg => plasticState(ph)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:) + del%rho_sgl_mob_edg_neg => plasticState(ph)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:) - stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:) - dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:) - del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:) + stt%rho_sgl_mob_scr_pos => plasticState(ph)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:) + dot%rho_sgl_mob_scr_pos => plasticState(ph)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:) + del%rho_sgl_mob_scr_pos => plasticState(ph)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:) - stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + stt%rho_sgl_mob_scr_neg => plasticState(ph)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - stt%rhoSglImmobile => plasticState(p)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - dot%rhoSglImmobile => plasticState(p)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - del%rhoSglImmobile => plasticState(p)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + stt%rhoSglImmobile => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + dot%rhoSglImmobile => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + del%rhoSglImmobile => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) - dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) - del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) + stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) + dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) + del%rho_sgl_imm_edg_pos => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) - stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:) - dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:) - del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:) + stt%rho_sgl_imm_edg_neg => plasticState(ph)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:) + dot%rho_sgl_imm_edg_neg => plasticState(ph)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:) + del%rho_sgl_imm_edg_neg => plasticState(ph)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:) - stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:) - dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:) - del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:) + stt%rho_sgl_imm_scr_pos => plasticState(ph)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:) + dot%rho_sgl_imm_scr_pos => plasticState(ph)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:) + del%rho_sgl_imm_scr_pos => plasticState(ph)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:) - stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + stt%rho_sgl_imm_scr_neg => plasticState(ph)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - stt%rhoDip => plasticState(p)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - dot%rhoDip => plasticState(p)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - del%rhoDip => plasticState(p)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + stt%rhoDip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + dot%rhoDip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + del%rhoDip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - stt%rho_dip_edg => plasticState(p)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) - dot%rho_dip_edg => plasticState(p)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) - del%rho_dip_edg => plasticState(p)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) + stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) + dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) + del%rho_dip_edg => plasticState(ph)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) - stt%rho_dip_scr => plasticState(p)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + stt%rho_dip_scr => plasticState(ph)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) - dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) - del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) - plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) - if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & + stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) + dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) + del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) + plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) + if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & extmsg = trim(extmsg)//' atol_gamma' - plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) + plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) - stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents) - stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) - stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents) - stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents) - stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents) - stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) + stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents) + stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) + stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents) + stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents) + stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents) + stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) end associate - if (Nconstituents > 0) call stateInit(ini,p,Nconstituents) - plasticState(p)%state0 = plasticState(p)%state + if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents) + plasticState(ph)%state0 = plasticState(ph)%state !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range @@ -526,35 +522,34 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0) allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0) - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length - if(.not. myPlasticity(p)) cycle - i = p + if(.not. myPlasticity(ph)) cycle - Nconstituents = count(material_phaseAt2 == p) + phase => phases%get(ph) + Nconstituents = count(material_phaseAt2 == ph) l = 0 do t = 1,4 - do s = 1,param(i)%sum_N_sl + do s = 1,param(ph)%sum_N_sl l = l + 1 - iRhoU(s,t,i) = l + iRhoU(s,t,ph) = l enddo enddo - l = l + (4+2+1+1)*param(i)%sum_N_sl ! immobile(4), dipole(2), shear, forest + l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest do t = 1,4 - do s = 1,param(i)%sum_N_sl + do s = 1,param(ph)%sum_N_sl l = l + 1 - iV(s,t,i) = l + iV(s,t,ph) = l enddo enddo do t = 1,2 - do s = 1,param(i)%sum_N_sl + do s = 1,param(ph)%sum_N_sl l = l + 1 - iD(s,t,i) = l + iD(s,t,ph) = l enddo enddo - if (iD(param(i)%sum_N_sl,2,i) /= plasticState(p)%sizeState) & - call IO_error(0, ext_msg = 'state indices not properly set (nonlocal)') + if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) & + error stop 'state indices not properly set (nonlocal)' enddo end function plastic_nonlocal_init diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 94214f14a..d769431b4 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -70,7 +70,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - p, i, & + ph, i, & Nconstituents, & sizeState, sizeDotState, & startIndex, endIndex @@ -101,14 +101,13 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) allocate(state(phases%length)) allocate(dotState(phases%length)) - do p = 1, phases%length - if(.not. myPlasticity(p)) cycle - phase => phases%get(p) + do ph = 1, phases%length + if(.not. myPlasticity(ph)) cycle + + associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) + + phase => phases%get(ph) mech => phase%get('mechanics') - i = p - associate(prm => param(i), & - dot => dotState(i), & - stt => state(i)) pl => mech%get('plasticity') !-------------------------------------------------------------------------------------------------- @@ -135,7 +134,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) xi_0_sl = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl)) prm%xi_inf_sl = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl)) prm%h_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), & - defaultVal=[(0.0_pReal,i=1,size(N_sl))]) + defaultVal=[(0.0_pReal,i=1,size(N_sl))]) prm%dot_gamma_0_sl = pl%get_asFloat('dot_gamma_0_sl') prm%n_sl = pl%get_asFloat('n_sl') @@ -224,49 +223,49 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == p) + Nconstituents = count(material_phaseAt2 == ph) sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState - call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:) stt%xi_slip = spread(xi_0_sl, 2, Nconstituents) - dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) - if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' + dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:) stt%xi_twin = spread(xi_0_tw, 2, Nconstituents) - dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) - if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' + dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) - dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + stt%gamma_slip => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) - dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally end associate