From d85ad0b554301469ee1577012aa210f9359a9725 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 18 Jul 2023 04:12:57 +0200 Subject: [PATCH 01/13] improved naming --- src/Marc/DAMASK_Marc.f90 | 4 +- src/phase_thermal.f90 | 75 ++++++++++--------- ...0 => phase_thermal_source_dissipation.f90} | 14 ++-- ... => phase_thermal_source_externalheat.f90} | 37 ++++----- 4 files changed, 61 insertions(+), 69 deletions(-) rename src/{phase_thermal_dissipation.f90 => phase_thermal_source_dissipation.f90} (90%) rename src/{phase_thermal_externalheat.f90 => phase_thermal_source_externalheat.f90} (78%) diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index 032c77394..b1d9d91ba 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -176,8 +176,8 @@ end module DAMASK_interface #include "../phase_mechanical_eigen_cleavageopening.f90" #include "../phase_mechanical_eigen_thermalexpansion.f90" #include "../phase_thermal.f90" -#include "../phase_thermal_dissipation.f90" -#include "../phase_thermal_externalheat.f90" +#include "../phase_thermal_source_dissipation.f90" +#include "../phase_thermal_source_externalheat.f90" #include "../phase_damage.f90" #include "../phase_damage_isobrittle.f90" #include "../phase_damage_anisobrittle.f90" diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 449e08ab8..8f59dc874 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -21,13 +21,13 @@ submodule(phase) thermal THERMAL_EXTERNALHEAT_ID end enum - type :: tDataContainer ! ?? not very telling name. Better: "fieldQuantities" ?? + type :: tFieldQuantities real(pREAL), dimension(:), allocatable :: T, dot_T - end type tDataContainer + end type tFieldQuantities integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & thermal_source - type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(en) reads quite good + type(tFieldQuantities), dimension(:), allocatable :: current type(tThermalParameters), dimension(:), allocatable :: param @@ -36,36 +36,36 @@ submodule(phase) thermal interface - module function dissipation_init(source_length) result(mySources) + module function source_dissipation_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources - end function dissipation_init + end function source_dissipation_init - module function externalheat_init(source_length) result(mySources) + module function source_externalheat_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources - end function externalheat_init + end function source_externalheat_init - module subroutine externalheat_dotState(ph, en) + module subroutine source_externalheat_dotState(ph, en) integer, intent(in) :: & ph, & en - end subroutine externalheat_dotState + end subroutine source_externalheat_dotState - module function dissipation_f_T(ph,en) result(f_T) + module function source_dissipation_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & en real(pREAL) :: f_T - end function dissipation_f_T + end function source_dissipation_f_T - module function externalheat_f_T(ph,en) result(f_T) + module function source_externalheat_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & en real(pREAL) :: f_T - end function externalheat_f_T + end function source_externalheat_f_T end interface @@ -132,8 +132,8 @@ module subroutine thermal_init(phases) allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) if (maxval(thermal_Nsources) /= 0) then - where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID - where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID + where(source_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID + where(source_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID end if thermal_source_maxSizeDotState = 0 @@ -151,7 +151,7 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- -!< @brief Calculate thermal source. +!< @brief Calculate thermal source (forcing term). !---------------------------------------------------------------------------------------------- module function phase_f_T(ph,en) result(f) @@ -168,10 +168,10 @@ module function phase_f_T(ph,en) result(f) select case(thermal_source(so,ph)) case (THERMAL_DISSIPATION_ID) - f = f + dissipation_f_T(ph,en) + f = f + source_dissipation_f_T(ph,en) case (THERMAL_EXTERNALHEAT_ID) - f = f + externalheat_f_T(ph,en) + f = f + source_externalheat_f_T(ph,en) end select @@ -183,22 +183,22 @@ end function phase_f_T !-------------------------------------------------------------------------------------------------- !> @brief tbd. !-------------------------------------------------------------------------------------------------- -function phase_thermal_collectDotState(ph,en) result(broken) +function phase_thermal_collectDotState(ph,en) result(ok) integer, intent(in) :: ph, en - logical :: broken + logical :: ok integer :: i - broken = .false. + ok = .true. SourceLoop: do i = 1, thermal_Nsources(ph) if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & - call externalheat_dotState(ph,en) + call source_externalheat_dotState(ph,en) - broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) + ok = ok .and. .not. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) end do SourceLoop @@ -241,34 +241,35 @@ module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) logical :: converged_ - converged_ = .not. integrateThermalState(Delta_t,ph,en) + converged_ = integrateThermalState(Delta_t,ph,en) end function phase_thermal_constitutive !-------------------------------------------------------------------------------------------------- -!> @brief integrate state with 1st order explicit Euler method +!> @brief Integrate state with 1st order explicit Euler method. !-------------------------------------------------------------------------------------------------- -function integrateThermalState(Delta_t, ph,en) result(broken) +function integrateThermalState(Delta_t, ph,en) result(converged) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: & - broken + logical :: converged integer :: & so, & sizeDotState - broken = phase_thermal_collectDotState(ph,en) - if (broken) return + converged = phase_thermal_collectDotState(ph,en) + if (converged) then - do so = 1, thermal_Nsources(ph) - sizeDotState = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & - + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t - end do + do so = 1, thermal_Nsources(ph) + sizeDotState = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & + + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t + end do + + end if end function integrateThermalState @@ -318,7 +319,7 @@ end subroutine thermal_forward !---------------------------------------------------------------------------------------------- -!< @brief Get temperature (for use by non-thermal physics) +!< @brief Get temperature (for use by non-thermal physics). !---------------------------------------------------------------------------------------------- pure module function thermal_T(ph,en) result(T) @@ -332,7 +333,7 @@ end function thermal_T !---------------------------------------------------------------------------------------------- -!< @brief Get rate of temperature (for use by non-thermal physics) +!< @brief Get rate of temperature (for use by non-thermal physics). !---------------------------------------------------------------------------------------------- module function thermal_dot_T(ph,en) result(dot_T) diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_source_dissipation.f90 similarity index 90% rename from src/phase_thermal_dissipation.f90 rename to src/phase_thermal_source_dissipation.f90 index 573921670..c42e4a817 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_source_dissipation.f90 @@ -5,7 +5,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:thermal) dissipation +submodule(phase:thermal) source_dissipation type :: tParameters !< container type for internal constitutive parameters real(pREAL) :: & @@ -22,7 +22,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function dissipation_init(source_length) result(mySources) +module function source_dissipation_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -41,7 +41,7 @@ module function dissipation_init(source_length) result(mySources) mySources = thermal_active('dissipation',source_length) if (count(mySources) == 0) return - print'(/,1x,a)', '<<<+- phase:thermal:dissipation init -+>>>' + print'(/,1x,a)', '<<<+- phase:thermal:source_dissipation init -+>>>' print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) @@ -71,13 +71,13 @@ module function dissipation_init(source_length) result(mySources) end do -end function dissipation_init +end function source_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module function dissipation_f_T(ph,en) result(f_T) +module function source_dissipation_f_T(ph,en) result(f_T) integer, intent(in) :: ph, en real(pREAL) :: & @@ -91,6 +91,6 @@ module function dissipation_f_T(ph,en) result(f_T) f_T = prm%kappa*sum(abs(Mp*mechanical_L_p(ph,en))) end associate -end function dissipation_f_T +end function source_dissipation_f_T -end submodule dissipation +end submodule source_dissipation diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_source_externalheat.f90 similarity index 78% rename from src/phase_thermal_externalheat.f90 rename to src/phase_thermal_source_externalheat.f90 index cdd037592..f86880e62 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_source_externalheat.f90 @@ -4,11 +4,11 @@ !> @author Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- -submodule(phase:thermal) externalheat +submodule(phase:thermal) source_externalheat integer, dimension(:), allocatable :: & - source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? + source_ID !< which source is my current thermal dissipation mechanism? type :: tParameters !< container type for internal constitutive parameters type(tTable) :: f @@ -24,7 +24,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function externalheat_init(source_length) result(mySources) +module function source_externalheat_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -43,13 +43,13 @@ module function externalheat_init(source_length) result(mySources) mySources = thermal_active('externalheat',source_length) if (count(mySources) == 0) return - print'(/,1x,a)', '<<<+- phase:thermal:externalheat init -+>>>' + print'(/,1x,a)', '<<<+- phase:thermal:source_externalheat init -+>>>' print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) phases => config_material%get_dict('phase') allocate(param(phases%length)) - allocate(source_thermal_externalheat_offset (phases%length), source=0) + allocate(source_ID(phases%length), source=0) do ph = 1, phases%length phase => phases%get_dict(ph) @@ -58,7 +58,7 @@ module function externalheat_init(source_length) result(mySources) sources => thermal%get_list('source') do so = 1, sources%length if (mySources(so,ph)) then - source_thermal_externalheat_offset(ph) = so + source_ID(ph) = so associate(prm => param(ph)) src => sources%get_dict(so) print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so @@ -74,33 +74,29 @@ module function externalheat_init(source_length) result(mySources) end do end do -end function externalheat_init +end function source_externalheat_init !-------------------------------------------------------------------------------------------------- !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- -module subroutine externalheat_dotState(ph, en) +module subroutine source_externalheat_dotState(ph, en) integer, intent(in) :: & ph, & en - integer :: & - so - so = source_thermal_externalheat_offset(ph) + thermalState(ph)%p(source_ID(ph))%dotState(1,en) = 1.0_pREAL ! state is current time - thermalState(ph)%p(so)%dotState(1,en) = 1.0_pREAL ! state is current time - -end subroutine externalheat_dotState +end subroutine source_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module function externalheat_f_T(ph,en) result(f_T) +module function source_externalheat_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & @@ -108,16 +104,11 @@ module function externalheat_f_T(ph,en) result(f_T) real(pREAL) :: & f_T - integer :: & - so - - - so = source_thermal_externalheat_offset(ph) associate(prm => param(ph)) - f_T = prm%f%at(thermalState(ph)%p(so)%state(1,en)) + f_T = prm%f%at(thermalState(ph)%p(source_ID(ph))%state(1,en)) end associate -end function externalheat_f_T +end function source_externalheat_f_T -end submodule externalheat +end submodule source_externalheat From 2ae1f1f829c96e0a418dd58f107e6370b4af9e36 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 18 Jul 2023 04:42:14 +0200 Subject: [PATCH 02/13] simplified --- src/Marc/DAMASK_Marc.f90 | 1 - src/phase_damage.f90 | 6 ++-- src/phase_mechanical_eigen.f90 | 6 ++-- ...phase_mechanical_eigen_cleavageopening.f90 | 30 ------------------- 4 files changed, 6 insertions(+), 37 deletions(-) delete mode 100644 src/phase_mechanical_eigen_cleavageopening.f90 diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index b1d9d91ba..69ee7091c 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -173,7 +173,6 @@ end module DAMASK_interface #include "../phase_mechanical_plastic_dislotungsten.f90" #include "../phase_mechanical_plastic_nonlocal.f90" #include "../phase_mechanical_eigen.f90" -#include "../phase_mechanical_eigen_cleavageopening.f90" #include "../phase_mechanical_eigen_thermalexpansion.f90" #include "../phase_thermal.f90" #include "../phase_thermal_source_dissipation.f90" diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 729309b82..f3afaec99 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -18,14 +18,14 @@ submodule(phase) damage integer :: phase_damage_maxSizeDotState - type :: tDataContainer + type :: tFieldQuantities real(pREAL), dimension(:), allocatable :: phi - end type tDataContainer + end type tFieldQuantities integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & phase_damage !< active sources mechanisms of each phase - type(tDataContainer), dimension(:), allocatable :: current + type(tFieldQuantities), dimension(:), allocatable :: current type(tDamageParameters), dimension(:), allocatable :: param diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index ec1bcfbbc..7b965bf23 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -68,7 +68,7 @@ module subroutine eigen_init(phases) allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID) - where(damage_anisobrittle_init()) model_damage = EIGEN_cleavage_opening_ID + where(kinematics_active2('anisobrittle')) model_damage = EIGEN_cleavage_opening_ID end subroutine eigen_init @@ -191,13 +191,13 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end select kinematicsType end do KinematicsLoop - select case (model_damage(ph)) + damageType: select case (model_damage(ph)) case (EIGEN_cleavage_opening_ID) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. - end select + end select damageType if (.not. active) return diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 deleted file mode 100644 index 780ed22b2..000000000 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ /dev/null @@ -1,30 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes -!> @details to be done -!-------------------------------------------------------------------------------------------------- -submodule(phase:eigen) cleavageopening - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -module function damage_anisobrittle_init() result(myKinematics) - - logical, dimension(:), allocatable :: myKinematics - - - myKinematics = kinematics_active2('anisobrittle') - if (count(myKinematics) == 0) return - - print'(/,1x,a)', '<<<+- phase:mechanical:eigen:cleavageopening init -+>>>' - print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) - -end function damage_anisobrittle_init - - -end submodule cleavageopening From 26014aec1f50c88c486a4b90ef051e718f2d74f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 18 Jul 2023 04:51:16 +0200 Subject: [PATCH 03/13] 'true' kinematics first --- src/phase_mechanical_eigen.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 7b965bf23..dc2d3d598 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -173,14 +173,6 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & dLi_dFi = 0.0_pREAL - plasticType: select case (phase_plasticity(ph)) - case (PLASTIC_isotropic_ID) plasticType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - active = .true. - end select plasticType - KinematicsLoop: do k = 1, Nmodels(ph) kinematicsType: select case (model(k,ph)) case (EIGEN_thermal_expansion_ID) kinematicsType @@ -191,6 +183,14 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end select kinematicsType end do KinematicsLoop + plasticType: select case (phase_plasticity(ph)) + case (PLASTIC_isotropic_ID) plasticType + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. + end select plasticType + damageType: select case (model_damage(ph)) case (EIGEN_cleavage_opening_ID) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) From dc9d4bb5a9d26fea87353609812ae86d8a8a8e11 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 18 Jul 2023 05:28:44 +0200 Subject: [PATCH 04/13] centralized ID handling to enable cross-talking --- src/phase.f90 | 22 ++++++++ src/phase_damage.f90 | 47 +++++++--------- src/phase_mechanical.f90 | 34 +++--------- src/phase_mechanical_eigen.f90 | 26 +++------ src/phase_mechanical_elastic.f90 | 4 +- src/phase_mechanical_plastic.f90 | 68 +++++++++++------------ src/phase_mechanical_plastic_nonlocal.f90 | 4 +- src/phase_thermal.f90 | 22 +++----- 8 files changed, 102 insertions(+), 125 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index f889a854f..c82cc9692 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -49,6 +49,28 @@ module phase type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase end type + enum, bind(c); enumerator :: & + UNDEFINED, & + MECHANICAL_PLASTICITY_NONE, & + MECHANICAL_PLASTICITY_ISOTROPIC, & + MECHANICAL_PLASTICITY_PHENOPOWERLAW, & + MECHANICAL_PLASTICITY_KINEHARDENING, & + MECHANICAL_PLASTICITY_DISLOTWIN, & + MECHANICAL_PLASTICITY_DISLOTUNGSTEN, & + MECHANICAL_PLASTICITY_NONLOCAL, & + MECHANICAL_EIGEN_THERMALEXPANSION, & + DAMAGE_ISOBRITTLE, & + DAMAGE_ANISOBRITTLE, & + THERMAL_SOURCE_DISSIPATION, & + THERMAL_SOURCE_EXTERNALHEAT + end enum + + + integer(kind(UNDEFINED)), dimension(:), allocatable :: & + mechanical_plasticity_type, & !< plasticity of each phase + damage_type !< active sources mechanisms of each phase + integer(kind(UNDEFINED)), dimension(:,:), allocatable :: & + thermal_source_type character(len=2), allocatable, dimension(:) :: phase_lattice real(pREAL), allocatable, dimension(:) :: phase_cOverA diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index f3afaec99..44c52bc80 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -9,21 +9,12 @@ submodule(phase) damage l_c = 0.0_pREAL !< characteristic length end type tDamageParameters - enum, bind(c); enumerator :: & - DAMAGE_UNDEFINED_ID, & - DAMAGE_ISOBRITTLE_ID, & - DAMAGE_ANISOBRITTLE_ID - end enum - integer :: phase_damage_maxSizeDotState - type :: tFieldQuantities real(pREAL), dimension(:), allocatable :: phi end type tFieldQuantities - integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & - phase_damage !< active sources mechanisms of each phase type(tFieldQuantities), dimension(:), allocatable :: current @@ -114,11 +105,11 @@ module subroutine damage_init() end do - allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) + allocate(damage_type(phases%length), source = UNDEFINED) if (damage_active) then - where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID - where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID + where(isobrittle_init() ) damage_type = DAMAGE_ISOBRITTLE + where(anisobrittle_init()) damage_type = DAMAGE_ANISOBRITTLE end if phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) @@ -159,8 +150,8 @@ module function phase_damage_C66(C66,ph,en) result(C66_degraded) real(pREAL), dimension(6,6) :: C66_degraded - damageType: select case (phase_damage(ph)) - case (DAMAGE_ISOBRITTLE_ID) damageType + damageType: select case (damage_type(ph)) + case (DAMAGE_ISOBRITTLE) damageType C66_degraded = C66 * damage_phi(ph,en)**2 case default damageType C66_degraded = C66 @@ -207,8 +198,8 @@ module function phase_f_phi(phi,co,ce) result(f) ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) + select case(damage_type(ph)) + case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) f = 1.0_pREAL & - 2.0_pREAL * phi*damageState(ph)%state(1,en) case default @@ -318,8 +309,8 @@ module subroutine damage_restartWrite(groupHandle,ph) integer, intent(in) :: ph - select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) + select case(damage_type(ph)) + case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) call HDF5_write(damageState(ph)%state,groupHandle,'omega_damage') end select @@ -332,8 +323,8 @@ module subroutine damage_restartRead(groupHandle,ph) integer, intent(in) :: ph - select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) + select case(damage_type(ph)) + case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) call HDF5_read(damageState(ph)%state0,groupHandle,'omega_damage') end select @@ -350,15 +341,15 @@ module subroutine damage_result(group,ph) integer, intent(in) :: ph - if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) & + if (damage_type(ph) /= UNDEFINED) & call result_closeGroup(result_addGroup(group//'damage')) - sourceType: select case (phase_damage(ph)) + sourceType: select case (damage_type(ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType + case (DAMAGE_ISOBRITTLE) sourceType call isobrittle_result(ph,group//'damage/') - case (DAMAGE_ANISOBRITTLE_ID) sourceType + case (DAMAGE_ANISOBRITTLE) sourceType call anisobrittle_result(ph,group//'damage/') end select sourceType @@ -381,9 +372,9 @@ function phase_damage_collectDotState(ph,en) result(broken) if (damageState(ph)%sizeState > 0) then - sourceType: select case (phase_damage(ph)) + sourceType: select case (damage_type(ph)) - case (DAMAGE_ANISOBRITTLE_ID) sourceType + case (DAMAGE_ANISOBRITTLE) sourceType call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! ToDo: use M_d end select sourceType @@ -446,9 +437,9 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) if (damageState(ph)%sizeState == 0) return - sourceType: select case (phase_damage(ph)) + sourceType: select case (damage_type(ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType + 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 diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 0f931517a..2ec8956aa 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -3,21 +3,6 @@ !---------------------------------------------------------------------------------------------------- submodule(phase) mechanical - - enum, bind(c); enumerator :: & - PLASTIC_UNDEFINED_ID, & - PLASTIC_NONE_ID, & - PLASTIC_ISOTROPIC_ID, & - PLASTIC_PHENOPOWERLAW_ID, & - PLASTIC_KINEHARDENING_ID, & - PLASTIC_DISLOTWIN_ID, & - PLASTIC_DISLOTUNGSTEN_ID, & - PLASTIC_NONLOCAL_ID, & - EIGEN_UNDEFINED_ID, & - EIGEN_CLEAVAGE_OPENING_ID, & - EIGEN_THERMAL_EXPANSION_ID - end enum - type(tTensorContainer), dimension(:), allocatable :: & ! current value phase_mechanical_Fe, & @@ -37,9 +22,6 @@ submodule(phase) mechanical phase_mechanical_S0 - integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: & - phase_plasticity !< plasticity of each phase - interface module subroutine eigen_init(phases) @@ -283,7 +265,7 @@ module subroutine mechanical_init(phases) call elastic_init(phases) allocate(plasticState(phases%length)) - allocate(phase_plasticity(phases%length),source = PLASTIC_UNDEFINED_ID) + allocate(mechanical_plasticity_type(phases%length),source = UNDEFINED) call plastic_init() do ph = 1,phases%length plasticState(ph)%state0 = plasticState(ph)%state @@ -327,24 +309,24 @@ module subroutine mechanical_result(group,ph) call results(group,ph) - select case(phase_plasticity(ph)) + select case(mechanical_plasticity_type(ph)) - case(PLASTIC_ISOTROPIC_ID) + case(MECHANICAL_PLASTICITY_ISOTROPIC) call plastic_isotropic_result(ph,group//'mechanical/') - case(PLASTIC_PHENOPOWERLAW_ID) + case(MECHANICAL_PLASTICITY_PHENOPOWERLAW) call plastic_phenopowerlaw_result(ph,group//'mechanical/') - case(PLASTIC_KINEHARDENING_ID) + case(MECHANICAL_PLASTICITY_KINEHARDENING) call plastic_kinehardening_result(ph,group//'mechanical/') - case(PLASTIC_DISLOTWIN_ID) + case(MECHANICAL_PLASTICITY_DISLOTWIN) call plastic_dislotwin_result(ph,group//'mechanical/') - case(PLASTIC_DISLOTUNGSTEN_ID) + case(MECHANICAL_PLASTICITY_DISLOTUNGSTEN) call plastic_dislotungsten_result(ph,group//'mechanical/') - case(PLASTIC_NONLOCAL_ID) + case(MECHANICAL_PLASTICITY_NONLOCAL) call plastic_nonlocal_result(ph,group//'mechanical/') end select diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index dc2d3d598..4d5771cca 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -3,15 +3,10 @@ submodule(phase:mechanical) eigen integer, dimension(:), allocatable :: & Nmodels - integer(kind(EIGEN_UNDEFINED_ID)), dimension(:,:), allocatable :: & + integer(kind(UNDEFINED)), dimension(:,:), allocatable :: & model - integer(kind(EIGEN_UNDEFINED_ID)), dimension(:), allocatable :: & - model_damage interface - module function damage_anisobrittle_init() result(myKinematics) - logical, dimension(:), allocatable :: myKinematics - end function damage_anisobrittle_init module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -60,17 +55,12 @@ module subroutine eigen_init(phases) Nmodels(ph) = kinematics%length end do - allocate(model(maxval(Nmodels),phases%length), source = EIGEN_undefined_ID) + allocate(model(maxval(Nmodels),phases%length), source = UNDEFINED) if (maxval(Nmodels) /= 0) then - where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID + where(thermalexpansion_init(maxval(Nmodels))) model = MECHANICAL_EIGEN_THERMALEXPANSION end if - allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID) - - where(kinematics_active2('anisobrittle')) model_damage = EIGEN_cleavage_opening_ID - - end subroutine eigen_init @@ -175,7 +165,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & KinematicsLoop: do k = 1, Nmodels(ph) kinematicsType: select case (model(k,ph)) - case (EIGEN_thermal_expansion_ID) kinematicsType + case (MECHANICAL_EIGEN_THERMALEXPANSION) kinematicsType call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS @@ -183,16 +173,16 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end select kinematicsType end do KinematicsLoop - plasticType: select case (phase_plasticity(ph)) - case (PLASTIC_isotropic_ID) plasticType + plasticType: select case (mechanical_plasticity_type(ph)) + case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. end select plasticType - damageType: select case (model_damage(ph)) - case (EIGEN_cleavage_opening_ID) + damageType: select case (damage_type(ph)) + case (DAMAGE_ANISOBRITTLE) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 75a8753a5..08055af98 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -199,8 +199,8 @@ module function phase_homogenizedC66(ph,en) result(C) integer, intent(in) :: ph, en - plasticType: select case (phase_plasticity(ph)) - case (PLASTIC_DISLOTWIN_ID) plasticType + plasticType: select case (mechanical_plasticity_type(ph)) + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType C = elastic_C66(ph,en) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 0c1959660..5d1a462e7 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -211,17 +211,17 @@ contains module subroutine plastic_init - print'(/,1x,a)', '<<<+- phase:mechanical:plastic init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:plasticity init -+>>>' - where(plastic_none_init()) phase_plasticity = PLASTIC_NONE_ID - where(plastic_isotropic_init()) phase_plasticity = PLASTIC_ISOTROPIC_ID - where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTIC_PHENOPOWERLAW_ID - where(plastic_kinehardening_init()) phase_plasticity = PLASTIC_KINEHARDENING_ID - where(plastic_dislotwin_init()) phase_plasticity = PLASTIC_DISLOTWIN_ID - where(plastic_dislotungsten_init()) phase_plasticity = PLASTIC_DISLOTUNGSTEN_ID - where(plastic_nonlocal_init()) phase_plasticity = PLASTIC_NONLOCAL_ID + where(plastic_none_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONE + where(plastic_isotropic_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_ISOTROPIC + where(plastic_phenopowerlaw_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_PHENOPOWERLAW + where(plastic_kinehardening_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_KINEHARDENING + where(plastic_dislotwin_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTWIN + where(plastic_dislotungsten_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTUNGSTEN + where(plastic_nonlocal_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONLOCAL - if (any(phase_plasticity == PLASTIC_undefined_ID)) call IO_error(201) + if (any(mechanical_plasticity_type == UNDEFINED)) call IO_error(201) end subroutine plastic_init @@ -251,7 +251,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & i, j - if (phase_plasticity(ph) == PLASTIC_NONE_ID) then + if (mechanical_plasticity_type(ph) == MECHANICAL_PLASTICITY_NONE) then Lp = 0.0_pREAL dLp_dFi = 0.0_pREAL dLp_dS = 0.0_pREAL @@ -259,24 +259,24 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & Mp = matmul(matmul(transpose(Fi),Fi),S) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_ISOTROPIC_ID) plasticType + case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_PHENOPOWERLAW_ID) plasticType + case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_KINEHARDENING_ID) plasticType + case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_DISLOTWIN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) end select plasticType @@ -308,28 +308,28 @@ module function plastic_dotState(subdt,ph,en) result(dotState) dotState - if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then + if (mechanical_plasticity_type(ph) /= MECHANICAL_PLASTICITY_NONE) then Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_ISOTROPIC_ID) plasticType + case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType dotState = isotropic_dotState(Mp,ph,en) - case (PLASTIC_PHENOPOWERLAW_ID) plasticType + case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType dotState = phenopowerlaw_dotState(Mp,ph,en) - case (PLASTIC_KINEHARDENING_ID) plasticType + case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType dotState = plastic_kinehardening_dotState(Mp,ph,en) - case (PLASTIC_DISLOTWIN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType dotState = dislotwin_dotState(Mp,ph,en) - case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType dotState = dislotungsten_dotState(Mp,ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_dotState(Mp,subdt,ph,en) dotState = plasticState(ph)%dotState(:,en) @@ -349,15 +349,15 @@ module subroutine plastic_dependentState(ph,en) en - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_DISLOTWIN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType call dislotwin_dependentState(ph,en) - case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType call dislotungsten_dependentState(ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_dependentState(ph,en) end select plasticType @@ -384,19 +384,19 @@ module function plastic_deltaState(ph, en) result(broken) broken = .false. - select case (phase_plasticity(ph)) - case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID) + select case (mechanical_plasticity_type(ph)) + case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING) Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_S(ph)%data(1:3,1:3,en)) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_KINEHARDENING_ID) plasticType + case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType call plastic_kinehardening_deltaState(Mp,ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call plastic_nonlocal_deltaState(Mp,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 4e45066b5..272755936 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1252,7 +1252,7 @@ function rhoDotFlux(timestep,ph,en) !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility if (neighbor_n > 0) then - if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. & + if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL .and. & any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then forall (s = 1:ns, t = 1:4) @@ -1298,7 +1298,7 @@ function rhoDotFlux(timestep,ph,en) !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. if (opposite_n > 0) then - if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then + if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL) then normal_me2neighbor_defConf = math_det33(Favg) & * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 8f59dc874..20b8560ac 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -15,17 +15,9 @@ submodule(phase) thermal type(tSourceState), allocatable, dimension(:) :: & thermalState - enum, bind(c); enumerator :: & - THERMAL_UNDEFINED_ID ,& - THERMAL_DISSIPATION_ID, & - THERMAL_EXTERNALHEAT_ID - end enum - type :: tFieldQuantities real(pREAL), dimension(:), allocatable :: T, dot_T end type tFieldQuantities - integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & - thermal_source type(tFieldQuantities), dimension(:), allocatable :: current @@ -129,11 +121,11 @@ module subroutine thermal_init(phases) end do - allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) + allocate(thermal_source_type(maxval(thermal_Nsources),phases%length), source = UNDEFINED) if (maxval(thermal_Nsources) /= 0) then - where(source_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID - where(source_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID + where(source_dissipation_init (maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_DISSIPATION + where(source_externalheat_init(maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_EXTERNALHEAT end if thermal_source_maxSizeDotState = 0 @@ -165,12 +157,12 @@ module function phase_f_T(ph,en) result(f) f = 0.0_pREAL do so = 1, thermal_Nsources(ph) - select case(thermal_source(so,ph)) + select case(thermal_source_type(so,ph)) - case (THERMAL_DISSIPATION_ID) + case (THERMAL_SOURCE_DISSIPATION) f = f + source_dissipation_f_T(ph,en) - case (THERMAL_EXTERNALHEAT_ID) + case (THERMAL_SOURCE_EXTERNALHEAT) f = f + source_externalheat_f_T(ph,en) end select @@ -195,7 +187,7 @@ function phase_thermal_collectDotState(ph,en) result(ok) SourceLoop: do i = 1, thermal_Nsources(ph) - if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & + 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))) From 5b56b13c64ca3e4cbe2578fdc43edb328304c930 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Jul 2023 02:30:07 +0200 Subject: [PATCH 05/13] report information on staggered iterations --- src/grid/DAMASK_grid.f90 | 10 ++++++---- src/grid/grid_damage_spectral.f90 | 3 --- src/grid/grid_thermal_spectral.f90 | 3 --- 3 files changed, 6 insertions(+), 10 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index d52a6f211..9dcae70a1 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -341,7 +341,7 @@ program DAMASK_grid if (worldrank == 0) then writeHeader: if (CLI_restartInc < 1) then open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') - write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file + write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded StagIterationsNeeded' ! statistics file else writeHeader open(newunit=statUnit,file=trim(getSolverJobName())//& '.sta',form='FORMATTED', position='APPEND', status='OLD') @@ -415,9 +415,11 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! solve fields - stagIter = 0 + stagIter = 1 stagIterate = .true. do while (stagIterate) + + if (nActiveFields > 1) print'(/,1x,a,i0)', 'Staggered Iteration ',stagIter do field = 1, nActiveFields select case(ID(field)) case(FIELD_MECH_ID) @@ -432,7 +434,7 @@ program DAMASK_grid end do stagIter = stagIter + 1 - stagIterate = stagIter < stagItMax & + stagIterate = stagIter <= stagItMax & .and. all(solres(:)%converged) & .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration end do @@ -448,7 +450,7 @@ program DAMASK_grid guess = .true. ! start guessing after first converged (sub)inc if (worldrank == 0) then write(statUnit,*) totalIncsCounter, t, cutBackLevel, & - solres(1)%converged, solres(1)%iterationsNeeded + solres(1)%converged, solres(1)%iterationsNeeded, StagIter flush(statUnit) end if elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 27295b76f..31aec3912 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -92,7 +92,6 @@ subroutine grid_damage_spectral_init(num_grid) petsc_options - print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' @@ -227,8 +226,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution) SNESConvergedReason :: reason - solution%converged = .false. - !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%Delta_t = Delta_t diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index f6b5c19a0..a03af881b 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -88,7 +88,6 @@ subroutine grid_thermal_spectral_init(num_grid) petsc_options - print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' @@ -204,8 +203,6 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) SNESConvergedReason :: reason - solution%converged = .false. - !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%Delta_t = Delta_t From ef5fe61ff65a0d9cfe0c30158268ff8fb8af3dcf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Jul 2023 02:58:08 +0200 Subject: [PATCH 06/13] limit scope length --- src/grid/DAMASK_grid.f90 | 241 ++++++++++++++++++++------------------- 1 file changed, 125 insertions(+), 116 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 9dcae70a1..3cf378d7c 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -75,7 +75,7 @@ program DAMASK_grid cutBack = .false.,& sig integer :: & - i, j, m, field, & + i, j, field, & errorID = 0, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ stepFraction = 0, & !< fraction of current time interval @@ -110,15 +110,8 @@ program DAMASK_grid load, & num_solver, & num_grid, & - load_step, & - solver, & - step_bc, & - step_mech, & - step_discretization + solver type(tList), pointer :: & -#ifdef __INTEL_LLVM_COMPILER - tensor, & -#endif load_steps character(len=:), allocatable :: & fileContent, fname @@ -210,113 +203,6 @@ program DAMASK_grid ID(field) = FIELD_DAMAGE_ID end if damageActive - -!-------------------------------------------------------------------------------------------------- - load_steps => load%get_list('loadstep') - allocate(loadCases(load_steps%length)) ! array of load cases - - do l = 1, load_steps%length - - load_step => load_steps%get_dict(l) - step_bc => load_step%get_dict('boundary_conditions') - step_mech => step_bc%get_dict('mechanical') - loadCases(l)%stress%myType='' - readMech: do m = 1, step_mech%length - select case (step_mech%key(m)) - case ('L','dot_F','F') ! assign values for the deformation BC matrix - loadCases(l)%deformation%myType = step_mech%key(m) -#ifdef __INTEL_LLVM_COMPILER - tensor => step_mech%get_list(m) - call getMaskedTensor(loadCases(l)%deformation%values,loadCases(l)%deformation%mask,tensor) -#else - call getMaskedTensor(loadCases(l)%deformation%values,loadCases(l)%deformation%mask,step_mech%get_list(m)) -#endif - case ('dot_P','P') - loadCases(l)%stress%myType = step_mech%key(m) -#ifdef __INTEL_LLVM_COMPILER - tensor => step_mech%get_list(m) - call getMaskedTensor(loadCases(l)%stress%values,loadCases(l)%stress%mask,tensor) -#else - call getMaskedTensor(loadCases(l)%stress%values,loadCases(l)%stress%mask,step_mech%get_list(m)) -#endif - end select - call loadCases(l)%rot%fromAxisAngle(step_mech%get_as1dReal('R',defaultVal = real([0.0,0.0,1.0,0.0],pREAL)),degrees=.true.) - end do readMech - if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing') - - step_discretization => load_step%get_dict('discretization') - loadCases(l)%t = step_discretization%get_asReal('t') - loadCases(l)%N = step_discretization%get_asInt ('N') - loadCases(l)%r = step_discretization%get_asReal('r',defaultVal= 1.0_pREAL) - - loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0)) - if (load_step%get_asStr('f_out',defaultVal='n/a') == 'none') then - loadCases(l)%f_out = huge(0) - else - loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1) - end if - loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) - - reportAndCheck: if (worldrank == 0) then - print'(/,1x,a,1x,i0)', 'load case:', l - print'(2x,a,1x,l1)', 'estimate_rate:', loadCases(l)%estimate_rate - if (loadCases(l)%deformation%myType == 'F') then - print'(2x,a)', 'F:' - else - print'(2x,a)', loadCases(l)%deformation%myType//' / 1/s:' - end if - do i = 1, 3; do j = 1, 3 - if (loadCases(l)%deformation%mask(i,j)) then - write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' - else - write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j) - end if - end do; write(IO_STDOUT,'(/)',advance='no') - end do - if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 - if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) & - errorID = 838 ! no rotation is allowed by stress BC - - if (loadCases(l)%stress%myType == 'P') print'(2x,a)', 'P / MPa:' - if (loadCases(l)%stress%myType == 'dot_P') print'(2x,a)', 'dot_P / MPa/s:' - - if (loadCases(l)%stress%myType /= '') then - do i = 1, 3; do j = 1, 3 - if (loadCases(l)%stress%mask(i,j)) then - write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' - else - write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pREAL - end if - end do; write(IO_STDOUT,'(/)',advance='no') - end do - end if - if (any(dNeq(loadCases(l)%rot%asMatrix(), math_I3))) & - write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',& - transpose(loadCases(l)%rot%asMatrix()) - - if (loadCases(l)%r <= 0.0_pREAL) errorID = 833 - if (loadCases(l)%t < 0.0_pREAL) errorID = 834 - if (loadCases(l)%N < 1) errorID = 835 - if (loadCases(l)%f_out < 1) errorID = 836 - if (loadCases(l)%f_restart < 1) errorID = 839 - - if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then - print'(2x,a)', 'r: 1 (constant step width)' - else - print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r - end if - print'(2x,a,1x,f0.3)', 't:', loadCases(l)%t - print'(2x,a,1x,i0)', 'N:', loadCases(l)%N - if (loadCases(l)%f_out < huge(0)) & - print'(2x,a,1x,i0)', 'f_out:', loadCases(l)%f_out - if (loadCases(l)%f_restart < huge(0)) & - print'(2x,a,1x,i0)', 'f_restart:', loadCases(l)%f_restart - - if (errorID > 0) call IO_error(errorID,label1='line',ID1=l) - - end if reportAndCheck - end do - !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers call spectral_utilities_init() @@ -354,6 +240,8 @@ program DAMASK_grid call materialpoint_result(0,0.0_pREAL) end if writeUndeformed + loadCases = parseLoadSteps(load%get_list('loadstep')) + loadCaseLooping: do l = 1, size(loadCases) t_0 = t ! load case start time guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc @@ -540,4 +428,125 @@ subroutine getMaskedTensor(values,mask,tensor) end subroutine getMaskedTensor + +function parseLoadsteps(load_steps) result(loadCases) + + type(tList), intent(in), target :: load_steps + type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases + + integer :: l,m + type(tDict), pointer :: & + load_step, & + step_bc, & + step_mech, & + step_discretization +#ifdef __INTEL_LLVM_COMPILER + type(tList), pointer :: & + tensor +#endif + + + allocate(loadCases(load_steps%length)) + do l = 1, load_steps%length + load_step => load_steps%get_dict(l) + step_bc => load_step%get_dict('boundary_conditions') + step_mech => step_bc%get_dict('mechanical') + loadCases(l)%stress%myType='' + readMech: do m = 1, step_mech%length + select case (step_mech%key(m)) + case ('L','dot_F','F') ! assign values for the deformation BC matrix + loadCases(l)%deformation%myType = step_mech%key(m) +#ifdef __INTEL_LLVM_COMPILER + tensor => step_mech%get_list(m) + call getMaskedTensor(loadCases(l)%deformation%values,loadCases(l)%deformation%mask,tensor) +#else + call getMaskedTensor(loadCases(l)%deformation%values,loadCases(l)%deformation%mask,step_mech%get_list(m)) +#endif + case ('dot_P','P') + loadCases(l)%stress%myType = step_mech%key(m) +#ifdef __INTEL_LLVM_COMPILER + tensor => step_mech%get_list(m) + call getMaskedTensor(loadCases(l)%stress%values,loadCases(l)%stress%mask,tensor) +#else + call getMaskedTensor(loadCases(l)%stress%values,loadCases(l)%stress%mask,step_mech%get_list(m)) +#endif + end select + call loadCases(l)%rot%fromAxisAngle(step_mech%get_as1dReal('R',defaultVal = real([0.0,0.0,1.0,0.0],pREAL)),degrees=.true.) + end do readMech + if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing') + + step_discretization => load_step%get_dict('discretization') + loadCases(l)%t = step_discretization%get_asReal('t') + loadCases(l)%N = step_discretization%get_asInt ('N') + loadCases(l)%r = step_discretization%get_asReal('r',defaultVal= 1.0_pREAL) + + loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0)) + if (load_step%get_asStr('f_out',defaultVal='n/a') == 'none') then + loadCases(l)%f_out = huge(0) + else + loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1) + end if + loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) + + reportAndCheck: if (worldrank == 0) then + print'(/,1x,a,1x,i0)', 'load case:', l + print'(2x,a,1x,l1)', 'estimate_rate:', loadCases(l)%estimate_rate + if (loadCases(l)%deformation%myType == 'F') then + print'(2x,a)', 'F:' + else + print'(2x,a)', loadCases(l)%deformation%myType//' / 1/s:' + end if + do i = 1, 3; do j = 1, 3 + if (loadCases(l)%deformation%mask(i,j)) then + write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' + else + write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j) + end if + end do; write(IO_STDOUT,'(/)',advance='no') + end do + if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 + if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) & + errorID = 838 ! no rotation is allowed by stress BC + + if (loadCases(l)%stress%myType == 'P') print'(2x,a)', 'P / MPa:' + if (loadCases(l)%stress%myType == 'dot_P') print'(2x,a)', 'dot_P / MPa/s:' + + if (loadCases(l)%stress%myType /= '') then + do i = 1, 3; do j = 1, 3 + if (loadCases(l)%stress%mask(i,j)) then + write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' + else + write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pREAL + end if + end do; write(IO_STDOUT,'(/)',advance='no') + end do + end if + if (any(dNeq(loadCases(l)%rot%asMatrix(), math_I3))) & + write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',& + transpose(loadCases(l)%rot%asMatrix()) + + if (loadCases(l)%r <= 0.0_pREAL) errorID = 833 + if (loadCases(l)%t < 0.0_pREAL) errorID = 834 + if (loadCases(l)%N < 1) errorID = 835 + if (loadCases(l)%f_out < 1) errorID = 836 + if (loadCases(l)%f_restart < 1) errorID = 839 + + if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then + print'(2x,a)', 'r: 1 (constant step width)' + else + print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r + end if + print'(2x,a,1x,f0.3)', 't:', loadCases(l)%t + print'(2x,a,1x,i0)', 'N:', loadCases(l)%N + if (loadCases(l)%f_out < huge(0)) & + print'(2x,a,1x,i0)', 'f_out:', loadCases(l)%f_out + if (loadCases(l)%f_restart < huge(0)) & + print'(2x,a,1x,i0)', 'f_restart:', loadCases(l)%f_restart + + if (errorID > 0) call IO_error(errorID,label1='line',ID1=l) + + end if reportAndCheck + end do +end function parseLoadsteps + end program DAMASK_grid From 535ccf13781c0d0548f1702a5a065227952ff665 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Jul 2023 04:19:57 +0200 Subject: [PATCH 07/13] clearer name --- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarization.f90 | 6 +++--- src/grid/spectral_utilities.f90 | 16 ++++++++-------- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 0db9db64b..ed1ff1e2f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -365,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ if (stress_BC%myType=='dot_P') P_aim = P_aim & + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t - F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average + F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3]) call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc) CHKERRQ(err_PETSc) diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 5342fa51a..c6c7c2fd5 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -408,11 +408,11 @@ subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,D if (stress_BC%myType=='dot_P') P_aim = P_aim & + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t - F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average - rotation_BC%rotate(F_aim,active=.true.)),& + F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average + rotation_BC%rotate(F_aim,active=.true.)),& [9,cells(1),cells(2),cells3]) if (guess) then - F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), & + F_tau = reshape(Utilities_forwardTensorField(Delta_t,F_tau_lastInc,F_taudot), & [9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition else do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 453723a9b..03d945bee 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -124,7 +124,7 @@ module spectral_utilities utilities_maskedCompliance, & utilities_constitutiveResponse, & utilities_calculateRate, & - utilities_forwardField, & + utilities_forwardTensorField, & utilities_updateCoords contains @@ -864,7 +864,7 @@ end function utilities_calculateRate !> @brief forwards a field with a pointwise given rate, if aim is given, !> ensures that the average matches the aim !-------------------------------------------------------------------------------------------------- -function utilities_forwardField(Delta_t,field_lastInc,rate,aim) +function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim) real(pREAL), intent(in) :: & Delta_t !< Delta_t of current step @@ -875,22 +875,22 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim) aim !< average field value aim real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & - utilities_forwardField + utilities_forwardTensorField real(pREAL), dimension(3,3) :: fieldDiff !< - aim integer(MPI_INTEGER_KIND) :: err_MPI - utilities_forwardField = field_lastInc + rate*Delta_t + utilities_forwardTensorField = field_lastInc + rate*Delta_t if (present(aim)) then !< correct to match average - fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt + fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' fieldDiff = fieldDiff - aim - utilities_forwardField = utilities_forwardField & - - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) + utilities_forwardTensorField = utilities_forwardTensorField & + - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) end if -end function utilities_forwardField +end function utilities_forwardTensorField !-------------------------------------------------------------------------------------------------- From 1fdc5443b295589bb05ca383a70b132316f7cf75 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Jul 2023 23:33:43 +0200 Subject: [PATCH 08/13] avoid code (loop) duplication --- src/grid/grid_damage_spectral.f90 | 25 ++++--------------------- src/grid/grid_thermal_spectral.f90 | 26 +++++++------------------- src/homogenization.f90 | 11 ++++------- src/homogenization_damage.f90 | 19 +++++++++---------- src/homogenization_thermal.f90 | 15 +++++++++------ 5 files changed, 33 insertions(+), 63 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 31aec3912..48023cc54 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -75,7 +75,6 @@ subroutine grid_damage_spectral_init(num_grid) type(tDict), pointer, intent(in) :: num_grid integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global - integer :: i, j, k, ce DM :: DM_damage real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed Vec :: uBound, lBound @@ -139,7 +138,7 @@ subroutine grid_damage_spectral_init(num_grid) 1_pPETSCINT, 1_pPETSCINT, int(worldsize,pPETSCINT), & 1_pPETSCINT, 0_pPETSCINT, & ! #dof (phi, scalar), ghost boundary width (domain overlap) [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],int(cells3_global,pPETSCINT), & ! local cells - DM_damage,err_PETSc) ! handle, error + DM_damage,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(DM_damage,err_PETSc) CHKERRQ(err_PETSc) @@ -193,11 +192,7 @@ subroutine grid_damage_spectral_init(num_grid) phi_stagInc = phi_lastInc end if restartRead - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 - ce = ce + 1 - call homogenization_set_phi(phi(i,j,k),ce) - end do; end do; end do + call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3])) call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) CHKERRQ(err_PETSc) @@ -215,7 +210,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution) real(pREAL), intent(in) :: & Delta_t !< increment in time for current solution - integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm @@ -253,13 +247,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' phi_stagInc = phi -!-------------------------------------------------------------------------------------------------- -! updating damage state - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0,cells(1)-1 - ce = ce + 1 - call homogenization_set_phi(phi(i,j,k),ce) - end do; end do; end do + call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3])) call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) CHKERRQ(err_PETSc) @@ -280,7 +268,6 @@ subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, ce DM :: DM_damage real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed PetscErrorCode :: err_PETSc @@ -292,11 +279,7 @@ subroutine grid_damage_spectral_forward(cutBack) CHKERRQ(err_PETSc) if (cutBack) then - ce = 0 - do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) - ce = ce + 1 - call homogenization_set_phi(phi_lastInc(i,j,k),ce) - end do; end do; end do + call homogenization_set_phi(reshape(phi_lastInc,[product(cells(1:2))*cells3])) phi = phi_lastInc phi_stagInc = phi_lastInc else diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index a03af881b..86d8e04f3 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -74,7 +74,7 @@ subroutine grid_thermal_spectral_init(num_grid) type(tDict), pointer, intent(in) :: num_grid integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global - integer :: i, j, k, ce + integer :: ce DM :: DM_thermal real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed integer(MPI_INTEGER_KIND) :: err_MPI @@ -170,11 +170,8 @@ subroutine grid_thermal_spectral_init(num_grid) dotT_lastInc = 0.0_pREAL * T_lastInc end if restartRead - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 - ce = ce + 1 - call homogenization_thermal_setField(T(i,j,k),0.0_pREAL,ce) - end do; end do; end do + call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), & + [(0.0_pReal, ce = 1,product(cells(1:2))*cells3)]) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) CHKERRQ(err_PETSc) @@ -230,13 +227,8 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' T_stagInc = T -!-------------------------------------------------------------------------------------------------- -! updating thermal state - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 - ce = ce + 1 - call homogenization_thermal_setField(T(i,j,k),(T(i,j,k)-T_lastInc(i+1,j+1,k+1))/params%Delta_t,ce) - end do; end do; end do + call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), & + reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) CHKERRQ(err_PETSc) @@ -257,7 +249,6 @@ subroutine grid_thermal_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, ce DM :: DM_thermal real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed PetscErrorCode :: err_PETSc @@ -269,11 +260,8 @@ subroutine grid_thermal_spectral_forward(cutBack) CHKERRQ(err_PETSc) if (cutBack) then - ce = 0 - do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) - ce = ce + 1 - call homogenization_thermal_setField(T_lastInc(i,j,k),dotT_lastInc(i,j,k),ce) - end do; end do; end do + call homogenization_thermal_setField(reshape(T_lastInc,[product(cells(1:2))*cells3]), & + reshape(dotT_lastInc,[product(cells(1:2))*cells3])) T = T_lastInc T_stagInc = T_lastInc else diff --git a/src/homogenization.f90 b/src/homogenization.f90 index c85547917..63eacdf9d 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -145,9 +145,8 @@ module homogenization real(pREAL) :: f end function homogenization_f_T - module subroutine homogenization_thermal_setField(T,dot_T, ce) - integer, intent(in) :: ce - real(pREAL), intent(in) :: T, dot_T + module subroutine homogenization_thermal_setField(T,dot_T) + real(pREAL), dimension(:), intent(in) :: T, dot_T end subroutine homogenization_thermal_setField module function homogenization_damage_active() result(active) @@ -170,10 +169,8 @@ module homogenization real(pREAL) :: f end function homogenization_f_phi - module subroutine homogenization_set_phi(phi,ce) - integer, intent(in) :: ce - real(pREAL), intent(in) :: & - phi + module subroutine homogenization_set_phi(phi) + real(pREAL), dimension(:), intent(in) :: phi end subroutine homogenization_set_phi end interface diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 466b8b47b..233425ebe 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -151,20 +151,19 @@ end function homogenization_f_phi !-------------------------------------------------------------------------------------------------- !> @brief Set damage field. !-------------------------------------------------------------------------------------------------- -module subroutine homogenization_set_phi(phi,ce) +module subroutine homogenization_set_phi(phi) - integer, intent(in) :: ce - real(pREAL), intent(in) :: phi + real(pREAL), dimension(:), intent(in) :: phi - integer :: & - ho, & - en + integer :: ho, en, ce - ho = material_ID_homogenization(ce) - en = material_entry_homogenization(ce) - damagestate_h(ho)%state(1,en) = phi - current(ho)%phi(en) = phi + do ce=1, ubound(phi,1) + ho = material_ID_homogenization(ce) + en = material_entry_homogenization(ce) + damagestate_h(ho)%state(1,en) = phi(ce) + current(ho)%phi(en) = phi(ce) + end do end subroutine homogenization_set_phi diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 789ac994b..791286912 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -173,15 +173,18 @@ end function homogenization_f_T !-------------------------------------------------------------------------------------------------- !> @brief Set thermal field and its rate (T and dot_T). !-------------------------------------------------------------------------------------------------- -module subroutine homogenization_thermal_setField(T,dot_T, ce) +module subroutine homogenization_thermal_setField(T,dot_T) - integer, intent(in) :: ce - real(pREAL), intent(in) :: T, dot_T + real(pREAL), dimension(:), intent(in) :: T, dot_T + + integer :: ce - current(material_ID_homogenization(ce))%T(material_entry_homogenization(ce)) = T - current(material_ID_homogenization(ce))%dot_T(material_entry_homogenization(ce)) = dot_T - call thermal_partition(ce) + do ce=1, min(ubound(T,1),ubound(dot_T,1)) + current(material_ID_homogenization(ce))%T(material_entry_homogenization(ce)) = T(ce) + current(material_ID_homogenization(ce))%dot_T(material_entry_homogenization(ce)) = dot_T(ce) + call thermal_partition(ce) + end do end subroutine homogenization_thermal_setField From 4086f40a1661ef56fd0a792baa9ebb85ef15dbcc Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 22 Jul 2023 14:01:53 -0400 Subject: [PATCH 09/13] pull "constants" out of inner loops --- src/phase_thermal_source_dissipation.f90 | 4 ++-- src/phase_thermal_source_externalheat.f90 | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/phase_thermal_source_dissipation.f90 b/src/phase_thermal_source_dissipation.f90 index c42e4a817..4713fb809 100644 --- a/src/phase_thermal_source_dissipation.f90 +++ b/src/phase_thermal_source_dissipation.f90 @@ -49,8 +49,9 @@ module function source_dissipation_init(source_length) result(mySources) allocate(param(phases%length)) do ph = 1, phases%length - phase => phases%get_dict(ph) if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 + Nmembers = count(material_ID_phase == ph) + phase => phases%get_dict(ph) thermal => phase%get_dict('thermal') sources => thermal%get_list('source') do so = 1, sources%length @@ -62,7 +63,6 @@ module function source_dissipation_init(source_length) result(mySources) if (len(refs) > 0) print'(/,1x,a)', refs prm%kappa = src%get_asReal('kappa') - Nmembers = count(material_ID_phase == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) end associate diff --git a/src/phase_thermal_source_externalheat.f90 b/src/phase_thermal_source_externalheat.f90 index f86880e62..24a5ea7f8 100644 --- a/src/phase_thermal_source_externalheat.f90 +++ b/src/phase_thermal_source_externalheat.f90 @@ -52,8 +52,9 @@ module function source_externalheat_init(source_length) result(mySources) allocate(source_ID(phases%length), source=0) do ph = 1, phases%length - phase => phases%get_dict(ph) if (count(mySources(:,ph)) == 0) cycle + Nmembers = count(material_ID_phase == ph) + phase => phases%get_dict(ph) thermal => phase%get_dict('thermal') sources => thermal%get_list('source') do so = 1, sources%length @@ -66,8 +67,6 @@ module function source_externalheat_init(source_length) result(mySources) if (len(refs) > 0) print'(/,1x,a)', refs prm%f = table(src,'t','f') - - Nmembers = count(material_ID_phase == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate end if From 2d0b64e00dafedf5523355fcab662c0eb12bf47d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 22 Jul 2023 15:04:11 -0400 Subject: [PATCH 10/13] better variable names --- src/phase_thermal.f90 | 12 ++++++------ src/phase_thermal_source_dissipation.f90 | 18 +++++++++--------- src/phase_thermal_source_externalheat.f90 | 22 +++++++++++----------- 3 files changed, 26 insertions(+), 26 deletions(-) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 7f31cecb2..e21f173ad 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -28,14 +28,14 @@ submodule(phase) thermal interface - module function source_dissipation_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function source_dissipation_init(maxNsources) result(isMySource) + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource end function source_dissipation_init - module function source_externalheat_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function source_externalheat_init(maxNsources) result(isMySource) + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource end function source_externalheat_init diff --git a/src/phase_thermal_source_dissipation.f90 b/src/phase_thermal_source_dissipation.f90 index 4713fb809..aaaffad13 100644 --- a/src/phase_thermal_source_dissipation.f90 +++ b/src/phase_thermal_source_dissipation.f90 @@ -22,10 +22,10 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_dissipation_init(source_length) result(mySources) +module function source_dissipation_init(maxNsources) result(isMySource) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource type(tDict), pointer :: & phases, & @@ -35,27 +35,27 @@ module function source_dissipation_init(source_length) result(mySources) class(tList), pointer :: & sources character(len=:), allocatable :: refs - integer :: so,Nmembers,ph + integer :: ph,Nmembers,so - mySources = thermal_active('dissipation',source_length) - if (count(mySources) == 0) return + isMySource = thermal_active('dissipation',maxNsources) + if (count(isMySource) == 0) return print'(/,1x,a)', '<<<+- phase:thermal:source_dissipation init -+>>>' - print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT) phases => config_material%get_dict('phase') allocate(param(phases%length)) do ph = 1, phases%length - if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 + if (count(isMySource(:,ph)) == 0) cycle !ToDo: error if > 1 Nmembers = count(material_ID_phase == ph) phase => phases%get_dict(ph) thermal => phase%get_dict('thermal') sources => thermal%get_list('source') do so = 1, sources%length - if (mySources(so,ph)) then + if (isMySource(so,ph)) then associate(prm => param(ph)) src => sources%get_dict(so) print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so diff --git a/src/phase_thermal_source_externalheat.f90 b/src/phase_thermal_source_externalheat.f90 index 24a5ea7f8..76fd9ac03 100644 --- a/src/phase_thermal_source_externalheat.f90 +++ b/src/phase_thermal_source_externalheat.f90 @@ -8,10 +8,10 @@ submodule(phase:thermal) source_externalheat integer, dimension(:), allocatable :: & - source_ID !< which source is my current thermal dissipation mechanism? + source_ID !< index in phase source list corresponding to this source type :: tParameters !< container type for internal constitutive parameters - type(tTable) :: f + type(tTable) :: f !< external heat power as (tabulated) function of time end type tParameters type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) @@ -24,10 +24,10 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_externalheat_init(source_length) result(mySources) +module function source_externalheat_init(maxNsources) result(isMySource) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource type(tDict), pointer :: & phases, & @@ -37,14 +37,14 @@ module function source_externalheat_init(source_length) result(mySources) type(tList), pointer :: & sources character(len=:), allocatable :: refs - integer :: so,Nmembers,ph + integer :: ph,Nmembers,so - mySources = thermal_active('externalheat',source_length) - if (count(mySources) == 0) return + isMySource = thermal_active('externalheat',maxNsources) + if (count(isMySource) == 0) return print'(/,1x,a)', '<<<+- phase:thermal:source_externalheat init -+>>>' - print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT) phases => config_material%get_dict('phase') @@ -52,13 +52,13 @@ module function source_externalheat_init(source_length) result(mySources) allocate(source_ID(phases%length), source=0) do ph = 1, phases%length - if (count(mySources(:,ph)) == 0) cycle + if (count(isMySource(:,ph)) == 0) cycle Nmembers = count(material_ID_phase == ph) phase => phases%get_dict(ph) thermal => phase%get_dict('thermal') sources => thermal%get_list('source') do so = 1, sources%length - if (mySources(so,ph)) then + if (isMySource(so,ph)) then source_ID(ph) = so associate(prm => param(ph)) src => sources%get_dict(so) From dafd39566dddcb22bf8c7794f511c2c1d43ebf4b Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 22 Jul 2023 15:17:30 -0400 Subject: [PATCH 11/13] can only have one instance per thermal source --- src/IO.f90 | 2 ++ src/phase_thermal_source_dissipation.f90 | 8 +++++--- src/phase_thermal_source_externalheat.f90 | 7 +++++-- 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index ba33a8a5c..9e0b4c40d 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -548,6 +548,8 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2) !-------------------------------------------------------------------------------------------------- ! user errors + case (600) + msg = 'only one source entry allowed' case (603) msg = 'invalid data for table' case (610) diff --git a/src/phase_thermal_source_dissipation.f90 b/src/phase_thermal_source_dissipation.f90 index aaaffad13..760259286 100644 --- a/src/phase_thermal_source_dissipation.f90 +++ b/src/phase_thermal_source_dissipation.f90 @@ -35,7 +35,7 @@ module function source_dissipation_init(maxNsources) result(isMySource) class(tList), pointer :: & sources character(len=:), allocatable :: refs - integer :: ph,Nmembers,so + integer :: ph,Nmembers,so,Nsources isMySource = thermal_active('dissipation',maxNsources) @@ -49,7 +49,9 @@ module function source_dissipation_init(maxNsources) result(isMySource) allocate(param(phases%length)) do ph = 1, phases%length - if (count(isMySource(:,ph)) == 0) cycle !ToDo: error if > 1 + Nsources = count(isMySource(:,ph)) + if (Nsources == 0) cycle + if (Nsources > 1) call IO_error(600,ext_msg='dissipation') Nmembers = count(material_ID_phase == ph) phase => phases%get_dict(ph) thermal => phase%get_dict('thermal') @@ -64,8 +66,8 @@ module function source_dissipation_init(maxNsources) result(isMySource) prm%kappa = src%get_asReal('kappa') call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) - end associate + exit end if end do end do diff --git a/src/phase_thermal_source_externalheat.f90 b/src/phase_thermal_source_externalheat.f90 index 76fd9ac03..1a803549c 100644 --- a/src/phase_thermal_source_externalheat.f90 +++ b/src/phase_thermal_source_externalheat.f90 @@ -37,7 +37,7 @@ module function source_externalheat_init(maxNsources) result(isMySource) type(tList), pointer :: & sources character(len=:), allocatable :: refs - integer :: ph,Nmembers,so + integer :: ph,Nmembers,so,Nsources isMySource = thermal_active('externalheat',maxNsources) @@ -52,7 +52,9 @@ module function source_externalheat_init(maxNsources) result(isMySource) allocate(source_ID(phases%length), source=0) do ph = 1, phases%length - if (count(isMySource(:,ph)) == 0) cycle + Nsources = count(isMySource(:,ph)) + if (Nsources == 0) cycle + if (Nsources > 1) call IO_error(600,ext_msg='externalheat') Nmembers = count(material_ID_phase == ph) phase => phases%get_dict(ph) thermal => phase%get_dict('thermal') @@ -69,6 +71,7 @@ module function source_externalheat_init(maxNsources) result(isMySource) prm%f = table(src,'t','f') call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate + exit end if end do end do From d953e6bedf687a0f4d56079b48bb2a33519250b4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 24 Jul 2023 18:29:36 +0200 Subject: [PATCH 12/13] unified style --- src/homogenization_damage.f90 | 2 +- src/homogenization_thermal.f90 | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 233425ebe..7e83c34d8 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -158,7 +158,7 @@ module subroutine homogenization_set_phi(phi) integer :: ho, en, ce - do ce=1, ubound(phi,1) + do ce=lbound(phi,1), ubound(phi,1) ho = material_ID_homogenization(ce) en = material_entry_homogenization(ce) damagestate_h(ho)%state(1,en) = phi(ce) diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 791286912..069c402eb 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -177,12 +177,14 @@ module subroutine homogenization_thermal_setField(T,dot_T) real(pREAL), dimension(:), intent(in) :: T, dot_T - integer :: ce + integer :: ho, en, ce - do ce=1, min(ubound(T,1),ubound(dot_T,1)) - current(material_ID_homogenization(ce))%T(material_entry_homogenization(ce)) = T(ce) - current(material_ID_homogenization(ce))%dot_T(material_entry_homogenization(ce)) = dot_T(ce) + do ce=max(lbound(T,1),lbound(dot_T,1)), min(ubound(T,1),ubound(dot_T,1)) + ho = material_ID_homogenization(ce) + en = material_entry_homogenization(ce) + current(ho)%T(en) = T(ce) + current(ho)%dot_T(en) = dot_T(ce) call thermal_partition(ce) end do From 039d5e0fce7b408ff5907a721c8c994e9c33720e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 24 Jul 2023 18:30:02 +0200 Subject: [PATCH 13/13] improved comments and naming --- src/phase.f90 | 5 +++-- src/phase_damage.f90 | 3 ++- src/phase_mechanical_eigen.f90 | 9 +++------ 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 7dbaa9f7c..ab365dbd6 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -68,9 +68,10 @@ module phase integer(kind(UNDEFINED)), dimension(:), allocatable :: & mechanical_plasticity_type, & !< plasticity of each phase - damage_type !< active sources mechanisms of each phase + damage_type !< damage type of each phase integer(kind(UNDEFINED)), dimension(:,:), allocatable :: & - thermal_source_type + thermal_source_type, & + mechanical_eigen_kinematics_type character(len=2), allocatable, dimension(:) :: phase_lattice real(pREAL), allocatable, dimension(:) :: phase_cOverA diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 44c52bc80..2c605c963 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -195,13 +195,14 @@ module function phase_f_phi(phi,co,ce) result(f) ph, & en + ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) select case(damage_type(ph)) case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) f = 1.0_pREAL & - - 2.0_pREAL * phi*damageState(ph)%state(1,en) + - 2.0_pREAL * phi*damageState(ph)%state(1,en) ! ToDo: MD: seems to be phi**2 case default f = 0.0_pREAL end select diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 4d5771cca..5bdb57f57 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -3,9 +3,6 @@ submodule(phase:mechanical) eigen integer, dimension(:), allocatable :: & Nmodels - integer(kind(UNDEFINED)), dimension(:,:), allocatable :: & - model - interface module function thermalexpansion_init(kinematics_length) result(myKinematics) @@ -55,10 +52,10 @@ module subroutine eigen_init(phases) Nmodels(ph) = kinematics%length end do - allocate(model(maxval(Nmodels),phases%length), source = UNDEFINED) + allocate(mechanical_eigen_kinematics_type(maxval(Nmodels),phases%length), source = UNDEFINED) if (maxval(Nmodels) /= 0) then - where(thermalexpansion_init(maxval(Nmodels))) model = MECHANICAL_EIGEN_THERMALEXPANSION + where(thermalexpansion_init(maxval(Nmodels))) mechanical_eigen_kinematics_type = MECHANICAL_EIGEN_THERMALEXPANSION end if end subroutine eigen_init @@ -164,7 +161,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & KinematicsLoop: do k = 1, Nmodels(ph) - kinematicsType: select case (model(k,ph)) + kinematicsType: select case (mechanical_eigen_kinematics_type(k,ph)) case (MECHANICAL_EIGEN_THERMALEXPANSION) kinematicsType call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en) Li = Li + my_Li