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

Merge branch 'less-instance' into 'development' 

don't use instances in plastic models

See merge request damask/DAMASK!339
This commit is contained in:
Sharan Roongta 2021-02-19 13:32:12 +00:00
parent b6146a8cc6 53bab41b47
commit c9e4dc21f8
11 changed files with 613 additions and 677 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
src/phase.f90
View File

@ -59,8 +59,7 @@ module phase
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
phase_elasticityInstance, &
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_plasticInstance
phase_NstiffnessDegradations
logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler)
phase_localPlasticity !< flags phases with local constitutive law
@ -247,9 +246,9 @@ module phase
TDot
end subroutine phase_thermal_getRate
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
integer, intent(in) :: &
instance, &
ph, &
i, &
e
type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
@ -616,7 +615,7 @@ subroutine crystallite_orientations(co,ip,el)
if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticInstance(material_phaseAt(1,el)),ip,el)
material_phaseAt(1,el),ip,el)
end subroutine crystallite_orientations

65
src/phase_mechanical.f90
View File

@ -60,7 +60,7 @@ submodule(phase) mechanical
module subroutine plastic_init
end subroutine plastic_init
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
@ -68,7 +68,7 @@ submodule(phase) mechanical
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
end subroutine plastic_isotropic_LiAndItsTangent
@ -85,11 +85,8 @@ submodule(phase) mechanical
logical :: broken
end function plastic_dotState
module function plastic_deltaState(co, ip, el, ph, me) result(broken)
module function plastic_deltaState(ph, me) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
logical :: &
@ -114,11 +111,9 @@ submodule(phase) mechanical
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, co, ip, el)
S, Fi, ph,me)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
ph,me
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
@ -130,33 +125,33 @@ submodule(phase) mechanical
end subroutine plastic_LpAndItsTangents
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
module subroutine plastic_isotropic_results(ph,group)
integer, intent(in) :: ph
character(len=*), intent(in) :: group
end subroutine plastic_isotropic_results
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
module subroutine plastic_phenopowerlaw_results(ph,group)
integer, intent(in) :: ph
character(len=*), intent(in) :: group
end subroutine plastic_phenopowerlaw_results
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
module subroutine plastic_kinehardening_results(ph,group)
integer, intent(in) :: ph
character(len=*), intent(in) :: group
end subroutine plastic_kinehardening_results
module subroutine plastic_dislotwin_results(instance,group)
integer, intent(in) :: instance
module subroutine plastic_dislotwin_results(ph,group)
integer, intent(in) :: ph
character(len=*), intent(in) :: group
end subroutine plastic_dislotwin_results
module subroutine plastic_dislotungsten_results(instance,group)
integer, intent(in) :: instance
module subroutine plastic_dislotungsten_results(ph,group)
integer, intent(in) :: ph
character(len=*), intent(in) :: group
end subroutine plastic_dislotungsten_results
module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance
module subroutine plastic_nonlocal_results(ph,group)
integer, intent(in) :: ph
character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results
@ -303,14 +298,12 @@ module subroutine mechanical_init(phases)
! initialize plasticity
allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
allocate(phase_plasticInstance(phases%length),source = 0)
allocate(phase_localPlasticity(phases%length), source=.true.)
call plastic_init()
do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
phase_plasticInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph))
enddo
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -403,22 +396,22 @@ module subroutine mechanical_results(group,ph)
select case(phase_plasticity(ph))
case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticInstance(ph),group//'plastic/')
call plastic_isotropic_results(ph,group//'plastic/')
case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticInstance(ph),group//'plastic/')
call plastic_phenopowerlaw_results(ph,group//'plastic/')
case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticInstance(ph),group//'plastic/')
call plastic_kinehardening_results(ph,group//'plastic/')
case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticInstance(ph),group//'plastic/')
call plastic_dislotwin_results(ph,group//'plastic/')
case(PLASTICITY_DISLOTUNGSTEN_ID)
call plastic_dislotungsten_results(phase_plasticInstance(ph),group//'plastic/')
call plastic_dislotungsten_results(ph,group//'plastic/')
case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/')
call plastic_nonlocal_results(ph,group//'plastic/')
end select
@ -538,7 +531,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
Fe, Fi_new, co, ip, el)
call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
S, Fi_new, co, ip, el)
S, Fi_new, ph,me)
!* update current residuum and check for convergence of loop
atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -702,7 +695,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) &
- r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then
broken = plastic_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(ph,me)
exit iteration
endif
@ -765,7 +758,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
plasticState(ph)%state(1:sizeDotState,me) = subState0 &
+ plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t
broken = plastic_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(ph,me)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -807,7 +800,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
plasticState(ph)%state(1:sizeDotState,me) = subState0 &
+ plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t
broken = plastic_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(ph,me)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -956,7 +949,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
if(broken) return
broken = plastic_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(ph,me)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -1329,7 +1322,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
phase_mechanical_S(ph)%data(1:3,1:3,me), &
phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el)
phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!--------------------------------------------------------------------------------------------------

3
src/phase_mechanical_eigen.f90
View File

@ -180,13 +180,12 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_isotropic_ID) plasticType
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me)
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,me)
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
active = .true.
end select plasticType
KinematicsLoop: do k = 1, Nmodels(ph)
kinematicsType: select case (model(k,ph))
case (KINEMATICS_thermal_expansion_ID) kinematicsType

65
src/phase_mechanical_plastic.f90
View File

