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