From 4c3d5806886a3c7e59d72fa3aed557ee333aad48 Mon Sep 17 00:00:00 2001
From: Achal H P <achal_hp@hotmail.com>
Date: Sun, 21 Jan 2024 11:19:32 +0530
Subject: [PATCH] Deltastate allocation

---
 src/grid/DAMASK_grid.f90                       |  2 +-
 src/phase_mechanical_plastic.f90               |  6 +++---
 src/phase_mechanical_plastic_phenopowerlaw.f90 | 18 +++++++++++-------
 3 files changed, 15 insertions(+), 11 deletions(-)

diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90
index 4395c1581..fd17644f2 100644
--- a/src/grid/DAMASK_grid.f90
+++ b/src/grid/DAMASK_grid.f90
@@ -132,7 +132,7 @@ program DAMASK_grid
 ! reading field paramters from numerics file and do sanity checks
   num_grid => config_numerics%get('grid', defaultVal=emptyDict)
   stagItMax  = num_grid%get_asInt('maxStaggeredIter',defaultVal=10)
-  maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3)
+  maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3)                                      ! Achal. To increase this?
 
   if (stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
   if (maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90
index 17bea1311..0179f5791 100644
--- a/src/phase_mechanical_plastic.f90
+++ b/src/phase_mechanical_plastic.f90
@@ -404,7 +404,7 @@ module function plastic_deltaState(ph, en) result(broken)
   broken = .false.
 
   select case (phase_plasticity(ph))
-    case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID)!, PLASTIC_PHENOPOWERLAW_ID)       !> Achal  
+    case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID, PLASTIC_PHENOPOWERLAW_ID)       !> Achal  
 
       Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                          phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
@@ -418,8 +418,8 @@ module function plastic_deltaState(ph, en) result(broken)
         case (PLASTIC_NONLOCAL_ID) plasticType
           call plastic_nonlocal_deltaState(Mp,ph,en)
 
-        !case (PLASTIC_PHENOPOWERLAW_ID) plasticType                                  !> Achal
-        !  call plastic_phenopowerlaw_deltaState(ph,en)                               !> Achal
+        case (PLASTIC_PHENOPOWERLAW_ID) plasticType                                  !> Achal
+          call plastic_phenopowerlaw_deltaState(ph,en)                               !> Achal
 
       end select plasticType
 
diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90
index c6d71c772..0365532cf 100644
--- a/src/phase_mechanical_plastic_phenopowerlaw.f90
+++ b/src/phase_mechanical_plastic_phenopowerlaw.f90
@@ -86,7 +86,7 @@ logical, dimension(:), allocatable :: myPlasticity
 integer :: &
   ph, i, o, &
   Nmembers, &
-  sizeState, sizeDotState, &
+  sizeState, sizeDotState, sizeDeltaState, &
   startIndex, endIndex
 integer,     dimension(:), allocatable :: &
   N_sl, N_tw
@@ -235,13 +235,17 @@ do ph = 1, phases%length
 ! allocate state arrays
   Nmembers = count(material_phaseID == ph)
   sizeDotState = size(['xi_sl   ','gamma_sl']) * prm%sum_N_sl &
-               + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw !&
-               !+ size(['xi_tw   ','f_twin  ']) * prm%sum_N_tw                !Achal
+               + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw &
+               + size(['xi_tw   ','f_twin  ']) * prm%sum_N_tw                !Achal
   sizeState = sizeDotState
 
-  !write(6,*)"size fn", size(['xi_sl   ','gamma_sl'])                          ! Achal Delete
+  write(6,*)"size fn", sizeDotState                          ! Achal Delete
 
-  call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
+  sizeDeltaState = size(['xi_sl   ','gamma_sl']) * prm%sum_N_sl &              !Achal
+                 + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw &
+                 + size(['xi_tw   ','f_twin  ']) * prm%sum_N_tw                !Achal
+
+  call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
   deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
 
 !--------------------------------------------------------------------------------------------------
@@ -385,7 +389,7 @@ associate(prm => param(ph), stt => state(ph), &
   !write(6,*)'deltaFp', deltaFp                                                                  ! delete this
   !write(6,*)'characteristicShearTwin', prm%gamma_char
   !write(6,*)'Schmid_twin',prm%P_sl
-  if(en==1) write(6,*)'maxF',maxval(stt%gamma_tw(:,en)/prm%gamma_char)
+  !if(en==1) write(6,*)'maxF',maxval(stt%gamma_tw(:,en)/prm%gamma_char)                           ! delete Achal
 
   sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
   xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
@@ -469,7 +473,7 @@ integer, intent(in)::&
   ph, &
   en
 
-deltastate(ph)%f_twin=0.0_pReal
+deltastate(ph)%f_twin=1.0_pReal
   
 
 end subroutine plastic_phenopowerlaw_deltaState