@ -104,7 +104,7 @@ submodule(phase:mechanical) plastic
end subroutine dislotungsten_LpAndItsTangent
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,ph,me,ip,el)
Mp,Temperature,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
@ -116,9 +116,7 @@ submodule(phase:mechanical) plastic
Temperature
integer, intent(in) :: &
ph, &
me, &
ip, & !< current integration point
el !< current element number
me
end subroutine nonlocal_LpAndItsTangent
@ -179,44 +177,42 @@ submodule(phase:mechanical) plastic
el !< current element number
end subroutine nonlocal_dotState
module subroutine dislotwin_dependentState(T,instance,me)
module subroutine dislotwin_dependentState(T,ph,me)
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), intent(in) :: &
T
end subroutine dislotwin_dependentState
module subroutine dislotungsten_dependentState(instance,me)
module subroutine dislotungsten_dependentState(ph,me)
integer, intent(in) :: &
instance, &
ph, &
me
end subroutine dislotungsten_dependentState
module subroutine nonlocal_dependentState(instance, me, ip, el)
module subroutine nonlocal_dependentState(ph, me, ip, el)
integer, intent(in) :: &
instance, &
ph, &
me, &
ip, & !< current integration point
el !< current element number
end subroutine nonlocal_dependentState
module subroutine plastic_kinehardening_deltaState(Mp,instance,me)
module subroutine plastic_kinehardening_deltaState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
end subroutine plastic_kinehardening_deltaState
module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
instance, &
me, &
ip, &
el
ph, &
me
end subroutine plastic_nonlocal_deltaState
end interface
@ -244,11 +240,9 @@ end subroutine plastic_init
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, co, ip, el)
S, Fi, ph,me)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
ph,me
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
@ -263,14 +257,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp
integer :: &
i, j, me, ph
i, j
Mp = matmul(matmul(transpose(Fi),Fi),S)
me = material_phasememberAt(co,ip,el)
ph = material_phaseAt(co,el)
plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_NONE_ID) plasticType
Lp = 0.0_pReal
@ -286,7 +279,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el)
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
@ -364,23 +357,22 @@ module subroutine plastic_dependentState(co, ip, el)
integer :: &
ph, &
instance, me
me
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
instance = phase_plasticInstance(ph)
plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_dependentState(thermal_T(ph,me),instance,me)
call dislotwin_dependentState(thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dependentState(instance,me)
call dislotungsten_dependentState(ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dependentState(instance,me,ip,el)
call nonlocal_dependentState(ph,me,ip,el)
end select plasticType
@ -391,12 +383,9 @@ end subroutine plastic_dependentState
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
module function plastic_deltaState(co, ip, el, ph, me) result(broken)
module function plastic_deltaState(ph, me) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
logical :: &
@ -405,23 +394,21 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken)
real(pReal), dimension(3,3) :: &
Mp
integer :: &
instance, &
myOffset, &
mySize
Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),&
phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me))
instance = phase_plasticInstance(ph)
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_deltaState(Mp,instance,me)
call plastic_kinehardening_deltaState(Mp,ph,me)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case (PLASTICITY_NONLOCAL_ID) plasticType
call plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
call plastic_nonlocal_deltaState(Mp,ph,me)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case default

105
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -78,8 +78,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, i, &
ph, i, &
Nconstituents, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -97,33 +96,31 @@ module function plastic_dislotungsten_init() result(myPlasticity)
mech, &
pl
myPlasticity = plastic_active('dislotungsten')
Ninstances = count(myPlasticity)
if(Ninstances == 0) return
if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotungsten init -+>>>'
print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
allocate(param(Ninstances))
allocate(state(Ninstances))
allocate(dotState(Ninstances))
allocate(dependentState(Ninstances))
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(dependentState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
if(.not. myPlasticity(p)) cycle
i = i + 1
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i), &
dst => dependentState(i))
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
@ -133,7 +130,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
#endif
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%mu = lattice_mu(ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -223,41 +220,41 @@ module function plastic_dislotungsten_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,Nconstituents)
dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,Nconstituents)
dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
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.0e-2_pReal
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
@ -290,16 +287,16 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
integer :: &
i,k,l,m,n
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos,dot_gamma_neg, &
ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticInstance(ph)))
associate(prm => param(ph))
call kinetics(Mp,T,phase_plasticInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
call kinetics(Mp,T,ph,me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl
Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -328,7 +325,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
real(pReal) :: &
VacancyDiffusion
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, gdot_neg,&
tau_pos,&
tau_neg, &
@ -337,10 +334,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
dot_rho_dip_climb, &
dip_distance
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),&
dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph),&
dot => dotState(ph), dst => dependentState(ph))
call kinetics(Mp,T,phase_plasticInstance(ph),me,&
call kinetics(Mp,T,ph,me,&
gdot_pos,gdot_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
@ -377,16 +374,16 @@ end subroutine dislotungsten_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Calculate derived quantities from state.
!--------------------------------------------------------------------------------------------------
module subroutine dislotungsten_dependentState(instance,me)
module subroutine dislotungsten_dependentState(ph,me)
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dislocationSpacing
associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph),dst => dependentState(ph))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
dst%threshold_stress(:,me) = prm%mu*prm%b_sl &
@ -402,14 +399,14 @@ end subroutine dislotungsten_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotungsten_results(instance,group)
module subroutine plastic_dislotungsten_results(ph,group)
integer, intent(in) :: instance
integer, intent(in) :: ph
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('rho_mob')
@ -441,7 +438,7 @@ end subroutine plastic_dislotungsten_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,T,instance,me, &
pure subroutine kinetics(Mp,T,ph,me, &
dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out)
real(pReal), dimension(3,3), intent(in) :: &
@ -449,18 +446,18 @@ pure subroutine kinetics(Mp,T,instance,me, &
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos, &
dot_gamma_neg
real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
ddot_gamma_dtau_pos, &
ddot_gamma_dtau_neg, &
tau_pos_out, &
tau_neg_out
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
StressRatio, &
StressRatio_p,StressRatio_pminus1, &
dvel, vel, &
@ -469,7 +466,7 @@ pure subroutine kinetics(Mp,T,instance,me, &
needsGoodName ! ToDo: @Karo: any idea?
integer :: j
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do j = 1, prm%sum_N_sl
tau_pos(j) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,j))

