Deltas updated in deltastate subroutine

src/phase_mechanical.f90

@@ -110,6 +110,15 @@ submodule(phase) mechanical
         dLp_dFi                                                                                         !< derivative of Lp with respect to Fi
     end subroutine plastic_LpAndItsTangents
 
+    module subroutine plastic_KinematicJump(ph, en, Jump_occurr,deltaFp)
+      integer, intent(in) :: &
+        ph, &
+        en
+      logical ,                     intent(out) :: &
+        Jump_occurr
+      real(pReal), dimension(3,3),  intent(out) :: &
+        deltaFp
+    end subroutine plastic_KinematicJump
 
     module subroutine plastic_isotropic_result(ph,group)
       integer,          intent(in) :: ph
@@ -376,6 +385,13 @@ end subroutine mechanical_result
 !--------------------------------------------------------------------------------------------------
 !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
 !> intermediate acceleration of the Newton-Raphson correction
+!> @modified by Satya and Achal
+!> check for detour i.e. if twinning is possible (we are not going ahead in Lp loop consistency) 
+!> checking by calling the deltaFp subroutine that should return 4 things
+!1) deltaFp= Cij (correspondance matrix) representing twinning shear and reorientation
+!2) -(twin volume fraction) for each twin system to make it harder for twinned material point to twin again by any twin system
+!3) jump in the last sampled volume fraction
+!4) logical true if twinning possible false if not occurring
 !--------------------------------------------------------------------------------------------------
 function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(status)
 
@@ -989,6 +1005,7 @@ end subroutine mechanical_forward
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
+!> @modified by Satya and Achal
 !--------------------------------------------------------------------------------------------------
 module function phase_mechanical_constitutive(Delta_t,co,ce) result(status)
 
@@ -1001,8 +1018,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(status)
   real(pREAL) :: &
     formerStep
   integer :: &
-    ph, en, sizeDotState
+    ph, en, sizeDotState, o, sd
-  logical :: todo
+  logical :: todo, FpJumped
   real(pREAL) :: stepFrac,step
   real(pREAL), dimension(3,3) :: &
     Fp0, &
@@ -1010,7 +1027,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(status)
     Lp0, &
     Li0, &
     F0, &
-    F
+    F, &
+    deltaFp
   real(pREAL), dimension(plasticState(material_ID_phase(co,ce))%sizeState) :: state0
 
 
@@ -1031,6 +1049,28 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(status)
   todo = .true.
   cutbackLooping: do while (todo)
 
+  ! achal calling Kinematic DeltaFp here
+  !** Starting to implement changes for accommodating large shear and reorientation caused by twinning**
+  !if(.not. FpJumped .and. NiterationStressLp>1) then                !Achal: Reason for this if statement?
+    call plastic_KinematicJump(ph, en, FpJumped,deltaFp)
+    !if(FpJumped) write(6,*) 'element jumped', en 
+    !if(FpJumped) then
+      !Fp0 = matmul(deltaFp,phase_mechanical_Fp0(ph)%data(1:3,1:3,en))
+      
+      o = plasticState(ph)%offsetDeltaState
+      sd = plasticState(ph)%sizeDeltaState
+      
+      !update current state by jump
+      !plasticState(ph)%state(o+1:o+sd,en) = plasticState(ph)%state(o+1:o+sd,en) & 
+                                             !+ plasticState(ph)%deltaState(o+1:o+sd,en)
+      
+      !store jumped state as initial value for next iteration
+      !plasticState(ph)%state0(o+1:o+sd,en) = plasticState(ph)%state(o+1:o+sd,en)
+
+      !store jumped state as initial value for for substate, partitioned state as well
+
+    !endif
+
     if (status == STATUS_OK) then
       formerStep = step
       stepFrac = stepFrac + step

src/phase_mechanical_plastic.f90

@@ -407,7 +407,8 @@ module function plastic_deltaState(ph, en) result(status)
   status = STATUS_OK
 
   select case (mechanical_plasticity_type(ph))
-    case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING)
+    case (MECHANICAL_PLASTICITY_PHENOPOWERLAW, MECHANICAL_PLASTICITY_NONLOCAL,& 
+            MECHANICAL_PLASTICITY_KINEHARDENING)
 
       Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                          phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -312,7 +312,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     startIndex =  endIndex + 1
     endIndex   =  endIndex + 1
     stt%frozen => plasticState(ph)%state(startIndex,:)
-    stt%frozen = 0.0_pReal-1.0_pReal 
+    !stt%frozen = 0.0_pReal-1.0_pReal 
     dlt%frozen => plasticState(ph)%deltaState(startIndex-o,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asReal('atol_gamma',defaultVal=1.0e-6_pReal)
   
@@ -424,6 +424,7 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
 
     call kinetics_sl(Mp,ph,en, dot_gamma_sl)
     call kinetics_tw(Mp,ph,en, dot_gamma_tw, fdot_twin)
+    
     dot_gamma_sl = abs(dot_gamma_sl)
     sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
 
@@ -444,6 +445,48 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
 
 end function phenopowerlaw_dotState
 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates (instantaneous) incremental change of microstructure
+!> Satya, Achal
+!--------------------------------------------------------------------------------------------------
+module subroutine plastic_phenopowerlaw_deltaState(ph,en)
+  implicit none
+  
+  integer, intent(in)::&
+    ph, &
+    en
+
+  logical :: &
+    twinJump
+
+  integer :: &
+    twin_var
+
+  real(pREAL), dimension(3,3) :: &
+    deltaFp
+  
+  ! These are updated at every strain increment. What should these initilizations be?
+  
+  associate(prm => param(ph), stt => state(ph), dlt => deltastate(ph))
+
+    twin_var = maxloc(stt%f_twin(:,en),dim=1)
+
+    !write(6,*)'twin_var',twin_var
+
+    call plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp)
+      !write(6,*)'twinJump',twinJump
+      if(twinJump) then
+        !write(6,*)'el',en
+        dlt%f_twin(:,en)     = 0.0_pReal - stt%f_twin(:,en)
+        dlt%fmc_twin(:,en)   = 0.0_pReal - stt%fmc_twin(:,en)
+        dlt%frozen(en)       = 1.0_pReal - stt%frozen(en)
+        write(6,*)'frozen',en,dlt%frozen(en),stt%frozen(en)
+        dlt%variant_twin(en) = twin_var !- stt%variant_twin(en)
+      endif
+  
+  end associate
+  
+  end subroutine plastic_phenopowerlaw_deltaState
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates instantaneous incremental change of kinematics and associated jump state
@@ -487,14 +530,16 @@ module subroutine plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp)
 
     Ability_Nucleation: if(stt%f_twin(twin_var,en)>(stt%fmc_twin(twin_var,en)+prm%checkstep(twin_var))) then
       stt%fmc_twin(twin_var,en) = stt%fmc_twin(twin_var,en)+prm%checkstep(twin_var)
