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

more readable

This commit is contained in:
Sharan Roongta 2023-07-25 23:13:00 +02:00
parent 2547d4a25c
commit af2a4c6e72
2 changed files with 84 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
examples/config/numerics.yaml
View File

@ -73,23 +73,23 @@ mesh:
phase: phase:
mechanical: mechanical:
step_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation r_cutback_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_cutback: 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_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) 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_state_max: 10 # state loop limit
plastic: plastic:
r_cutback_step_Lp: 0.5 # factor to decrease the step when cutback in Lp calculation 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 phase stress loop (Lp residuum) eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum
eps_abs_Lp: 1.0e-8 # absolute tolerance in phase stress loop (Lp residuum) eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum
N_iter_Lp_max: 40 # stress loop limit for Lp N_iter_Lp_max: 40 # stress loop limit for Lp
f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in 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) 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: eigen:
r_cutback_step_Li: 0.5 # factor to decrease the step when cutback in Li calculation 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 phase stress loop (Li residuum) eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum
eps_abs_Li: 1.0e-8 # absolute tolerance in phase stress loop (Li residuum) eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum
N_iter_Li_max: 40 # stress loop limit for Li N_iter_Li_max: 40 # stress loop limit for Li
f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li

150
src/phase_mechanical.f90
View File

