no need to store dot state

saves memory and flow is easier to understand

src/phase_mechanical_plastic.f90

@@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic
     end subroutine nonlocal_LpAndItsTangent
 
 
-    module subroutine isotropic_dotState(Mp,ph,en)
+    module function isotropic_dotState(Mp,ph,en) result(dotState)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
         en
-    end subroutine isotropic_dotState
+      real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+        dotState
+    end function isotropic_dotState
 
-    module subroutine phenopowerlaw_dotState(Mp,ph,en)
+    module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
         en
-    end subroutine phenopowerlaw_dotState
+      real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+        dotState
+    end function phenopowerlaw_dotState
 
-    module subroutine plastic_kinehardening_dotState(Mp,ph,en)
+    module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
         en
-    end subroutine plastic_kinehardening_dotState
+      real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+        dotState
+    end function plastic_kinehardening_dotState
 
     module subroutine dislotwin_dotState(Mp,T,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
@@ -144,15 +150,15 @@ submodule(phase:mechanical) plastic
         en
     end subroutine dislotwin_dotState
 
-    module subroutine dislotungsten_dotState(Mp,T,ph,en)
+    module function dislotungsten_dotState(Mp,ph,en) result(dotState)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
-      real(pReal),                  intent(in) :: &
-        T
       integer,                      intent(in) :: &
         ph, &
         en
-    end subroutine dislotungsten_dotState
+      real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+        dotState
+    end function dislotungsten_dotState
 
     module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el)
       real(pReal), dimension(3,3), intent(in) :: &
@@ -318,27 +324,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
     plasticType: select case (phase_plasticity(ph))
 
       case (PLASTIC_ISOTROPIC_ID) plasticType
-        call isotropic_dotState(Mp,ph,en)
+        dotState = isotropic_dotState(Mp,ph,en)
 
       case (PLASTIC_PHENOPOWERLAW_ID) plasticType
-        call phenopowerlaw_dotState(Mp,ph,en)
+        dotState = phenopowerlaw_dotState(Mp,ph,en)
 
       case (PLASTIC_KINEHARDENING_ID) plasticType
-        call plastic_kinehardening_dotState(Mp,ph,en)
+        dotState = plastic_kinehardening_dotState(Mp,ph,en)
 
       case (PLASTIC_DISLOTWIN_ID) plasticType
         call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
+        dotState = plasticState(ph)%dotState(:,en)
 
       case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
-        call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
+        dotState = dislotungsten_dotState(Mp,ph,en)
 
       case (PLASTIC_NONLOCAL_ID) plasticType
         call nonlocal_dotState(Mp,subdt,ph,en,ip,el)
+        dotState = plasticState(ph)%dotState(:,en)
+
     end select plasticType
   end if
 
-  dotState = plasticState(ph)%dotState(:,en)
-
 end function plastic_dotState
 
 

src/phase_mechanical_plastic_dislotungsten.f90

@@ -43,6 +43,13 @@ submodule(phase:plastic) dislotungsten
       systems_sl
   end type tParameters                                                                              !< container type for internal constitutive parameters
 
+  type :: tIndexDotState
+    integer, dimension(2) :: &
+      rho_mob, &
+      rho_dip, &
+      gamma_sl
+  end type tIndexDotState
+
   type :: tDislotungstenState
     real(pReal), dimension(:,:), pointer :: &
       rho_mob, &
@@ -59,9 +66,8 @@ submodule(phase:plastic) dislotungsten
 !--------------------------------------------------------------------------------------------------
 ! containers for parameters and state
   type(tParameters),                  allocatable, dimension(:) :: param
-  type(tDisloTungstenState),          allocatable, dimension(:) :: &
+  type(tIndexDotState),               allocatable, dimension(:) :: indexDotState
-    dotState, &
+  type(tDisloTungstenState),          allocatable, dimension(:) :: state
-    state
   type(tDisloTungstenDependentState), allocatable, dimension(:) :: dependentState
 
 contains
@@ -103,18 +109,17 @@ module function plastic_dislotungsten_init() result(myPlasticity)
   print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
   print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
 
-
   phases => config_material%get('phase')
   allocate(param(phases%length))
+  allocate(indexDotState(phases%length))
   allocate(state(phases%length))
-  allocate(dotState(phases%length))
   allocate(dependentState(phases%length))
 
-
   do ph = 1, phases%length
     if (.not. myPlasticity(ph)) cycle
 
