notation as in paper

src/phase_mechanical_plastic_kinehardening.f90

@@ -18,6 +18,8 @@ submodule(phase:plastic) kinehardening
       h_inf_b, &                                                                                    !< asymptotic hardening rate of back stress for each slip
       xi_inf_f, &
       xi_inf_b
+    real(pReal),              allocatable, dimension(:,:) :: &
+      h_sl_sl                                                                                       !< slip resistance from slip activity
     real(pReal),              allocatable, dimension(:,:,:) :: &
       P, &
       nonSchmid_pos, &
@@ -122,7 +124,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
         prm%nonSchmid_pos  = prm%P
         prm%nonSchmid_neg  = prm%P
       endif
-      prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, &
+      prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, &
                                                    pl%get_as1dFloat('h_sl-sl'), &
                                                    phase%get_asString('lattice'))
 
@@ -158,7 +160,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
     else slipActive
       xi_0 = emptyRealArray
       allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
-      allocate(prm%interaction_SlipSlip(0,0))
+      allocate(prm%h_sl_sl(0,0))
     endif slipActive
 
 !--------------------------------------------------------------------------------------------------
@@ -288,7 +290,7 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
   sumGamma = sum(stt%accshear(:,en))
 
 
-  dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%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) &

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -46,10 +46,10 @@ submodule(phase:plastic) phenopowerlaw
 
   type :: tPhenopowerlawState
     real(pReal), pointer, dimension(:,:) :: &
-      xi_slip, &
+      xi_sl, &
-      xi_twin, &
+      xi_tw, &
-      gamma_slip, &
+      gamma_sl, &
-      gamma_twin
+      gamma_tw
   end type tPhenopowerlawState
 
 !--------------------------------------------------------------------------------------------------
@@ -235,30 +235,30 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
-    stt%xi_slip => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%xi_sl => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_slip =  spread(xi_0_sl, 2, Nmembers)
+    stt%xi_sl =  spread(xi_0_sl, 2, Nmembers)
-    dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
+    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_twin => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%xi_tw => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_twin =  spread(xi_0_tw, 2, Nmembers)
+    stt%xi_tw =  spread(xi_0_tw, 2, Nmembers)
-    dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
+    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_slip => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%gamma_sl => plasticState(ph)%state   (startIndex:endIndex,:)
-    dot%gamma_slip => plasticState(ph)%dotState(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'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw
-    stt%gamma_twin => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%gamma_tw => plasticState(ph)%state   (startIndex:endIndex,:)
-    dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
 
     end associate
@@ -293,28 +293,28 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
   integer :: &
     i,k,l,m,n
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    gdot_slip_pos,gdot_slip_neg, &
+    gdot_sl_pos,gdot_sl_neg, &
     dgdot_dtauslip_pos,dgdot_dtauslip_neg
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    gdot_twin,dgdot_dtautwin
+    gdot_tw,dgdot_dtautwin
 
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
   associate(prm => param(ph))
 
-  call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
+  call kinetics_sl(Mp,ph,en,gdot_sl_pos,gdot_sl_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
   slipSystems: do i = 1, prm%sum_N_sl
-    Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
+    Lp = Lp + (gdot_sl_pos(i)+gdot_sl_neg(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) &
                        + dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) &
                        + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
   enddo slipSystems
 
-  call kinetics_twin(Mp,ph,en,gdot_twin,dgdot_dtautwin)
+  call kinetics_tw(Mp,ph,en,gdot_tw,dgdot_dtautwin)
   twinSystems: do i = 1, prm%sum_N_tw
-    Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
+    Lp = Lp + gdot_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) &
                        + dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
@@ -337,46 +337,33 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
     en
 
   real(pReal) :: &
-    c_SlipSlip,c_TwinSlip,c_TwinTwin, &
+    xi_sl_sat_offset,&
-    xi_slip_sat_offset,&
+    sumF
-    sumGamma,sumF
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    left_SlipSlip,right_SlipSlip, &
+    gdot_sl_pos,gdot_sl_neg, &
-    gdot_slip_pos,gdot_slip_neg
+    right_SlipSlip
 
-  associate(prm => param(ph), stt => state(ph), &
-  dot => dotState(ph))
 
-  sumGamma = sum(stt%gamma_slip(:,en))
+  associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
-  sumF     = sum(stt%gamma_twin(:,en)/prm%gamma_char)
 
-!--------------------------------------------------------------------------------------------------
+    call kinetics_sl(Mp,ph,en,gdot_sl_pos,gdot_sl_neg)
-! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
+    dot%gamma_sl(:,en) = abs(gdot_sl_pos+gdot_sl_neg)
-  c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2)
+    call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
-  c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3
-  c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4
 
-!--------------------------------------------------------------------------------------------------
-!  calculate left and right vectors
-  left_SlipSlip  = 1.0_pReal + prm%h_int
-  xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
-  right_SlipSlip = sign(abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl, &
-                        1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset))
 
-!--------------------------------------------------------------------------------------------------
+    sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
-! shear rates
+    xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
-  call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg)
+    right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset)) **prm%a_sl, &
-  dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg)
+                          1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
-  call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en))
 
-!--------------------------------------------------------------------------------------------------
+    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) &
-! hardening
+                      * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
-  dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * &
+                    + matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
-                      matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) &
+
-                    + matmul(prm%h_sl_tw,dot%gamma_twin(:,en))
+    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_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) &
-                    + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en))
   end associate
 
 end subroutine phenopowerlaw_dotState
@@ -397,17 +384,17 @@ module subroutine plastic_phenopowerlaw_results(ph,group)
     select case(trim(prm%output(o)))
 
       case('xi_sl')