193
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -48,7 +48,7 @@ submodule(phase:plastic) dislotwin
dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system
t_tw, & !< twin thickness [m] for each twin system
i_sl, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
t_tr, & !< martensite lamellar thickness [m] for each trans system and instance
t_tr, & !< martensite lamellar thickness [m] for each trans system
p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity
r, & !< r-exponent in twin nucleation rate
@ -126,8 +126,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, i, &
ph, i, &
Nconstituents, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -144,12 +143,12 @@ module function plastic_dislotwin_init() result(myPlasticity)
mech, &
pl
myPlasticity = plastic_active('dislotwin')
Ninstances = count(myPlasticity)
if(Ninstances == 0) return
if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotwin init -+>>>'
print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
@ -160,22 +159,21 @@ module function plastic_dislotwin_init() result(myPlasticity)
print*, 'Wong et al., Acta Materialia 118:140–151, 2016'
print*, 'https://doi.org/10.1016/j.actamat.2016.07.032'
allocate(param(Ninstances))
allocate(state(Ninstances))
allocate(dotState(Ninstances))
allocate(dependentState(Ninstances))
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(dependentState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
if(.not. myPlasticity(p)) cycle
i = i + 1
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i), &
dst => dependentState(i))
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
@ -185,9 +183,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
#endif
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%nu = lattice_nu(p)
prm%C66 = lattice_C66(1:6,1:6,p)
prm%mu = lattice_mu(ph)
prm%nu = lattice_nu(ph)
prm%C66 = lattice_C66(1:6,1:6,ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -204,8 +202,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == lattice_FCC_ID) &
.and. (N_sl(1) == 12)
prm%fccTwinTransNucleation = lattice_structure(ph) == lattice_FCC_ID .and. (N_sl(1) == 12)
if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl))
@ -234,7 +231,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
prm%omega = pl%get_asFloat('omega', defaultVal = 1000.0_pReal) &
* merge(12.0_pReal,8.0_pReal,any(lattice_structure(p) == [lattice_FCC_ID,lattice_HEX_ID]))
* merge(12.0_pReal,8.0_pReal,any(lattice_structure(ph) == [lattice_FCC_ID,lattice_HEX_ID]))
! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl)
@ -342,7 +339,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
pl%get_asFloat('a_cI', defaultVal=0.0_pReal), &
pl%get_asFloat('a_cF', defaultVal=0.0_pReal))
if (lattice_structure(p) /= lattice_FCC_ID) then
if (lattice_structure(ph) /= lattice_FCC_ID) then
prm%dot_N_0_tr = pl%get_asFloats('dot_N_0_tr')
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr)
endif
@ -357,7 +354,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr'
if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr'
if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr'
if (lattice_structure(p) /= lattice_FCC_ID) then
if (lattice_structure(ph) /= lattice_FCC_ID) then
if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr'
endif
else transActive
@ -408,53 +405,53 @@ module function plastic_dislotwin_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol
startIndex = 1
endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob= spread(rho_mob_0,2,Nconstituents)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip= spread(rho_dip_0,2,Nconstituents)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
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.0e-2_pReal
stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin'
stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tr
stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
@ -469,7 +466,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
@ -496,8 +493,8 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
real(pReal) :: f_unrotated
associate(prm => param(phase_plasticInstance(ph)),&
stt => state(phase_plasticInstance(ph)))
associate(prm => param(ph),&
stt => state(ph))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) &
@ -535,11 +532,11 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
BoltzmannRatio, &
ddot_gamma_dtau, &
tau
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_twin,ddot_gamma_dtau_twin
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: &
real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr,ddot_gamma_dtau_trans
real(pReal):: dot_gamma_sb
real(pReal), dimension(3,3) :: eigVectors, P_sb
@ -564,7 +561,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
0, 1, 1 &
],pReal),[ 3,6])
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) &
@ -573,7 +570,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip)
call kinetics_slip(Mp,T,ph,me,dot_gamma_sl,ddot_gamma_dtau_slip)
slipContribution: do i = 1, prm%sum_N_sl
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -581,7 +578,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin)
call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin,ddot_gamma_dtau_twin)
twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -589,7 +586,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans)
call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_trans)
transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -653,24 +650,24 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
tau, &
sigma_cl, & !< climb stress
b_d !< ratio of Burgers vector to stacking fault width
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_rho_dip_formation, &
dot_rho_dip_climb, &
rho_dip_distance_min, &
dot_gamma_sl
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_twin
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: &
real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph), dst => dependentState(ph))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me))
call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl)
call kinetics_slip(Mp,T,ph,me,dot_gamma_sl)
dot%gamma_sl(:,me) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%D_a*prm%b_sl
@ -721,10 +718,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
- dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin)
call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin)
dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr)
call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr)
dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
end associate
@ -735,33 +732,33 @@ end subroutine dislotwin_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Calculate derived quantities from state.
!--------------------------------------------------------------------------------------------------
module subroutine dislotwin_dependentState(T,instance,me)
module subroutine dislotwin_dependentState(T,ph,me)
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), intent(in) :: &
T
real(pReal) :: &
sumf_twin,Gamma,sumf_trans
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation
inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
real(pReal), dimension(param(instance)%sum_N_tw) :: &
real(pReal), dimension(param(ph)%sum_N_tw) :: &
inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
f_over_t_tw
real(pReal), dimension(param(instance)%sum_N_tr) :: &
real(pReal), dimension(param(ph)%sum_N_tr) :: &
inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite
f_over_t_tr
real(pReal), dimension(:), allocatable :: &
x0
associate(prm => param(instance),&
stt => state(instance),&
dst => dependentState(instance))
associate(prm => param(ph),&
stt => state(ph),&
dst => dependentState(ph))
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me))
@ -827,14 +824,14 @@ end subroutine dislotwin_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotwin_results(instance,group)
module subroutine plastic_dislotwin_results(ph,group)
integer, intent(in) :: instance
integer, intent(in) :: ph
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
@ -882,7 +879,7 @@ end subroutine plastic_dislotwin_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,T,instance,me, &
pure subroutine kinetics_slip(Mp,T,ph,me, &
dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip)
real(pReal), dimension(3,3), intent(in) :: &
@ -890,18 +887,18 @@ pure subroutine kinetics_slip(Mp,T,instance,me, &
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
dot_gamma_sl
real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: &
ddot_gamma_dtau_slip, &
tau_slip
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
ddot_gamma_dtau
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, &
stressRatio, &
StressRatio_p, &
@ -914,7 +911,7 @@ pure subroutine kinetics_slip(Mp,T,instance,me, &
tau_eff !< effective resolved stress
integer :: i
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_sl
tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
@ -959,7 +956,7 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,&
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_twin,ddot_gamma_dtau_twin)
real(pReal), dimension(3,3), intent(in) :: &
@ -967,17 +964,17 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,&
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: &
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
dot_gamma_sl
real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
dot_gamma_twin
real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
ddot_gamma_dtau_twin
real, dimension(param(instance)%sum_N_tw) :: &
real, dimension(param(ph)%sum_N_tw) :: &
tau, &
Ndot0, &
stressRatio_r, &
@ -985,7 +982,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,&
integer :: i,s1,s2
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_tw
tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
@ -1028,7 +1025,7 @@ end subroutine kinetics_twin
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,me,&
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_tr,ddot_gamma_dtau_trans)
real(pReal), dimension(3,3), intent(in) :: &
@ -1036,24 +1033,24 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,me,&
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: &
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
dot_gamma_sl
real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: &
dot_gamma_tr
real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
ddot_gamma_dtau_trans
real, dimension(param(instance)%sum_N_tr) :: &
real, dimension(param(ph)%sum_N_tr) :: &
tau, &
Ndot0, &
stressRatio_s, &
ddot_gamma_dtau
integer :: i,s1,s2
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_tr
tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))