-
       Success_Nucleation: if (random <= stt%f_twin(twin_var,en)) then          ! Instead of sum take max
+        write(6,*)'element twinned',en,'random',random,'variant',twin_var,'volume fraction', stt%f_twin(twin_var,en)
         twinJump = .true.
         deltaFp  = prm%CorrespondenceMatrix(:,:,twin_var)
-        dlt%f_twin(:,en)     = 0.0_pReal - stt%f_twin(:,en)
+        !write(6,*) twinJump
-        dlt%fmc_twin(:,en)   = 0.0_pReal - stt%fmc_twin(:,en)
+        !dlt%f_twin(:,en)     = 0.0_pReal - stt%f_twin(:,en)
-        dlt%frozen(en)       = 1.0_pReal - stt%frozen(en)
+        !dlt%fmc_twin(:,en)   = 0.0_pReal - stt%fmc_twin(:,en)
-        dlt%variant_twin(en) = twin_var !- stt%variant_twin(en)                ! Achal LHS is real, RHS integer   ! why this equation?
+        !dlt%frozen(en)       = 1.0_pReal - stt%frozen(en)
+        !write(6,*)'frozen',en,dlt%frozen(en),stt%frozen(en)
+        !dlt%variant_twin(en) = twin_var !- stt%variant_twin(en)                ! Achal LHS is real, RHS integer   ! why this equation?
 
       endif Success_Nucleation
 
@@ -504,30 +549,6 @@ module subroutine plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp)
 
 end subroutine plastic_kinematic_deltaFp
 
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates (instantaneous) incremental change of microstructure
-!> Satya, Achal
-!--------------------------------------------------------------------------------------------------
-module subroutine plastic_phenopowerlaw_deltaState(ph,en)
-  implicit none
-  
-  integer, intent(in)::&
-    ph, &
-    en
-  
-  ! These are updated at every strain increment. What should these initilizations be?
-  
-  associate(dlt => deltastate(ph))
-  
-    dlt%f_twin(:,en)   = 0.0_pReal
-    dlt%fmc_twin(:,en) = 0.0_pReal
-    !dlt%variant_twin(en) = 0.0_pReal                                  
-    !dlt%frozen(en) = 0.0_pReal
-  
-  end associate
-  
-  end subroutine plastic_phenopowerlaw_deltaState
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Write results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
@@ -559,6 +580,17 @@ module subroutine plastic_phenopowerlaw_result(ph,group)
           call result_writeDataset(stt%gamma_tw,group,trim(prm%output(ou)), &
                                    'twinning shear','1',prm%systems_tw)
 
+        case('f_twin')
+          call result_writeDataset(stt%f_twin,group,trim(prm%output(ou)), &
+                                   'volume fraction','1',prm%systems_tw)                          !Achal
+
+        case('variant_twin')
+          call result_writeDataset(stt%variant_twin,group,trim(prm%output(ou)), &
+                                    'twin variant','1')                                            !Achal
+
+        case('fbinary_twin')
+          call result_writeDataset(stt%frozen,group,trim(prm%output(ou)), &
+                                    'binary twin flag','1')                                       !Achal
       end select
 
     end do
@@ -643,14 +675,14 @@ pure subroutine kinetics_tw(Mp,ph,en,&
   integer :: i
 
 
-  associate(prm => param(ph), stt => state(ph))
+  associate(prm => param(ph), stt => state(ph), dlt => state(ph))
 
     tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)]
 
-    where(tau_tw > 0.0_pREAL .and. stt%frozen(en) < 0.9_pReal)                                      ! Achal
+    where(tau_tw > 0.0_pREAL .and. stt%frozen(en) < 0.9_pReal)                                      ! Achal .and. stt%frozen(en) < 0.9_pReal
-      dot_gamma_tw = (1.0_pREAL-sum(stt%gamma_tw(:,en)/prm%gamma_char)) &                           ! only twin in untwinned volume fraction
+      dot_gamma_tw = 0.0_pREAL !(1.0_pREAL-sum(stt%gamma_tw(:,en)/prm%gamma_char)) &                           ! only twin in untwinned volume fraction
-                   * prm%dot_gamma_0_tw*(tau_tw/stt%xi_tw(:,en))**prm%n_tw
+                   !* prm%dot_gamma_0_tw*(tau_tw/stt%xi_tw(:,en))**prm%n_tw
-      fdot_twin = (0.05_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char               ! Achal 0.05 is a constant
+      fdot_twin = (0.005_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char               ! Achal 0.005 is reference slip rate
     else where
       dot_gamma_tw = 0.0_pREAL
       fdot_twin = 0.0_pREAL