From f5af3526441f6131d2e12d29d828a0786859b665 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 00:09:51 +0200 Subject: [PATCH] whitespace adjustments --- ...phase_mechanical_plastic_dislotungsten.f90 | 9 +- src/phase_mechanical_plastic_dislotwin.f90 | 76 +++++++-------- src/phase_mechanical_plastic_isotropic.f90 | 56 +++++------ ...phase_mechanical_plastic_kinehardening.f90 | 96 +++++++++---------- ...phase_mechanical_plastic_phenopowerlaw.f90 | 91 +++++++++--------- 5 files changed, 163 insertions(+), 165 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index a47337f4d..71776ba47 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -378,13 +378,14 @@ module subroutine dislotungsten_dependentState(ph,en) real(pReal), dimension(param(ph)%sum_N_sl) :: & dislocationSpacing + associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) - dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%threshold_stress(:,en) = prm%mu*prm%b_sl & - * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dst%threshold_stress(:,en) = prm%mu*prm%b_sl & + * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) end associate diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 33f148036..6ecf9ee7a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -598,8 +598,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & - * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & - * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) + * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & + * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) Lp = Lp + dot_gamma_sb * P_sb forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -734,54 +734,52 @@ module subroutine dislotwin_dependentState(T,ph,en) x0 - associate(prm => param(ph),& - stt => state(ph),& - dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) - sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) + sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) + sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) - Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T + Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T - !* rescaled volume fraction for topology - f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... - f_over_t_tr = sumf_tr/prm%t_tr ! but this not + !* rescaled volume fraction for topology + f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... + f_over_t_tr = sumf_tr/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl - if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) - if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) - dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) + inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl + if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & + inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) + if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & + inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) + dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) - inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) - dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) + dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) - inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) - dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) + inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) + dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) - !* threshold stress for dislocation motion - dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + !* threshold stress for dislocation motion + dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - !* threshold stress for growing twin/martensite - if(prm%sum_N_tw == prm%sum_N_sl) & - dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & - + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? - if(prm%sum_N_tr == prm%sum_N_sl) & - dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & - + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? - + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) + !* threshold stress for growing twin/martensite + if(prm%sum_N_tw == prm%sum_N_sl) & + dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & + + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? + if(prm%sum_N_tr == prm%sum_N_sl) & + dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & + + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? + + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) - dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal - dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal + dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal + dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal - x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans - dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans + dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) - x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip - dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) + x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip + dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) end associate @@ -880,9 +878,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & 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)) - enddo + tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] tau_eff = abs(tau)-dst%tau_pass(:,en) diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 39a619897..8b10ed5ef 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -170,27 +170,28 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: & k, l, m, n + associate(prm => param(ph), stt => state(ph)) - Mp_dev = math_deviatoric33(Mp) - squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) - norm_Mp_dev = sqrt(squarenorm_Mp_dev) + Mp_dev = math_deviatoric33(Mp) + squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) + 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 + 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 - Lp = dot_gamma * Mp_dev/norm_Mp_dev - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev - forall (k=1:3,l=1:3) & - dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal - forall (k=1:3,m=1:3) & - dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal - dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev - else - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - end if + Lp = dot_gamma * Mp_dev/norm_Mp_dev + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev + forall (k=1:3,l=1:3) & + dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal + forall (k=1:3,m=1:3) & + dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal + dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev + else + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + end if end associate @@ -218,19 +219,20 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) integer :: & k, l, m, n + associate(prm => param(ph), stt => state(ph)) - tr=math_trace33(math_spherical33(Mi)) + 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 - Li = math_I3 & - * prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-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 - endif + 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 * (3.0_pReal*prm%M*stt%xi(en))**(-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 + endif end associate diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index b68229c12..9195cba45 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -280,26 +280,25 @@ 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 => dotState(ph)) - call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) - dot%accshear(:,en) = abs(dot_gamma_pos+dot_gamma_neg) - sumGamma = sum(stt%accshear(:,en)) + call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) + dot%accshear(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + sumGamma = sum(stt%accshear(:,en)) - dot%crss(:,en) = matmul(prm%h_sl_sl,dot%accshear(:,en)) & - * ( 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%crss(:,en) = matmul(prm%h_sl_sl,dot%accshear(:,en)) & + * ( 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%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * & - ( 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%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))& - ) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) & - ) + dot%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * & + ( 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%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))& + ) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) & + ) end associate @@ -325,9 +324,8 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,en) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) sense = merge(state(ph)%sense(:,en), & ! keep existing... - sign(1.0_pReal,dot_gamma_pos+dot_gamma_neg), & ! ...or have a defined - dEq0(dot_gamma_pos+dot_gamma_neg,1e-10_pReal)) ! current sense of shear direction - + sign(1.0_pReal,dot_gamma_pos+dot_gamma_neg), & ! ...or have a defined + dEq0(dot_gamma_pos+dot_gamma_neg,1e-10_pReal)) ! current sense of shear direction !-------------------------------------------------------------------------------------------------- ! switch in sense of shear? @@ -412,42 +410,44 @@ pure subroutine kinetics(Mp,ph,en, & tau_neg integer :: i + associate(prm => param(ph), stt => state(ph)) - do i = 1, prm%sum_N_sl - tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%crss_back(i,en) - tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%crss_back(i,en), & - 0.0_pReal, prm%nonSchmidActive) - enddo + do i = 1, prm%sum_N_sl + tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%crss_back(i,en) + tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%crss_back(i,en), & + 0.0_pReal, prm%nonSchmidActive) + enddo - where(dNeq0(tau_pos)) - dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos) - else where - dot_gamma_pos = 0.0_pReal - end where - - where(dNeq0(tau_neg)) - dot_gamma_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg) - else where - dot_gamma_neg = 0.0_pReal - end where - - if (present(ddot_gamma_dtau_pos)) then - where(dNeq0(dot_gamma_pos)) - ddot_gamma_dtau_pos = dot_gamma_pos*prm%n/tau_pos + where(dNeq0(tau_pos)) + dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos) else where - ddot_gamma_dtau_pos = 0.0_pReal + dot_gamma_pos = 0.0_pReal end where - endif - if (present(ddot_gamma_dtau_neg)) then - where(dNeq0(dot_gamma_neg)) - ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg + + where(dNeq0(tau_neg)) + dot_gamma_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg) else where - ddot_gamma_dtau_neg = 0.0_pReal + dot_gamma_neg = 0.0_pReal end where - endif + + if (present(ddot_gamma_dtau_pos)) then + where(dNeq0(dot_gamma_pos)) + ddot_gamma_dtau_pos = dot_gamma_pos*prm%n/tau_pos + else where + ddot_gamma_dtau_pos = 0.0_pReal + end where + endif + if (present(ddot_gamma_dtau_neg)) then + where(dNeq0(dot_gamma_neg)) + ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg + else where + ddot_gamma_dtau_neg = 0.0_pReal + end where + endif + end associate end subroutine kinetics diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 2a1acb228..3bfb360b4 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -344,10 +344,9 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en) dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en)) - 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, & + 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) & @@ -428,40 +427,41 @@ pure subroutine kinetics_sl(Mp,ph,en, & associate(prm => param(ph), stt => state(ph)) - do i = 1, prm%sum_N_sl - tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), & - 0.0_pReal, prm%nonSchmidActive) - enddo + do i = 1, prm%sum_N_sl + tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) + tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), & + 0.0_pReal, prm%nonSchmidActive) + enddo - where(dNeq0(tau_sl_pos)) - dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos) - else where - dot_gamma_sl_pos = 0.0_pReal - end where - - where(dNeq0(tau_sl_neg)) - dot_gamma_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg) - else where - dot_gamma_sl_neg = 0.0_pReal - end where - - if (present(ddot_gamma_dtau_sl_pos)) then - where(dNeq0(dot_gamma_sl_pos)) - ddot_gamma_dtau_sl_pos = dot_gamma_sl_pos*prm%n_sl/tau_sl_pos + where(dNeq0(tau_sl_pos)) + dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos) else where - ddot_gamma_dtau_sl_pos = 0.0_pReal + dot_gamma_sl_pos = 0.0_pReal end where - endif - if (present(ddot_gamma_dtau_sl_neg)) then - where(dNeq0(dot_gamma_sl_neg)) - ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg + + where(dNeq0(tau_sl_neg)) + dot_gamma_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg) else where - ddot_gamma_dtau_sl_neg = 0.0_pReal + dot_gamma_sl_neg = 0.0_pReal end where - endif + + if (present(ddot_gamma_dtau_sl_pos)) then + where(dNeq0(dot_gamma_sl_pos)) + ddot_gamma_dtau_sl_pos = dot_gamma_sl_pos*prm%n_sl/tau_sl_pos + else where + ddot_gamma_dtau_sl_pos = 0.0_pReal + end where + endif + if (present(ddot_gamma_dtau_sl_neg)) then + where(dNeq0(dot_gamma_sl_neg)) + ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg + else where + ddot_gamma_dtau_sl_neg = 0.0_pReal + end where + endif + end associate end subroutine kinetics_sl @@ -492,26 +492,25 @@ pure subroutine kinetics_tw(Mp,ph,en,& tau_tw integer :: i + associate(prm => param(ph), stt => state(ph)) - do i = 1, prm%sum_N_tw - tau_tw(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - enddo + tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] - where(tau_tw > 0.0_pReal) - dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction - * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw - else where - dot_gamma_tw = 0.0_pReal - end where - - if (present(ddot_gamma_dtau_tw)) then - where(dNeq0(dot_gamma_tw)) - ddot_gamma_dtau_tw = dot_gamma_tw*prm%n_tw/tau_tw + where(tau_tw > 0.0_pReal) + dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction + * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw else where - ddot_gamma_dtau_tw = 0.0_pReal + dot_gamma_tw = 0.0_pReal end where - endif + + if (present(ddot_gamma_dtau_tw)) then + where(dNeq0(dot_gamma_tw)) + ddot_gamma_dtau_tw = dot_gamma_tw*prm%n_tw/tau_tw + else where + ddot_gamma_dtau_tw = 0.0_pReal + end where + endif end associate