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

no need for outputID 

just adds overhead, one string comparison per output and increment is
computationally not an issue

also unified to PEP recommendation of function description
This commit is contained in:
Martin Diehl 2020-02-14 07:17:30 +01:00
parent 6adb116712
commit 486385978c
3 changed files with 264 additions and 393 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

153
src/constitutive_plastic_isotropic.f90
View File

@ -8,14 +8,7 @@
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_isotropic
enum, bind(c)
enumerator :: &
undefined_ID, &
xi_ID, &
dot_gamma_ID
end enum
type :: tParameters
real(pReal) :: &
M, & !< Taylor factor
@ -34,12 +27,12 @@ submodule(constitutive) plastic_isotropic
aTol_gamma
integer :: &
of_debug = 0
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
logical :: &
dilatation
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tIsotropicState
real(pReal), pointer, dimension(:) :: &
xi, &
@ -56,50 +49,45 @@ submodule(constitutive) plastic_isotropic
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_init
integer :: &
Ninstance, &
p, i, &
p, &
NipcMyPhase, &
sizeState, sizeDotState
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
#ifdef DEBUG
if (p==material_phaseAt(debug_g,debug_e)) &
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
#endif
prm%xi_0 = config%getFloat('tau0')
prm%xi_inf = config%getFloat('tausat')
prm%dot_gamma_0 = config%getFloat('gdot0')
@ -114,9 +102,9 @@ module subroutine plastic_isotropic_init
prm%a = config%getFloat('a')
prm%aTol_xi = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
prm%aTol_gamma = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%dilatation = config%keyExists('/dilatation/')
!--------------------------------------------------------------------------------------------------
! sanity checks
extmsg = ''
@ -128,79 +116,62 @@ module subroutine plastic_isotropic_init
if (prm%M <= 0.0_pReal) extmsg = trim(extmsg)//' m'
if (prm%aTol_xi <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('flowstress')
outputID = xi_ID
case ('strainrate')
outputID = dot_gamma_ID
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['xi ','accumulated_shear'])
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
stt%xi => plasticState(p)%state (1,:)
stt%xi = prm%xi_0
dot%xi => plasticState(p)%dotState(1,:)
plasticState(p)%aTolState(1) = prm%aTol_xi
stt%gamma => plasticState(p)%state (2,:)
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%aTolState(2) = prm%aTol_gamma
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(3,3) :: &
Mp_dev !< deviatoric part of the Mandel stress
real(pReal) :: &
@ -209,16 +180,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n
Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
@ -240,42 +211,42 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
end if
end associate
end subroutine plastic_isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate inelastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
tr !< trace of spherical part of Mandel stress (= 3 x pressure)
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
tr=math_trace33(math_spherical33(Mi))
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
@ -292,38 +263,38 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
Li = 0.0_pReal
dLi_dMi = 0.0_pReal
endif
end associate
end subroutine plastic_isotropic_LiAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
if (prm%dilatation) then
norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
else
norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n
if (dot_gamma > 1e-12_pReal) then
if (dEq0(prm%c_1)) then
xi_inf_star = prm%xi_inf
@ -339,28 +310,28 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
else
dot%xi(of) = 0.0_pReal
endif
dot%gamma(of) = dot_gamma ! ToDo: not really used
end associate
end subroutine plastic_isotropic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (xi_ID)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('flowstress') ! ToDo: should be 'xi'
call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa')
end select
enddo outputsLoop

255
src/constitutive_plastic_kinehardening.f90
View File

