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

centralized ID handling to enable cross-talking

This commit is contained in:
Martin Diehl 2023-07-18 05:28:44 +02:00
parent 26014aec1f
commit dc9d4bb5a9
8 changed files with 102 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/phase.f90
View File

@ -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

47
src/phase_damage.f90
View File

@ -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

34
src/phase_mechanical.f90
View File

@ -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

26
src/phase_mechanical_eigen.f90
View File

@ -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

4
src/phase_mechanical_elastic.f90
View File

@ -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)

68
src/phase_mechanical_plastic.f90
View File

@ -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

4
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -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)

22
src/phase_thermal.f90
View File

@ -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)))