Twin variant selection works!

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -62,8 +62,7 @@ submodule(phase:plastic) phenopowerlaw
     xi_sl, &
     xi_tw, &
     gamma_sl, &
-      gamma_tw, &
+    gamma_tw
-      f_twin                                                         !Achal
 end type tPhenopowerlawState
 
 !--------------------------------------------------------------------------------------------------
@@ -269,8 +268,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   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
 
 !--------------------------------------------------------------------------------------------------
@@ -323,6 +320,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 end do slipSystems
 
 call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
+
 twinSystems: 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) &
@@ -354,8 +352,9 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
 real(pReal), dimension(param(ph)%sum_N_sl) :: &
   dot_gamma_sl_pos,dot_gamma_sl_neg, &
   right_SlipSlip
-  real(pReal), dimension(param(ph)%sum_N_tw) :: &
+
-    fdot_twin
+logical :: twinJump                                                                               !delete this
+real(pReal), dimension(3,3) :: deltaFp                                                            !delete this
 
 associate(prm => param(ph), stt => state(ph), &
           dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
@@ -365,11 +364,9 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
 
   call kinetics_sl(Mp,ph,en,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,fdot_twin)
+  call kinetics_tw(Mp,ph,en,dot_gamma_tw)
-    !write(6,*)'fdot_twin',fdot_twin
+  call plastic_kinematic_deltaFp(Mp,ph,en,twinJump,deltaFp)                                      ! delete this
-    if ( en==6 ) then
+  !write(6,*)'param(ph)%sum_N_sl', param(ph)%sum_N_sl                                        ! delete this
-      write(6,*)'dot_gamma_tw',dot_gamma_tw
-    end if
 
   sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
   xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
@@ -389,35 +386,39 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
 end function phenopowerlaw_dotState
 
 
-----------------------------------------------------------f_twin_nucl----------------------------------------
+!--------------------------------------------------------------------------------------------------
-> @brief calculates instantaneous incremental change of kinematics and associated jump state
+!> @brief calculates instantaneous incremental change of kinematics and associated jump state
---------------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-module subroutine plastic_kinematic_deltaFp(twinJump,deltaFp,Mp,ph,en)
+module subroutine plastic_kinematic_deltaFp(Mp,ph,en,twinJump,deltaFp)
 use math, only: &
   math_I3
-  logical ,                     intent(out) :: &
-    twinJump
-  real(pReal), dimension(3,3),  intent(out) :: &
-    deltaFp
 real(pReal), dimension(3,3), intent(in) :: &
   Mp  
 integer, intent(in) :: &
   ph, &
   en
+logical ,                     intent(out) :: &
+  twinJump
+real(pReal), dimension(3,3),  intent(out) :: &
+  deltaFp
+
 real(pReal)   :: &
   random
 integer :: &
   twin_var
 real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    dot_gamma_tw, fdot_twin
+  fdot_twin
+real(pReal), dimension(param(ph)%sum_N_tw) :: &
+  tau_tw
+integer :: i
 twinJump = .false.
 deltaFp  = math_I3
-  call kinetics_tw(Mp,ph,en,dot_gamma_tw,fdot_twin)
 
+tau_tw = [(math_tensordot(Mp,param(ph)%P_tw(1:3,1:3,i)),i=1,param(ph)%sum_N_tw)]
+twin_var = maxloc((0.05_pReal*(abs(tau_tw)/state(ph)%xi_tw(:,en))**param(ph)%n_tw)/param(ph)%gamma_char,dim=1)          ! This prints values from 1 to 6
+fdot_twin = (0.05_pReal*(abs(tau_tw)/state(ph)%xi_tw(:,en))**param(ph)%n_tw)/param(ph)%gamma_char                       ! This is sometimes >1
 
-  twin_var = maxloc(fdot_twin(:),dim=1)                            ! Correct this
+write(6,*) 'twin_var', twin_var
-
-  !write(6,*)'fdot_twin',fdot_twin
 
 call RANDOM_NUMBER(random)
 
@@ -592,7 +593,7 @@ end subroutine kinetics_sl
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_tw(Mp,ph,en,&
-                            dot_gamma_tw,fdot_twin,ddot_gamma_dtau_tw)
+                          dot_gamma_tw,ddot_gamma_dtau_tw)
 
 real(pReal), dimension(3,3),  intent(in) :: &
   Mp                                                                                              !< Mandel stress
@@ -601,7 +602,7 @@ pure subroutine kinetics_tw(Mp,ph,en,&
   en
 
 real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
-    dot_gamma_tw, fdot_twin
+  dot_gamma_tw
 real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
   ddot_gamma_dtau_tw
 
@@ -617,10 +618,8 @@ pure subroutine kinetics_tw(Mp,ph,en,&
   where(tau_tw > 0.0_pReal)
     dot_gamma_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_tw)/stt%xi_tw(:,en))**prm%n_tw
-      fdot_twin = (prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char
   else where
     dot_gamma_tw = 0.0_pReal
-      fdot_twin = 0.0_pReal
   end where
 
   if (present(ddot_gamma_dtau_tw)) then