@ -7,26 +7,13 @@
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_kinehardening
enum, bind(c)
enumerator :: &
undefined_ID, &
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, &
shearrate_ID, &
resolvedstress_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip
n, & !< stress exponent for slip
aTolResistance, &
aTolShear
real(pReal), allocatable, dimension(:) :: &
real(pReal), allocatable, dimension(:) :: &
crss0, & !< initial critical shear stress for slip
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip
@ -35,21 +22,21 @@ submodule(constitutive) plastic_kinehardening
tau1, &
tau1_b, &
nonSchmidCoeff
real(pReal), allocatable, dimension(:,:) :: &
real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: &
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
totalNslip, & !< total number of active slip system
of_debug = 0
integer, allocatable, dimension(:) :: &
integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
@ -59,7 +46,7 @@ submodule(constitutive) plastic_kinehardening
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
end type tKinehardeningState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
@ -72,37 +59,32 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_init
integer :: &
Ninstance, &
p, i, o, &
p, o, &
NipcMyPhase, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'; flush(6)
Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
do p = 1, size(phase_plasticityInstance)
if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
@ -110,22 +92,22 @@ module subroutine plastic_kinehardening_init
dlt => deltaState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)),&
config => config_phase(p))
#ifdef DEBUG
if (p==material_phaseAt(debug_g,debug_e)) then
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
endif
#endif
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
@ -133,7 +115,7 @@ module subroutine plastic_kinehardening_init
slipActive: if (prm%totalNslip > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,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)
@ -146,7 +128,7 @@ module subroutine plastic_kinehardening_init
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
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))
@ -154,10 +136,10 @@ module subroutine plastic_kinehardening_init
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%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
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)
@ -166,9 +148,9 @@ module subroutine plastic_kinehardening_init
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)
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
@ -176,58 +158,29 @@ module subroutine plastic_kinehardening_init
if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
!ToDo: Any sensible checks for theta?
endif slipActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance')
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0)
case ('accumulatedshear')
outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0)
case ('shearrate')
outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0)
case ('resolvedstress')
outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0)
case ('backstress')
outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0)
case ('sense')
outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0)
case ('chi0')
outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0)
case ('gamma0')
outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0)
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID , outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip
sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
@ -236,13 +189,13 @@ module subroutine plastic_kinehardening_init
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
@ -250,34 +203,34 @@ module subroutine plastic_kinehardening_init
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
@ -285,26 +238,26 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%totalNslip
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -312,14 +265,14 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
+ dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
end subroutine plastic_kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,instance,of)
@ -328,40 +281,40 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
sumGamma
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
* exp(-sumGamma*prm%theta0/prm%tau1) &
)
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + &
(prm%theta0_b - prm%theta1_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
)
end associate
end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates (instantaneous) incremental change of microstructure
!> @brief Calculate (instantaneous) incremental change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
@ -370,18 +323,18 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
sense
associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,of), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug &
@ -390,7 +343,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
write(6,*) sense,state(instance)%sense(:,of)
endif
#endif
!--------------------------------------------------------------------------------------------------
! switch in sense of shear?
where(dNeq(sense,stt%sense(:,of),0.1_pReal))
@ -402,45 +355,46 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
dlt%chi0 (:,of) = 0.0_pReal
dlt%gamma0(:,of) = 0.0_pReal
end where
end associate
end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('resistance')
if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
case(crss_back_ID)
call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case (sense_ID)
call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
case('backstress') ! ToDo: should be 'tau_back'
if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case ('sense')
if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
case ('chi0')
if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
case ('gamma0')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
case ('accumulatedshear')
if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
case (chi0_ID)
call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
case (gamma0_ID)
call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
case (accshear_ID)
call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
end select
enddo outputsLoop
end associate
@ -449,10 +403,11 @@ end subroutine plastic_kinehardening_results
!--------------------------------------------------------------------------------------------------
!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
!> @details: Shear rates are calculated only optionally.
!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved
! stress.
!> @details: Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,instance,of, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
@ -462,44 +417,44 @@ pure subroutine kinetics(Mp,instance,of, &
integer, intent(in) :: &
instance, &
of
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_pos, &
gdot_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_pos, &
dgdot_dtau_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_pos, &
tau_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
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), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_pos))
gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
else where
gdot_pos = 0.0_pReal
end where
where(dNeq0(tau_neg))
gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
else where
gdot_neg = 0.0_pReal
end where
if (present(dgdot_dtau_pos)) then
where(dNeq0(gdot_pos))
dgdot_dtau_pos = gdot_pos*prm%n/tau_pos

