investigating the reason for the poor performance

src/plastic_phenopowerlaw.f90

@@ -145,7 +145,7 @@ subroutine plastic_phenopowerlaw_init
 
  integer(pInt) :: &
    maxNinstance, &
-   instance,p,j,k, o, i,&
+   instance,p, i,&
    NipcMyPhase, outputSize, &
    sizeState,sizeDotState, &
    startIndex, endIndex
@@ -445,7 +445,10 @@ end subroutine plastic_phenopowerlaw_init
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
-
+ use prec, only: &
+   dNeq0
+ use math, only: &
+   math_mul33xx33
  implicit none
  real(pReal), dimension(3,3), intent(out) :: &
    Lp                                                                                               !< plastic velocity gradient
@@ -460,9 +463,10 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
 
  integer(pInt) :: &
    i,k,l,m,n
- real(pReal), dimension(param(instance)%totalNslip) :: &
-   dgdot_dtauslip_pos,dgdot_dtauslip_neg, &
-   gdot_slip_pos,gdot_slip_neg
+ real(pReal) :: &
+   tau, &
+   gdot_slip, &
+   dgdot_dtau_slip
  real(pReal), dimension(param(instance)%totalNtwin) :: &
    gdot_twin,dgdot_dtautwin
 
@@ -474,14 +478,42 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
  Lp = 0.0_pReal
  dLp_dMp = 0.0_pReal
 
- call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
- slipSystems: do i = 1_pInt, prm%totalNslip
-   Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
-   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-     dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                      + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
-                      + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
- enddo slipSystems
+ do i = 1_pInt, prm%totalNslip
+
+   tau = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
+   if (tau > 1.0e-8_pReal) then
+
+     gdot_slip = prm%gdot0_slip &
+               * sign(abs(tau/stt%xi_slip(i,of))**prm%n_slip,  tau)
+     if (size(prm%nonSchmidCoeff) > 0_pInt)  gdot_slip = 0.5 * gdot_slip
+
+     dgdot_dtau_slip = gdot_slip*prm%n_slip/tau
+
+     Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip)*prm%Schmid_slip(1:3,1:3,i)
+
+     forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                        + dgdot_dtau_slip * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i)
+   endif
+
+   if (size(prm%nonSchmidCoeff) > 0_pInt) then
+     
+     tau = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i))
+     if (tau > 1.0e-8_pReal) then
+
+       gdot_slip = 0.5_pReal * prm%gdot0_slip &
+                 * sign(abs(tau/stt%xi_slip(i,of))**prm%n_slip,  tau)
+       
+       dgdot_dtau_slip = gdot_slip*prm%n_slip/tau
+
+       Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip)*prm%Schmid_slip(1:3,1:3,i)
+
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                          + dgdot_dtau_slip * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
+     endif
+   endif
+ enddo
 
  call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin)
  twinSystems: do i = 1_pInt, prm%totalNtwin
@@ -500,7 +532,8 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent
 !> @brief calculates the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
-
+ use math, only: &
+   math_mul33xx33
  implicit none
  real(pReal), dimension(3,3),  intent(in) :: &
    Mp                                                                                               !< Mandel stress
@@ -512,11 +545,10 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
    i
  real(pReal) :: &
    c_SlipSlip,c_TwinSlip,c_TwinTwin, &
-   xi_slip_sat_offset
+   xi_slip_sat_offset,tau
 
  real(pReal), dimension(param(instance)%totalNslip) :: &
-   left_SlipSlip,right_SlipSlip, &
-   gdot_slip_pos,gdot_slip_neg
+   left_SlipSlip,right_SlipSlip
 
  type(tParameters)         :: prm
  type(tPhenopowerlawState) :: dot,stt
@@ -540,9 +572,23 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
 
 !--------------------------------------------------------------------------------------------------
 ! shear rates
- call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg)
- dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
- dot%sumGamma(of) = sum(dot%gamma_slip(:,of))
+ do i = 1_pInt, prm%totalNslip
+   tau = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
+   if (tau > 1.0e-8_pReal) then
+     dot%gamma_slip(i,of) = abs(prm%gdot0_slip &
+                          * sign(abs(tau/stt%xi_slip(i,of))**prm%n_slip,  tau))
+     if (size(prm%nonSchmidCoeff) > 0_pInt)  dot%gamma_slip(i,of) = 0.5 * dot%gamma_slip(i,of)
+   endif
+
+   if (size(prm%nonSchmidCoeff) > 0_pInt) then
+     
+     tau = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i))
+     if (tau > 1.0e-8_pReal) &
+       dot%gamma_slip(i,of) = abs(prm%gdot0_slip *0.5_pReal &
+                            * sign(abs(tau/stt%xi_slip(i,of))**prm%n_slip,  tau))
+   endif
+ enddo
+
  call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of))
  if (stt%sumF(of) < 0.98_pReal) dot%sumF(of) = sum(dot%gamma_twin(:,of)/prm%gamma_twin_char)
 
@@ -701,7 +747,7 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
    postResults
 
  integer(pInt) :: &
-   o,c,i,j
+   o,c,i
  real(pReal), dimension(param(instance)%totalNslip) :: &
    gdot_slip_pos,gdot_slip_neg