-    associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
+    associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
+              idx_dot => indexDotState(ph))
 
     phase => phases%get(ph)
     mech  => phase%get('mechanical')
@@ -214,28 +219,29 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     sizeState = sizeDotState
 
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
+    deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
+    idx_dot%rho_mob = [startIndex,endIndex]
     stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
     stt%rho_mob =  spread(rho_mob_0,2,Nmembers)
-    dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
     if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
+    idx_dot%rho_dip = [startIndex,endIndex]
     stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
     stt%rho_dip =  spread(rho_dip_0,2,Nmembers)
-    dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
+    idx_dot%gamma_sl = [startIndex,endIndex]
     stt%gamma_sl => plasticState(ph)%state(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'
 
@@ -300,15 +306,15 @@ end subroutine dislotungsten_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotungsten_dotState(Mp,T,ph,en)
+module function dislotungsten_dotState(Mp,ph,en) result(dotState)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
-  real(pReal),                  intent(in) :: &
-    T                                                                                               !< temperature
   integer,                      intent(in) :: &
     ph, &
     en
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+    dotState
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_pos, dot_gamma_neg,&
@@ -319,17 +325,22 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
     dot_rho_dip_climb, &
     d_hat
   real(pReal) :: &
-    mu
+    mu, T
 
-  associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
+
+  associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
+            dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), &
+            dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
+            dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)))
 
     mu = elastic_mu(ph,en)
+    T = thermal_T(ph,en)
 
     call kinetics(Mp,T,ph,en,&
                   dot_gamma_pos,dot_gamma_neg, &
                   tau_pos_out = tau_pos,tau_neg_out = tau_neg)
 
-    dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
+    dot_gamma_sl = abs(dot_gamma_pos+dot_gamma_neg)
 
     where(dEq0((tau_pos+tau_neg)*0.5_pReal))
       dot_rho_dip_formation = 0.0_pReal
@@ -338,7 +349,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
       d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
                         prm%d_caron, &                                                              ! lower limit
                         dst%Lambda_sl(:,en))                                                        ! upper limit
-      dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
+      dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot_gamma_sl/prm%b_sl, &
                                     0.0_pReal, &
                                     prm%dipoleformation)
       v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) &
@@ -346,16 +357,16 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
       dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron)                    ! ToDo: Discuss with Franz: Stress dependency?
     end where
 
-    dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) &                         ! multiplication
+    dot_rho_mob = dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) &                                     ! multiplication
                       - dot_rho_dip_formation &
-                      - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en)       ! Spontaneous annihilation of 2 edges
+                - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot_gamma_sl                   ! Spontaneous annihilation of 2 edges
-    dot%rho_dip(:,en) = dot_rho_dip_formation &
+    dot_rho_dip = dot_rho_dip_formation &
-                      - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) &     ! Spontaneous annihilation of an edge with a dipole
+                - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot_gamma_sl &                 ! Spontaneous annihilation of an edge with a dipole
                       - dot_rho_dip_climb
 
   end associate
 
-end subroutine dislotungsten_dotState
+end function dislotungsten_dotState
 
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_isotropic.f90

@@ -37,9 +37,7 @@ submodule(phase:plastic) isotropic
 !--------------------------------------------------------------------------------------------------
 ! containers for parameters and state
   type(tParameters),     allocatable, dimension(:) :: param
-  type(tIsotropicState), allocatable, dimension(:) :: &
+  type(tIsotropicState), allocatable, dimension(:) :: state
-    dotState, &
-    state
 
 contains
 
@@ -77,12 +75,11 @@ module function plastic_isotropic_init() result(myPlasticity)
   phases => config_material%get('phase')
   allocate(param(phases%length))
   allocate(state(phases%length))
-  allocate(dotState(phases%length))
 
   do ph = 1, phases%length
     if(.not. myPlasticity(ph)) cycle
 
-    associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
+    associate(prm => param(ph), stt => state(ph))
 
     phase => phases%get(ph)
     mech => phase%get('mechanical')
@@ -125,12 +122,12 @@ module function plastic_isotropic_init() result(myPlasticity)
     sizeState = sizeDotState
 
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
+    deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
-    stt%xi  => plasticState(ph)%state   (1,:)
+    stt%xi => plasticState(ph)%state(1,:)
     stt%xi = xi_0
-    dot%xi  => plasticState(ph)%dotState(1,:)
     plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
 
