From c09c2a6c8e123d42ee33527675a5096275cdf236 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 00:50:42 +0100 Subject: [PATCH] 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))