223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

notation from DAMASK paper

This commit is contained in:
Martin Diehl 2020-03-15 21:16:28 +01:00
parent 16e23b113a
commit c7f3c2cb56
3 changed files with 91 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/constitutive_plastic_dislotwin.f90
View File

@ -71,9 +71,9 @@ submodule(constitutive) plastic_dislotwin
forestProjection, & forestProjection, &
C66 C66
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
P_tr, &
P_sl, & P_sl, &
P_tw, & P_tw, &
P_tr, &
C66_tw, & C66_tw, &
C66_tr C66_tr
integer :: & integer :: &

64
src/constitutive_plastic_kinehardening.f90
View File

@ -23,11 +23,11 @@ submodule(constitutive) plastic_kinehardening
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, & P, &
nonSchmid_pos, & nonSchmid_pos, &
nonSchmid_neg nonSchmid_neg
integer :: & integer :: &
totalNslip, & !< total number of active slip system sum_N_sl, & !< total number of active slip system
of_debug = 0 of_debug = 0
integer, allocatable, dimension(:) :: & integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family Nslip !< number of active slip systems for each family
@ -102,9 +102,9 @@ module subroutine plastic_kinehardening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%sum_N_sl = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& prm%P = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then if(trim(config%getString('lattice_structure')) == 'bcc') then
@ -113,8 +113,8 @@ module subroutine plastic_kinehardening_init
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else else
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_pos = prm%P
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%P
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
@ -157,8 +157,8 @@ module subroutine plastic_kinehardening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
@ -166,7 +166,7 @@ module subroutine plastic_kinehardening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
startIndex = 1 startIndex = 1
endIndex = prm%totalNslip endIndex = prm%sum_N_sl
stt%crss => plasticState(p)%state (startIndex:endIndex,:) stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(prm%crss0, 2, NipcMyPhase) stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
@ -175,13 +175,13 @@ module subroutine plastic_kinehardening_init
extmsg = trim(extmsg)//' atol_crss' extmsg = trim(extmsg)//' atol_crss'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%sum_N_sl
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%sum_N_sl
stt%accshear => plasticState(p)%state (startIndex:endIndex,:) stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:) dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
@ -192,17 +192,17 @@ module subroutine plastic_kinehardening_init
o = plasticState(p)%offsetDeltaState o = plasticState(p)%offsetDeltaState
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%sum_N_sl
stt%sense => plasticState(p)%state (startIndex :endIndex ,:) stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%sum_N_sl
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%sum_N_sl
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
@ -238,7 +238,7 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg dgdot_dtau_pos,dgdot_dtau_neg
@ -249,12 +249,12 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%totalNslip do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtau_pos(i) * prm%P(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) + dgdot_dtau_neg(i) * prm%P(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo enddo
end associate end associate
@ -275,7 +275,7 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of)
real(pReal) :: & real(pReal) :: &
sumGamma sumGamma
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
gdot_pos,gdot_neg gdot_pos,gdot_neg
@ -315,7 +315,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
instance, & instance, &
of of
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
sense sense
@ -366,22 +366,22 @@ module subroutine plastic_kinehardening_results(instance,group)
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case('resistance') case('resistance')
if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa') 'resistance against plastic slip','Pa')
case('backstress') ! ToDo: should be 'tau_back' case('backstress') ! ToDo: should be 'tau_back'
if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa') 'back stress against plastic slip','Pa')
case ('sense') case ('sense')
if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%sense,'sense_of_shear', &
'tbd','1') 'tbd','1')
case ('chi0') case ('chi0')
if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%chi0,'chi0', &
'tbd','Pa') 'tbd','Pa')
case ('gamma0') case ('gamma0')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma0,'gamma0', &
'tbd','1') 'tbd','1')
case ('accumulatedshear') case ('accumulatedshear')
if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1') 'plastic shear','1')
end select end select
enddo outputsLoop enddo outputsLoop
@ -406,14 +406,14 @@ pure subroutine kinetics(Mp,instance,of, &
instance, & instance, &
of of
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
gdot_pos, & gdot_pos, &
gdot_neg gdot_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
dgdot_dtau_pos, & dgdot_dtau_pos, &
dgdot_dtau_neg dgdot_dtau_neg
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
tau_pos, & tau_pos, &
tau_neg tau_neg
integer :: i integer :: i
@ -423,7 +423,7 @@ pure subroutine kinetics(Mp,instance,of, &
nonSchmidActive = size(prm%nonSchmidCoeff) > 0 nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip do i = 1, prm%sum_N_sl
tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
0.0_pReal, nonSchmidActive) 0.0_pReal, nonSchmidActive)

