From 00230d482f30ce8155852d1e201613bbf6b37a06 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 22:07:23 +0100 Subject: [PATCH] use data from other physics directly more clear code, simplified interfaces --- src/phase.f90 | 2 +- src/phase_mechanical_plastic.f90 | 22 +--- ...phase_mechanical_plastic_dislotungsten.f90 | 8 +- src/phase_mechanical_plastic_dislotwin.f90 | 112 +++++++++--------- src/phase_mechanical_plastic_nonlocal.f90 | 21 ++-- src/phase_thermal.f90 | 2 +- 6 files changed, 81 insertions(+), 86 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 214b0f5fa..6035b4491 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -155,7 +155,7 @@ module phase real(pReal), dimension(3,3) :: P end function phase_P - module function thermal_T(ph,en) result(T) + pure module function thermal_T(ph,en) result(T) integer, intent(in) :: ph,en real(pReal) :: T end function thermal_T diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 4951712fd..0a24b0cbf 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -73,47 +73,37 @@ submodule(phase:mechanical) plastic en end subroutine kinehardening_LpAndItsTangent - module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp - real(pReal), dimension(3,3), intent(in) :: & Mp - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en end subroutine dislotwin_LpAndItsTangent - pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp - real(pReal), dimension(3,3), intent(in) :: & Mp - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en end subroutine dislotungsten_LpAndItsTangent - module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,en) + module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp - real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - Temperature integer, intent(in) :: & ph, & en @@ -282,13 +272,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_DISLOTWIN_ID) plasticType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index c2a584ea1..cd71a7fd7 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -257,27 +257,27 @@ end function plastic_dislotungsten_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,ph,en) + Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature integer, intent(in) :: & ph, & en integer :: & i,k,l,m,n + real(pReal) :: & + T !< temperature real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg + T = thermal_T(ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index c5e043e3f..cfc9a002f 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -513,20 +513,20 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) +module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp real(pReal), dimension(3,3), intent(in) :: Mp integer, intent(in) :: ph,en - real(pReal), intent(in) :: T integer :: i,k,l,m,n real(pReal) :: & - f_unrotated,StressRatio_p,& - E_kB_T, & - ddot_gamma_dtau, & - tau + f_unrotated,StressRatio_p,& + E_kB_T, & + ddot_gamma_dtau, & + tau, & + T real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & @@ -556,69 +556,71 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) 0, 1, 1 & ],pReal),[ 3,6]) - associate(prm => param(ph), stt => state(ph)) - - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + T = thermal_T(ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) - 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) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) - end do slipContribution + associate(prm => param(ph), stt => state(ph)) - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) - twinContibution: do i = 1, prm%sum_N_tw - Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) - end do twinContibution + f_unrotated = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) - 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) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) - end do transContibution + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) + 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) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + end do slipContribution - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) + twinContibution: do i = 1, prm%sum_N_tw + Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + end do twinContibution - shearBandingContribution: if (dNeq0(prm%v_sb)) then + call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) + 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) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + end do transContibution - E_kB_T = prm%E_sb/(K_B*T) - call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? + Lp = Lp * f_unrotated + dLp_dMp = dLp_dMp * f_unrotated - do i = 1,6 - P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& - matmul(eigVectors,sb_mComposition(1:3,i))) - tau = math_tensordot(Mp,P_sb) + shearBandingContribution: if (dNeq0(prm%v_sb)) then - significantShearBandStress: if (abs(tau) > tol_math_check) then - StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb - dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*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) + E_kB_T = prm%E_sb/(K_B*T) + call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? - Lp = Lp + dot_gamma_sb * P_sb - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) - end if significantShearBandStress - end do + do i = 1,6 + P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& + matmul(eigVectors,sb_mComposition(1:3,i))) + tau = math_tensordot(Mp,P_sb) - end if shearBandingContribution + significantShearBandStress: if (abs(tau) > tol_math_check) then + StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb + dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*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) - end associate + Lp = Lp + dot_gamma_sb * P_sb + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) + end if significantShearBandStress + end do + + end if shearBandingContribution + + end associate end subroutine dislotwin_LpAndItsTangent diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 7ff016514..7118055fe 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -741,7 +741,7 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,en) + Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -749,9 +749,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & integer, intent(in) :: & ph, & en - real(pReal), intent(in) :: & - Temperature !< temperature - real(pReal), dimension(3,3), intent(in) :: & Mp !< derivative of Lp with respect to Mp @@ -771,8 +768,14 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms dot_gamma !< shear rate + real(pReal) :: & + Temperature !< temperature + Temperature = thermal_T(ph,en) + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph)) !*** shortcut to state variables @@ -821,8 +824,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal do s = 1,prm%sum_N_sl Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & @@ -1649,10 +1650,12 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_P, & !< maximum obstacle strength criticalStress_S !< maximum obstacle strength + + v = 0.0_pReal + dv_dtau = 0.0_pReal + dv_dtauNS = 0.0_pReal + associate(prm => param(ph)) - v = 0.0_pReal - dv_dtau = 0.0_pReal - dv_dtauNS = 0.0_pReal do s = 1,prm%sum_N_sl if (abs(tau(s)) > tauThreshold(s)) then diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 4ea21d039..a3e4dd628 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -271,7 +271,7 @@ end subroutine thermal_forward !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- -module function thermal_T(ph,en) result(T) +pure module function thermal_T(ph,en) result(T) integer, intent(in) :: ph, en real(pReal) :: T