From 5343aecd6f5e33ad993329881dec61c7ea5ce483 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jan 2024 18:02:28 +0100 Subject: [PATCH 01/15] 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) From 0ad1a14e3f394214f32f94bca776067fc147015b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 7 Jan 2024 08:44:11 +0100 Subject: [PATCH 02/15] sophisticasted error handling --- src/phase_mechanical.f90 | 91 +++++++++++++++++++++------------------- 1 file changed, 48 insertions(+), 43 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 9f4c02863..b87c12b27 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -583,7 +583,7 @@ 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,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(:) :: state0 @@ -591,8 +591,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) integer, intent(in) :: & ph, & en - logical :: & - broken + integer(kind(STATUS_OK)) :: status integer :: & 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 - broken = .true. + status = STATUS_FAILED_PHASE_STATE dotState = plastic_dotState(Delta_t,ph,en) 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,1) = dotState - broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) - if (broken) exit iteration + status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + if (status /= STATUS_OK) exit iteration dotState = plastic_dotState(Delta_t,ph,en) 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 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 end if @@ -671,7 +670,7 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @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(:) :: state0 @@ -679,8 +678,8 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) integer, intent(in) :: & ph, & en !< grain index in grain loop - logical :: & - broken + integer(kind(STATUS_OK)) :: & + status real(pREAL), dimension(plasticState(ph)%sizeDotState) :: & dotState @@ -688,7 +687,7 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) sizeDotState - broken = .true. + status = STATUS_FAILED_PHASE_STATE dotState = plastic_dotState(Delta_t,ph,en) 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 plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = STATUS_OK /= plastic_deltaState(ph,en) - if (broken) return + status = plastic_deltaState(ph,en) + 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 @@ -707,7 +706,7 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @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(:) :: state0 @@ -715,8 +714,8 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( integer, intent(in) :: & ph, & en - logical :: & - broken + integer(kind(STATUS_OK)) :: & + status integer :: & sizeDotState @@ -725,7 +724,7 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( dotState - broken = .true. + status = STATUS_FAILED_PHASE_STATE dotState = plastic_dotState(Delta_t,ph,en) 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 plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = STATUS_OK /= plastic_deltaState(ph,en) - if (broken) return + status = plastic_deltaState(ph,en) + if (status /= STATUS_OK) return - broken = STATUS_OK /= integrateStress(F,Fp0,Fi0,Delta_t,ph,en) - if (broken) return + status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + if (status /= STATUS_OK) return dotState = plastic_dotState(Delta_t,ph,en) + ! ToDo: MD: need to set status to failed 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)%atol(1:sizeDotState)) + plasticState(ph)%atol(1:sizeDotState))) end function integrateStateAdaptiveEuler @@ -754,13 +756,13 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @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(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: broken + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(3,3), parameter :: & 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) - 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 @@ -782,13 +784,13 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @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(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: broken + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(5,5), parameter :: & 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] - 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 @@ -818,7 +820,7 @@ 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,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(:) :: 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) :: & ph, & en - logical :: broken + integer(kind(STATUS_OK)) :: status integer :: & 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 - broken = .true. + status = STATUS_FAILED_PHASE_STATE dotState = plastic_dotState(Delta_t,ph,en) 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 - broken = STATUS_OK /= integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) - if (broken) exit + status = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) + if (status /= STATUS_OK) return dotState = plastic_dotState(Delta_t*C(stage), ph,en) + ! ToDo: MD: need to set status to failed if (any(IEEE_is_NaN(dotState))) exit end do - if (broken) return + if (status /= STATUS_OK) return 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 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)%atol(1:sizeDotState)) + plasticState(ph)%atol(1:sizeDotState))) - if (broken) return + if (status /= STATUS_OK) return - broken = STATUS_OK /= plastic_deltaState(ph,en) - if (broken) return + status = plastic_deltaState(ph,en) + 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 @@ -1062,7 +1067,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) sizeDotState = plasticState(ph)%sizeDotState 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(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 do cutbackLooping From 70e720c032bd1c2ba0da3acfc71f08c968e88264 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 9 Jan 2024 07:26:28 +0100 Subject: [PATCH 03/15] better to understand and more flexible --- src/constants.f90 | 4 +++- src/phase_damage.f90 | 40 +++++++++++++++++++--------------------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 450c51f0f..d775d92b9 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -24,7 +24,9 @@ module constants STATUS_OK, & STATUS_FAILED_PHASE_STATE, & STATUS_FAILED_PHASE_DELTASTATE, & - STATUS_FAILED_PHASE_STRESS + STATUS_FAILED_PHASE_STRESS, & + STATUS_FAILED_DAMAGE_STATE, & + STATUS_FAILED_DAMAGE_DELTASTATE end enum end module constants diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 2c605c963..cb3cb6efc 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -135,7 +135,7 @@ module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - converged_ = .not. integrateDamageState(Delta_t,ph,en) + converged_ = STATUS_OK == integrateDamageState(Delta_t,ph,en) end function phase_damage_constitutive @@ -214,13 +214,13 @@ end function phase_f_phi !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateDamageState(Delta_t,ph,en) result(broken) +function integrateDamageState(Delta_t,ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & ph, & en - logical :: broken + integer(kind(STATUS_OK)) :: status integer :: & NiterationState, & !< number of iterations in state loop @@ -235,13 +235,13 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (damageState(ph)%sizeState == 0) then - broken = .false. + status = STATUS_OK return end if converged_ = .true. - broken = phase_damage_collectDotState(ph,en) - if (broken) return + status = phase_damage_collectDotState(ph,en) + if (status /= STATUS_OK) return size_so = damageState(ph)%sizeDotState damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & @@ -253,8 +253,8 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) - broken = phase_damage_collectDotState(ph,en) - if (broken) exit iteration + status = phase_damage_collectDotState(ph,en) + if (status /= STATUS_OK) exit iteration zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) @@ -270,14 +270,13 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (converged_) then - broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) + status = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) exit iteration end if end do iteration - broken = broken .or. .not. converged_ - + if (.not. converged_) status = STATUS_FAILED_DAMAGE_STATE contains !-------------------------------------------------------------------------------------------------- @@ -361,15 +360,15 @@ end subroutine damage_result !-------------------------------------------------------------------------------------------------- !> @brief Constitutive equation for calculating the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(ph,en) result(broken) +function phase_damage_collectDotState(ph,en) result(status) integer, intent(in) :: & ph, & en !< counter in source loop - logical :: broken + integer(kind(STATUS_OK)) :: status - broken = .false. + status = STATUS_OK if (damageState(ph)%sizeState > 0) then @@ -380,7 +379,7 @@ function phase_damage_collectDotState(ph,en) result(broken) end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) + if (any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))) status = STATUS_FAILED_DAMAGE_STATE end if @@ -419,7 +418,7 @@ end function phase_K_phi !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function phase_damage_deltaState(Fe, ph, en) result(broken) +function phase_damage_deltaState(Fe, ph, en) result(status) integer, intent(in) :: & ph, & @@ -430,11 +429,10 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) integer :: & myOffset, & mySize - logical :: & - broken + integer(kind(STATUS_OK)) :: status - broken = .false. + status = STATUS_OK if (damageState(ph)%sizeState == 0) return @@ -442,8 +440,8 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) case (DAMAGE_ISOBRITTLE) sourceType call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) - broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) - if (.not. broken) then + if (any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))) status = STATUS_FAILED_DAMAGE_DELTASTATE + if (status == STATUS_OK) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & From 845d3eed3319f1c807f21941527e60a98dd87872 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 9 Jan 2024 12:03:25 +0100 Subject: [PATCH 04/15] harmonize status reporting --- src/Marc/materialpoint_Marc.f90 | 17 +++++++++-------- src/constants.f90 | 4 +++- src/grid/grid_mech_utilities.f90 | 6 ++++-- src/homogenization.f90 | 18 +++++++++--------- src/mesh/FEM_utilities.f90 | 5 ++++- 5 files changed, 29 insertions(+), 21 deletions(-) diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index 18ba90a61..f0040cba2 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -20,6 +20,7 @@ module materialpoint_Marc use material use phase use homogenization + use constants use discretization use discretization_Marc @@ -27,18 +28,18 @@ module materialpoint_Marc implicit none(type,external) private - real(pREAL), dimension (:,:,:), allocatable, private :: & + real(pREAL), dimension (:,:,:), allocatable :: & materialpoint_cs !< Cauchy stress - real(pREAL), dimension (:,:,:,:), allocatable, private :: & + real(pREAL), dimension (:,:,:,:), allocatable :: & materialpoint_dcsdE, & !< Cauchy stress tangent materialpoint_F !< deformation gradient - real(pREAL), dimension (:,:,:,:), allocatable, private :: & + real(pREAL), dimension (:,:,:,:), allocatable :: & materialpoint_dcsdE_knownGood !< known good tangent integer, public :: & cycleCounter = 0 !< needs description - logical, public :: & - broken = .false. !< needs description + integer(kind(STATUS_OK)) :: & + status integer, parameter, public :: & materialpoint_CALCRESULTS = 2**0, & @@ -150,17 +151,17 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, if (iand(mode, materialpoint_CALCRESULTS) /= 0) then - validCalculation: if (broken) then + validCalculation: if (status /= STATUS_OK) then call random_number(rnd) if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL materialpoint_cs(1:6,ip,elCP) = ODD_STRESS * rnd materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - call homogenization_mechanical_response(broken, & + call homogenization_mechanical_response(status, & dt,(elCP-1)*discretization_nIPs + ip, (elCP-1)*discretization_nIPs + ip) - terminalIllness: if (broken) then + terminalIllness: if (status /= STATUS_OK) then call random_number(rnd) if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL diff --git a/src/constants.f90 b/src/constants.f90 index d775d92b9..7d33d38d8 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -26,7 +26,9 @@ module constants STATUS_FAILED_PHASE_DELTASTATE, & STATUS_FAILED_PHASE_STRESS, & STATUS_FAILED_DAMAGE_STATE, & - STATUS_FAILED_DAMAGE_DELTASTATE + STATUS_FAILED_DAMAGE_DELTASTATE, & + STATUS_FAILED_DAMAGE, & + STATUS_FAILED_MECHANICAL end enum end module constants diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index a767c8623..d05d4a51f 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -18,7 +18,7 @@ module grid_mech_utilities use discretization use spectral_utilities use homogenization - + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -129,6 +129,7 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min real(pREAL) :: dPdF_norm_max, dPdF_norm_min real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF + integer(kind(STATUS_OK)) :: status print'(/,1x,a)', '... evaluating constitutive response ......................................' @@ -136,7 +137,8 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(broken,Delta_t,1,product(cells(1:2))*cells3) ! calculate P field + call homogenization_mechanical_response(status,Delta_t,1,product(cells(1:2))*cells3) ! calculate P field + broken = STATUS_OK /= status P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ab2b8a8b1..553387ea0 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -210,9 +210,9 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end) +subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status real(pREAL), intent(in) :: Delta_t !< time increment integer, intent(in) :: & cell_start, cell_end @@ -224,7 +224,7 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end doneAndHappy - broken = .false. + status = STATUS_OK !$OMP PARALLEL DO PRIVATE(en,ho,co,converged,doneAndHappy) do ce = cell_start, cell_end @@ -239,7 +239,7 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end doneAndHappy = [.false.,.true.] - convergenceLooping: do while (.not. (broken .or. doneAndHappy(1))) + convergenceLooping: do while (.not. (status /= STATUS_OK .or. doneAndHappy(1))) call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) @@ -251,19 +251,19 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end end if end do convergenceLooping if (.not. converged) then - if (.not. broken) print*, ' Cell ', ce, ' failed (mechanics)' - broken = .true. + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (mechanics)' + status = STATUS_FAILED_MECHANICAL end if converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then - if (.not. broken) print*, ' Cell ', ce, ' failed (damage)' - broken = .true. + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)' + status = STATUS_FAILED_DAMAGE end if end do !$OMP END PARALLEL DO - if (broken) return + if (status /= STATUS_OK) return !$OMP PARALLEL DO PRIVATE(ho) do ce = cell_start, cell_end diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 9f1584795..0ab8f76b2 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -22,6 +22,7 @@ module FEM_utilities use discretization_mesh use homogenization use FEM_quadrature + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -128,11 +129,13 @@ subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData) real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress integer(MPI_INTEGER_KIND) :: err_MPI + integer(kind(STATUS_OK)) :: status print'(/,1x,a)', '... evaluating constitutive response ......................................' - call homogenization_mechanical_response(broken,Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field + call homogenization_mechanical_response(status,Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field + broken = STATUS_OK /= status cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt From 62bce163cd10d6ee1be65e96f62129b4a31aee8d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 9 Jan 2024 12:55:02 +0100 Subject: [PATCH 05/15] better status propagation --- src/constants.f90 | 3 ++- src/grid/grid_thermal_spectral.f90 | 5 ++++- src/homogenization.f90 | 15 +++++++-------- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 7d33d38d8..a7b0a1877 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -28,7 +28,8 @@ module constants STATUS_FAILED_DAMAGE_STATE, & STATUS_FAILED_DAMAGE_DELTASTATE, & STATUS_FAILED_DAMAGE, & - STATUS_FAILED_MECHANICAL + STATUS_FAILED_MECHANICAL, & + STATUS_PHASE_THERMAL end enum end module constants diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 35f2b47b4..f86280b2c 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -25,6 +25,7 @@ module grid_thermal_spectral use homogenization use YAML_types use config + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -322,9 +323,11 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc) integer :: i, j, k, ce real(pREAL), dimension(3,cells(1),cells(2),cells3) :: vectorField + integer(kind(STATUS_OK)) :: status - call homogenization_thermal_response(broken,Delta_t_,1,product(cells(1:2))*cells3) + call homogenization_thermal_response(status,Delta_t_,1,product(cells(1:2))*cells3) + broken = STATUS_OK /= status associate(T => x_scal) vectorField = utilities_ScalarGradient(T) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 553387ea0..d02162650 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -257,7 +257,7 @@ subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then - if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)' + if (STATUS_OK == status) print*, ' Cell ', ce, ' failed (damage)' status = STATUS_FAILED_DAMAGE end if end do @@ -281,10 +281,10 @@ end subroutine homogenization_mechanical_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_thermal_response(broken, & +subroutine homogenization_thermal_response(status, & Delta_t,cell_start,cell_end) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status real(pREAL), intent(in) :: Delta_t !< time increment integer, intent(in) :: & cell_start, cell_end @@ -293,20 +293,19 @@ subroutine homogenization_thermal_response(broken, & co, ce, ho - broken = .false. + status = STATUS_OK !$OMP PARALLEL DO PRIVATE(ho) do ce = cell_start, cell_end - if (broken) continue + if (STATUS_OK /= status) continue ho = material_ID_homogenization(ce) do co = 1, homogenization_Nconstituents(ho) if (.not. phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce))) then - if (.not. broken) print*, ' Cell ', ce, ' failed (thermal)' - broken = .true. + if (STATUS_OK == status) print*, ' Cell ', ce, ' failed (thermal)' + status = STATUS_PHASE_THERMAL end if end do end do !$OMP END PARALLEL DO - broken = broken end subroutine homogenization_thermal_response From 5e7565f99f135e871fe7a58cc40d80284e2a3bdc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 9 Jan 2024 14:25:02 +0100 Subject: [PATCH 06/15] status needs to be passed --- src/Marc/DAMASK_Marc.f90 | 2 +- src/parallelization.f90 | 21 +++++++++++++-------- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index 539bc49d2..6f34f5354 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -16,8 +16,8 @@ #endif #include "../prec.f90" -#include "../parallelization.f90" #include "../constants.f90" +#include "../parallelization.f90" #include "../misc.f90" #include "../IO.f90" #include "../YAML_types.f90" diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 20ab75336..775a68f6c 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -7,6 +7,8 @@ module parallelization OUTPUT_UNIT, & ERROR_UNIT + use constants + #ifdef PETSC #include use PETScSys @@ -63,6 +65,7 @@ subroutine parallelization_init() !$ integer(pI32) :: OMP_NUM_THREADS !$ character(len=6) NumThreadsString PetscErrorCode :: err_PETSc + integer(kind(STATUS_OK)) :: status #ifdef _OPENMP @@ -116,28 +119,30 @@ subroutine parallelization_init() #endif call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) & - error stop 'Could not determine worldsize' + call parallelization_chkerr(err_MPI) if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) & - error stop 'Could not determine size of MPI_INTEGER' + call parallelization_chkerr(err_MPI) if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0),MPI_INTEGER_KIND)) & error stop 'Mismatch between MPI_INTEGER and DAMASK default integer' call MPI_Type_size(MPI_INTEGER8,typeSize,err_MPI) - if (err_MPI /= 0) & - error stop 'Could not determine size of MPI_INTEGER8' + call parallelization_chkerr(err_MPI) if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0_pI64),MPI_INTEGER_KIND)) & error stop 'Mismatch between MPI_INTEGER8 and DAMASK pI64' call MPI_Type_size(MPI_DOUBLE,typeSize,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) & - error stop 'Could not determine size of MPI_DOUBLE' + call parallelization_chkerr(err_MPI) if (typeSize*8_MPI_INTEGER_KIND /= int(storage_size(0.0_pREAL),MPI_INTEGER_KIND)) & error stop 'Mismatch between MPI_DOUBLE and DAMASK pREAL' + call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) + call parallelization_chkerr(err_MPI) + if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(status),MPI_INTEGER_KIND)) & + error stop 'Mismatch between MPI_INTEGER and DAMASK status' + + !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env) !$ if (got_env /= 0) then !$ print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default' From 48a38f3cf516ed5596dfb7052d461cb87b4f0948 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 9 Jan 2024 14:43:13 +0100 Subject: [PATCH 07/15] harmonize status reporting --- src/grid/grid_mech_FEM.f90 | 14 +++++++------- src/grid/grid_mech_spectral_basic.f90 | 13 +++++++------ src/grid/grid_mech_spectral_polarization.f90 | 15 ++++++++------- src/grid/grid_mech_utilities.f90 | 6 ++---- 4 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 140b22817..258a1b511 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -28,7 +28,7 @@ module grid_mechanical_FEM use homogenization use discretization use discretization_grid - + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -84,7 +84,7 @@ module grid_mechanical_FEM err_BC !< deviation from stress BC integer :: totalIter = 0 !< total iteration in current increment - logical :: broken + integer(kind(STATUS_OK)) :: status public :: & grid_mechanical_FEM_init, & @@ -271,7 +271,7 @@ subroutine grid_mechanical_FEM_init(num_grid) end if restartRead call utilities_updateCoords(F) - call utilities_constitutiveResponse(broken,P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_constitutiveResponse(status,P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F 0.0_pREAL) ! time increment call DMDAVecRestoreArrayF90(DM_mech,u_PETSc,u,err_PETSc) @@ -325,7 +325,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution) solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = broken + solution%termIll = STATUS_OK /= status P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_FEM_solution @@ -492,7 +492,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,e BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) & - .or. broken) then + .or. STATUS_OK /= status) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -570,10 +570,10 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(broken,P_current,& + call utilities_constitutiveResponse(status,P_current,& P_av,C_volAvg,devNull, & F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index cc92d4e13..2653be0d6 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -26,6 +26,7 @@ module grid_mechanical_spectral_basic use grid_mech_utilities use homogenization use discretization_grid + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -84,7 +85,7 @@ module grid_mechanical_spectral_basic err_div !< RMS of div of P integer :: totalIter = 0 !< total iteration in current increment - logical :: broken + integer(kind(STATUS_OK)) :: status public :: & grid_mechanical_spectral_basic_init, & @@ -228,7 +229,7 @@ subroutine grid_mechanical_spectral_basic_init(num_grid) end if restartRead call utilities_updateCoords(reshape(F,shape(F_lastInc))) - call utilities_constitutiveResponse(broken,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_constitutiveResponse(status,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pREAL) ! time increment call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc) ! deassociate pointer @@ -287,7 +288,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution) solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = broken + solution%termIll = STATUS_OK /= status P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_basic_solution @@ -452,7 +453,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) & - .or. broken) then + .or. STATUS_OK /= status) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -514,10 +515,10 @@ subroutine formResidual(residual_subdomain, F, & end if newIteration associate (P => r) - call utilities_constitutiveResponse(broken,P, & + call utilities_constitutiveResponse(status,P, & P_av,C_volAvg,C_minMaxAvg, & F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) err_div = utilities_divergenceRMS(P) end associate diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 4328e591f..93e3cfc4a 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -27,6 +27,7 @@ module grid_mechanical_spectral_polarization use config use homogenization use discretization_grid + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -96,7 +97,7 @@ module grid_mechanical_spectral_polarization err_div !< RMS of div of P integer :: totalIter = 0 !< total iteration in current increment - logical :: broken + integer(kind(STATUS_OK)) :: status public :: & grid_mechanical_spectral_polarization_init, & @@ -257,7 +258,7 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid) end if restartRead call utilities_updateCoords(reshape(F,shape(F_lastInc))) - call utilities_constitutiveResponse(broken,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_constitutiveResponse(status,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pREAL) ! time increment call DMDAVecRestoreArrayF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc) ! deassociate pointer @@ -322,7 +323,7 @@ function grid_mechanical_spectral_polarization_solution(incInfoIn) result(soluti solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = broken + solution%termIll = STATUS_OK /= status P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_polarization_solution @@ -516,7 +517,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pREAL)) & - .or. broken) then + .or. STATUS_OK /= status) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -604,14 +605,14 @@ subroutine formResidual(residual_subdomain, FandF_tau, & err_curl = utilities_curlRMS(F) #ifdef __GFORTRAN__ - call utilities_constitutiveResponse(broken,r_F, & + call utilities_constitutiveResponse(status,r_F, & #else associate (P => r_F) - call utilities_constitutiveResponse(broken, P, & + call utilities_constitutiveResponse(status, P, & #endif P_av,C_volAvg,C_minMaxAvg, & F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) #ifdef __GFORTRAN__ err_div = utilities_divergenceRMS(r_F) #else diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index d05d4a51f..a4080efc2 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -113,10 +113,10 @@ end function utilities_maskedCompliance !-------------------------------------------------------------------------------------------------- !> @brief Calculate constitutive response. !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& +subroutine utilities_constitutiveResponse(status, P,P_av,C_volAvg,C_minmaxAvg,& F,Delta_t,rotation_BC) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status real(pREAL), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pREAL), intent(out), dimension(3,3) :: P_av !< average PK stress real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress @@ -129,7 +129,6 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min real(pREAL) :: dPdF_norm_max, dPdF_norm_min real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF - integer(kind(STATUS_OK)) :: status print'(/,1x,a)', '... evaluating constitutive response ......................................' @@ -138,7 +137,6 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field call homogenization_mechanical_response(status,Delta_t,1,product(cells(1:2))*cells3) ! calculate P field - broken = STATUS_OK /= status P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt From 89ebad19d51a536d7f4df037e92f08cfed449e86 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 9 Jan 2024 10:00:42 -0500 Subject: [PATCH 08/15] consistent ordering of status logical tests --- src/grid/grid_mech_FEM.f90 | 4 ++-- src/grid/grid_mech_spectral_basic.f90 | 4 ++-- src/grid/grid_mech_spectral_polarization.f90 | 4 ++-- src/grid/grid_thermal_spectral.f90 | 2 +- src/homogenization.f90 | 8 ++++---- src/mesh/FEM_utilities.f90 | 2 +- src/phase_damage.f90 | 2 +- 7 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 258a1b511..d18155033 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -325,7 +325,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution) solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = STATUS_OK /= status + solution%termIll = status /= STATUS_OK P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_FEM_solution @@ -492,7 +492,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,e BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) & - .or. STATUS_OK /= status) then + .or. status /= STATUS_OK) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 2653be0d6..591299ad1 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -288,7 +288,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution) solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = STATUS_OK /= status + solution%termIll = status /= STATUS_OK P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_basic_solution @@ -453,7 +453,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) & - .or. STATUS_OK /= status) then + .or. status /= STATUS_OK) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 93e3cfc4a..4710c7e1e 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -323,7 +323,7 @@ function grid_mechanical_spectral_polarization_solution(incInfoIn) result(soluti solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = STATUS_OK /= status + solution%termIll = status /= STATUS_OK P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_polarization_solution @@ -517,7 +517,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pREAL)) & - .or. STATUS_OK /= status) then + .or. status /= STATUS_OK) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index f86280b2c..193d85c23 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -327,7 +327,7 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc) call homogenization_thermal_response(status,Delta_t_,1,product(cells(1:2))*cells3) - broken = STATUS_OK /= status + broken = status /= STATUS_OK associate(T => x_scal) vectorField = utilities_ScalarGradient(T) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index d02162650..a80517b2b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -239,7 +239,7 @@ subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end doneAndHappy = [.false.,.true.] - convergenceLooping: do while (.not. (status /= STATUS_OK .or. doneAndHappy(1))) + convergenceLooping: do while (status == STATUS_OK .and. .not. doneAndHappy(1)) call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) @@ -257,7 +257,7 @@ subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then - if (STATUS_OK == status) print*, ' Cell ', ce, ' failed (damage)' + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)' status = STATUS_FAILED_DAMAGE end if end do @@ -296,11 +296,11 @@ subroutine homogenization_thermal_response(status, & status = STATUS_OK !$OMP PARALLEL DO PRIVATE(ho) do ce = cell_start, cell_end - if (STATUS_OK /= status) continue + if (status /= STATUS_OK) continue ho = material_ID_homogenization(ce) do co = 1, homogenization_Nconstituents(ho) if (.not. phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce))) then - if (STATUS_OK == status) print*, ' Cell ', ce, ' failed (thermal)' + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (thermal)' status = STATUS_PHASE_THERMAL end if end do diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 0ab8f76b2..64509ed87 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -135,7 +135,7 @@ subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData) print'(/,1x,a)', '... evaluating constitutive response ......................................' call homogenization_mechanical_response(status,Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field - broken = STATUS_OK /= status + broken = status /= STATUS_OK cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index cb3cb6efc..2ca6506c1 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -135,7 +135,7 @@ module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - converged_ = STATUS_OK == integrateDamageState(Delta_t,ph,en) + converged_ = integrateDamageState(Delta_t,ph,en) == STATUS_OK end function phase_damage_constitutive From 0384f93d466c409bc1a623d70be1a9062dfd95cd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 10 Jan 2024 00:36:22 +0100 Subject: [PATCH 09/15] propagate non-binary status information --- src/grid/grid_thermal_spectral.f90 | 6 ++---- src/mesh/FEM_utilities.f90 | 6 ++---- src/mesh/mesh_mech_FEM.f90 | 12 +++++++----- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 193d85c23..cc4b89afb 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -57,7 +57,7 @@ module grid_thermal_spectral integer :: totalIter = 0 !< total iteration in current increment real(pREAL), dimension(3,3) :: K_ref real(pREAL) :: mu_ref, Delta_t_ - logical :: broken + integer(kind(STATUS_OK)) :: status public :: & grid_thermal_spectral_init, & @@ -208,7 +208,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc) CHKERRQ(err_PETSc) - solution%converged = reason > 0 .and. .not. broken + solution%converged = reason > 0 .and. status == STATUS_OK solution%iterationsNeeded = merge(totalIter,num%itmax,solution%converged) call SNESGetDM(SNES_thermal,DM_thermal,err_PETSc) @@ -323,11 +323,9 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc) integer :: i, j, k, ce real(pREAL), dimension(3,cells(1),cells(2),cells3) :: vectorField - integer(kind(STATUS_OK)) :: status call homogenization_thermal_response(status,Delta_t_,1,product(cells(1:2))*cells3) - broken = status /= STATUS_OK associate(T => x_scal) vectorField = utilities_ScalarGradient(T) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 64509ed87..ea4e4b87f 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -121,21 +121,19 @@ end subroutine FEM_utilities_init !-------------------------------------------------------------------------------------------------- !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData) +subroutine utilities_constitutiveResponse(status, Delta_t,P_av,forwardData) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status real(pREAL), intent(in) :: Delta_t !< loading time logical, intent(in) :: forwardData !< age results real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress integer(MPI_INTEGER_KIND) :: err_MPI - integer(kind(STATUS_OK)) :: status print'(/,1x,a)', '... evaluating constitutive response ......................................' call homogenization_mechanical_response(status,Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field - broken = status /= STATUS_OK cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 0baebe6c3..6eafe66e3 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -25,6 +25,7 @@ module mesh_mechanical_FEM use FEM_quadrature use homogenization use math + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -68,7 +69,8 @@ module mesh_mechanical_FEM character(len=pSTRLEN) :: incInfo real(pREAL), dimension(3,3) :: & P_av = 0.0_pREAL - logical :: ForwardData, broken + logical :: ForwardData + integer(kind(STATUS_OK)) :: status real(pREAL), parameter :: eps = 1.0e-18_pREAL external :: & ! ToDo: write interfaces @@ -311,7 +313,7 @@ subroutine FEM_mechanical_init(mechBC,num_mesh) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,err_PETSc) CHKERRQ(err_PETSc) end do - call utilities_constitutiveResponse(broken,0.0_pREAL,devNull,.true.) + call utilities_constitutiveResponse(status,0.0_pREAL,devNull,.true.) end subroutine FEM_mechanical_init @@ -458,8 +460,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(broken,params%Delta_t,P_av,ForwardData) - call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call utilities_constitutiveResponse(status,params%Delta_t,P_av,ForwardData) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) ForwardData = .false. @@ -747,7 +749,7 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc) CHKERRQ(err_PETSc) - if (broken) reason = SNES_DIVERGED_FUNCTION_DOMAIN + if (status /= STATUS_OK) reason = SNES_DIVERGED_FUNCTION_DOMAIN print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), & ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & From 7ca02b3f3df014bf15411b767eb1fc39037134bd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jan 2024 06:56:38 +0100 Subject: [PATCH 10/15] consistently use status code as return value --- src/constants.f90 | 4 +++- src/homogenization.f90 | 6 +++--- src/phase.f90 | 14 ++++++-------- src/phase_damage.f90 | 6 +++--- src/phase_mechanical.f90 | 12 ++++++------ src/phase_thermal.f90 | 22 +++++++++++----------- 6 files changed, 32 insertions(+), 32 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index a7b0a1877..6f3e21d0d 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -29,7 +29,9 @@ module constants STATUS_FAILED_DAMAGE_DELTASTATE, & STATUS_FAILED_DAMAGE, & STATUS_FAILED_MECHANICAL, & - STATUS_PHASE_THERMAL + STATUS_PHASE_THERMAL, & + STATUS_PHASE_THERMAL_DOTSTATE, & + STATUS_ITERATING end enum end module constants diff --git a/src/homogenization.f90 b/src/homogenization.f90 index a80517b2b..2c0134b42 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -242,7 +242,7 @@ subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end convergenceLooping: do while (status == STATUS_OK .and. .not. doneAndHappy(1)) call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) - converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) + converged = all([(phase_mechanical_constitutive(Delta_t,co,ce) == STATUS_OK,co=1,homogenization_Nconstituents(ho))]) if (converged) then doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) @@ -254,7 +254,7 @@ subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (mechanics)' status = STATUS_FAILED_MECHANICAL end if - converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce)==STATUS_OK,co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)' @@ -299,7 +299,7 @@ subroutine homogenization_thermal_response(status, & if (status /= STATUS_OK) continue ho = material_ID_homogenization(ce) do co = 1, homogenization_Nconstituents(ho) - if (.not. phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce))) then + if (phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce)) /= STATUS_OK) then if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (thermal)' status = STATUS_PHASE_THERMAL end if diff --git a/src/phase.f90 b/src/phase.f90 index 39a8679be..ed2f6fa02 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -282,24 +282,22 @@ module phase ! == cleaned:end =================================================================================== - module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) - + module function phase_thermal_constitutive(Delta_t,ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: converged_ - + integer(kind(STATUS_OK)) :: status end function phase_thermal_constitutive - module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) + module function phase_damage_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: co, ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status end function phase_damage_constitutive - module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) + module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: co, ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status end function phase_mechanical_constitutive !ToDo: Merge all the stiffness functions diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 2ca6506c1..6b5785a5c 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -120,13 +120,13 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) +module function phase_damage_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & co, & ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status integer :: & ph, en @@ -135,7 +135,7 @@ module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - converged_ = integrateDamageState(Delta_t,ph,en) == STATUS_OK + status = integrateDamageState(Delta_t,ph,en) end function phase_damage_constitutive diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index b87c12b27..6ddd43c5e 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -990,13 +990,13 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) +module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & co, & ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status real(pREAL) :: & formerStep @@ -1025,13 +1025,13 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) F0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) stepFrac = 0.0_pREAL todo = .true. - step = 1.0_pREAL/num%stepSizeCryst - converged_ = .false. ! pretend failed step of 1/stepSizeCryst + step = 1.0_pREAL/num%stepSizeCryst ! pretend failed step of 1/stepSizeCryst + status = STATUS_ITERATING todo = .true. cutbackLooping: do while (todo) - if (converged_) then + if (status == STATUS_OK) then formerStep = step stepFrac = stepFrac + step step = min(1.0_pREAL - stepFrac, num%stepIncreaseCryst * step) @@ -1067,7 +1067,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) sizeDotState = plasticState(ph)%sizeDotState F = F0 & + step * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = STATUS_OK == integrateState(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en) + status = integrateState(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en) end if end do cutbackLooping diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index f63eb5e19..2598381ff 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -187,22 +187,22 @@ end function phase_f_T !-------------------------------------------------------------------------------------------------- !> @brief tbd. !-------------------------------------------------------------------------------------------------- -function phase_thermal_collectDotState(ph,en) result(ok) +function phase_thermal_collectDotState(ph,en) result(status) integer, intent(in) :: ph, en - logical :: ok + integer(kind(STATUS_OK)) :: status integer :: i - ok = .true. + status = STATUS_OK SourceLoop: do i = 1, thermal_Nsources(ph) if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) & call source_externalheat_dotState(ph,en) - ok = ok .and. .not. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) + if (any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))) status = STATUS_PHASE_THERMAL_DOTSTATE end do SourceLoop @@ -238,14 +238,14 @@ module function phase_K_T(co,ce) result(K) end function phase_K_T -module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) +module function phase_thermal_constitutive(Delta_t,ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: converged_ + integer(kind(STATUS_OK)) :: status - converged_ = integrateThermalState(Delta_t,ph,en) + status = integrateThermalState(Delta_t,ph,en) end function phase_thermal_constitutive @@ -253,19 +253,19 @@ end function phase_thermal_constitutive !-------------------------------------------------------------------------------------------------- !> @brief Integrate state with 1st order explicit Euler method. !-------------------------------------------------------------------------------------------------- -function integrateThermalState(Delta_t, ph,en) result(converged) +function integrateThermalState(Delta_t, ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: converged + integer(kind(STATUS_OK)) :: status integer :: & so, & sizeDotState - converged = phase_thermal_collectDotState(ph,en) - if (converged) then + status = phase_thermal_collectDotState(ph,en) + if (status == STATUS_OK) then do so = 1, thermal_Nsources(ph) sizeDotState = thermalState(ph)%p(so)%sizeDotState From 9d3b4a270f6af3245704faf9f77c13886ad66c9f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Jan 2024 07:30:37 +0100 Subject: [PATCH 11/15] reasonable order --- src/constants.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 6f3e21d0d..f212ffb7d 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -22,6 +22,7 @@ module constants enum, bind(c); enumerator :: & STATUS_OK, & + STATUS_ITERATING, & STATUS_FAILED_PHASE_STATE, & STATUS_FAILED_PHASE_DELTASTATE, & STATUS_FAILED_PHASE_STRESS, & @@ -30,8 +31,7 @@ module constants STATUS_FAILED_DAMAGE, & STATUS_FAILED_MECHANICAL, & STATUS_PHASE_THERMAL, & - STATUS_PHASE_THERMAL_DOTSTATE, & - STATUS_ITERATING + STATUS_PHASE_THERMAL_DOTSTATE end enum end module constants From 1d60ce29c663a87efd9e0f9f0879a56432123ab0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Jan 2024 08:27:17 +0100 Subject: [PATCH 12/15] systematic naming of status code enums --- src/constants.f90 | 18 +++++++++--------- src/homogenization.f90 | 6 +++--- src/phase_damage.f90 | 6 +++--- src/phase_mechanical.f90 | 14 +++++++------- src/phase_mechanical_plastic.f90 | 2 +- src/phase_thermal.f90 | 2 +- 6 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index f212ffb7d..38cdcb195 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -23,15 +23,15 @@ module constants enum, bind(c); enumerator :: & STATUS_OK, & STATUS_ITERATING, & - STATUS_FAILED_PHASE_STATE, & - STATUS_FAILED_PHASE_DELTASTATE, & - STATUS_FAILED_PHASE_STRESS, & - STATUS_FAILED_DAMAGE_STATE, & - STATUS_FAILED_DAMAGE_DELTASTATE, & - STATUS_FAILED_DAMAGE, & - STATUS_FAILED_MECHANICAL, & - STATUS_PHASE_THERMAL, & - STATUS_PHASE_THERMAL_DOTSTATE + STATUS_FAIL_PHASE_MECHANICAL, & + STATUS_FAIL_PHASE_MECHANICAL_STATE, & + STATUS_FAIL_PHASE_MECHANICAL_DELTASTATE, & + STATUS_FAIL_PHASE_MECHANICAL_STRESS, & + STATUS_FAIL_PHASE_DAMAGE, & + STATUS_FAIL_PHASE_DAMAGE_STATE, & + STATUS_FAIL_PHASE_DAMAGE_DELTASTATE, & + STATUS_FAIL_PHASE_THERMAL, & + STATUS_FAIL_PHASE_THERMAL_DOTSTATE end enum end module constants diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 2c0134b42..81ca6cb93 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -252,13 +252,13 @@ subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end end do convergenceLooping if (.not. converged) then if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (mechanics)' - status = STATUS_FAILED_MECHANICAL + status = STATUS_FAIL_PHASE_MECHANICAL end if converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce)==STATUS_OK,co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)' - status = STATUS_FAILED_DAMAGE + status = STATUS_FAIL_PHASE_DAMAGE end if end do !$OMP END PARALLEL DO @@ -301,7 +301,7 @@ subroutine homogenization_thermal_response(status, & do co = 1, homogenization_Nconstituents(ho) if (phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce)) /= STATUS_OK) then if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (thermal)' - status = STATUS_PHASE_THERMAL + status = STATUS_FAIL_PHASE_THERMAL end if end do end do diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 6b5785a5c..ade9e9930 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -276,7 +276,7 @@ function integrateDamageState(Delta_t,ph,en) result(status) end do iteration - if (.not. converged_) status = STATUS_FAILED_DAMAGE_STATE + if (.not. converged_) status = STATUS_FAIL_PHASE_DAMAGE_STATE contains !-------------------------------------------------------------------------------------------------- @@ -379,7 +379,7 @@ function phase_damage_collectDotState(ph,en) result(status) end select sourceType - if (any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))) status = STATUS_FAILED_DAMAGE_STATE + if (any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))) status = STATUS_FAIL_PHASE_DAMAGE_STATE end if @@ -440,7 +440,7 @@ function phase_damage_deltaState(Fe, ph, en) result(status) case (DAMAGE_ISOBRITTLE) sourceType call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) - if (any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))) status = STATUS_FAILED_DAMAGE_DELTASTATE + if (any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))) status = STATUS_FAIL_PHASE_DAMAGE_DELTASTATE if (status == STATUS_OK) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 6ddd43c5e..ef4914db4 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -434,7 +434,7 @@ function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(status) logical :: error - status = STATUS_FAILED_PHASE_STRESS + status = STATUS_FAIL_PHASE_MECHANICAL_STRESS call plastic_dependentState(ph,en) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess @@ -605,7 +605,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) dotState_last - status = STATUS_FAILED_PHASE_STATE + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -687,7 +687,7 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) sizeDotState - status = STATUS_FAILED_PHASE_STATE + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -724,7 +724,7 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( dotState - status = STATUS_FAILED_PHASE_STATE + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -745,7 +745,7 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( if (any(IEEE_is_NaN(dotState))) return status = merge(STATUS_OK, & - STATUS_FAILED_PHASE_STATE, & + STATUS_FAIL_PHASE_MECHANICAL_STATE, & converged(r + 0.5_pREAL * dotState * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState))) @@ -843,7 +843,7 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(st plastic_RKdotState - status = STATUS_FAILED_PHASE_STATE + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -878,7 +878,7 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(st if (present(DB)) & status = merge(STATUS_OK, & - STATUS_FAILED_PHASE_STATE, & + STATUS_FAIL_PHASE_MECHANICAL_STATE, & converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState))) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 616d8bcf9..548957764 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -407,7 +407,7 @@ module function plastic_deltaState(ph, en) result(status) end select plasticType - if (any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))) status = STATUS_FAILED_PHASE_DELTASTATE + if (any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))) status = STATUS_FAIL_PHASE_MECHANICAL_DELTASTATE if (status == STATUS_OK) then mySize = plasticState(ph)%sizeDeltaState plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 2598381ff..aad7ef470 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -202,7 +202,7 @@ function phase_thermal_collectDotState(ph,en) result(status) if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) & call source_externalheat_dotState(ph,en) - if (any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))) status = STATUS_PHASE_THERMAL_DOTSTATE + if (any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))) status = STATUS_FAIL_PHASE_THERMAL_DOTSTATE end do SourceLoop From 712b22b678026864e2da2061d1f50a93e08d8573 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Jan 2024 06:57:15 +0100 Subject: [PATCH 13/15] only needed for older (broken) ifort --- src/homogenization.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 81ca6cb93..13e6910ab 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -50,9 +50,9 @@ module homogenization ! General variables for the homogenization at a material point real(pREAL), dimension(:,:,:), allocatable, public :: & homogenization_F !< def grad of IP to be reached at end of FE increment - real(pREAL), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort + real(pREAL), dimension(:,:,:), allocatable, public, protected :: & homogenization_P !< first P--K stress of IP - real(pREAL), dimension(:,:,:,:,:), allocatable, public :: & !, protected :: & + real(pREAL), dimension(:,:,:,:,:), allocatable, public, protected :: & homogenization_dPdF !< tangent of first P--K stress at IP !-------------------------------------------------------------------------------------------------- From 27a8b63a3d7e0bb0b1f8093b4c4271febccd33f0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jan 2024 23:59:02 +0100 Subject: [PATCH 14/15] safer reduction MPI_SUM can lead to invalid/undefined values for enum --- src/grid/grid_mech_FEM.f90 | 2 +- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarization.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index d18155033..f303e71de 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -573,7 +573,7 @@ subroutine formResidual(da_local,x_local, & call utilities_constitutiveResponse(status,P_current,& P_av,C_volAvg,devNull, & F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 591299ad1..d71181479 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -518,7 +518,7 @@ subroutine formResidual(residual_subdomain, F, & call utilities_constitutiveResponse(status,P, & P_av,C_volAvg,C_minMaxAvg, & F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) err_div = utilities_divergenceRMS(P) end associate diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 4710c7e1e..88261303f 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -612,7 +612,7 @@ subroutine formResidual(residual_subdomain, FandF_tau, & #endif P_av,C_volAvg,C_minMaxAvg, & F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) #ifdef __GFORTRAN__ err_div = utilities_divergenceRMS(r_F) #else diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 6eafe66e3..f409c0385 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -461,7 +461,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response call utilities_constitutiveResponse(status,params%Delta_t,P_av,ForwardData) - call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) ForwardData = .false. From db263c3af7fdf82a8ae74f2231a633f99080281a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jan 2024 23:37:54 +0100 Subject: [PATCH 15/15] consistent naming --- src/grid/grid_mech_FEM.f90 | 12 ++++++------ src/grid/grid_mech_spectral_polarization.f90 | 8 ++++---- src/grid/grid_mech_utilities.f90 | 16 ++++++++-------- src/mesh/mesh_mech_FEM.f90 | 10 +++++----- 4 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index f303e71de..00653a508 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -525,7 +525,7 @@ subroutine formResidual(da_local,x_local, & real(pREAL), pointer,dimension(:,:,:,:) :: x_scal, r real(pREAL), dimension(8,3) :: x_elem, f_elem - PetscInt :: i, ii, j, jj, k, kk, ctr, ele + PetscInt :: i, ii, j, jj, k, kk, ctr, ce PetscInt :: & PETScIter, & nfuncs @@ -587,7 +587,7 @@ subroutine formResidual(da_local,x_local, & CHKERRQ(err_PETSc) call DMDAVecGetArrayReadF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) - ele = 0 + ce = 0 r = 0.0_pREAL do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 @@ -595,11 +595,11 @@ subroutine formResidual(da_local,x_local, & ctr = ctr + 1 x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) end do; end do; end do - ele = ele + 1 + ce = ce + 1 f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + & - matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + & - homogenization_dPdF(2,2,2,2,ele) + & - homogenization_dPdF(3,3,3,3,ele))/3.0_pREAL + matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ce) + & + homogenization_dPdF(2,2,2,2,ce) + & + homogenization_dPdF(3,3,3,3,ce))/3.0_pREAL ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 88261303f..7304d3640 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -563,7 +563,7 @@ subroutine formResidual(residual_subdomain, FandF_tau, & nfuncs integer(MPI_INTEGER_KIND) :: err_MPI integer :: & - i, j, k, e + i, j, k, ce F => FandF_tau(1:3,1:3,1,1:cells(1),1:cells(2),1:cells3) @@ -618,11 +618,11 @@ subroutine formResidual(residual_subdomain, FandF_tau, & #else err_div = utilities_divergenceRMS(P) #endif - e = 0 + ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) - e = e + 1 + ce = ce + 1 r_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & + math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,ce) + C_scale), & #ifdef __GFORTRAN__ r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & #else diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index a4080efc2..be426497e 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -124,7 +124,7 @@ subroutine utilities_constitutiveResponse(status, P,P_av,C_volAvg,C_minmaxAvg,& real(pREAL), intent(in) :: Delta_t !< loading time type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame - integer :: i + integer :: ce integer(MPI_INTEGER_KIND) :: err_MPI real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min real(pREAL) :: dPdF_norm_max, dPdF_norm_min @@ -156,14 +156,14 @@ subroutine utilities_constitutiveResponse(status, P,P_av,C_volAvg,C_minmaxAvg,& dPdF_norm_max = 0.0_pREAL dPdF_min = huge(1.0_pREAL) dPdF_norm_min = huge(1.0_pREAL) - do i = 1, product(cells(1:2))*cells3 - if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then - dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) - dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) + do ce = 1, product(cells(1:2))*cells3 + if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2)) then + dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) + dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2) end if - if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then - dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i) - dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) + if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2)) then + dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) + dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2) end if end do diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index f409c0385..17475ba48 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -531,7 +531,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P real(pREAL),dimension(cellDOF,cellDOF) :: K_eA, K_eB PetscInt :: cellStart, cellEnd, cell, component, face, & - qPt, basis, comp, cidx,bcSize, m, i + qPt, basis, comp, cidx,bcSize, ce, i IS :: bcPoints @@ -583,7 +583,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P FAvg = 0.0_pREAL BMatAvg = 0.0_pREAL do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT - m = cell*nQuadrature + qPt + 1_pPETSCINT + ce = cell*nQuadrature + qPt + 1_pPETSCINT BMat = 0.0_pREAL do basis = 0_pPETSCINT, nBasis-1_pPETSCINT do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT @@ -593,7 +593,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) end do end do - MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), & + MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,ce), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT) if (num%BBarStabilization) then @@ -601,11 +601,11 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P FInv = math_inv33(F) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pREAL/real(dimPlex,pREAL)) K_eB = K_eB - & - matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), & + matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,ce),shape=[dimPlex**2,1_pPETSCINT]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & shape=[1_pPETSCINT,dimPlex**2],order=[2,1]),BMat))),MatA) MatB = MatB & - + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1_pPETSCINT,dimPlex**2]),MatA) + + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,ce),shape=[1_pPETSCINT,dimPlex**2]),MatA) FAvg = FAvg + F BMatAvg = BMatAvg + BMat else