88
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -34,13 +34,13 @@ submodule(constitutive) plastic_phenopowerlaw
interaction_TwinSlip, & !< twin resistance from slip activity interaction_TwinSlip, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity interaction_TwinTwin !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
Schmid_slip, & P_sl, &
Schmid_twin, & P_tw, &
nonSchmid_pos, & nonSchmid_pos, &
nonSchmid_neg nonSchmid_neg
integer :: & integer :: &
totalNslip, & !< total number of active slip system sum_N_sl, & !< total number of active slip system
totalNtwin !< total number of active twin systems sum_N_tw !< total number of active twin systems
integer, allocatable, dimension(:) :: & integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family Ntwin !< number of active twin systems for each family
@ -109,9 +109,9 @@ module subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%sum_N_sl = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& prm%P_sl = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then if(trim(config%getString('lattice_structure')) == 'bcc') then
@ -121,8 +121,8 @@ module subroutine plastic_phenopowerlaw_init
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else else
allocate(prm%nonSchmidCoeff(0)) allocate(prm%nonSchmidCoeff(0))
prm%nonSchmid_pos = prm%Schmid_slip prm%nonSchmid_pos = prm%P_sl
prm%nonSchmid_neg = prm%Schmid_slip prm%nonSchmid_neg = prm%P_sl
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
@ -157,9 +157,9 @@ module subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! twin related parameters ! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin) prm%sum_N_tw = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),& prm%P_tw = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,& prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), & config%getFloats('interaction_twintwin'), &
@ -188,7 +188,7 @@ module subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip-twin related parameters ! slip-twin related parameters
slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h0_TwinSlip = config%getFloat('h0_twinslip') prm%h0_TwinSlip = config%getFloat('h0_twinslip')
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,& prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), & config%getFloats('interaction_sliptwin'), &
@ -197,8 +197,8 @@ module subroutine plastic_phenopowerlaw_init
config%getFloats('interaction_twinslip'), & config%getFloats('interaction_twinslip'), &
config%getString('lattice_structure')) config%getString('lattice_structure'))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0 allocate(prm%interaction_SlipTwin(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 allocate(prm%interaction_TwinSlip(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive endif slipAndTwinActive
@ -209,8 +209,8 @@ module subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & sizeDotState = size(['tau_slip ','gamma_slip']) * prm%sum_N_sl &
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin + size(['tau_twin ','gamma_twin']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
@ -218,7 +218,7 @@ module subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
startIndex = 1 startIndex = 1
endIndex = prm%totalNslip endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
@ -227,14 +227,14 @@ module subroutine plastic_phenopowerlaw_init
extmsg = trim(extmsg)//' atol_xi' extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%sum_N_sl
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
@ -244,7 +244,7 @@ module subroutine plastic_phenopowerlaw_init
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin endIndex = endIndex + prm%sum_N_tw
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_twinfrac',defaultVal=1.0e-6_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_twinfrac',defaultVal=1.0e-6_pReal)
@ -284,10 +284,10 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
gdot_slip_pos,gdot_slip_neg, & gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(instance)%totalNtwin) :: & real(pReal), dimension(param(instance)%sum_N_tw) :: &
gdot_twin,dgdot_dtautwin gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal Lp = 0.0_pReal
@ -296,20 +296,20 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
associate(prm => param(instance)) associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%totalNslip slipSystems: do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%totalNtwin twinSystems: do i = 1, prm%sum_N_tw
Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) + dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinSystems enddo twinSystems
end associate end associate
@ -332,7 +332,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
c_SlipSlip,c_TwinSlip,c_TwinTwin, & c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,& xi_slip_sat_offset,&
sumGamma,sumF sumGamma,sumF
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
left_SlipSlip,right_SlipSlip, & left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg gdot_slip_pos,gdot_slip_neg
@ -388,17 +388,17 @@ module subroutine plastic_phenopowerlaw_results(instance,group)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case('resistance_slip') case('resistance_slip')
if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa') 'resistance against plastic slip','Pa')
case('accumulatedshear_slip') case('accumulatedshear_slip')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1') 'plastic shear','1')
case('resistance_twin') case('resistance_twin')
if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', & if(prm%sum_N_tw>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa') 'resistance against twinning','Pa')
case('accumulatedshear_twin') case('accumulatedshear_twin')
if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & if(prm%sum_N_tw>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1') 'twinning shear','1')
end select end select
@ -424,14 +424,14 @@ pure subroutine kinetics_slip(Mp,instance,of, &
instance, & instance, &
of of
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
gdot_slip_pos, & gdot_slip_pos, &
gdot_slip_neg gdot_slip_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
dgdot_dtau_slip_pos, & dgdot_dtau_slip_pos, &
dgdot_dtau_slip_neg dgdot_dtau_slip_neg
real(pReal), dimension(param(instance)%totalNslip) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
tau_slip_pos, & tau_slip_pos, &
tau_slip_neg tau_slip_neg
integer :: i integer :: i
@ -441,7 +441,7 @@ pure subroutine kinetics_slip(Mp,instance,of, &
nonSchmidActive = size(prm%nonSchmidCoeff) > 0 nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip do i = 1, prm%sum_N_sl
tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
0.0_pReal, nonSchmidActive) 0.0_pReal, nonSchmidActive)
@ -496,19 +496,19 @@ pure subroutine kinetics_twin(Mp,instance,of,&
instance, & instance, &
of of
real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: &
gdot_twin gdot_twin
real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & real(pReal), dimension(param(instance)%sum_N_tw), intent(out), optional :: &
dgdot_dtau_twin dgdot_dtau_twin
real(pReal), dimension(param(instance)%totalNtwin) :: & real(pReal), dimension(param(instance)%sum_N_tw) :: &
tau_twin tau_twin
integer :: i integer :: i
associate(prm => param(instance), stt => state(instance)) associate(prm => param(instance), stt => state(instance))
do i = 1, prm%totalNtwin do i = 1, prm%sum_N_tw
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) tau_twin(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i))
enddo enddo
where(tau_twin > 0.0_pReal) where(tau_twin > 0.0_pReal)