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

sophisticasted error handling

This commit is contained in:
Martin Diehl 2024-01-07 08:44:11 +01:00
parent 5343aecd6f
commit 0ad1a14e3f
No known key found for this signature in database
GPG Key ID: 1FD50837275A0A9B
1 changed files with 48 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

91
src/phase_mechanical.f90
View File

@ -583,7 +583,7 @@ 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,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status)
real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in),dimension(:) :: state0
@ -591,8 +591,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: & integer(kind(STATUS_OK)) :: status
broken
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
@ -606,7 +605,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
dotState_last dotState_last
broken = .true. status = STATUS_FAILED_PHASE_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return if (any(IEEE_is_NaN(dotState))) return
@ -619,8 +618,8 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
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 = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
if (broken) exit iteration if (status /= STATUS_OK) exit iteration
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) exit iteration if (any(IEEE_is_NaN(dotState))) exit iteration
@ -634,7 +633,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
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
if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
broken = STATUS_OK /= plastic_deltaState(ph,en) status = plastic_deltaState(ph,en)
exit iteration exit iteration
end if end if
@ -671,7 +670,7 @@ 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,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status)
real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in),dimension(:) :: state0
@ -679,8 +678,8 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en !< grain index in grain loop en !< grain index in grain loop
logical :: & integer(kind(STATUS_OK)) :: &
broken status
real(pREAL), dimension(plasticState(ph)%sizeDotState) :: & real(pREAL), dimension(plasticState(ph)%sizeDotState) :: &
dotState dotState
@ -688,7 +687,7 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
sizeDotState sizeDotState
broken = .true. status = STATUS_FAILED_PHASE_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return if (any(IEEE_is_NaN(dotState))) return
@ -696,10 +695,10 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
broken = STATUS_OK /= plastic_deltaState(ph,en) status = plastic_deltaState(ph,en)
if (broken) return if (status /= STATUS_OK) return
broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
end function integrateStateEuler end function integrateStateEuler
@ -707,7 +706,7 @@ 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,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status)
real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in),dimension(:) :: state0
@ -715,8 +714,8 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: & integer(kind(STATUS_OK)) :: &
broken status
integer :: & integer :: &
sizeDotState sizeDotState
@ -725,7 +724,7 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(
dotState dotState
broken = .true. status = STATUS_FAILED_PHASE_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return if (any(IEEE_is_NaN(dotState))) return
@ -735,18 +734,21 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(
r = - dotState * 0.5_pREAL * Delta_t r = - dotState * 0.5_pREAL * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
broken = STATUS_OK /= plastic_deltaState(ph,en) status = plastic_deltaState(ph,en)
if (broken) return if (status /= STATUS_OK) return
broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
if (broken) return if (status /= STATUS_OK) return
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
! ToDo: MD: need to set status to failed
if (any(IEEE_is_NaN(dotState))) return if (any(IEEE_is_NaN(dotState))) return
broken = .not. converged(r + 0.5_pREAL * dotState * Delta_t, & status = merge(STATUS_OK, &
STATUS_FAILED_PHASE_STATE, &
converged(r + 0.5_pREAL * dotState * Delta_t, &
plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState)) plasticState(ph)%atol(1:sizeDotState)))
end function integrateStateAdaptiveEuler end function integrateStateAdaptiveEuler
@ -754,13 +756,13 @@ 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,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) function integrateStateRK4(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status)
real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 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 integer(kind(STATUS_OK)) :: status
real(pREAL), dimension(3,3), parameter :: & real(pREAL), dimension(3,3), parameter :: &
A = reshape([& A = reshape([&
@ -774,7 +776,7 @@ function integrateStateRK4(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
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,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C) status = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C)
end function integrateStateRK4 end function integrateStateRK4
@ -782,13 +784,13 @@ 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,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) function integrateStateRKCK45(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status)
real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 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 integer(kind(STATUS_OK)) :: status
real(pREAL), dimension(5,5), parameter :: & real(pREAL), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
@ -809,7 +811,7 @@ function integrateStateRKCK45(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
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,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) status = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB)
end function integrateStateRKCK45 end function integrateStateRKCK45
@ -818,7 +820,7 @@ 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,Fp0,Fi0,state0,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(status)
real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in),dimension(:) :: state0
@ -829,7 +831,7 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: broken integer(kind(STATUS_OK)) :: status
integer :: & integer :: &
stage, & ! stage index in integration stage loop stage, & ! stage index in integration stage loop
@ -841,7 +843,7 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br
plastic_RKdotState plastic_RKdotState
broken = .true. status = STATUS_FAILED_PHASE_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return if (any(IEEE_is_NaN(dotState))) return
@ -859,14 +861,15 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br
plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
broken = STATUS_OK /= integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) status = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en)
if (broken) exit if (status /= STATUS_OK) return
dotState = plastic_dotState(Delta_t*C(stage), ph,en) dotState = plastic_dotState(Delta_t*C(stage), ph,en)
! ToDo: MD: need to set status to failed
if (any(IEEE_is_NaN(dotState))) exit if (any(IEEE_is_NaN(dotState))) exit
end do end do
if (broken) return if (status /= STATUS_OK) return
plastic_RKdotState(1:sizeDotState,size(B)) = dotState plastic_RKdotState(1:sizeDotState,size(B)) = dotState
@ -874,16 +877,18 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br
plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
if (present(DB)) & if (present(DB)) &
broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & status = merge(STATUS_OK, &
STATUS_FAILED_PHASE_STATE, &
converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState)) plasticState(ph)%atol(1:sizeDotState)))
if (broken) return if (status /= STATUS_OK) return
broken = STATUS_OK /= plastic_deltaState(ph,en) status = plastic_deltaState(ph,en)
if (broken) return if (status /= STATUS_OK) return
broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
end function integrateStateRK end function integrateStateRK
@ -1062,7 +1067,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
F = F0 & 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(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en) converged_ = STATUS_OK == integrateState(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en)
end if end if
end do cutbackLooping end do cutbackLooping