From 66edad1cf8a0813272f3ee15f11598a469848088 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Sep 2018 22:07:59 +0200 Subject: [PATCH] avoid code doubling shear rates needed multiple times. For now, just introduced but not used. Will become active once the test passes again --- src/plastic_dislotwin.f90 | 107 ++++++++++++++++++++++++++++++++------ 1 file changed, 91 insertions(+), 16 deletions(-) diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index d840e19c4..85234ae4f 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -1227,8 +1227,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature Ndot0_twin=prm%Ndot0_twin(i) endif isFCCtwin - gdot_twin = (1.0_pReal-sumf_twin-sumf_trans)* prm%shear_twin(i) * mse%twinVolume(i,of) & - * Ndot0_twin*exp(-StressRatio_r) + gdot_twin = prm%shear_twin(i) * mse%twinVolume(i,of) * Ndot0_twin*exp(-StressRatio_r) + gdot_twin = (1.0_pReal-sumf_twin-sumf_trans) * gdot_twin dgdot_dtau = ((gdot_twin*prm%r(i))/tau)*StressRatio_r Lp = Lp + gdot_twin*prm%Schmid_twin(1:3,1:3,i) @@ -1261,8 +1261,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature Ndot0_trans=prm%Ndot0_trans(i) endif isFCCtrans - gdot_trans = (1.0_pReal-sumf_twin-sumf_trans)* mse%martensiteVolume(i,of) & - * Ndot0_trans*exp(-StressRatio_s) + gdot_trans = mse%martensiteVolume(i,of) * Ndot0_trans*exp(-StressRatio_s) + gdot_trans = (1.0_pReal-sumf_twin-sumf_trans)* gdot_trans dgdot_dtau = ((gdot_trans*prm%s(i))/tau)*StressRatio_s Lp = Lp + gdot_trans*prm%Schmid_trans(1:3,1:3,i) @@ -1482,40 +1482,115 @@ subroutine kinetics_slip(prm,stt,mse,of,S,temperature,gdot_slip,dgdot_dtau_slip) gdot_slip real, dimension(prm%totalNslip), optional, intent(out) :: & dgdot_dtau_slip + real, dimension(prm%totalNslip) :: & + dgdot_dtau real(pReal), dimension(3,3), intent(in) :: & S real(pReal), intent(in) :: & temperature real, dimension(prm%totalNslip) :: & - tau_slip, & + tau, & stressRatio, & StressRatio_p, & BoltzmannRatio integer(pInt) :: i do i = 1_pInt, prm%totalNslip - tau_slip = math_mul33xx33(S,prm%Schmid_slip(1:3,1:3,i)) + tau = math_mul33xx33(S,prm%Schmid_slip(1:3,1:3,i)) enddo - where((abs(tau_slip)-mse%threshold_stress_slip(:,of)) > tol_math_check) - stressRatio = ((abs(tau_slip)- mse%threshold_stress_slip(:,of))/& + significantStress: where((abs(tau)-mse%threshold_stress_slip(:,of)) > tol_math_check) + stressRatio = ((abs(tau)- mse%threshold_stress_slip(:,of))/& (prm%SolidSolutionStrength+prm%tau_peierls(:))) StressRatio_p = stressRatio** prm%p BoltzmannRatio = prm%Qedge/(kB*Temperature) gdot_slip = stt%rhoEdge(:,of)*prm%burgers_slip* prm%v0 & - * sign(exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q), tau_slip) - dgdot_dtau_slip = abs(gdot_slip)*BoltzmannRatio*prm%p * prm%q & - / (prm%SolidSolutionStrength+prm%tau_peierls) & - * stressRatio**(prm%p-1.0_pReal)*(1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) - else where + * sign(exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q), tau) + dgdot_dtau = abs(gdot_slip)*BoltzmannRatio*prm%p * prm%q & + / (prm%SolidSolutionStrength+prm%tau_peierls) & + * stressRatio**(prm%p-1.0_pReal)*(1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) + else where significantStress gdot_slip = 0.0_pReal - dgdot_dtau_slip = 0.0_pReal - end where - + dgdot_dtau = 0.0_pReal + end where significantStress + + if(present(dgdot_dtau_slip)) dgdot_dtau_slip = dgdot_dtau end subroutine +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on slip systems +!-------------------------------------------------------------------------------------------------- +subroutine kinetics_twin(prm,stt,mse,of,S,temperature,gdot_slip,gdot_twin,dgdot_dtau_twin) + use prec, only: & + tol_math_check, & + dNeq0 + use math, only: & + math_mul33xx33 + + implicit none + type(tParameters), intent(in) :: & + prm + type(tDislotwinState), intent(in) :: & + stt + integer(pInt), intent(in) :: & + of + type(tDislotwinMicrostructure) :: & + mse + real, dimension(prm%totalNslip), intent(out) :: & + gdot_slip + real, dimension(prm%totalNtwin), intent(out) :: & + gdot_twin + real, dimension(prm%totalNtwin), optional, intent(out) :: & + dgdot_dtau_twin + real(pReal), dimension(3,3), intent(in) :: & + S + real(pReal), intent(in) :: & + temperature + + real, dimension(prm%totalNtwin) :: & + tau, & + Ndot0_twin, & + stressRatio_r, & + dgdot_dtau + + integer(pInt) :: i,s1,s2 + + do i = 1_pInt, prm%totalNtwin + tau = math_mul33xx33(S,prm%Schmid_twin(1:3,1:3,i)) + isFCC: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < mse%tau_r_twin(i,of)) then + Ndot0_twin=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& !!!!! correct? + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_twin*prm%burgers_slip(i))*& + (1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)*& + (mse%tau_r_twin(i,of)-tau))) + else + Ndot0_twin=0.0_pReal + end if + else isFCC + Ndot0_twin=prm%Ndot0_twin(i) + endif isFCC + enddo + + + significantStress: where(tau > tol_math_check) + StressRatio_r = (mse%threshold_stress_twin(:,of)/tau)**prm%r + gdot_twin = prm%shear_twin * mse%twinVolume(:,of) * Ndot0_twin*exp(-StressRatio_r) + dgdot_dtau = ((gdot_twin*prm%r)/tau)*StressRatio_r + else where significantStress + gdot_twin = 0.0_pReal + dgdot_dtau = 0.0_pReal + end where significantStress + + if(present(dgdot_dtau_twin)) dgdot_dtau_twin = dgdot_dtau + +end subroutine + + !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !--------------------------------------------------------------------------------------------------