@@ -178,7 +175,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
     norm_Mp_dev = sqrt(squarenorm_Mp_dev)
 
     if (norm_Mp_dev > 0.0_pReal) then
-      dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
+      dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en)))**prm%n
 
       Lp = dot_gamma * Mp_dev/norm_Mp_dev
       forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -242,27 +239,26 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine isotropic_dotState(Mp,ph,en)
+module function isotropic_dotState(Mp,ph,en) result(dotState)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
     en
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+    dotState
 
   real(pReal) :: &
     dot_gamma, &                                                                                    !< strainrate
     xi_inf_star, &                                                                                  !< saturation xi
     norm_Mp                                                                                         !< norm of the (deviatoric) Mandel stress
 
-  associate(prm => param(ph), stt => state(ph), &
+  associate(prm => param(ph), stt => state(ph), dot_xi => dotState(1))
-   dot => dotState(ph))
 
-  if (prm%dilatation) then
+    norm_Mp = merge(sqrt(math_tensordot(Mp,Mp)), &
-    norm_Mp = sqrt(math_tensordot(Mp,Mp))
+                    sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))), &
-  else
+                    prm%dilatation)
-    norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
-  end if
 
   dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
 
@@ -274,16 +270,16 @@ module subroutine isotropic_dotState(Mp,ph,en)
                   + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
                   / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
     end if
-    dot%xi(en) = dot_gamma &
+    dot_xi = dot_gamma &
            * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
            * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
   else
-    dot%xi(en) = 0.0_pReal
+    dot_xi = 0.0_pReal
   end if
 
   end associate
 
-end subroutine isotropic_dotState
+end function isotropic_dotState
 
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_kinehardening.f90

@@ -34,6 +34,13 @@ submodule(phase:plastic) kinehardening
       systems_sl
   end type tParameters
 
+  type :: tIndexDotState
+    integer, dimension(2) :: &
+      xi, &
+      chi, &
+      gamma
+  end type tIndexDotState
+
   type :: tKinehardeningState
     real(pReal), pointer, dimension(:,:) :: &
       xi, &                                                                                         !< resistance against plastic slip
@@ -47,10 +54,8 @@ submodule(phase:plastic) kinehardening
 !--------------------------------------------------------------------------------------------------
 ! containers for parameters and state
   type(tParameters),         allocatable, dimension(:) :: param
-  type(tKinehardeningState), allocatable, dimension(:) :: &
+  type(tIndexDotState),      allocatable, dimension(:) :: indexDotState
-    dotState, &
+  type(tKinehardeningState), allocatable, dimension(:) :: state, deltaState
-    deltaState, &
-    state
 
 contains
 
@@ -91,15 +96,16 @@ module function plastic_kinehardening_init() result(myPlasticity)
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
+  allocate(indexDotState(phases%length))
   allocate(state(phases%length))
-  allocate(dotState(phases%length))
   allocate(deltaState(phases%length))
 
 
   do ph = 1, phases%length
     if (.not. myPlasticity(ph)) cycle
 
-    associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
+    associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), &
+              idx_dot => indexDotState(ph))
 
     phase => phases%get(ph)
     mech => phase%get('mechanical')
@@ -173,27 +179,28 @@ module function plastic_kinehardening_init() result(myPlasticity)
     sizeState = sizeDotState + sizeDeltaState
 
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
+    deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
-    stt%xi => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%xi = [startIndex,endIndex]
+    stt%xi => plasticState(ph)%state(startIndex:endIndex,:)
     stt%xi = spread(xi_0, 2, Nmembers)
-    dot%xi => 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_sl
-    stt%chi => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%chi = [startIndex,endIndex]
-    dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:)
+    stt%chi => plasticState(ph)%state(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 => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%gamma = [startIndex,endIndex]
-    dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:)
+    stt%gamma => plasticState(ph)%state(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'
 
@@ -270,13 +277,15 @@ end subroutine kinehardening_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_kinehardening_dotState(Mp,ph,en)
+module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
     en
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+    dotState
 
   real(pReal) :: &
     sumGamma
@@ -284,29 +293,32 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
     dot_gamma_pos,dot_gamma_neg
 
 
-  associate(prm => param(ph), stt => state(ph),dot => dotState(ph))
+  associate(prm => param(ph), stt => state(ph), &
+            dot_xi => dotState(IndexDotState(ph)%xi(1):IndexDotState(ph)%xi(2)),&
+            dot_chi => dotState(IndexDotState(ph)%chi(1):IndexDotState(ph)%chi(2)),&
+            dot_gamma => dotState(IndexDotState(ph)%gamma(1):IndexDotState(ph)%gamma(2)))
 
     call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
-    dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
+    dot_gamma = abs(dot_gamma_pos+dot_gamma_neg)
     sumGamma = sum(stt%gamma(:,en))
 
 
-    dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) &
+    dot_xi = matmul(prm%h_sl_sl,dot_gamma) &
                    * (  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) &
                      )
 