-        if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_slip,group,trim(prm%output(o)), &
+        if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_sl,group,trim(prm%output(o)), &
                                                      'resistance against plastic slip','Pa')
       case('gamma_sl')
-        if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_slip,group,trim(prm%output(o)), &
+        if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_sl,group,trim(prm%output(o)), &
                                                      'plastic shear','1')
 
       case('xi_tw')
-        if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_twin,group,trim(prm%output(o)), &
+        if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_tw,group,trim(prm%output(o)), &
                                                      'resistance against twinning','Pa')
       case('gamma_tw')
-        if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_twin,group,trim(prm%output(o)), &
+        if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_tw,group,trim(prm%output(o)), &
                                                      'twinning shear','1')
 
     end select
@@ -424,8 +411,8 @@ end subroutine plastic_phenopowerlaw_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_slip(Mp,ph,en, &
+pure subroutine kinetics_sl(Mp,ph,en, &
-                              gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
+                            gdot_sl_pos,gdot_sl_neg,dgdot_dtau_sl_pos,dgdot_dtau_sl_neg)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -434,56 +421,56 @@ pure subroutine kinetics_slip(Mp,ph,en, &
     en
 
   real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl) :: &
-    gdot_slip_pos, &
+    gdot_sl_pos, &
-    gdot_slip_neg
+    gdot_sl_neg
   real(pReal),                  intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
-    dgdot_dtau_slip_pos, &
+    dgdot_dtau_sl_pos, &
-    dgdot_dtau_slip_neg
+    dgdot_dtau_sl_neg
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    tau_slip_pos, &
+    tau_sl_pos, &
-    tau_slip_neg
+    tau_sl_neg
   integer :: i
 
   associate(prm => param(ph), stt => state(ph))
 
   do i = 1, prm%sum_N_sl
-    tau_slip_pos(i) =       math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i))
+    tau_sl_pos(i) =       math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i))
-    tau_slip_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
+    tau_sl_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
                           0.0_pReal, prm%nonSchmidActive)
   enddo
 
-  where(dNeq0(tau_slip_pos))
+  where(dNeq0(tau_sl_pos))
-    gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &          ! 1/2 if non-Schmid active
+    gdot_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_slip_pos/stt%xi_slip(:,en))**prm%n_sl,  tau_slip_pos)
+                * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl,  tau_sl_pos)
   else where
-    gdot_slip_pos = 0.0_pReal
+    gdot_sl_pos = 0.0_pReal
   end where
 
-  where(dNeq0(tau_slip_neg))
+  where(dNeq0(tau_sl_neg))
-    gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal &                                                ! only used if non-Schmid active, always 1/2
+    gdot_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal &                                                  ! only used if non-Schmid active, always 1/2
-                  * sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl,  tau_slip_neg)
+                * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl,  tau_sl_neg)
   else where
-    gdot_slip_neg = 0.0_pReal
+    gdot_sl_neg = 0.0_pReal
   end where
 
-  if (present(dgdot_dtau_slip_pos)) then
+  if (present(dgdot_dtau_sl_pos)) then
-    where(dNeq0(gdot_slip_pos))
+    where(dNeq0(gdot_sl_pos))
-      dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos
+      dgdot_dtau_sl_pos = gdot_sl_pos*prm%n_sl/tau_sl_pos
     else where
-      dgdot_dtau_slip_pos = 0.0_pReal
+      dgdot_dtau_sl_pos = 0.0_pReal
     end where
   endif
-  if (present(dgdot_dtau_slip_neg)) then
+  if (present(dgdot_dtau_sl_neg)) then
-    where(dNeq0(gdot_slip_neg))
+    where(dNeq0(gdot_sl_neg))
-      dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg
+      dgdot_dtau_sl_neg = gdot_sl_neg*prm%n_sl/tau_sl_neg
     else where
-      dgdot_dtau_slip_neg = 0.0_pReal
+      dgdot_dtau_sl_neg = 0.0_pReal
     end where
   endif
   end associate
 
-end subroutine kinetics_slip
+end subroutine kinetics_sl
 
 
 !--------------------------------------------------------------------------------------------------
@@ -493,8 +480,8 @@ end subroutine kinetics_slip
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_twin(Mp,ph,en,&
+pure subroutine kinetics_tw(Mp,ph,en,&
-                              gdot_twin,dgdot_dtau_twin)
+                            gdot_tw,dgdot_dtau_tw)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -503,37 +490,37 @@ pure subroutine kinetics_twin(Mp,ph,en,&
     en
 
   real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
-    gdot_twin
+    gdot_tw
   real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
-    dgdot_dtau_twin
+    dgdot_dtau_tw
 
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    tau_twin
+    tau_tw
   integer :: i
 
   associate(prm => param(ph), stt => state(ph))
 
   do i = 1, prm%sum_N_tw
-    tau_twin(i)  = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
+    tau_tw(i)  = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
   enddo
 
-  where(tau_twin > 0.0_pReal)
+  where(tau_tw > 0.0_pReal)
-    gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) &                              ! only twin in untwinned volume fraction
+    gdot_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_twin)/stt%xi_twin(:,en))**prm%n_tw
+              * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw
   else where
-    gdot_twin = 0.0_pReal
+    gdot_tw = 0.0_pReal
   end where
 
-  if (present(dgdot_dtau_twin)) then
+  if (present(dgdot_dtau_tw)) then
-    where(dNeq0(gdot_twin))
+    where(dNeq0(gdot_tw))
-      dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin
+      dgdot_dtau_tw = gdot_tw*prm%n_tw/tau_tw
     else where
-      dgdot_dtau_twin = 0.0_pReal
+      dgdot_dtau_tw = 0.0_pReal
     end where
   endif
 
   end associate
 
-end subroutine kinetics_twin
+end subroutine kinetics_tw
 
 end submodule phenopowerlaw