87
src/phase_mechanical_plastic_isotropic.f90
View File

@ -22,8 +22,6 @@ submodule(phase:plastic) isotropic
c_4, &
c_3, &
c_2
integer :: &
of_debug = 0
logical :: &
dilatation
character(len=pStringLen), allocatable, dimension(:) :: &
@ -53,9 +51,7 @@ module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, &
i, &
ph, &
Nconstituents, &
sizeState, sizeDotState
real(pReal) :: &
@ -68,33 +64,29 @@ module function plastic_isotropic_init() result(myPlasticity)
mech, &
pl
myPlasticity = plastic_active('isotropic')
Ninstances = count(myPlasticity)
if(Ninstances == 0) return
if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:isotropic init -+>>>'
print'(a,i0)', ' # phses: ',Ninstances; flush(IO_STDOUT)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
allocate(param(Ninstances))
allocate(state(Ninstances))
allocate(dotState(Ninstances))
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
mech => phase%get('mechanics')
if(.not. myPlasticity(p)) cycle
i = i + 1
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i))
pl => mech%get('plasticity')
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(pl)
@ -102,11 +94,6 @@ module function plastic_isotropic_init() result(myPlasticity)
prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray)
#endif
#ifdef DEBUG
if (p==material_phaseAt(debugConstitutive%grain,debugConstitutive%element)) &
prm%of_debug = material_phasememberAt(debugConstitutive%grain,debugConstitutive%ip,debugConstitutive%element)
#endif
xi_0 = pl%get_asFloat('xi_0')
prm%xi_inf = pl%get_asFloat('xi_inf')
prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0')
@ -132,28 +119,28 @@ module function plastic_isotropic_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
stt%xi => plasticState(p)%state (1,:)
stt%xi => plasticState(ph)%state (1,:)
stt%xi = xi_0
dot%xi => plasticState(p)%dotState(1,:)
plasticState(p)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(p)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
dot%xi => plasticState(ph)%dotState(1,:)
plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
stt%gamma => plasticState(p)%state (2,:)
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (plasticState(p)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
stt%gamma => plasticState(ph)%state (2,:)
dot%gamma => plasticState(ph)%dotState(2,:)
plasticState(ph)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (plasticState(ph)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
@ -191,7 +178,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: &
k, l, m, n
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph))
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
@ -221,7 +208,7 @@ end subroutine isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate inelastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
@ -231,7 +218,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal) :: &
@ -239,7 +226,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
tr=math_trace33(math_spherical33(Mi))
@ -247,9 +234,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
Li = math_I3 &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
else
Li = 0.0_pReal
dLi_dMi = 0.0_pReal
@ -276,8 +261,8 @@ module subroutine isotropic_dotState(Mp,ph,me)
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph))
if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp))
@ -313,14 +298,14 @@ end subroutine isotropic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_results(instance,group)
module subroutine plastic_isotropic_results(ph,group)
integer, intent(in) :: instance
integer, intent(in) :: ph
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('xi')