249
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -6,19 +6,6 @@
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_phenopowerlaw
enum, bind(c)
enumerator :: &
undefined_ID, &
resistance_slip_ID, &
accumulatedshear_slip_ID, &
shearrate_slip_ID, &
resolvedstress_slip_ID, &
resistance_twin_ID, &
accumulatedshear_twin_ID, &
shearrate_twin_ID, &
resolvedstress_twin_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip
@ -37,19 +24,19 @@ submodule(constitutive) plastic_phenopowerlaw
aTolResistance, & !< absolute tolerance for integration of xi
aTolShear, & !< absolute tolerance for integration of gamma
aTolTwinfrac !< absolute tolerance for integration of f
real(pReal), allocatable, dimension(:) :: &
real(pReal), allocatable, dimension(:) :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
xi_slip_sat, & !< maximum critical shear stress for slip
nonSchmidCoeff, &
H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: &
real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity
interaction_SlipTwin, & !< slip resistance from twin activity
interaction_TwinSlip, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: &
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid_slip, &
Schmid_twin, &
nonSchmid_pos, &
@ -57,13 +44,13 @@ submodule(constitutive) plastic_phenopowerlaw
integer :: &
totalNslip, & !< total number of active slip system
totalNtwin !< total number of active twin systems
integer, allocatable, dimension(:) :: &
integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
@ -71,7 +58,7 @@ submodule(constitutive) plastic_phenopowerlaw
gamma_slip, &
gamma_twin
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
@ -83,7 +70,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_init
@ -91,51 +78,46 @@ module subroutine plastic_phenopowerlaw_init
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, outputSize, &
NipcMyPhase, &
sizeState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'; flush(6)
Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal)
prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal)
prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal)
prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal)
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
@ -143,7 +125,7 @@ module subroutine plastic_phenopowerlaw_init
slipActive: if (prm%totalNslip > 0) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,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)
@ -157,22 +139,22 @@ module subroutine plastic_phenopowerlaw_init
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
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%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
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)
! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip'
@ -183,7 +165,7 @@ module subroutine plastic_phenopowerlaw_init
allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0))
endif slipActive
!--------------------------------------------------------------------------------------------------
! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
@ -196,17 +178,17 @@ module subroutine plastic_phenopowerlaw_init
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
prm%gdot0_twin = config%getFloat('gdot0_twin')
prm%n_twin = config%getFloat('n_twin')
prm%spr = config%getFloat('s_pr')
prm%h0_TwinTwin = config%getFloat('h0_twintwin')
! expand: family => system
prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin)
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin'
@ -215,7 +197,7 @@ module subroutine plastic_phenopowerlaw_init
allocate(prm%xi_twin_0(0))
allocate(prm%gamma_twin_char(0))
endif twinActive
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
@ -231,63 +213,25 @@ module subroutine plastic_phenopowerlaw_init
allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance_slip')
outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('accumulatedshear_slip')
outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('shearrate_slip')
outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resolvedstress_slip')
outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resistance_twin')
outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('accumulatedshear_twin')
outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('shearrate_twin')
outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('resolvedstress_twin')
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip &
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
@ -296,14 +240,14 @@ module subroutine plastic_phenopowerlaw_init
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
@ -317,18 +261,18 @@ module subroutine plastic_phenopowerlaw_init
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume
!--------------------------------------------------------------------------------------------------
@ -338,13 +282,13 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
@ -352,12 +296,12 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(instance)%totalNtwin) :: &
gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%totalNslip
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
@ -366,7 +310,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(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)
enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%totalNtwin
Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
@ -374,14 +318,14 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
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)
enddo twinSystems
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
@ -390,7 +334,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
@ -398,9 +342,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
real(pReal), dimension(param(instance)%totalNslip) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
@ -409,26 +353,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int
xi_slip_sat_offset = prm%spr*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset))
!--------------------------------------------------------------------------------------------------
! shear rates
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of))
!--------------------------------------------------------------------------------------------------
! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
end associate
@ -437,33 +381,33 @@ end subroutine plastic_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('resistance_slip')
if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case('accumulatedshear_slip')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case('resistance_twin')
if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case('accumulatedshear_twin')
if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
case (resistance_slip_ID)
call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case (accumulatedshear_slip_ID)
call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case (resistance_twin_ID)
call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case (accumulatedshear_twin_ID)
call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
end select
enddo outputsLoop
end associate
@ -472,10 +416,11 @@ end subroutine plastic_phenopowerlaw_results
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved.
! stress.
!> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,of, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
@ -485,44 +430,44 @@ pure subroutine kinetics_slip(Mp,instance,of, &
integer, intent(in) :: &
instance, &
of
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos, &
gdot_slip_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_slip_pos, &
dgdot_dtau_slip_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_slip_pos, &
tau_slip_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
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)), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
else where
gdot_slip_pos = 0.0_pReal
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg)
else where
gdot_slip_neg = 0.0_pReal
end where
if (present(dgdot_dtau_slip_pos)) then
where(dNeq0(gdot_slip_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
@ -543,9 +488,9 @@ end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
! twinning is assumed to take place only in untwinned volume.
!> @details Derivates are calculated only optionally.
!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved
! stress. Twinning is assumed to take place only in untwinned volume.
!> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
@ -557,29 +502,29 @@ pure subroutine kinetics_twin(Mp,instance,of,&
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
gdot_twin
real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
dgdot_dtau_twin
real(pReal), dimension(param(instance)%totalNtwin) :: &
tau_twin
integer :: i
associate(prm => param(instance), stt => state(instance))
do i = 1, prm%totalNtwin
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
else where
gdot_twin = 0.0_pReal
end where
if (present(dgdot_dtau_twin)) then
where(dNeq0(gdot_twin))
dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
@ -587,7 +532,7 @@ pure subroutine kinetics_twin(Mp,instance,of,&
dgdot_dtau_twin = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_twin