use data from other physics directly

more clear code, simplified interfaces

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

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, &
+    module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
-                                                       Mp,Temperature,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
 

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
 

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,&
+    f_unrotated,StressRatio_p,&
-     E_kB_T, &
+    E_kB_T, &
-     ddot_gamma_dtau, &
+    ddot_gamma_dtau, &
-     tau
+    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)
+  associate(prm => param(ph), stt => state(ph))
-  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
 
-  call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
+    f_unrotated = 1.0_pReal &
-  twinContibution: do i = 1, prm%sum_N_tw
+                - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-    Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
+                - sum(stt%f_tr(1:prm%sum_N_tr,en))
-    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
 
-  call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
+    call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
-  transContibution: do i = 1, prm%sum_N_tr
+    slipContribution: do i = 1, prm%sum_N_sl
-    Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
+      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) &
+      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) &
+        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)
+                         + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
-  end do transContibution
+    end do slipContribution
 
-  Lp      = Lp      * f_unrotated
+    call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
-  dLp_dMp = dLp_dMp * f_unrotated
+    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)
+    Lp      = Lp      * f_unrotated
-    call math_eigh33(eigValues,eigVectors,Mp)                                                       ! is Mp symmetric by design?
+    dLp_dMp = dLp_dMp * f_unrotated
 
-    do i = 1,6
+    shearBandingContribution: if (dNeq0(prm%v_sb)) then
-      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)
 
-      significantShearBandStress: if (abs(tau) > tol_math_check) then
+      E_kB_T = prm%E_sb/(K_B*T)
-        StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
+      call math_eigh33(eigValues,eigVectors,Mp)                                                       ! is Mp symmetric by design?
-        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)
 
-        Lp = Lp + dot_gamma_sb * P_sb
+      do i = 1,6
-        forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
-          dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                                      matmul(eigVectors,sb_mComposition(1:3,i)))
-                           + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
+        tau = math_tensordot(Mp,P_sb)
-      end if significantShearBandStress
-    end do
 
-  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
 

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

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