@ -277,28 +277,28 @@ module subroutine mechanical_init(phases, num_mech)
num_mech_plastic => num_mech%get_dict('plastic', defaultVal=emptyDict) num_mech_plastic => num_mech%get_dict('plastic', defaultVal=emptyDict)
num_mech_eigen => num_mech%get_dict('eigen', 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%stepMinCryst = num_mech%get_asReal ('r_cutback_min', defaultVal=1.0e-3_pREAL)
num%stepSizeCryst = num_mech%get_asReal ('r_cutback_step', defaultVal=0.25_pREAL) num%stepSizeCryst = num_mech%get_asReal ('r_cutback', defaultVal=0.25_pREAL)
num%stepIncreaseCryst = num_mech%get_asReal ('r_increase_step', defaultVal=1.5_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%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%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%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%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%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%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%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%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%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%iJacoLiresiduum = num_mech_eigen%get_asInt ('f_update_jacobi_Li', defaultVal=num%iJacoLpresiduum)
extmsg = '' extmsg = ''
if (num%stepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' sub_step_min' if (num%stepMinCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_min'
if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step' if (num%stepSizeCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback'
if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_increase_step' if (num%stepIncreaseCryst <= 0.0_pREAL) extmsg = trim(extmsg)//' r_increase'
if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step_Lp' if (num%stepSizeLp <= 0.0_pREAL) extmsg = trim(extmsg)//' r_linesearch_Lp'
if (num%stepSizeLi <= 0.0_pREAL) extmsg = trim(extmsg)//' r_cutback_step_Li' 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%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%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%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 !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction !> 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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en 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 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 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 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 if (error) return ! error
A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp 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 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> 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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: subState0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: &
ph, & 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 if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState 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 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,2) = merge(dotState_last(1:sizeDotState,1),0.0_pREAL, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState 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 if (broken) exit iteration
dotState = plastic_dotState(Delta_t,ph,en) 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 = dotState * zeta &
+ dotState_last(1:sizeDotState,1) * (1.0_pREAL - zeta) + dotState_last(1:sizeDotState,1) * (1.0_pREAL - zeta)
r = plasticState(ph)%state(1:sizeDotState,en) & r = plasticState(ph)%state(1:sizeDotState,en) &
- subState0 & - state0 &
- dotState * Delta_t - dotState * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r 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 !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: subState0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
@ -694,15 +694,15 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
#ifndef __INTEL_LLVM_COMPILER #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 #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 #endif
broken = plastic_deltaState(ph,en) broken = plastic_deltaState(ph,en)
if (broken) return if (broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
end function integrateStateEuler end function integrateStateEuler
@ -710,10 +710,10 @@ end function integrateStateEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: subState0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
@ -737,15 +737,15 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en
r = - dotState * 0.5_pREAL * Delta_t r = - dotState * 0.5_pREAL * Delta_t
#ifndef __INTEL_LLVM_COMPILER #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 #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 #endif
broken = plastic_deltaState(ph,en) broken = plastic_deltaState(ph,en)
if (broken) return if (broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
if (broken) return if (broken) return
dotState = plastic_dotState(Delta_t,ph,en) 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 !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: subState0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: broken 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) 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 end function integrateStateRK4
@ -789,10 +789,10 @@ end function integrateStateRK4
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: subState0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: broken 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] 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 end function integrateStateRKCK45
@ -825,10 +825,10 @@ end function integrateStateRKCK45
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method !! 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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: subState0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
real(pREAL), dimension(:,:), intent(in) :: A real(pREAL), dimension(:,:), intent(in) :: A
real(pREAL), dimension(:), intent(in) :: B, C 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 end do
#ifndef __INTEL_LLVM_COMPILER #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 #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 #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 if (broken) exit
dotState = plastic_dotState(Delta_t*C(stage), ph,en) 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 plastic_RKdotState(1:sizeDotState,size(B)) = dotState
dotState = matmul(plastic_RKdotState,B) dotState = matmul(plastic_RKdotState,B)
#ifndef __INTEL_LLVM_COMPILER #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 #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 #endif
if (present(DB)) & 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) broken = plastic_deltaState(ph,en)
if (broken) return if (broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
end function integrateStateRK end function integrateStateRK
@ -1019,24 +1019,24 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
logical :: todo logical :: todo
real(pREAL) :: stepFrac,step real(pREAL) :: stepFrac,step
real(pREAL), dimension(3,3) :: & real(pREAL), dimension(3,3) :: &
subFp0, & Fp0, &
subFi0, & Fi0, &
subLp0, & Lp0, &
subLi0, & Li0, &
subF0, & F0, &
subF F
real(pREAL), dimension(plasticState(material_ID_phase(co,ce))%sizeState) :: subState0 real(pREAL), dimension(plasticState(material_ID_phase(co,ce))%sizeState) :: state0
ph = material_ID_phase(co,ce) ph = material_ID_phase(co,ce)
en = material_entry_phase(co,ce) en = material_entry_phase(co,ce)
subState0 = plasticState(ph)%state0(:,en) state0 = plasticState(ph)%state0(:,en)
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) Li0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) Lp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) Fp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) Fi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) F0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
stepFrac = 0.0_pREAL stepFrac = 0.0_pREAL
todo = .true. todo = .true.
step = 1.0_pREAL/num%stepSizeCryst 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? todo = step > 0.0_pREAL ! still time left to integrate on?
if (todo) then if (todo) then
subF0 = subF F0 = F
subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,en) Lp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,en)
subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,en) Li0 = phase_mechanical_Li(ph)%data(1:3,1:3,en)
subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) Fp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) Fi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%state(:,en) state0 = plasticState(ph)%state(:,en)
end if end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore) ! cut back (reduced time and restore)
else else
step = num%stepSizeCryst * step step = num%stepSizeCryst * step
phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp0
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 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) 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 if (step < 1.0_pREAL) then ! actual (not initial) cutback
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Lp(ph)%data(1:3,1:3,en) = Lp0
phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = Li0
end if end if
plasticState(ph)%state(:,en) = subState0 plasticState(ph)%state(:,en) = state0
todo = step > num%stepMinCryst ! still on track or already done (beyond repair) todo = step > num%stepMinCryst ! still on track or already done (beyond repair)
end if end if
@ -1079,9 +1079,9 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
! prepare for integration ! prepare for integration
if (todo) then if (todo) then
sizeDotState = plasticState(ph)%sizeDotState 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)) + 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 if
end do cutbackLooping end do cutbackLooping