diff --git a/PRIVATE b/PRIVATE index 162106e63..5c93094df 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 162106e6379d484ee101981c66e3f159d2f8821a +Subproject commit 5c93094df289da74588a3b5dd644d31923c9c651 diff --git a/src/phase.f90 b/src/phase.f90 index 573e71069..e264f54ac 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -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') diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 0f931517a..f131c552b 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -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)