142
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -25,8 +25,7 @@ submodule(phase:plastic) kinehardening
nonSchmid_pos, &
nonSchmid_neg
integer :: &
sum_N_sl, & !< total number of active slip system
of_debug = 0
sum_N_sl
logical :: &
nonSchmidActive = .false.
character(len=pStringLen), allocatable, dimension(:) :: &
@ -62,8 +61,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, i, o, &
ph, o, &
Nconstituents, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
@ -81,29 +79,26 @@ module function plastic_kinehardening_init() result(myPlasticity)
pl
myPlasticity = plastic_active('kinehardening')
Ninstances = count(myPlasticity)
if(Ninstances == 0) return
if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:kinehardening init -+>>>'
print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
allocate(param(Ninstances))
allocate(state(Ninstances))
allocate(dotState(Ninstances))
allocate(deltaState(Ninstances))
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(deltaState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
if(.not. myPlasticity(p)) cycle
i = i + 1
associate(prm => param(i), &
dot => dotState(i), &
dlt => deltaState(i), &
stt => state(i))
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
@ -112,12 +107,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray)
#endif
#ifdef DEBUG
if (p==material_phaseAt(debugConstitutive%grain,debugConstitutive%element)) then
prm%of_debug = material_phasememberAt(debugConstitutive%grain,debugConstitutive%ip,debugConstitutive%element)
endif
#endif
!--------------------------------------------------------------------------------------------------
! slip related parameters
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
@ -176,55 +165,55 @@ module function plastic_kinehardening_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss => plasticState(ph)%state (startIndex:endIndex,:)
stt%crss = spread(xi_0, 2, Nconstituents)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%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) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
stt%crss_back => plasticState(ph)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('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) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
stt%accshear => plasticState(ph)%state (startIndex:endIndex,:)
dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
o = plasticState(ph)%offsetDeltaState
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
stt%sense => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
stt%chi0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
@ -256,16 +245,16 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: &
i,k,l,m,n
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticInstance(ph)))
associate(prm => param(ph))
call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
call kinetics(Mp,ph,me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
@ -293,14 +282,14 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,me)
real(pReal) :: &
sumGamma
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),&
dot => dotState(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph),&
dot => dotState(ph))
call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg)
call kinetics(Mp,ph,me,gdot_pos,gdot_neg)
dot%accshear(:,me) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,me))
@ -326,32 +315,25 @@ end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Calculate (instantaneous) incremental change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_deltaState(Mp,instance,me)
module subroutine plastic_kinehardening_deltaState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg, &
sense
associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
call kinetics(Mp,instance,me,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,me), & ! keep existing...
call kinetics(Mp,ph,me,gdot_pos,gdot_neg)
sense = merge(state(ph)%sense(:,me), & ! 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 (debugConstitutive%extensive &
.and. (me == prm%of_debug .or. .not. debugConstitutive%selective)) then
print*, '======= kinehardening delta state ======='
print*, sense,state(instance)%sense(:,me)
endif
#endif
!--------------------------------------------------------------------------------------------------
! switch in sense me shear?
@ -373,14 +355,14 @@ end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_results(instance,group)
module subroutine plastic_kinehardening_results(ph,group)
integer, intent(in) :: instance
integer, intent(in) :: ph
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('xi')
@ -415,28 +397,28 @@ end subroutine plastic_kinehardening_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,instance,me, &
pure subroutine kinetics(Mp,ph,me, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, &
gdot_neg
real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
dgdot_dtau_pos, &
dgdot_dtau_neg
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau_pos, &
tau_neg
integer :: i
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
do i = 1, prm%sum_N_sl
tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,me)

10
src/phase_mechanical_plastic_none.f90
View File

@ -16,8 +16,7 @@ module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, &
Nconstituents
ph
class(tNode), pointer :: &
phases
@ -29,10 +28,9 @@ module function plastic_none_init() result(myPlasticity)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get('phase')
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
Nconstituents = count(material_phaseAt2 == p)
call phase_allocateState(plasticState(p),Nconstituents,0,0,0)
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0)
enddo
end function plastic_none_init

