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