numerical parameters related to phase state and stress integration

This commit is contained in:
Sharan Roongta 2023-07-21 23:41:12 +02:00
parent d104f4810e
commit 59ae287564
3 changed files with 64 additions and 51 deletions

@ -1 +1 @@
Subproject commit 162106e6379d484ee101981c66e3f159d2f8821a
Subproject commit 5c93094df289da74588a3b5dd644d31923c9c651

View File

@ -61,6 +61,7 @@ module phase
type :: tNumerics
integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
iJacoLiresiduum, & !< frequency of Jacobian update of residuum in Li
nState, & !< state loop limit
nStress !< stress loop limit
real(pREAL) :: &
@ -69,9 +70,11 @@ module phase
subStepSizeLp, & !< size of first substep when cutback in Lp calculation
subStepSizeLi, & !< size of first substep when cutback in Li calculation
stepIncreaseCryst, & !< increase of next substep size when previous substep converged
rtol_crystalliteState, & !< relative tolerance in state loop
rtol_crystalliteStress, & !< relative tolerance in stress loop
atol_crystalliteStress !< absolute tolerance in stress loop
rtol_crystalliteState, &
rtol_Lp, & !< relative tolerance in stress loop for Lp
atol_Lp, & !< absolute tolerance in stress loop for Lp
rtol_Li, & !< relative tolerance in stress loop for Li
atol_Li !< absolute tolerance in stress loop for Li
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
@ -85,8 +88,8 @@ module phase
interface
! == cleaned:begin =================================================================================
module subroutine mechanical_init(phases)
type(tDict), pointer :: phases
module subroutine mechanical_init(phases,num_mech)
type(tDict), pointer :: phases, num_mech
end subroutine mechanical_init
module subroutine damage_init
@ -381,7 +384,9 @@ subroutine phase_init
ph, ce, co, ma
type(tDict), pointer :: &
phases, &
phase
phase, &
num_phase, &
num_mech
character(len=:), allocatable :: refs
@ -420,7 +425,10 @@ subroutine phase_init
phase_O(ph)%data = phase_O_0(ph)%data
end do
call mechanical_init(phases)
num_phase => config_numerics%get_dict('phase',defaultVal=emptyDict)
num_mech => num_phase%get_dict('mechanical', defaultVal=emptyDict)
call mechanical_init(phases,num_mech)
call damage_init()
call thermal_init(phases)
@ -531,39 +539,8 @@ subroutine crystallite_init()
el, & !< counter in element loop
en, ph
type(tDict), pointer :: &
num_crystallite, &
num_phase, &
phases
character(len=:), allocatable :: extmsg
num_crystallite => config_numerics%get_dict('crystallite',defaultVal=emptyDict)
num%subStepMinCryst = num_crystallite%get_asReal ('subStepMin', defaultVal=1.0e-3_pREAL)
num%subStepSizeCryst = num_crystallite%get_asReal ('subStepSize', defaultVal=0.25_pREAL)
num%stepIncreaseCryst = num_crystallite%get_asReal ('stepIncrease', defaultVal=1.5_pREAL)
num%subStepSizeLp = num_crystallite%get_asReal ('subStepSizeLp', defaultVal=0.5_pREAL)
num%subStepSizeLi = num_crystallite%get_asReal ('subStepSizeLi', defaultVal=0.5_pREAL)
num%rtol_crystalliteState = num_crystallite%get_asReal ('rtol_State', defaultVal=1.0e-6_pREAL)
num%rtol_crystalliteStress = num_crystallite%get_asReal ('rtol_Stress', defaultVal=1.0e-6_pREAL)
num%atol_crystalliteStress = num_crystallite%get_asReal ('atol_Stress', defaultVal=1.0e-8_pREAL)
num%iJacoLpresiduum = num_crystallite%get_asInt ('iJacoLpresiduum', defaultVal=1)
num%nState = num_crystallite%get_asInt ('nState', defaultVal=20)
num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40)
extmsg = ''
if (num%subStepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' subStepMinCryst'
if (num%subStepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' subStepSizeCryst'
if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' stepIncreaseCryst'
if (num%subStepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' subStepSizeLp'
if (num%subStepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' subStepSizeLi'
if (num%rtol_crystalliteState <= 0.0_pREAL) extmsg = trim(extmsg)//' rtol_crystalliteState'
if (num%rtol_crystalliteStress <= 0.0_pREAL) extmsg = trim(extmsg)//' rtol_crystalliteStress'
if (num%atol_crystalliteStress <= 0.0_pREAL) extmsg = trim(extmsg)//' atol_crystalliteStress'
if (num%iJacoLpresiduum < 1) extmsg = trim(extmsg)//' iJacoLpresiduum'
if (num%nState < 1) extmsg = trim(extmsg)//' nState'
if (num%nStress < 1) extmsg = trim(extmsg)//' nStress'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
phases => config_material%get_dict('phase')

View File

@ -198,10 +198,11 @@ contains
!> @brief Initialize mechanical field related constitutive models
!> @details Initialize elasticity, plasticity and stiffness degradation models.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_init(phases)
module subroutine mechanical_init(phases, num_mech)
type(tDict), pointer :: &
phases
phases, &
num_mech
integer :: &
ce, &
@ -211,9 +212,11 @@ module subroutine mechanical_init(phases)
en, &
Nmembers
type(tDict), pointer :: &
num_crystallite, &
phase, &
mech
mech, &
num_mech_plastic, &
num_mech_eigen
character(len=:), allocatable :: extmsg
print'(/,1x,a)', '<<<+- phase:mechanical init -+>>>'
@ -289,9 +292,42 @@ module subroutine mechanical_init(phases)
plasticState(ph)%state0 = plasticState(ph)%state
end do
num_crystallite => config_numerics%get_dict('crystallite',defaultVal=emptyDict)
num_mech_plastic => num_mech%get_dict('plastic', defaultVal=emptyDict)
num_mech_eigen => num_mech%get_dict('eigen', defaultVal=emptyDict)
select case(num_crystallite%get_asStr('integrator',defaultVal='FPI'))
num%subStepMinCryst = num_mech%get_asReal ('sub_step_min', defaultVal=1.0e-3_pREAL)
num%subStepSizeCryst = num_mech%get_asReal ('sub_step_size', defaultVal=0.25_pREAL)
num%stepIncreaseCryst = num_mech%get_asReal ('step_increase', defaultVal=1.5_pREAL)
num%rtol_crystalliteState = num_mech%get_asReal ('eps_rel_state', defaultVal=1.0e-6_pREAL)
num%nState = num_mech%get_asInt ('N_iter_state_max', defaultVal=20)
num%nStress = num_mech%get_asInt ('N_iter_stress_max', defaultVal=40)
num%subStepSizeLp = num_mech_plastic%get_asReal ('sub_step_size_Lp', defaultVal=0.5_pREAL)
num%rtol_Lp = num_mech_plastic%get_asReal ('eps_rel_Lp', defaultVal=1.0e-6_pREAL)
num%atol_Lp = num_mech_plastic%get_asReal ('eps_abs_Lp', defaultVal=1.0e-8_pREAL)
num%iJacoLpresiduum = num_mech_plastic%get_asInt ('f_update_jacobi_Lp', defaultVal=1)
num%subStepSizeLi = num_mech_eigen%get_asReal ('sub_step_size_Li', defaultVal=0.5_pREAL)
num%rtol_Li = num_mech_eigen%get_asReal ('eps_rel_Li', defaultVal=num%rtol_Lp)
num%atol_Li = num_mech_eigen%get_asReal ('eps_abs_Li', defaultVal=num%atol_Lp)
num%iJacoLiresiduum = num_mech_eigen%get_asInt ('f_update_jacobi_Li', defaultVal=num%iJacoLpresiduum)
extmsg = ''
if (num%subStepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_min'
if (num%subStepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size'
if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' step_increase'
if (num%subStepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size_Lp'
if (num%subStepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size_Li'
if (num%rtol_Lp <= 0.0_pREAL) extmsg = trim(extmsg)//' epl_rel_Lp'
if (num%atol_Lp <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_Lp'
if (num%rtol_Li <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_Li'
if (num%atol_Li <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_Li'
if (num%iJacoLpresiduum < 1) extmsg = trim(extmsg)//' f_jacobi_residuum_update_Lp'
if (num%iJacoLiresiduum < 1) extmsg = trim(extmsg)//' f_jacobi_residuum_update_Li'
if (num%nState < 1) extmsg = trim(extmsg)//' N_iter_state_max'
if (num%nStress < 1) extmsg = trim(extmsg)//' N_iter_stress_max'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
select case(num_mech_plastic%get_asStr('integrator_state',defaultVal='FPI'))
case('FPI')
integrateState => integrateStateFPI
@ -458,8 +494,8 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
S, Fi_new, ph,en)
!* 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
num%atol_crystalliteStress) ! minimum lower cutoff
atol_Lp = max(num%rtol_Lp * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%atol_Lp) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
if (any(IEEE_is_NaN(residuumLp))) then
@ -499,8 +535,8 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
S, Fi_new, ph,en)
!* update current residuum and check for convergence of loop
atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%atol_crystalliteStress) ! minimum lower cutoff
atol_Li = max(num%rtol_Li * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%atol_Li) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive
if (any(IEEE_is_NaN(residuumLi))) then
return ! error
@ -517,7 +553,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
cycle LiLoop
end if
calculateJacobiLi: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
calculateJacobiLi: if (mod(jacoCounterLi, num%iJacoLiresiduum) == 0) then
jacoCounterLi = jacoCounterLi + 1
temp_33 = matmul(matmul(A,B),invFi_current)