no need to store # sytems/family

This commit is contained in:
Martin Diehl 2020-03-15 23:07:23 +01:00
parent c7f3c2cb56
commit 0735d2da7c
5 changed files with 108 additions and 119 deletions

View File

@ -13,7 +13,7 @@ submodule(constitutive) plastic_disloUCLA
type :: tParameters
real(pReal) :: &
D, & !< grain size
mu, &
mu, & !< equivalent shear modulus
D_0, & !< prefactor for self-diffusion coefficient
Q_cl !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: &
@ -43,8 +43,6 @@ submodule(constitutive) plastic_disloUCLA
nonSchmid_neg
integer :: &
sum_N_sl !< total number of active slip system
integer, allocatable, dimension(:) :: &
N_sl !< number of active slip systems for each family
character(len=pStringLen), allocatable, dimension(:) :: &
output
logical :: &
@ -87,7 +85,8 @@ module subroutine plastic_disloUCLA_init
NipcMyPhase, &
sizeState, sizeDotState, &
startIndex, endIndex
integer, dimension(:), allocatable :: &
N_sl
character(len=pStringLen) :: &
extmsg = ''
@ -121,46 +120,46 @@ module subroutine plastic_disloUCLA_init
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(prm%N_sl)
N_sl = config%getInts('nslip',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(N_sl)
slipActive: if (prm%sum_N_sl > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),&
prm%Schmid = lattice_SchmidMatrix_slip(N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%N_sl,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%N_sl,prm%nonSchmidCoeff,-1)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1)
else
allocate(prm%nonSchmidCoeff(0))
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, &
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%forestProjectionEdge = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),&
prm%forestProjectionEdge = lattice_forestProjection_edge(N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%forestProjectionEdge = transpose(prm%forestProjectionEdge)
prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl))
prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl))
prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl))
prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl))
prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl))
prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(N_sl))
prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(N_sl))
prm%v0 = config%getFloats('v0', requiredSize=size(N_sl))
prm%b_sl = config%getFloats('slipburgers', requiredSize=size(N_sl))
prm%delta_F = config%getFloats('qedge', requiredSize=size(N_sl))
prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl))
prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl))
prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), &
defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))])
prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl), &
defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))])
prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%N_sl))
prm%w = config%getFloats('kink_width', requiredSize=size(prm%N_sl))
prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl))
prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl))
prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(N_sl))
prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(N_sl))
prm%p = config%getFloats('p_slip', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%q = config%getFloats('q_slip', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%kink_height = config%getFloats('kink_height', requiredSize=size(N_sl))
prm%w = config%getFloats('kink_width', requiredSize=size(N_sl))
prm%omega = config%getFloats('omega', requiredSize=size(N_sl))
prm%B = config%getFloats('friction_coeff', requiredSize=size(N_sl))
prm%D = config%getFloat('grainsize')
prm%D_0 = config%getFloat('d0')
@ -170,21 +169,21 @@ module subroutine plastic_disloUCLA_init
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
! expand: family => system
prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl)
prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl)
prm%q = math_expand(prm%q, prm%N_sl)
prm%p = math_expand(prm%p, prm%N_sl)
prm%delta_F = math_expand(prm%delta_F, prm%N_sl)
prm%b_sl = math_expand(prm%b_sl, prm%N_sl)
prm%kink_height = math_expand(prm%kink_height, prm%N_sl)
prm%w = math_expand(prm%w, prm%N_sl)
prm%omega = math_expand(prm%omega, prm%N_sl)
prm%tau_0 = math_expand(prm%tau_0, prm%N_sl)
prm%v0 = math_expand(prm%v0, prm%N_sl)
prm%B = math_expand(prm%B, prm%N_sl)
prm%i_sl = math_expand(prm%i_sl, prm%N_sl)
prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl)
prm%D_a = math_expand(prm%D_a, prm%N_sl)
prm%rho_mob_0 = math_expand(prm%rho_mob_0, N_sl)
prm%rho_dip_0 = math_expand(prm%rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%delta_F = math_expand(prm%delta_F, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%kink_height = math_expand(prm%kink_height, N_sl)
prm%w = math_expand(prm%w, N_sl)
prm%omega = math_expand(prm%omega, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl)
prm%v0 = math_expand(prm%v0, N_sl)
prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%atomicVolume = math_expand(prm%atomicVolume, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl)
! sanity checks
@ -213,15 +212,14 @@ module subroutine plastic_disloUCLA_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho')
if (any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_rho'
if (any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
@ -234,7 +232,7 @@ module subroutine plastic_disloUCLA_init
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = 1.0e6_pReal ! Don't use for convergence check
plasticState(p)%atol(startIndex:endIndex) = 1.0e6_pReal ! ARRG
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)

View File

@ -120,7 +120,7 @@ module subroutine plastic_isotropic_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol
! state aliases and initialization
stt%xi => plasticState(p)%state (1,:)
stt%xi = prm%xi_0
dot%xi => plasticState(p)%dotState(1,:)
@ -131,7 +131,6 @@ module subroutine plastic_isotropic_init
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%atol(2) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
if (plasticState(p)%atol(2) <= 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)

View File

@ -29,8 +29,6 @@ submodule(constitutive) plastic_kinehardening
integer :: &
sum_N_sl, & !< total number of active slip system
of_debug = 0
integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
@ -68,7 +66,8 @@ module subroutine plastic_kinehardening_init
NipcMyPhase, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer, dimension(:), allocatable :: &
N_sl
character(len=pStringLen) :: &
extmsg = ''
@ -101,44 +100,44 @@ module subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(prm%Nslip)
N_sl = config%getInts('nslip',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(N_sl)
slipActive: if (prm%sum_N_sl > 0) then
prm%P = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
prm%P = lattice_SchmidMatrix_slip(N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1)
else
prm%nonSchmid_pos = prm%P
prm%nonSchmid_neg = prm%P
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%crss0 = config%getFloats('crss0', requiredSize=size(N_sl))
prm%tau1 = config%getFloats('tau1', requiredSize=size(N_sl))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(N_sl))
prm%theta0 = config%getFloats('theta0', requiredSize=size(N_sl))
prm%theta1 = config%getFloats('theta1', requiredSize=size(N_sl))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(N_sl))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(N_sl))
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
! expand: family => system
prm%crss0 = math_expand(prm%crss0, prm%Nslip)
prm%tau1 = math_expand(prm%tau1, prm%Nslip)
prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip)
prm%theta0 = math_expand(prm%theta0, prm%Nslip)
prm%theta1 = math_expand(prm%theta1, prm%Nslip)
prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip)
prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip)
prm%crss0 = math_expand(prm%crss0, N_sl)
prm%tau1 = math_expand(prm%tau1, N_sl)
prm%tau1_b = math_expand(prm%tau1_b, N_sl)
prm%theta0 = math_expand(prm%theta0, N_sl)
prm%theta1 = math_expand(prm%theta1, N_sl)
prm%theta0_b = math_expand(prm%theta0_b,N_sl)
prm%theta1_b = math_expand(prm%theta1_b,N_sl)
@ -164,29 +163,27 @@ module subroutine plastic_kinehardening_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_crss'
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%crss_back => plasticState(p)%state (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_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)