-    dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * &
+    dot_chi = stt%sgn_gamma(:,en)*dot_gamma &
-             ( prm%h_inf_b + &
+            * ( prm%h_inf_b &
-               (prm%h_0_b - prm%h_inf_b &
+               + (prm%h_0_b - prm%h_inf_b &
                  + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))&
                ) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) &
              )
 
   end associate
 
-end subroutine plastic_kinehardening_dotState
+end function plastic_kinehardening_dotState
 
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -47,6 +47,14 @@ submodule(phase:plastic) phenopowerlaw
       systems_tw
   end type tParameters
 
+  type :: tIndexDotState
+    integer, dimension(2) :: &
+      xi_sl, &
+      xi_tw, &
+      gamma_sl, &
+      gamma_tw
+  end type tIndexDotState
+
   type :: tPhenopowerlawState
     real(pReal), pointer, dimension(:,:) :: &
       xi_sl, &
@@ -56,11 +64,10 @@ submodule(phase:plastic) phenopowerlaw
   end type tPhenopowerlawState
 
 !--------------------------------------------------------------------------------------------------
-! containers for parameters and state
+! containers for parameters, dot state index,  and state
   type(tParameters),         allocatable, dimension(:) :: param
-  type(tPhenopowerlawState), allocatable, dimension(:) :: &
+  type(tIndexDotState),      allocatable, dimension(:) :: indexDotState
-    dotState, &
+  type(tPhenopowerlawState), allocatable, dimension(:) :: state
-    state
 
 contains
 
@@ -101,13 +108,14 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
+  allocate(indexDotState(phases%length))
   allocate(state(phases%length))
-  allocate(dotState(phases%length))
 
   do ph = 1, phases%length
     if (.not. myPlasticity(ph)) cycle
 
-    associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
+    associate(prm => param(ph), stt => state(ph), &
+              idx_dot => indexDotState(ph))
 
     phase => phases%get(ph)
     mech => phase%get('mechanical')
@@ -224,37 +232,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
                  + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw
     sizeState = sizeDotState
 
-
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
+    deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
-    stt%xi_sl => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%xi_sl = [startIndex,endIndex]
+    stt%xi_sl => plasticState(ph)%state(startIndex:endIndex,:)
     stt%xi_sl =  spread(xi_0_sl, 2, Nmembers)
-    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_tw => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%xi_tw = [startIndex,endIndex]
+    stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:)
     stt%xi_tw =  spread(xi_0_tw, 2, Nmembers)
-    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_sl => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%gamma_sl = [startIndex,endIndex]
-    dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
+    stt%gamma_sl => plasticState(ph)%state(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_tw => plasticState(ph)%state   (startIndex:endIndex,:)
+    idx_dot%gamma_tw = [startIndex,endIndex]
-    dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
+    stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
 
     end associate
@@ -324,13 +332,15 @@ end subroutine phenopowerlaw_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine phenopowerlaw_dotState(Mp,ph,en)
+module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
     en
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
+    dotState
 
   real(pReal) :: &
     xi_sl_sat_offset,&
@@ -340,28 +350,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
     right_SlipSlip
 
 
-  associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
+  associate(prm => param(ph), stt => state(ph), &
+            dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
+            dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), &
+            dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
+            dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)))
 
     call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
-    dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
+    dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
-    call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
+    call kinetics_tw(Mp,ph,en,dot_gamma_tw)
 
     sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
     xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
     right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
                           1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
 
-    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) &
+    dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
-                      * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
+                * matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) &
-                    + matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
+              + matmul(prm%h_sl_tw,dot_gamma_tw)
 
-    dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
+    dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
-                      * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) &
+                * matmul(prm%h_tw_sl,dot_gamma_sl) &
-                    + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en))
+              + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw)
 
   end associate
 
-end subroutine phenopowerlaw_dotState
+end function phenopowerlaw_dotState
 
 
 !--------------------------------------------------------------------------------------------------