From 5343aecd6f5e33ad993329881dec61c7ea5ce483 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jan 2024 18:02:28 +0100 Subject: [PATCH] non-binary status allows to transfer more information --- src/constants.f90 | 9 ++++++++- src/phase_mechanical.f90 | 33 ++++++++++++++++---------------- src/phase_mechanical_plastic.f90 | 10 +++++----- 3 files changed, 30 insertions(+), 22 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 1402154c7..450c51f0f 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -1,6 +1,6 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven -!> @brief physical constants +!> @brief Constants. !-------------------------------------------------------------------------------------------------- module constants use prec @@ -20,4 +20,11 @@ module constants character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' character(len=len(LOWER)), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + enum, bind(c); enumerator :: & + STATUS_OK, & + STATUS_FAILED_PHASE_STATE, & + STATUS_FAILED_PHASE_DELTASTATE, & + STATUS_FAILED_PHASE_STRESS + end enum + end module constants diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index beb968875..9f4c02863 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -71,12 +71,12 @@ submodule(phase) mechanical dotState end function plastic_dotState - module function plastic_deltaState(ph, en) result(broken) + module function plastic_deltaState(ph, en) result(status) integer, intent(in) :: & ph, & en - logical :: & - broken + integer(kind(STATUS_OK)) :: & + status end function plastic_deltaState module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & @@ -377,11 +377,12 @@ 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,Fp0,Fi0,Delta_t,ph,en) result(broken) +function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(status) real(pREAL), dimension(3,3), intent(in) :: F,Fp0,Fi0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new @@ -430,10 +431,10 @@ function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken) p, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update - logical :: error,broken + logical :: error - broken = .true. + status = STATUS_FAILED_PHASE_STRESS call plastic_dependentState(ph,en) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess @@ -573,7 +574,7 @@ function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken) phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp_new / math_det33(Fp_new)**(1.0_pREAL/3.0_pREAL) ! regularize phase_mechanical_Fi(ph)%data(1:3,1:3,en) = Fi_new phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),invFi_new) - broken = .false. + status = STATUS_OK end function integrateStress @@ -618,7 +619,7 @@ 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,1) = dotState - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) if (broken) exit iteration dotState = plastic_dotState(Delta_t,ph,en) @@ -633,7 +634,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 if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then - broken = plastic_deltaState(ph,en) + broken = STATUS_OK /= plastic_deltaState(ph,en) exit iteration end if @@ -695,10 +696,10 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) sizeDotState = plasticState(ph)%sizeDotState plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = plastic_deltaState(ph,en) + broken = STATUS_OK /= plastic_deltaState(ph,en) if (broken) return - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) end function integrateStateEuler @@ -734,10 +735,10 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( r = - dotState * 0.5_pREAL * Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = plastic_deltaState(ph,en) + broken = STATUS_OK /= plastic_deltaState(ph,en) if (broken) return - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) if (broken) return dotState = plastic_dotState(Delta_t,ph,en) @@ -858,7 +859,7 @@ 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 - broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) + broken = STATUS_OK /= 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) @@ -879,10 +880,10 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br if (broken) return - broken = plastic_deltaState(ph,en) + broken = STATUS_OK /= plastic_deltaState(ph,en) if (broken) return - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) end function integrateStateRK diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 2d0324742..616d8bcf9 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -375,12 +375,12 @@ end subroutine plastic_dependentState !> @brief for constitutive models that have an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -module function plastic_deltaState(ph, en) result(broken) +module function plastic_deltaState(ph, en) result(status) integer, intent(in) :: & ph, & en - logical :: broken + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(3,3) :: & Mp @@ -388,7 +388,7 @@ module function plastic_deltaState(ph, en) result(broken) mySize - broken = .false. + status = STATUS_OK select case (mechanical_plasticity_type(ph)) case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING) @@ -407,8 +407,8 @@ module function plastic_deltaState(ph, en) result(broken) end select plasticType - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) - if (.not. broken) then + if (any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))) status = STATUS_FAILED_PHASE_DELTASTATE + if (status == STATUS_OK) then mySize = plasticState(ph)%sizeDeltaState plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & + plasticState(ph)%deltaState(1:mySize,en)