View File

@ -13,7 +13,7 @@ submodule(constitutive) plastic_nonlocal
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0
real(pReal), parameter :: &
KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
! storage order of dislocation types
integer, dimension(8), parameter :: &
@ -215,7 +215,7 @@ module subroutine plastic_nonlocal_init
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
prm%atol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%atol_rho = config%getFloat('atol_rho',defaultVal=1.0e4_pReal)
structure = config%getString('lattice_structure')
@ -350,7 +350,7 @@ module subroutine plastic_nonlocal_init
if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' significantN'
if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' significantrho'
if (prm%atol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' CFLfactor'
if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p'
@ -462,8 +462,8 @@ module subroutine plastic_nonlocal_init
stt%gamma => plasticState(p)%state (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
dot%gamma => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
del%gamma => plasticState(p)%deltaState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
plasticState(p)%atol(10*prm%totalNslip+1:11*prm%totalNslip ) = config%getFloat('atol_shear', defaultVal=0.0_pReal)
if(any(plasticState(p)%atol(10*prm%totalNslip+1:11*prm%totalNslip)<=0.0_pReal)) &
plasticState(p)%atol(10*prm%totalNslip+1:11*prm%totalNslip ) = config%getFloat('atol_gamma', defaultVal = 1.0e-20_pReal)
if(any(plasticState(p)%atol(10*prm%totalNslip+1:11*prm%totalNslip) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)

View File

@ -41,9 +41,6 @@ submodule(constitutive) plastic_phenopowerlaw
integer :: &
sum_N_sl, & !< total number of active slip system
sum_N_tw !< total number of active twin systems
integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
@ -78,7 +75,8 @@ module subroutine plastic_phenopowerlaw_init
NipcMyPhase, &
sizeState, sizeDotState, &
startIndex, endIndex
integer, dimension(:), allocatable :: &
N_sl, N_tw
character(len=pStringLen) :: &
extmsg = ''
@ -108,30 +106,30 @@ module subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(prm%Nslip)
N_sl = config%getInts('nslip',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(N_sl)
slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1)
else
allocate(prm%nonSchmidCoeff(0))
prm%nonSchmid_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))])
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(N_sl))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(N_sl))
prm%H_int = config%getFloats('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal,i=1,size(N_sl))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
@ -139,9 +137,9 @@ module subroutine plastic_phenopowerlaw_init
prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip)
prm%H_int = math_expand(prm%H_int, prm%Nslip)
prm%xi_slip_0 = math_expand(prm%xi_slip_0, N_sl)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl)
prm%H_int = math_expand(prm%H_int, N_sl)
! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip'
@ -156,18 +154,18 @@ module subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(prm%Ntwin)
N_tw = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(N_tw)
twinActive: if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,&
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(N_tw,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(N_tw))
prm%gdot0_twin = config%getFloat('gdot0_twin')
prm%n_twin = config%getFloat('n_twin')
@ -175,7 +173,7 @@ module subroutine plastic_phenopowerlaw_init
prm%h0_TwinTwin = config%getFloat('h0_twintwin')
! expand: family => system
prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin)
prm%xi_twin_0 = math_expand(prm%xi_twin_0, N_tw)
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin'
@ -190,10 +188,10 @@ module subroutine plastic_phenopowerlaw_init
! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h0_TwinSlip = config%getFloat('h0_twinslip')
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(N_sl,N_tw,&
config%getFloats('interaction_sliptwin'), &
config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,&
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(N_tw,N_sl,&
config%getFloats('interaction_twinslip'), &
config%getString('lattice_structure'))
else slipAndTwinActive
@ -216,15 +214,14 @@ module subroutine plastic_phenopowerlaw_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) &
extmsg = trim(extmsg)//' atol_xi'
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
@ -238,8 +235,7 @@ module subroutine plastic_phenopowerlaw_init
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma_slip'
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma_slip'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
@ -247,9 +243,8 @@ module subroutine plastic_phenopowerlaw_init
endIndex = endIndex + prm%sum_N_tw
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_twinfrac',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma_twin'
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally