From a9a5c8fb7374d9235ceecefd907bbd39f4095321 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 3 Feb 2022 23:06:14 +0100 Subject: [PATCH 1/2] simplify acces with pointer naming will be adjusted once global deltaState is removed --- src/phase.f90 | 2 ++ src/phase_mechanical.f90 | 8 ++++---- src/phase_mechanical_plastic.f90 | 7 +++---- src/prec.f90 | 2 ++ 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index b7bace81a..098cdebe8 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -436,6 +436,8 @@ subroutine phase_allocateState(state, & allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal) allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal) + state%deltaState2 => state%state(state%offsetDeltaState+1: & + state%offsetDeltaState+state%sizeDeltaState,:) end subroutine phase_allocateState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 04cc4b946..733219d56 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -81,13 +81,13 @@ submodule(phase) mechanical module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element + co, & !< constituent + ip, & !< integration point + el, & !< element ph, & en real(pReal), intent(in) :: & - subdt !< timestep + subdt !< timestep real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState end function plastic_dotState diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 912aadc03..bb5a9819a 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -415,10 +415,9 @@ module function plastic_deltaState(ph, en) result(broken) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) if (.not. broken) then - myOffset = plasticState(ph)%offsetDeltaState - mySize = plasticState(ph)%sizeDeltaState - plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = & - plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en) + mySize = plasticState(ph)%sizeDeltaState + plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & + + plasticState(ph)%deltaState(1:mySize,en) end if end select diff --git a/src/prec.f90 b/src/prec.f90 index 8de82fee8..61fa141ba 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -45,6 +45,8 @@ module prec state, & !< state dotState, & !< rate of state change deltaState !< increment of state change + real(pReal), pointer, dimension(:,:) :: & + deltaState2 end type type, extends(tState) :: tPlasticState From ce1eb4f59e4f0bd259edd3ec599ec8fa4cfd81cb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 3 Feb 2022 23:40:25 +0100 Subject: [PATCH 2/2] no need to store dot state saves memory and flow is easier to understand --- src/phase_mechanical_plastic.f90 | 39 ++++++----- ...phase_mechanical_plastic_dislotungsten.f90 | 57 +++++++++------ src/phase_mechanical_plastic_isotropic.f90 | 44 ++++++------ ...phase_mechanical_plastic_kinehardening.f90 | 56 +++++++++------ src/phase_mechanical_plastic_nonlocal.f90 | 4 +- ...phase_mechanical_plastic_phenopowerlaw.f90 | 70 +++++++++++-------- 6 files changed, 155 insertions(+), 115 deletions(-) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index bb5a9819a..70181d724 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic end subroutine nonlocal_LpAndItsTangent - module subroutine isotropic_dotState(Mp,ph,en) + module function isotropic_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en - end subroutine isotropic_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function isotropic_dotState - module subroutine phenopowerlaw_dotState(Mp,ph,en) + module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en - end subroutine phenopowerlaw_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function phenopowerlaw_dotState - module subroutine plastic_kinehardening_dotState(Mp,ph,en) + module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en - end subroutine plastic_kinehardening_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function plastic_kinehardening_dotState module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(3,3), intent(in) :: & @@ -144,15 +150,15 @@ submodule(phase:mechanical) plastic en end subroutine dislotwin_dotState - module subroutine dislotungsten_dotState(Mp,T,ph,en) + module function dislotungsten_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en - end subroutine dislotungsten_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function dislotungsten_dotState module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & @@ -318,27 +324,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) plasticType: select case (phase_plasticity(ph)) case (PLASTIC_ISOTROPIC_ID) plasticType - call isotropic_dotState(Mp,ph,en) + dotState = isotropic_dotState(Mp,ph,en) case (PLASTIC_PHENOPOWERLAW_ID) plasticType - call phenopowerlaw_dotState(Mp,ph,en) + dotState = phenopowerlaw_dotState(Mp,ph,en) case (PLASTIC_KINEHARDENING_ID) plasticType - call plastic_kinehardening_dotState(Mp,ph,en) + dotState = plastic_kinehardening_dotState(Mp,ph,en) case (PLASTIC_DISLOTWIN_ID) plasticType call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) + dotState = plasticState(ph)%dotState(:,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) + dotState = dislotungsten_dotState(Mp,ph,en) case (PLASTIC_NONLOCAL_ID) plasticType call nonlocal_dotState(Mp,subdt,ph,en,ip,el) + dotState = plasticState(ph)%dotState(:,en) + end select plasticType end if - dotState = plasticState(ph)%dotState(:,en) - end function plastic_dotState diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 9f0465158..5a6ac8f5c 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -43,6 +43,13 @@ submodule(phase:plastic) dislotungsten systems_sl end type tParameters !< container type for internal constitutive parameters + type :: tIndexDotState + integer, dimension(2) :: & + rho_mob, & + rho_dip, & + gamma_sl + end type tIndexDotState + type :: tDislotungstenState real(pReal), dimension(:,:), pointer :: & rho_mob, & @@ -58,10 +65,9 @@ submodule(phase:plastic) dislotungsten !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tDisloTungstenState), allocatable, dimension(:) :: & - dotState, & - state + type(tParameters), allocatable, dimension(:) :: param + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tDisloTungstenState), allocatable, dimension(:) :: state type(tDisloTungstenDependentState), allocatable, dimension(:) :: dependentState contains @@ -103,18 +109,17 @@ module function plastic_dislotungsten_init() result(myPlasticity) print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002' - phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) allocate(dependentState(phases%length)) - do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) mech => phase%get('mechanical') @@ -214,28 +219,29 @@ module function plastic_dislotungsten_init() result(myPlasticity) sizeState = sizeDotState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl + idx_dot%rho_mob = [startIndex,endIndex] stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob = spread(rho_mob_0,2,Nmembers) - 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 + idx_dot%rho_dip = [startIndex,endIndex] stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip = spread(rho_dip_0,2,Nmembers) - 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 + idx_dot%gamma_sl = [startIndex,endIndex] stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) - dot%gamma_sl => 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' @@ -300,15 +306,15 @@ end subroutine dislotungsten_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotungsten_dotState(Mp,T,ph,en) +module function dislotungsten_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos, dot_gamma_neg,& @@ -319,17 +325,22 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_climb, & d_hat real(pReal) :: & - mu + mu, T - associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) + + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), & + dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), & + dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), & + dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2))) mu = elastic_mu(ph,en) + T = thermal_T(ph,en) call kinetics(Mp,T,ph,en,& dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + dot_gamma_sl = abs(dot_gamma_pos+dot_gamma_neg) where(dEq0((tau_pos+tau_neg)*0.5_pReal)) dot_rho_dip_formation = 0.0_pReal @@ -338,7 +349,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), & prm%d_caron, & ! lower limit dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & + dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot_gamma_sl/prm%b_sl, & 0.0_pReal, & prm%dipoleformation) v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) & @@ -346,16 +357,16 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication + dot_rho_mob = dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges - dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot_gamma_sl ! Spontaneous annihilation of 2 edges + dot_rho_dip = dot_rho_dip_formation & + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot_gamma_sl & ! Spontaneous annihilation of an edge with a dipole - dot_rho_dip_climb end associate -end subroutine dislotungsten_dotState +end function dislotungsten_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index bea5339c7..49efc9fae 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -37,9 +37,7 @@ submodule(phase:plastic) isotropic !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param - type(tIsotropicState), allocatable, dimension(:) :: & - dotState, & - state + type(tIsotropicState), allocatable, dimension(:) :: state contains @@ -77,16 +75,15 @@ module function plastic_isotropic_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) do ph = 1, phases%length if(.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') #if defined (__GFORTRAN__) prm%output = output_as1dString(pl) @@ -125,12 +122,12 @@ module function plastic_isotropic_init() result(myPlasticity) sizeState = sizeDotState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization - stt%xi => plasticState(ph)%state (1,:) - stt%xi = xi_0 - dot%xi => plasticState(ph)%dotState(1,:) + stt%xi => plasticState(ph)%state(1,:) + stt%xi = xi_0 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' @@ -178,7 +175,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) norm_Mp_dev = sqrt(squarenorm_Mp_dev) if (norm_Mp_dev > 0.0_pReal) then - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en)))**prm%n Lp = dot_gamma * Mp_dev/norm_Mp_dev forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -242,27 +239,26 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine isotropic_dotState(Mp,ph,en) +module function isotropic_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal) :: & dot_gamma, & !< strainrate xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - associate(prm => param(ph), stt => state(ph), & - dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph), dot_xi => dotState(1)) - if (prm%dilatation) then - norm_Mp = sqrt(math_tensordot(Mp,Mp)) - else - norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))) - end if + norm_Mp = merge(sqrt(math_tensordot(Mp,Mp)), & + sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))), & + prm%dilatation) dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n @@ -274,16 +270,16 @@ module subroutine isotropic_dotState(Mp,ph,en) + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) end if - dot%xi(en) = dot_gamma & - * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & - * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star) + dot_xi = dot_gamma & + * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & + * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star) else - dot%xi(en) = 0.0_pReal + dot_xi = 0.0_pReal end if end associate -end subroutine isotropic_dotState +end function isotropic_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 03eb27f31..0f6bd53d4 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -34,6 +34,13 @@ submodule(phase:plastic) kinehardening systems_sl end type tParameters + type :: tIndexDotState + integer, dimension(2) :: & + xi, & + chi, & + gamma + end type tIndexDotState + type :: tKinehardeningState real(pReal), pointer, dimension(:,:) :: & xi, & !< resistance against plastic slip @@ -47,10 +54,8 @@ submodule(phase:plastic) kinehardening !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param - type(tKinehardeningState), allocatable, dimension(:) :: & - dotState, & - deltaState, & - state + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tKinehardeningState), allocatable, dimension(:) :: state, deltaState contains @@ -91,19 +96,20 @@ module function plastic_kinehardening_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) allocate(deltaState(phases%length)) do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') #if defined (__GFORTRAN__) prm%output = output_as1dString(pl) @@ -173,27 +179,28 @@ module function plastic_kinehardening_init() result(myPlasticity) sizeState = sizeDotState + sizeDeltaState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi => plasticState(ph)%state (startIndex:endIndex,:) + idx_dot%xi = [startIndex,endIndex] + stt%xi => plasticState(ph)%state(startIndex:endIndex,:) stt%xi = spread(xi_0, 2, Nmembers) - dot%xi => 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%chi => plasticState(ph)%state (startIndex:endIndex,:) - dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%chi = [startIndex,endIndex] + stt%chi => plasticState(ph)%state(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%gamma => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%gamma = [startIndex,endIndex] + stt%gamma => plasticState(ph)%state(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' @@ -270,13 +277,15 @@ end subroutine kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_dotState(Mp,ph,en) +module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal) :: & sumGamma @@ -284,29 +293,32 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en) dot_gamma_pos,dot_gamma_neg - associate(prm => param(ph), stt => state(ph),dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph), & + dot_xi => dotState(IndexDotState(ph)%xi(1):IndexDotState(ph)%xi(2)),& + dot_chi => dotState(IndexDotState(ph)%chi(1):IndexDotState(ph)%chi(2)),& + dot_gamma => dotState(IndexDotState(ph)%gamma(1):IndexDotState(ph)%gamma(2))) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) - dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + dot_gamma = abs(dot_gamma_pos+dot_gamma_neg) sumGamma = sum(stt%gamma(:,en)) - dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) & + dot_xi = matmul(prm%h_sl_sl,dot_gamma) & * ( prm%h_inf_f & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & ) - dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * & - ( prm%h_inf_b + & - (prm%h_0_b - prm%h_inf_b & + dot_chi = stt%sgn_gamma(:,en)*dot_gamma & + * ( prm%h_inf_b & + + (prm%h_0_b - prm%h_inf_b & + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))& ) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) & ) end associate -end subroutine plastic_kinehardening_dotState +end function plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index a8f56b5c5..d668e40fe 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -227,8 +227,8 @@ module function plastic_nonlocal_init() result(myPlasticity) st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.) #if defined (__GFORTRAN__) diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 6e79968f3..0b018562c 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -47,6 +47,14 @@ submodule(phase:plastic) phenopowerlaw systems_tw end type tParameters + type :: tIndexDotState + integer, dimension(2) :: & + xi_sl, & + xi_tw, & + gamma_sl, & + gamma_tw + end type tIndexDotState + type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & xi_sl, & @@ -56,11 +64,10 @@ submodule(phase:plastic) phenopowerlaw end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- -! containers for parameters and state +! containers for parameters, dot state index, and state type(tParameters), allocatable, dimension(:) :: param - type(tPhenopowerlawState), allocatable, dimension(:) :: & - dotState, & - state + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tPhenopowerlawState), allocatable, dimension(:) :: state contains @@ -101,17 +108,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -224,37 +232,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:) + idx_dot%xi_sl = [startIndex,endIndex] + stt%xi_sl => plasticState(ph)%state(startIndex:endIndex,:) stt%xi_sl = spread(xi_0_sl, 2, Nmembers) - dot%xi_sl => 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_tw => plasticState(ph)%state (startIndex:endIndex,:) + idx_dot%xi_tw = [startIndex,endIndex] + stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:) stt%xi_tw = spread(xi_0_tw, 2, Nmembers) - dot%xi_tw => 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%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%gamma_sl = [startIndex,endIndex] + stt%gamma_sl => plasticState(ph)%state(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' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%gamma_tw = [startIndex,endIndex] + stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) end associate @@ -324,13 +332,15 @@ end subroutine phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine phenopowerlaw_dotState(Mp,ph,en) +module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal) :: & xi_sl_sat_offset,& @@ -340,28 +350,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en) right_SlipSlip - associate(prm => param(ph), stt => state(ph), dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph), & + dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), & + dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), & + dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), & + dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2))) call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg) - dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) - call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en)) + dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) + call kinetics_tw(Mp,ph,en,dot_gamma_tw) sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char) xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, & 1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset)) - dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & - * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) & - + matmul(prm%h_sl_tw,dot%gamma_tw(:,en)) + dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & + * matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot_gamma_tw) - dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & - * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) & - + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en)) + dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + * matmul(prm%h_tw_sl,dot_gamma_sl) & + + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw) end associate -end subroutine phenopowerlaw_dotState +end function phenopowerlaw_dotState !--------------------------------------------------------------------------------------------------