diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 94c7a61c6..0179f5791 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -204,12 +204,24 @@ submodule(phase:mechanical) plastic en end subroutine plastic_nonlocal_deltaState - !module subroutine plastic_kinematic_deltaFp(Mp,ph,en,twinJump) - ! implicit none - ! - !contains + module subroutine plastic_phenopowerlaw_deltaState(ph,en) !< Achal + integer, intent(in)::& + ph, & + en + end subroutine plastic_phenopowerlaw_deltaState + + module subroutine plastic_kinematic_deltaFp(Mp,ph,en,twinJump,deltaFp) !< Achal + 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 - !end module subroutine + end subroutine plastic_kinematic_deltaFp end interface @@ -238,7 +250,7 @@ end subroutine plastic_init ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ph,en) + S, Fi, ph, en) integer, intent(in) :: & ph,en real(pReal), intent(in), dimension(3,3) :: & @@ -392,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) + 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)),& @@ -406,6 +418,9 @@ 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 + end select plasticType broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 7a21fea25..f74344a0b 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -36,7 +36,7 @@ type :: tParameters P_tw, & P_nS_pos, & P_nS_neg, & - CorrespondanceMatrix + CorrespondanceMatrix !< Achal integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw !< total number of active twin systems @@ -55,7 +55,7 @@ type :: tIndexDotState xi_tw, & gamma_sl, & gamma_tw, & - f_twin + f_twin !< Achal end type tIndexDotState type :: tPhenopowerlawState @@ -64,14 +64,14 @@ type :: tPhenopowerlawState xi_tw, & gamma_sl, & gamma_tw, & - f_twin + f_twin !< Achal end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- ! containers for parameters, dot state index, and state type(tParameters), allocatable, dimension(:) :: param type(tIndexDotState), allocatable, dimension(:) :: indexDotState -type(tPhenopowerlawState), allocatable, dimension(:) :: state, dot_state, deltastate +type(tPhenopowerlawState), allocatable, dimension(:) :: state, dot_state, deltastate !< Achal. added deltastate contains @@ -84,7 +84,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - ph, i, & + ph, i, o, & Nmembers, & sizeState, sizeDotState, & startIndex, endIndex @@ -114,6 +114,7 @@ phases => config_material%get('phase') allocate(param(phases%length)) allocate(indexDotState(phases%length)) allocate(state(phases%length)) +allocate(deltastate(phases%length)) !< Achal do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle @@ -181,7 +182,7 @@ do ph = 1, phases%length prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'),phase_lattice(ph)) prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) - prm%CorrespondanceMatrix = lattice_CorrespondanceMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) + prm%CorrespondanceMatrix = lattice_CorrespondanceMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) !< Achal xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw)) @@ -210,7 +211,7 @@ do ph = 1, phases%length !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then - prm%h_0_tw_sl = pl%get_asFloat('h_0_tw-sl') + prm%h_0_tw_sl = pl%get_asFloat('h_0_tw-sl') ! sanity checks prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), & phase_lattice(ph)) prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dFloat('h_tw-sl'), & @@ -234,10 +235,12 @@ 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 + + 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 + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely @@ -271,12 +274,14 @@ do ph = 1, phases%length stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_tw - idx_dot%f_twin = [startIndex,endIndex] - deltastate(ph)%f_twin => plasticState(ph)%state(startIndex:endIndex,:) - deltastate(ph)%f_twin = spread(xi_0_tw, 2, Nmembers) - + o = plasticState(ph)%offsetDeltaState + startIndex = endIndex + 1 ! Achal + endIndex = endIndex + prm%sum_N_tw ! Achal + idx_dot%f_twin = [startIndex,endIndex] ! Achal + stt%f_twin => plasticState(ph)%state(startIndex:endIndex,:) ! Achal + deltastate(ph)%f_twin => plasticState(ph)%state(startIndex-o:endIndex-o,:) ! Achal + + !write(6,*)"delta state", deltastate(ph)%f_twin ! Achal Delete end associate @@ -401,6 +406,7 @@ end function phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates instantaneous incremental change of kinematics and associated jump state +!> Satya, Achal !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinematic_deltaFp(Mp,ph,en,twinJump,deltaFp) use math, only: & @@ -433,20 +439,16 @@ fdot_twin = (0.05_pReal*(abs(tau_tw)/state(ph)%xi_tw(:,en))**param(ph)%n_tw)/par !write(6,*) 'twin_var', twin_var !delete this -if (en==1) write(6,*)'correspondanceMatrix1', param(ph)%CorrespondanceMatrix(:,:,1) !delete this -if (en==1) write(6,*)'correspondanceMatrix2', param(ph)%CorrespondanceMatrix(:,:,2) !delete this -if (en==1) write(6,*)'correspondanceMatrix3', param(ph)%CorrespondanceMatrix(:,:,3) !delete this -if (en==1) write(6,*)'correspondanceMatrix4', param(ph)%CorrespondanceMatrix(:,:,4) !delete this -if (en==1) write(6,*)'correspondanceMatrix5', param(ph)%CorrespondanceMatrix(:,:,5) !delete this -if (en==1) write(6,*)'correspondanceMatrix6', param(ph)%CorrespondanceMatrix(:,:,6) !delete this +!if (en==1) write(6,*)'correspondanceMatrix1', param(ph)%CorrespondanceMatrix(:,:,1) !delete this !delete this call RANDOM_NUMBER(random) !write(6,*)'random',random !delete this -Success_Growth: if (random*0.0000000000000000000000001_pReal <= sum(state(ph)%gamma_tw(:,en)/param(ph)%gamma_char)) then +Success_Growth: if (random*0.0000000000000000000000001_pReal <= sum(state(ph)%gamma_tw(:,en)/param(ph)%gamma_char)) then ! Instead of sum take max twinJump = .true. deltaFp = param(ph)%CorrespondanceMatrix(:,:,twin_var) + deltastate(ph)%f_twin(:,en) = 0.0_pReal - state(ph)%f_twin(:,en) end if Success_Growth @@ -458,14 +460,14 @@ end subroutine plastic_kinematic_deltaFp !> @brief calculates (instantaneous) incremental change of microstructure !> Satya, Achal !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_deltaState(ph,en) +module subroutine plastic_phenopowerlaw_deltaState(ph,en) implicit none integer, intent(in)::& ph, & en -deltastate(ph)%gamma_sl=0.0_pReal +deltastate(ph)%f_twin=0.0_pReal end subroutine plastic_phenopowerlaw_deltaState