From 59ae2875640fbb6278a14b7057c0353e6ac169a5 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 21 Jul 2023 23:41:12 +0200 Subject: [PATCH 1/6] numerical parameters related to phase state and stress integration --- PRIVATE | 2 +- src/phase.f90 | 55 +++++++++++-------------------------- src/phase_mechanical.f90 | 58 ++++++++++++++++++++++++++++++++-------- 3 files changed, 64 insertions(+), 51 deletions(-) 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) From 33e0048010d4a7b9cf67aae777d78cff6d6eec0a Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Mon, 24 Jul 2023 21:03:35 +0200 Subject: [PATCH 2/6] subStep ---> step --- src/phase.f90 | 8 ++-- src/phase_mechanical.f90 | 86 ++++++++++++++++++++-------------------- 2 files changed, 47 insertions(+), 47 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index e264f54ac..46b5c052e 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -65,10 +65,10 @@ module phase nState, & !< state loop limit nStress !< stress loop limit real(pREAL) :: & - subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback - subStepSizeCryst, & !< size of first substep when cutback - subStepSizeLp, & !< size of first substep when cutback in Lp calculation - subStepSizeLi, & !< size of first substep when cutback in Li calculation + stepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback + stepSizeCryst, & !< size of first substep when cutback + stepSizeLp, & !< size of first substep when cutback in Lp calculation + stepSizeLi, & !< size of first substep when cutback in Li calculation stepIncreaseCryst, & !< increase of next substep size when previous substep converged rtol_crystalliteState, & rtol_Lp, & !< relative tolerance in stress loop for Lp diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index f131c552b..91b373a0e 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -295,35 +295,35 @@ module subroutine mechanical_init(phases, num_mech) num_mech_plastic => num_mech%get_dict('plastic', defaultVal=emptyDict) num_mech_eigen => num_mech%get_dict('eigen', defaultVal=emptyDict) - 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) + num%stepMinCryst = num_mech%get_asReal ('step_min', defaultVal=1.0e-3_pREAL) + num%stepSizeCryst = num_mech%get_asReal ('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%stepSizeLp = num_mech_plastic%get_asReal ('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%stepSizeLi = num_mech_eigen%get_asReal ('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 (num%stepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_min' + if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size' + if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' step_increase' + if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size_Lp' + if (num%stepSizeLi <= 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)) @@ -507,7 +507,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) Lpguess_old = Lpguess steplengthLp = 1.0_pREAL ! ...proceed with normal step length (calculate new search direction) else ! not converged and residuum not improved... - steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction + steplengthLp = num%stepSizeLp * steplengthLp ! ...try with smaller step length in same direction Lpguess = Lpguess_old & + deltaLp * stepLengthLp cycle LpLoop @@ -547,7 +547,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) Liguess_old = Liguess steplengthLi = 1.0_pREAL ! ...proceed with normal step length (calculate new search direction) else ! not converged and residuum not improved... - steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction + steplengthLi = num%stepSizeLi * steplengthLi ! ...try with smaller step length in same direction Liguess = Liguess_old & + deltaLi * steplengthLi cycle LiLoop @@ -1029,11 +1029,11 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) logical :: converged_ real(pREAL) :: & - formerSubStep + formerStep integer :: & ph, en, sizeDotState logical :: todo - real(pREAL) :: subFrac,subStep + real(pREAL) :: stepFrac,step real(pREAL), dimension(3,3) :: & subFp0, & subFi0, & @@ -1053,20 +1053,20 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) - subFrac = 0.0_pREAL + stepFrac = 0.0_pREAL todo = .true. - subStep = 1.0_pREAL/num%subStepSizeCryst - converged_ = .false. ! pretend failed step of 1/subStepSizeCryst + step = 1.0_pREAL/num%stepSizeCryst + converged_ = .false. ! pretend failed step of 1/stepSizeCryst todo = .true. cutbackLooping: do while (todo) if (converged_) then - formerSubStep = subStep - subFrac = subFrac + subStep - subStep = min(1.0_pREAL - subFrac, num%stepIncreaseCryst * subStep) + formerStep = step + stepFrac = stepFrac + step + step = min(1.0_pREAL - stepFrac, num%stepIncreaseCryst * step) - todo = subStep > 0.0_pREAL ! still time left to integrate on? + todo = step > 0.0_pREAL ! still time left to integrate on? if (todo) then subF0 = subF @@ -1079,16 +1079,16 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) else - subStep = num%subStepSizeCryst * subStep + step = num%stepSizeCryst * step phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) - if (subStep < 1.0_pREAL) then ! actual (not initial) cutback + if (step < 1.0_pREAL) then ! actual (not initial) cutback phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 end if plasticState(ph)%state(:,en) = subState0 - todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) + todo = step > num%stepMinCryst ! still on track or already done (beyond repair) end if !-------------------------------------------------------------------------------------------------- @@ -1096,8 +1096,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) if (todo) then sizeDotState = plasticState(ph)%sizeDotState subF = subF0 & - + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en) + + step * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),step * Delta_t,ph,en) end if end do cutbackLooping From c4134f2358bcd0a7ba9b079bfdbb6648b936bffd Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 25 Jul 2023 11:04:38 +0200 Subject: [PATCH 3/6] updated examples --- examples/config/numerics.yaml | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index a0c7cb6b0..dcf330073 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -71,19 +71,26 @@ mesh: eps_struct_atol: 1.0e-10 # absolute tolerance for mechanical equilibrium eps_struct_rtol: 1.0e-4 # relative tolerance for mechanical equilibrium -crystallite: - subStepMin: 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in crystallite - subStepSize: 0.25 # size of substep when cutback introduced in crystallite (value between 0 and 1) - stepIncrease: 1.5 # increase of next substep size when previous substep converged in crystallite (value higher than 1) - subStepSizeLp: 0.5 # size of first substep when cutback in Lp calculation - subStepSizeLi: 0.5 # size of first substep when cutback in Li calculation - nState: 10 # state loop limit - nStress: 40 # stress loop limit - rtol_State: 1.0e-6 # relative tolerance in crystallite state loop (abs tol provided by constitutive law) - rtol_Stress: 1.0e-6 # relative tolerance in crystallite stress loop (Lp residuum) - atol_Stress: 1.0e-8 # absolute tolerance in crystallite stress loop (Lp residuum!) - integrator: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) - iJacoLpresiduum: 1 # frequency of Jacobian update of residuum in Lp +phase: + mechanical: + step_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in crystallite + step_size: 0.25 # size of step when cutback introduced in crystallite (value between 0 and 1) + step_increase: 1.5 # increase of next step size when previous step converged in crystallite (value higher than 1) + eps_rel_state: 1.0e-6 # relative tolerance in crystallite state loop (abs tol provided by constitutive law) + N_iter_state_max: 10 # state loop limit + N_iter_stress_max: 40 # stress loop limit + + plastic: + step_size_Lp: 0.5 # size of first substep when cutback in Lp calculation + eps_rel_Lp: 1.0e-6 # relative tolerance in crystallite stress loop (Lp residuum) + eps_abs_Lp: 1.0e-8 # absolute tolerance in crystallite stress loop (Lp residuum!) + f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp + integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) + eigen: + step_size_Li: 0.5 # size of first substep when cutback in Li calculation + eps_rel_Li: 1.0e-6 # relative tolerance in crystallite stress loop (Li residuum) + eps_abs_Li: 1.0e-8 # absolute tolerance in crystallite stress loop (Li residuum!) + f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li commercialFEM: unitlength: 1 # physical length of one computational length unit From 2547d4a25cfc2a9082687e04e23de3fd7d95eda7 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 25 Jul 2023 23:04:29 +0200 Subject: [PATCH 4/6] clear names --- examples/config/numerics.yaml | 23 ++++++++++++----------- src/phase.f90 | 3 ++- src/phase_mechanical.f90 | 30 ++++++++++++++++-------------- 3 files changed, 30 insertions(+), 26 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index dcf330073..ac37e7e08 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -73,23 +73,24 @@ mesh: phase: mechanical: - step_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in crystallite - step_size: 0.25 # size of step when cutback introduced in crystallite (value between 0 and 1) - step_increase: 1.5 # increase of next step size when previous step converged in crystallite (value higher than 1) - eps_rel_state: 1.0e-6 # relative tolerance in crystallite state loop (abs tol provided by constitutive law) + step_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation + r_cutback_step: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) + r_increase_step: 1.5 # factor to increase size of next step when previous step converged in phase state calculation + eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) N_iter_state_max: 10 # state loop limit - N_iter_stress_max: 40 # stress loop limit plastic: - step_size_Lp: 0.5 # size of first substep when cutback in Lp calculation - eps_rel_Lp: 1.0e-6 # relative tolerance in crystallite stress loop (Lp residuum) - eps_abs_Lp: 1.0e-8 # absolute tolerance in crystallite stress loop (Lp residuum!) + r_cutback_step_Lp: 0.5 # factor to decrease the step when cutback in Lp calculation + eps_rel_Lp: 1.0e-6 # relative tolerance in phase stress loop (Lp residuum) + eps_abs_Lp: 1.0e-8 # absolute tolerance in phase stress loop (Lp residuum) + N_iter_Lp_max: 40 # stress loop limit for Lp f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) eigen: - step_size_Li: 0.5 # size of first substep when cutback in Li calculation - eps_rel_Li: 1.0e-6 # relative tolerance in crystallite stress loop (Li residuum) - eps_abs_Li: 1.0e-8 # absolute tolerance in crystallite stress loop (Li residuum!) + r_cutback_step_Li: 0.5 # factor to decrease the step when cutback in Li calculation + eps_rel_Li: 1.0e-6 # relative tolerance in phase stress loop (Li residuum) + eps_abs_Li: 1.0e-8 # absolute tolerance in phase stress loop (Li residuum) + N_iter_Li_max: 40 # stress loop limit for Li f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li commercialFEM: diff --git a/src/phase.f90 b/src/phase.f90 index c233049f5..4b01c6ef0 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -86,7 +86,8 @@ module phase 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 + nStress_Lp, & !< stress loop limit for Lp + nStress_Li !< stress loop limit for Li real(pREAL) :: & stepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback stepSizeCryst, & !< size of first substep when cutback diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index e65745121..ddd31cbcc 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -278,34 +278,36 @@ module subroutine mechanical_init(phases, num_mech) num_mech_eigen => num_mech%get_dict('eigen', defaultVal=emptyDict) num%stepMinCryst = num_mech%get_asReal ('step_min', defaultVal=1.0e-3_pREAL) - num%stepSizeCryst = num_mech%get_asReal ('step_size', defaultVal=0.25_pREAL) - num%stepIncreaseCryst = num_mech%get_asReal ('step_increase', defaultVal=1.5_pREAL) + num%stepSizeCryst = num_mech%get_asReal ('r_cutback_step', defaultVal=0.25_pREAL) + num%stepIncreaseCryst = num_mech%get_asReal ('r_increase_step', 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%stepSizeLp = num_mech_plastic%get_asReal ('step_size_Lp', defaultVal=0.5_pREAL) + num%nStress_Lp = num_mech_plastic%get_asInt ('N_iter_Lp_max', defaultVal=40) + num%stepSizeLp = num_mech_plastic%get_asReal ('r_cutback_step_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%stepSizeLi = num_mech_eigen%get_asReal ('step_size_Li', defaultVal=0.5_pREAL) + num%nStress_Li = num_mech_eigen%get_asInt ('N_iter_Li_max', defaultVal=40) + num%stepSizeLi = num_mech_eigen%get_asReal ('r_cutback_step_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%stepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_min' - if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size' - if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' step_increase' - if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size_Lp' - if (num%stepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_size_Li' + if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step' + if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_increase_step' + if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step_Lp' + if (num%stepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step_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%iJacoLpresiduum < 1) extmsg = trim(extmsg)//' f_update_jacobi_Lp' + if (num%iJacoLiresiduum < 1) extmsg = trim(extmsg)//' f_update_jacobi_Li' if (num%nState < 1) extmsg = trim(extmsg)//' N_iter_state_max' - if (num%nStress < 1) extmsg = trim(extmsg)//' N_iter_stress_max' + if (num%nStress_Lp < 1) extmsg = trim(extmsg)//' N_iter_Lp_max' + if (num%nStress_Li < 1) extmsg = trim(extmsg)//' N_iter_Li_max' if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg)) @@ -452,7 +454,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) NiterationStressLi = 0 LiLoop: do NiterationStressLi = NiterationStressLi + 1 - if (NiterationStressLi>num%nStress) return ! error + if (NiterationStressLi>num%nStress_Li) return ! error invFi_new = matmul(invFi_current,math_I3 - Delta_t*Liguess) Fi_new = math_inv33(invFi_new) @@ -465,7 +467,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) NiterationStressLp = 0 LpLoop: do NiterationStressLp = NiterationStressLp + 1 - if (NiterationStressLp>num%nStress) return ! error + if (NiterationStressLp>num%nStress_Lp) return ! error B = math_I3 - Delta_t*Lpguess Fe = matmul(matmul(A,B), invFi_new) From af2a4c6e7290853cc2d0e5c2ca6b6aba9ca81123 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 25 Jul 2023 23:13:00 +0200 Subject: [PATCH 5/6] more readable --- examples/config/numerics.yaml | 18 ++-- src/phase_mechanical.f90 | 150 +++++++++++++++++----------------- 2 files changed, 84 insertions(+), 84 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index ac37e7e08..97885311d 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -73,23 +73,23 @@ mesh: phase: mechanical: - step_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation - r_cutback_step: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) - r_increase_step: 1.5 # factor to increase size of next step when previous step converged in phase state calculation + r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation + r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) + r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) N_iter_state_max: 10 # state loop limit plastic: - r_cutback_step_Lp: 0.5 # factor to decrease the step when cutback in Lp calculation - eps_rel_Lp: 1.0e-6 # relative tolerance in phase stress loop (Lp residuum) - eps_abs_Lp: 1.0e-8 # absolute tolerance in phase stress loop (Lp residuum) + r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation + eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum + eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum N_iter_Lp_max: 40 # stress loop limit for Lp f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) eigen: - r_cutback_step_Li: 0.5 # factor to decrease the step when cutback in Li calculation - eps_rel_Li: 1.0e-6 # relative tolerance in phase stress loop (Li residuum) - eps_abs_Li: 1.0e-8 # absolute tolerance in phase stress loop (Li residuum) + r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation + eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum + eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum N_iter_Li_max: 40 # stress loop limit for Li f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index ddd31cbcc..e7103bedc 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -277,28 +277,28 @@ module subroutine mechanical_init(phases, num_mech) num_mech_plastic => num_mech%get_dict('plastic', defaultVal=emptyDict) num_mech_eigen => num_mech%get_dict('eigen', defaultVal=emptyDict) - num%stepMinCryst = num_mech%get_asReal ('step_min', defaultVal=1.0e-3_pREAL) - num%stepSizeCryst = num_mech%get_asReal ('r_cutback_step', defaultVal=0.25_pREAL) - num%stepIncreaseCryst = num_mech%get_asReal ('r_increase_step', defaultVal=1.5_pREAL) + num%stepMinCryst = num_mech%get_asReal ('r_cutback_min', defaultVal=1.0e-3_pREAL) + num%stepSizeCryst = num_mech%get_asReal ('r_cutback', defaultVal=0.25_pREAL) + num%stepIncreaseCryst = num_mech%get_asReal ('r_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_Lp = num_mech_plastic%get_asInt ('N_iter_Lp_max', defaultVal=40) - num%stepSizeLp = num_mech_plastic%get_asReal ('r_cutback_step_Lp', defaultVal=0.5_pREAL) + num%stepSizeLp = num_mech_plastic%get_asReal ('r_linesearch_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%nStress_Li = num_mech_eigen%get_asInt ('N_iter_Li_max', defaultVal=40) - num%stepSizeLi = num_mech_eigen%get_asReal ('r_cutback_step_Li', defaultVal=0.5_pREAL) + num%stepSizeLi = num_mech_eigen%get_asReal ('r_linesearch_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%stepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_min' - if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step' - if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_increase_step' - if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step_Lp' - if (num%stepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step_Li' + if (num%stepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_min' + if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback' + if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_increase' + if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' r_linesearch_Lp' + if (num%stepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' r_linesearch_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' @@ -377,9 +377,9 @@ end subroutine mechanical_result !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) +function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken) - real(pREAL), dimension(3,3), intent(in) :: F,subFp0,subFi0 + real(pREAL), dimension(3,3), intent(in) :: F,Fp0,Fi0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en @@ -439,9 +439,9 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess - call math_invert33(invFp_current,error=error,A=subFp0) + call math_invert33(invFp_current,error=error,A=Fp0) if (error) return ! error - call math_invert33(invFi_current,error=error,A=subFi0) + call math_invert33(invFi_current,error=error,A=Fi0) if (error) return ! error A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp @@ -582,10 +582,10 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) +function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) - real(pREAL), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 - real(pREAL), intent(in),dimension(:) :: subState0 + real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 + real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & ph, & @@ -611,14 +611,14 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState * Delta_t iteration: do NiterationState = 1, num%nState dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0_pREAL, nIterationState > 1) dotState_last(1:sizeDotState,1) = dotState - broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) + broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) if (broken) exit iteration dotState = plastic_dotState(Delta_t,ph,en) @@ -628,7 +628,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b dotState = dotState * zeta & + dotState_last(1:sizeDotState,1) * (1.0_pREAL - zeta) r = plasticState(ph)%state(1:sizeDotState,en) & - - subState0 & + - state0 & - dotState * Delta_t plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r @@ -670,10 +670,10 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) +function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) - real(pREAL), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 - real(pREAL), intent(in),dimension(:) :: subState0 + real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 + real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & ph, & @@ -694,15 +694,15 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result sizeDotState = plasticState(ph)%sizeDotState #ifndef __INTEL_LLVM_COMPILER - plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t + plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t #else - plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,state0) #endif broken = plastic_deltaState(ph,en) if (broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) + broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) end function integrateStateEuler @@ -710,10 +710,10 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) +function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) - real(pREAL), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 - real(pREAL), intent(in),dimension(:) :: subState0 + real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 + real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & ph, & @@ -737,15 +737,15 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en r = - dotState * 0.5_pREAL * Delta_t #ifndef __INTEL_LLVM_COMPILER - plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t + plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t #else - plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,state0) #endif broken = plastic_deltaState(ph,en) if (broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) + broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) if (broken) return dotState = plastic_dotState(Delta_t,ph,en) @@ -761,10 +761,10 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) +function integrateStateRK4(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) - real(pREAL), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 - real(pREAL), intent(in),dimension(:) :: subState0 + real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 + real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en logical :: broken @@ -781,7 +781,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b B = [6.0_pREAL, 3.0_pREAL, 3.0_pREAL, 6.0_pREAL]**(-1) - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C) + broken = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C) end function integrateStateRK4 @@ -789,10 +789,10 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) +function integrateStateRKCK45(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) - real(pREAL), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 - real(pREAL), intent(in),dimension(:) :: subState0 + real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 + real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en logical :: broken @@ -816,7 +816,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) resul 13525.0_pREAL/55296.0_pREAL, 277.0_pREAL/14336.0_pREAL, 1._pREAL/4._pREAL] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) + broken = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) end function integrateStateRKCK45 @@ -825,10 +825,10 @@ end function integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) result(broken) +function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(broken) - real(pREAL), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 - real(pREAL), intent(in),dimension(:) :: subState0 + real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 + real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t real(pREAL), dimension(:,:), intent(in) :: A real(pREAL), dimension(:), intent(in) :: B, C @@ -869,12 +869,12 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) end do #ifndef __INTEL_LLVM_COMPILER - plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t + plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t #else - plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,state0) #endif - broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en) + broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) if (broken) exit dotState = plastic_dotState(Delta_t*C(stage), ph,en) @@ -887,9 +887,9 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) plastic_RKdotState(1:sizeDotState,size(B)) = dotState dotState = matmul(plastic_RKdotState,B) #ifndef __INTEL_LLVM_COMPILER - plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t + plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t #else - plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,state0) #endif if (present(DB)) & @@ -902,7 +902,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) broken = plastic_deltaState(ph,en) if (broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) + broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) end function integrateStateRK @@ -1019,24 +1019,24 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) logical :: todo real(pREAL) :: stepFrac,step real(pREAL), dimension(3,3) :: & - subFp0, & - subFi0, & - subLp0, & - subLi0, & - subF0, & - subF - real(pREAL), dimension(plasticState(material_ID_phase(co,ce))%sizeState) :: subState0 + Fp0, & + Fi0, & + Lp0, & + Li0, & + F0, & + F + real(pREAL), dimension(plasticState(material_ID_phase(co,ce))%sizeState) :: state0 ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - subState0 = plasticState(ph)%state0(:,en) - subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) - subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) - subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) - subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) - subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) + state0 = plasticState(ph)%state0(:,en) + Li0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) + Lp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) + Fp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) + Fi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) + F0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) stepFrac = 0.0_pREAL todo = .true. step = 1.0_pREAL/num%stepSizeCryst @@ -1053,25 +1053,25 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) todo = step > 0.0_pREAL ! still time left to integrate on? if (todo) then - subF0 = subF - subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,en) - subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,en) - subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) - subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) - subState0 = plasticState(ph)%state(:,en) + F0 = F + Lp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,en) + Li0 = phase_mechanical_Li(ph)%data(1:3,1:3,en) + Fp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) + Fi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) + state0 = plasticState(ph)%state(:,en) end if !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) else step = num%stepSizeCryst * step - phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 - phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp0 + phase_mechanical_Fi(ph)%data(1:3,1:3,en) = Fi0 phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) if (step < 1.0_pREAL) then ! actual (not initial) cutback - phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 - phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 + phase_mechanical_Lp(ph)%data(1:3,1:3,en) = Lp0 + phase_mechanical_Li(ph)%data(1:3,1:3,en) = Li0 end if - plasticState(ph)%state(:,en) = subState0 + plasticState(ph)%state(:,en) = state0 todo = step > num%stepMinCryst ! still on track or already done (beyond repair) end if @@ -1079,9 +1079,9 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) ! prepare for integration if (todo) then sizeDotState = plasticState(ph)%sizeDotState - subF = subF0 & + F = F0 & + step * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),step * Delta_t,ph,en) + converged_ = .not. integrateState(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en) end if end do cutbackLooping From 776d75379cfe4007bbf0307ef8a9bc3b662b0986 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 26 Jul 2023 12:59:33 -0400 Subject: [PATCH 6/6] use modified numerics in test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 5c93094df..8c05965ef 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5c93094df289da74588a3b5dd644d31923c9c651 +Subproject commit 8c05965ef4437598898f467a213ffa88938e860a