diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index f4e400c38..6e892c2fb 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -473,15 +473,15 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) !** 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(en==15) write(6,*)'deltaFp',deltaFp !Achal Delete - + !if(en==15) write(6,*)'deltaFp',deltaFp !Achal Delete + if(en==15) write(6,*)'FpJumped',FpJumped !converged = .true. means no more iteration if(FpJumped) then !crystallite_converged(ipc,ip,el) = .true. !> See "phase_mechanical_constitutive" and "homogenization_mechanical_response" !crystallite_todo(ipc,ip,el) = .false. !> Can't find this ! _converged = .not. broken - !subFp0 = matmul(deltaFp,subFp0) + Fp_new = matmul(deltaFp,subFp0) ! subFp0 is input need to change "phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal)" !plasticState(ph)%state() @@ -1031,7 +1031,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) formerSubStep integer :: & ph, en, sizeDotState - logical :: todo + logical :: todo, FpJumped !Achal real(pReal) :: subFrac,subStep real(pReal), dimension(3,3) :: & subFp0, & @@ -1039,7 +1039,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) subLp0, & subLi0, & subF0, & - subF + subF, & + deltaFp !Achal real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0 @@ -1060,6 +1061,18 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) todo = .true. cutbackLooping: do while (todo) + !matmul(deltaFp,subFp0) + + ! 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) then + subFp0 = matmul(deltaFp,phase_mechanical_Fp0(ph)%data(1:3,1:3,en)) + endif + + !endif + if (converged_) then formerSubStep = subStep subFrac = subFrac + subStep diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index d58d0d812..a1931daf7 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -82,8 +82,11 @@ type :: tPhenopowerlawState xi_tw, & gamma_sl, & gamma_tw, & - f_twin, & - fmc_twin !< Achal, To control sampling frequency + f_twin, & !< Twin volume fraction + fmc_twin !< Achal, To control sampling frequency + real(pReal), pointer, dimension(:) :: & + variant_twin, & + frozen end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- @@ -261,11 +264,13 @@ do ph = 1, phases%length sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw','f_twin ']) * prm%sum_N_tw !Achal - sizeDeltaState = size(['f_twin ','fmc_twin']) * prm%sum_N_tw !Achal + sizeDeltaState = size(['f_twin ','fmc_twin']) * prm%sum_N_tw & !Achal + + size(['variant_twin','frozen ']) * prm%sum_N_tw sizeState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw & - + size(['f_twin ','fmc_twin']) * prm%sum_N_tw !Achal + + size(['f_twin ','fmc_twin']) * prm%sum_N_tw & + + size(['variant_twin','frozen ']) * prm%sum_N_tw !Achal @@ -326,11 +331,26 @@ do ph = 1, phases%length dlt%f_twin => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) ! Achal plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_tw + stt%frozen => plasticState(ph)%state(startIndex,:) + stt%frozen = 0.0_pReal-1.0_pReal + dlt%frozen => plasticState(ph)%deltaState(startIndex-o,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%fmc_twin => plasticState(ph)%state(startIndex:endIndex,:) dlt%fmc_twin => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_tw + stt%variant_twin => plasticState(ph)%state(startIndex,:) + dlt%variant_twin => plasticState(ph)%deltaState(startIndex-o,:) + plasticState(ph)%atol(startIndex:endIndex) = 0.0_pReal + + write(6,*)"offset", o ! Achal Delete @@ -531,7 +551,9 @@ associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dlt => deltas twinJump = .true. deltaFp = prm%CorrespondanceMatrix(:,:,twin_var) dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en) - dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en) + dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en) + dlt%frozen(en) = 1.0_pReal - stt%frozen(en) + dlt%variant_twin(en) = twin_var - stt%variant_twin(en) end if Success_Nucleation endif Ability_Nucleation @@ -553,6 +575,8 @@ integer, intent(in)::& deltastate(ph)%f_twin = 0.0_pReal deltastate(ph)%fmc_twin = 0.0_pReal +deltastate(ph)%variant_twin = 0 +deltastate(ph)%frozen = 0.0_pReal end subroutine plastic_phenopowerlaw_deltaState @@ -589,6 +613,17 @@ associate(prm => param(ph), stt => state(ph)) call results_writeDataset(stt%gamma_tw,group,trim(prm%output(ou)), & 'twinning shear','1',prm%systems_tw) + case('f_twin') + call results_writeDataset(stt%f_twin,group,trim(prm%output(ou)), & + 'volume fraction','1',prm%systems_tw) + + case('variant_twin') + call results_writeDataset(stt%variant_twin,group,trim(prm%output(ou)), & + 'twin variant','1',prm%systems_tw) + + case('fbinary_twin') + call results_writeDataset(stt%frozen,group,trim(prm%output(ou)), & + 'binary twin flag','1',prm%systems_tw) end select end do