465
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -13,6 +13,12 @@ submodule(phase:plastic) nonlocal
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, &
geometry_plastic_nonlocal_disable
type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0
end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
@ -154,7 +160,7 @@ submodule(phase:plastic) nonlocal
state, &
state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
@ -170,7 +176,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, i, &
ph, &
Nconstituents, &
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
s1, s2, &
@ -203,29 +209,28 @@ module function plastic_nonlocal_init() result(myPlasticity)
print*, 'Kords, Dissertation RWTH Aachen, 2014'
print*, 'http://publications.rwth-aachen.de/record/229993'
allocate(param(Ninstances))
allocate(state(Ninstances))
allocate(state0(Ninstances))
allocate(dotState(Ninstances))
allocate(deltaState(Ninstances))
allocate(microstructure(Ninstances))
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
allocate(geom(phases%length))
allocate(param(phases%length))
allocate(state(phases%length))
allocate(state0(phases%length))
allocate(dotState(phases%length))
allocate(deltaState(phases%length))
allocate(microstructure(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), &
st0 => state0(ph), del => deltaState(ph), dst => microstructure(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
if(.not. myPlasticity(p)) cycle
i = i + 1
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i), &
st0 => state0(i), &
del => deltaState(i), &
dst => microstructure(i))
pl => mech%get('plasticity')
phase_localPlasticity(p) = .not. pl%contains('nonlocal')
phase_localPlasticity(ph) = .not. pl%contains('nonlocal')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(pl)
@ -236,8 +241,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0e4_pReal)
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%nu = lattice_nu(p)
prm%mu = lattice_mu(ph)
prm%nu = lattice_nu(ph)
ini%N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl))
@ -393,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -407,98 +412,101 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
allocate(geom(ph)%V_0(Nconstituents))
call storeGeometry(ph)
plasticState(ph)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
call IO_error(212,ext_msg='IPneighborhood does not exist')
plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(p)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho => plasticState(p)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
plasticState(p)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
stt%rhoSgl => plasticState(p)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSgl => plasticState(p)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSgl => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSglMobile => plasticState(p)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rhoSglMobile => plasticState(p)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rhoSglMobile => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rhoSglMobile => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rhoSglMobile => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rhoSglMobile => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_pos => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_neg => plasticState(ph)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_neg => plasticState(ph)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_neg => plasticState(ph)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_pos => plasticState(ph)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_pos => plasticState(ph)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_pos => plasticState(ph)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_neg => plasticState(ph)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rhoSglImmobile => plasticState(p)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSglImmobile => plasticState(p)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSglImmobile => plasticState(p)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSglImmobile => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSglImmobile => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSglImmobile => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_pos => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_neg => plasticState(ph)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_neg => plasticState(ph)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_neg => plasticState(ph)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_pos => plasticState(ph)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_pos => plasticState(ph)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_pos => plasticState(ph)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_neg => plasticState(ph)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoDip => plasticState(p)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rhoDip => plasticState(p)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rhoDip => plasticState(p)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rhoDip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rhoDip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rhoDip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho_dip_edg => plasticState(p)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
dot%rho_dip_edg => plasticState(p)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
del%rho_dip_edg => plasticState(p)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
del%rho_dip_edg => plasticState(ph)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
stt%rho_dip_scr => plasticState(p)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho_dip_scr => plasticState(ph)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
end associate
if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i)
plasticState(p)%state0 = plasticState(p)%state
if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents)
plasticState(ph)%state0 = plasticState(ph)%state
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
@ -510,40 +518,38 @@ module function plastic_nonlocal_init() result(myPlasticity)
discretization_nIPs,discretization_Nelems), source=0.0_pReal)
! BEGIN DEPRECATED----------------------------------------------------------------------------------
allocate(iRhoU(maxval(param%sum_N_sl),4,Ninstances), source=0)
allocate(iV(maxval(param%sum_N_sl),4,Ninstances), source=0)
allocate(iD(maxval(param%sum_N_sl),2,Ninstances), source=0)
allocate(iRhoU(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0)
i = 0
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
if(.not. myPlasticity(p)) cycle
i = i + 1
if(.not. myPlasticity(ph)) cycle
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
phase => phases%get(ph)
Nconstituents = count(material_phaseAt2 == ph)
l = 0
do t = 1,4
do s = 1,param(i)%sum_N_sl
do s = 1,param(ph)%sum_N_sl
l = l + 1
iRhoU(s,t,i) = l
iRhoU(s,t,ph) = l
enddo
enddo
l = l + (4+2+1+1)*param(i)%sum_N_sl ! immobile(4), dipole(2), shear, forest
l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest
do t = 1,4
do s = 1,param(i)%sum_N_sl
do s = 1,param(ph)%sum_N_sl
l = l + 1
iV(s,t,i) = l
iV(s,t,ph) = l
enddo
enddo
do t = 1,2
do s = 1,param(i)%sum_N_sl
do s = 1,param(ph)%sum_N_sl
l = l + 1
iD(s,t,i) = l
iD(s,t,ph) = l
enddo
enddo
if (iD(param(i)%sum_N_sl,2,i) /= plasticState(p)%sizeState) &
call IO_error(0, ext_msg = 'state indices not properly set (nonlocal)')
if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) &
error stop 'state indices not properly set (nonlocal)'
enddo
end function plastic_nonlocal_init
@ -552,20 +558,18 @@ end function plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
module subroutine nonlocal_dependentState(instance, me, ip, el)
module subroutine nonlocal_dependentState(ph, me, ip, el)
integer, intent(in) :: &
instance, &
ph, &
me, &
ip, &
el
integer :: &
ph, &
no, & !< neighbor offset
neighbor_el, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point
neighbor_instance, & ! instance of this plasticity of neighboring material point
c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index
dir, &
@ -589,29 +593,29 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
invConnections
real(pReal), dimension(3,nIPneighbors) :: &
connection_latticeConf
real(pReal), dimension(2,param(instance)%sum_N_sl) :: &
real(pReal), dimension(2,param(ph)%sum_N_sl) :: &
rhoExcess
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
rho_edg_delta, &
rho_scr_delta
real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho, &
rho0, &
rho_neighbor0
real(pReal), dimension(param(instance)%sum_N_sl,param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: &
myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(param(instance)%sum_N_sl,nIPneighbors) :: &
real(pReal), dimension(param(ph)%sum_N_sl,nIPneighbors) :: &
rho_edg_delta_neighbor, &
rho_scr_delta_neighbor
real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: &
neighbor_rhoExcess, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,param(instance)%sum_N_sl,2) :: &
real(pReal), dimension(3,param(ph)%sum_N_sl,2) :: &
m ! direction of dislocation motion
associate(prm => param(instance),dst => microstructure(instance), stt => state(instance))
associate(prm => param(ph),dst => microstructure(ph), stt => state(ph))
rho = getRho(instance,me,ip,el)
rho = getRho(ph,me)
stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+ matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
@ -639,9 +643,8 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
! ToDo: MD: this is most likely only correct for F_i = I
!#################################################################################################
rho0 = getRho0(instance,me,ip,el)
rho0 = getRho0(ph,me)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
ph = material_phaseAt(1,el)
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me))
@ -651,7 +654,7 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
rhoExcess(1,:) = rho_edg_delta
rhoExcess(2,:) = rho_scr_delta
FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal)
FVsize = geom(ph)%V_0(me) ** (1.0_pReal/3.0_pReal)
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
@ -662,11 +665,10 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
neighbor_ip = IPneighborhood(2,n,ip,el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el))
if (neighbor_instance == instance) then
if (material_phaseAt(1,neighbor_el) == ph) then
nRealNeighbors = nRealNeighbors + 1.0_pReal
rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el)
rho_neighbor0 = getRho0(ph,no)
rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
@ -758,16 +760,14 @@ end subroutine nonlocal_dependentState
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,ph,me,ip,el)
Mp,Temperature,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
integer, intent(in) :: &
ph, &
me, &
ip, & !< current integration point
el !< current element number
me
real(pReal), intent(in) :: &
Temperature !< temperature
@ -782,25 +782,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
l, &
t, & !< dislocation type
s !< index of my current slip system
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: &
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: &
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: &
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate
associate(prm => param(phase_plasticInstance(ph)),dst=>microstructure(phase_plasticInstance(ph)),&
stt=>state(phase_plasticInstance(ph)))
associate(prm => param(ph),dst=>microstructure(ph),&
stt=>state(ph))
ns = prm%sum_N_sl
!*** shortcut to state variables
rho = getRho(phase_plasticInstance(ph),me,ip,el)
rho = getRho(ph,me)
rhoSgl = rho(:,sgl)
do s = 1,ns
@ -820,7 +820,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
! edges
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticInstance(ph))
tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, ph)
v(:,2) = v(:,1)
dv_dtau(:,2) = dv_dtau(:,1)
dv_dtauNS(:,2) = dv_dtauNS(:,1)
@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
else
do t = 3,4
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticInstance(ph))
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, ph)
enddo
endif
@ -866,47 +866,42 @@ end subroutine nonlocal_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
integer, intent(in) :: &
instance, & ! current instance of this plasticity
me, & !< offset
ip, &
el
ph, &
me
integer :: &
ph, & !< phase
ns, & ! short notation for the total number of active slip systems
c, & ! character of dislocation
t, & ! type of dislocation
s ! index of my current slip system
real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho ! current dislocation densities
real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v ! dislocation glide velocity
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau ! current resolved shear stress
real(pReal), dimension(param(instance)%sum_N_sl,2) :: &
real(pReal), dimension(param(ph)%sum_N_sl,2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
dUpper, & ! current maximum stable dipole distance for edges and screws
dUpperOld, & ! old maximum stable dipole distance for edges and screws
deltaDUpper ! change in maximum stable dipole distance for edges and screws
ph = material_phaseAt(1,el)
associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance))
associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph))
ns = prm%sum_N_sl
!*** shortcut to state variables
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me)
forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),me)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me)
forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),me)
rho = getRho(instance,me,ip,el)
rho = getRho(ph,me)
rhoDip = rho(:,dip)
!****************************************************************************
@ -951,20 +946,11 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
/ (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),me) = dUpper(s,c)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),me) = dUpper(s,c)
plasticState(ph)%deltaState(:,me) = 0.0_pReal
del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
#ifdef DEBUG
if (debugConstitutive%extensive &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then
print'(a,/,8(12x,12(e12.5,1x),/))', '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8)
print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
endif
#endif
end associate
end subroutine plastic_nonlocal_deltaState
@ -992,7 +978,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
c, & !< character of dislocation
t, & !< type of dislocation
s !< index of my current slip system
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: &
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho, &
rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution
@ -1000,17 +986,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: &
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: &
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity
v0, &
gdot !< shear rates
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,2) :: &
real(pReal), dimension(param(ph)%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws
@ -1022,22 +1008,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
return
endif
associate(prm => param(phase_plasticInstance(ph)), &
dst => microstructure(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticInstance(ph)), &
stt => state(phase_plasticInstance(ph)))
associate(prm => param(ph), &
dst => microstructure(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
tau = 0.0_pReal
gdot = 0.0_pReal
rho = getRho(phase_plasticInstance(ph),me,ip,el)
rho = getRho(ph,me)
rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip)
rho0 = getRho0(phase_plasticInstance(ph),me,ip,el)
rho0 = getRho0(ph,me)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticInstance(ph)),me)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG
@ -1086,7 +1072,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
* sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticInstance(ph)),me)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me)
!****************************************************************************
@ -1142,7 +1128,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, phase_plasticInstance(ph),me,ip,el) &
rhoDot = rhoDotFlux(timestep, ph,me,ip,el) &
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
@ -1171,19 +1157,18 @@ end subroutine nonlocal_dotState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
function rhoDotFlux(timestep,instance,me,ip,el)
function rhoDotFlux(timestep,ph,me,ip,el)
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
instance, &
ph, &
me, &
ip, & !< current integration point
el !< current element number
integer :: &
ph, &
neighbor_instance, & !< instance of my neighbor's plasticity
neighbor_ph, & !< phase of my neighbor's plasticity
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
n, & !< index of my current neighbor
@ -1199,20 +1184,20 @@ function rhoDotFlux(timestep,instance,me,ip,el)
np,& !< neighbor phase shortcut
topp, & !< type of dislocation with opposite sign to t
s !< index of my current slip system
real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho, &
rho0, & !< dislocation density at beginning of time step
rhoDotFlux !< density evolution by flux
real(pReal), dimension(param(instance)%sum_N_sl,8) :: &
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity
v0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates
real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: &
real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: &
m !< direction of dislocation motion
real(pReal), dimension(3,3) :: &
my_F, & !< my total deformation gradient
@ -1230,26 +1215,25 @@ function rhoDotFlux(timestep,instance,me,ip,el)
transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point
lineLength !< dislocation line length leaving the current interface
ph = material_phaseAt(1,el)
associate(prm => param(instance), &
dst => microstructure(instance), &
dot => dotState(instance), &
stt => state(instance))
associate(prm => param(ph), &
dst => microstructure(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
gdot = 0.0_pReal
rho = getRho(instance,me,ip,el)
rho = getRho(ph,me)
rhoSgl = rho(:,sgl)
rho0 = getRho0(instance,me,ip,el)
rho0 = getRho0(ph,me)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) !ToDo: MD: I think we should use state0 here
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),me)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me)
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity)
@ -1301,7 +1285,7 @@ function rhoDotFlux(timestep,instance,me,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el))
neighbor_ph = material_phaseAt(1,neighbor_el)
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F)
@ -1324,8 +1308,8 @@ function rhoDotFlux(timestep,instance,me,ip,el)
any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then
forall (s = 1:ns, t = 1:4)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal)
endforall
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
@ -1408,12 +1392,12 @@ end function rhoDotFlux
! plane normals and signed cosine of the angle between the slip directions. Only the largest values
! that sum up to a total of 1 are considered, all others are set to zero.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
orientation ! crystal orientation
integer, intent(in) :: &
instance, &
ph, &
i, &
e
@ -1421,24 +1405,21 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
n, & ! neighbor index
neighbor_e, & ! element index of my neighbor
neighbor_i, & ! integration point index of my neighbor
ph, &
neighbor_phase, &
ns, & ! number of active slip systems
s1, & ! slip system index (me)
s2 ! slip system index (my neighbor)
real(pReal), dimension(2,param(instance)%sum_N_sl,param(instance)%sum_N_sl,nIPneighbors) :: &
real(pReal), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: &
my_compatibility ! my_compatibility for current element and ip
real(pReal) :: &
my_compatibilitySum, &
thresholdValue, &
nThresholdValues
logical, dimension(param(instance)%sum_N_sl) :: &
logical, dimension(param(ph)%sum_N_sl) :: &
belowThreshold
type(rotation) :: mis
ph = material_phaseAt(1,e)
associate(prm => param(instance))
associate(prm => param(ph))
ns = prm%sum_N_sl
!*** start out fully compatible
@ -1524,14 +1505,14 @@ end subroutine plastic_nonlocal_updateCompatibility
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_results(instance,group)
module subroutine plastic_nonlocal_results(ph,group)
integer, intent(in) :: instance
integer, intent(in) :: ph
character(len=*),intent(in) :: group
integer :: o
associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance))
associate(prm => param(ph),dst => microstructure(ph),stt=>state(ph))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('rho_u_ed_pos')
@ -1595,17 +1576,16 @@ end subroutine plastic_nonlocal_results
!--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density
!--------------------------------------------------------------------------------------------------
subroutine stateInit(ini,phase,Nconstituents,instance)
subroutine stateInit(ini,phase,Nconstituents)
type(tInitialParameters) :: &
ini
integer,intent(in) :: &
phase, &
Nconstituents, &
instance
Nconstituents
integer :: &
e, &
i, &
e, &
f, &
from, &
upto, &
@ -1623,7 +1603,7 @@ subroutine stateInit(ini,phase,Nconstituents,instance)
volume
associate(stt => state(instance))
associate(stt => state(phase))
if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
do e = 1,discretization_Nelems
@ -1671,18 +1651,18 @@ end subroutine stateInit
!--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, ph)
integer, intent(in) :: &
c, & !< dislocation character (1:edge, 2:screw)
instance
ph
real(pReal), intent(in) :: &
Temperature !< temperature
real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: &
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including non Schmid effects)
tauThreshold !< threshold shear stress
real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
v, & !< velocity
dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
@ -1713,7 +1693,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
criticalStress_S, & !< maximum obstacle strength
mobility !< dislocation mobility
associate(prm => param(instance))
associate(prm => param(ph))
ns = prm%sum_N_sl
v = 0.0_pReal
dv_dtau = 0.0_pReal
@ -1786,20 +1766,21 @@ end subroutine kinetics
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho(instance,me,ip,el)
pure function getRho(ph,me)
integer, intent(in) :: instance, me,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho
integer, intent(in) :: ph, me
real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho
associate(prm => param(instance))
getRho = reshape(state(instance)%rho(:,me),[prm%sum_N_sl,10])
associate(prm => param(ph))
getRho = reshape(state(ph)%rho(:,me),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign)
getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
getRho(:,dip) = max(getRho(:,dip),0.0_pReal)
where(abs(getRho) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho = 0.0_pReal
end associate
@ -1811,24 +1792,46 @@ end function getRho
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho0(instance,me,ip,el)
pure function getRho0(ph,me)
integer, intent(in) :: instance, me,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0
integer, intent(in) :: ph, me
real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho0
associate(prm => param(instance))
getRho0 = reshape(state0(instance)%rho(:,me),[prm%sum_N_sl,10])
associate(prm => param(ph))
getRho0 = reshape(state0(ph)%rho(:,me),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
where(abs(getRho0) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho0 = 0.0_pReal
end associate
end function getRho0
subroutine storeGeometry(ph)
integer, intent(in) :: ph
integer :: ip, el, ce, co
ce = 0
do el = 1, size(material_homogenizationMemberAt,2)
do ip = 1, size(material_homogenizationMemberAt,1)
ce = ce + 1
do co = 1, homogenization_maxNconstituents
if(material_phaseAt2(co,ce) == ph) then
geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = IPvolume(ip,el)
endif
enddo
enddo
enddo
end subroutine
end submodule nonlocal

120
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -70,8 +70,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, i, &
ph, i, &
Nconstituents, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -91,27 +90,24 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
myPlasticity = plastic_active('phenopowerlaw')
Ninstances = count(myPlasticity)
if(Ninstances == 0) return
if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>'
print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
allocate(param(Ninstances))
allocate(state(Ninstances))
allocate(dotState(Ninstances))
phases => config_material%get('phase')
i = 0
do p = 1, phases%length
phase => phases%get(p)
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
if(.not. myPlasticity(p)) cycle
i = i + 1
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i))
pl => mech%get('plasticity')
!--------------------------------------------------------------------------------------------------
@ -227,49 +223,49 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_0_sl, 2, Nconstituents)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_0_tw, 2, Nconstituents)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
stt%gamma_slip => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1
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) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
@ -302,18 +298,18 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: &
i,k,l,m,n
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
real(pReal), dimension(param(ph)%sum_N_tw) :: &
gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticInstance(ph)))
associate(prm => param(ph))
call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%sum_N_sl
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) &
@ -322,7 +318,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+ dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
call kinetics_twin(Mp,phase_plasticInstance(ph),me,gdot_twin,dgdot_dtautwin)
call kinetics_twin(Mp,ph,me,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%sum_N_tw
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) &
@ -350,12 +346,12 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
sumGamma,sumF
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticInstance(ph)))
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph))
sumGamma = sum(stt%gamma_slip(:,me))
sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char)
@ -375,9 +371,9 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
!--------------------------------------------------------------------------------------------------
! shear rates
call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg)
call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,phase_plasticInstance(ph),me,dot%gamma_twin(:,me))
call kinetics_twin(Mp,ph,me,dot%gamma_twin(:,me))
!--------------------------------------------------------------------------------------------------
! hardening
@ -395,14 +391,14 @@ end subroutine phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_results(instance,group)
module subroutine plastic_phenopowerlaw_results(ph,group)
integer, intent(in) :: instance
integer, intent(in) :: ph
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
@ -434,28 +430,28 @@ end subroutine plastic_phenopowerlaw_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,me, &
pure subroutine kinetics_slip(Mp,ph,me, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
gdot_slip_pos, &
gdot_slip_neg
real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
dgdot_dtau_slip_pos, &
dgdot_dtau_slip_neg
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau_slip_pos, &
tau_slip_neg
integer :: i
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
do i = 1, prm%sum_N_sl
tau_slip_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i))
@ -503,25 +499,25 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,instance,me,&
pure subroutine kinetics_twin(Mp,ph,me,&
gdot_twin,dgdot_dtau_twin)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
ph, &
me
real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: &
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
gdot_twin
real(pReal), dimension(param(instance)%sum_N_tw), intent(out), optional :: &
real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
dgdot_dtau_twin
real(pReal), dimension(param(instance)%sum_N_tw) :: &
real(pReal), dimension(param(ph)%sum_N_tw) :: &
tau_twin
integer :: i
associate(prm => param(instance), stt => state(instance))
associate(prm => param(ph), stt => state(ph))
do i = 1, prm%sum_N_tw
tau_twin(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))