Merge remote-tracking branch 'origin/development' into parallel-mesh

This commit is contained in:
Martin Diehl 2022-02-22 14:56:25 +01:00
commit 947f228225
29 changed files with 200 additions and 111 deletions

@ -1 +1 @@
Subproject commit 142be5919d4a4e61e9ad909b6ad7a1ca334fc652
Subproject commit 4c8116ba3b9e9fbb325a580705028e8310139117

View File

@ -0,0 +1 @@
N_constituents: 8

View File

@ -1,4 +0,0 @@
[Parallel3]
mech isostrain
nconstituents 3
mapping sum # or 'parallel'

View File

@ -0,0 +1 @@
N_constituents: 2

View File

@ -0,0 +1,3 @@
# For single point calculations, requires N_constituents = 1
type: pass
output: ['T']

View File

@ -0,0 +1 @@
N_constituents: 1

View File

@ -1,10 +1,8 @@
8Grains:
N_constituents: 8
mechanical:
type: RGC
D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0
cluster_size: [2, 2, 2]
output: [M, Delta_V, avg_dot_a, max_dot_a]
xi_alpha: 10.0
# For Relaxed Grain Cluster homogenization, requires N_constituents = 8
type: RGC
D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0
cluster_size: [2, 2, 2]
output: [M, Delta_V, avg_dot_a, max_dot_a]
xi_alpha: 10.0

View File

@ -1,3 +0,0 @@
Taylor2:
N_constituents: 2
mechanical: {type: isostrain}

View File

@ -0,0 +1,3 @@
# For Taylor homogenization with N_constituents > 1
type: isostrain
output: ['F', 'P']

View File

@ -0,0 +1,3 @@
# For single point calculations, requires N_constituents = 1
type: pass
output: ['F', 'P']

View File

@ -0,0 +1,3 @@
# For homogenization with N_constituents > 1
type: isotemperature
output: ['T']

View File

@ -0,0 +1,3 @@
# For single point calculations, requires N_constituents = 1
type: pass
output: ['T']

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 3.4.1 (RRR=1000, T_min=200K, T_max=900K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 2.380e+2
K_11,T: 2.068e-3
K_11,T^2: -7.765e-5

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 2.4.1 (RRR=1000, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/copper.htm
output: [T]
K_11: 4.039e+2
K_11,T: -8.119e-2
K_11,T^2: 1.454e-5

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 4.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 8.055e+1
K_11,T: -1.051e-1
K_11,T^2: 5.464e-5

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 35R (T_min=150K, T_max=500K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 9.132e+1
K_11,T: -1.525e-1
K_11,T^2: 3.053e-4

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 61R (T_min=100K, T_max=400K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 7.414e+1
K_11,T: -6.465e-2
K_11,T^2: 2.066e-4

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 5.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/tungsten.htm
output: [T]
K_11: 1.758e+2
K_11,T: -1.605e-1
K_11,T^2: 1.160e-4

View File

@ -1 +1 @@
v3.0.0-alpha6-4-gca6a3e786
v3.0.0-alpha6-14-g3657b2316

View File

@ -6,9 +6,10 @@
!--------------------------------------------------------------------------------------------------
module homogenization
use prec
use math
use constants
use IO
use config
use math
use material
use phase
use discretization
@ -32,7 +33,7 @@ module homogenization
HOMOGENIZATION_RGC_ID
end enum
type(tState), allocatable, dimension(:), public :: &
type(tState), allocatable, dimension(:), public :: &
homogState, &
damageState_h
@ -66,15 +67,13 @@ module homogenization
!--------------------------------------------------------------------------------------------------
interface
module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: &
num_homog !< pointer to mechanical homogenization numerics data
module subroutine mechanical_init()
end subroutine mechanical_init
module subroutine thermal_init
module subroutine thermal_init()
end subroutine thermal_init
module subroutine damage_init
module subroutine damage_init()
end subroutine damage_init
module subroutine mechanical_partition(subF,ce)
@ -204,15 +203,15 @@ subroutine homogenization_init()
allocate(homogState (size(material_name_homogenization)))
allocate(damageState_h (size(material_name_homogenization)))
call material_parseHomogenization()
call parseHomogenization()
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
call mechanical_init(num_homog)
call mechanical_init()
call thermal_init()
call damage_init()
@ -323,13 +322,13 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
enddo
call mechanical_homogenize(Delta_t,ce)
enddo IpLooping3
enddo elementLooping3
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
end do
call mechanical_homogenize(Delta_t,ce)
end do IpLooping3
end do elementLooping3
!$OMP END PARALLEL DO
@ -447,7 +446,7 @@ end subroutine homogenization_restartRead
!--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration
!--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization
subroutine parseHomogenization
class(tNode), pointer :: &
material_homogenization, &
@ -459,8 +458,8 @@ subroutine material_parseHomogenization
material_homogenization => config_material%get('homogenization')
allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID)
allocate(thermal_type(size(material_name_homogenization)),source=THERMAL_isothermal_ID)
allocate(damage_type (size(material_name_homogenization)),source=DAMAGE_none_ID)
do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h)
@ -468,7 +467,7 @@ subroutine material_parseHomogenization
if (homog%contains('thermal')) then
homogThermal => homog%get('thermal')
select case (homogThermal%get_asString('type'))
case('pass')
case('pass','isotemperature')
thermal_type(h) = THERMAL_conduction_ID
case default
call IO_error(500,ext_msg=homogThermal%get_asString('type'))
@ -486,7 +485,7 @@ subroutine material_parseHomogenization
endif
enddo
end subroutine material_parseHomogenization
end subroutine parseHomogenization
end module homogenization

View File

@ -7,15 +7,13 @@ submodule(homogenization) mechanical
interface
module subroutine pass_init
module subroutine pass_init()
end subroutine pass_init
module subroutine isostrain_init
module subroutine isostrain_init()
end subroutine isostrain_init
module subroutine RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data
module subroutine RGC_init()
end subroutine RGC_init
@ -52,6 +50,12 @@ submodule(homogenization) mechanical
end interface
type :: tOutput !< requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:) :: output_mechanical
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: &
homogenization_type !< type of each homogenization
@ -60,27 +64,20 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: &
num_homog
class(tNode), pointer :: &
num_homogMech
module subroutine mechanical_init()
print'(/,1x,a)', '<<<+- homogenization:mechanical init -+>>>'
call material_parseHomogenization2()
call parseMechanical()
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity
homogenization_F = homogenization_F0 ! initialize to identity
allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
allocate(homogenization_dPdF(3,3,3,3,discretization_Ncells), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_Ncells)
homogenization_F = homogenization_F0
allocate(homogenization_P(3,3,discretization_Ncells),source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init(num_homogMech)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init()
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init()
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init()
end subroutine mechanical_init
@ -185,8 +182,10 @@ module subroutine mechanical_results(group_base,ho)
character(len=*), intent(in) :: group_base
integer, intent(in) :: ho
integer :: ou
character(len=:), allocatable :: group
group = trim(group_base)//'/mechanical'
call results_closeGroup(results_addGroup(group))
@ -197,12 +196,17 @@ module subroutine mechanical_results(group_base,ho)
end select
!temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1')
!temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa')
do ou = 1, size(output_mechanical(1)%label)
select case (output_mechanical(ho)%label(ou))
case('F')
call results_writeDataset(reshape(homogenization_F,[3,3,discretization_nCells]),group,'F', &
'deformation gradient','1')
case('P')
call results_writeDataset(reshape(homogenization_P,[3,3,discretization_nCells]),group,'P', &
'deformation gradient','1')
end select
end do
end subroutine mechanical_results
@ -210,35 +214,42 @@ end subroutine mechanical_results
!--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration
!--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization2()
subroutine parseMechanical()
class(tNode), pointer :: &
material_homogenization, &
homog, &
homogMech
mechanical
integer :: ho
integer :: h
material_homogenization => config_material%get('homogenization')
allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(output_mechanical(size(material_name_homogenization)))
do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h)
homogMech => homog%get('mechanical')
select case (homogMech%get_asString('type'))
do ho=1, size(material_name_homogenization)
homog => material_homogenization%get(ho)
mechanical => homog%get('mechanical')
#if defined(__GFORTRAN__)
output_mechanical(ho)%label = output_as1dString(mechanical)
#else
output_mechanical(ho)%label = mechanical%get_as1dString('output',defaultVal=emptyStringArray)
#endif
select case (mechanical%get_asString('type'))
case('pass')
homogenization_type(h) = HOMOGENIZATION_NONE_ID
homogenization_type(ho) = HOMOGENIZATION_NONE_ID
case('isostrain')
homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
homogenization_type(ho) = HOMOGENIZATION_ISOSTRAIN_ID
case('RGC')
homogenization_type(h) = HOMOGENIZATION_RGC_ID
homogenization_type(ho) = HOMOGENIZATION_RGC_ID
case default
call IO_error(500,ext_msg=homogMech%get_asString('type'))
call IO_error(500,ext_msg=mechanical%get_asString('type'))
end select
end do
end subroutine material_parseHomogenization2
end subroutine parseMechanical
end submodule mechanical

View File

@ -71,10 +71,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
module subroutine RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data
module subroutine RGC_init()
integer :: &
ho, &
@ -82,6 +79,8 @@ module subroutine RGC_init(num_homogMech)
sizeState, nIntFaceTot
class (tNode), pointer :: &
num_homogenization, &
num_mechanical, &
num_RGC, & ! pointer to RGC numerics data
material_homogenization, &
homog, &
@ -105,7 +104,9 @@ module subroutine RGC_init(num_homogMech)
allocate(state0(material_homogenization%length))
allocate(dependentState(material_homogenization%length))
num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
num_homogenization => config_numerics%get('homogenization',defaultVal=emptyDict)
num_mechanical => num_homogenization%get('mechanical',defaultVal=emptyDict)
num_RGC => num_mechanical%get('RGC',defaultVal=emptyDict)
num%atol = num_RGC%get_asFloat('atol', defaultVal=1.0e+4_pReal)
num%rtol = num_RGC%get_asFloat('rtol', defaultVal=1.0e-3_pReal)
@ -171,8 +172,8 @@ module subroutine RGC_init(num_homogMech)
allocate(homogState(ho)%state0 (sizeState,Nmembers), source=0.0_pReal)
allocate(homogState(ho)%state (sizeState,Nmembers), source=0.0_pReal)
stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:)
stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:)
allocate(dst%volumeDiscrepancy( Nmembers), source=0.0_pReal)
allocate(dst%relaxationRate_avg( Nmembers), source=0.0_pReal)

View File

@ -52,10 +52,11 @@ module subroutine thermal_init()
allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length
allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal)
allocate(current(ho)%T(count(material_homogenizationID==ho)), source=T_ROOM)
allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal)
configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho))
if (configHomogenization%contains('thermal')) then
configHomogenizationThermal => configHomogenization%get('thermal')
#if defined (__GFORTRAN__)
@ -63,13 +64,22 @@ module subroutine thermal_init()
#else
prm%output = configHomogenizationThermal%get_as1dString('output',defaultVal=emptyStringArray)
#endif
select case (configHomogenizationThermal%get_asString('type'))
case ('pass')
call pass_init()
case ('isothermal')
call isotemperature_init()
end select
else
prm%output = emptyStringArray
end if
end associate
end do
call pass_init()
end subroutine thermal_init

View File

@ -1,13 +1,14 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Dummy homogenization scheme for 1 constituent per material point
!> @brief Isotemperature homogenization
!--------------------------------------------------------------------------------------------------
submodule(homogenization:thermal) isotemperature
contains
module subroutine isotemperature_init
module subroutine isotemperature_init()
print'(/,1x,a)', '<<<+- homogenization:thermal:isotemperature init -+>>>'
end subroutine isotemperature_init

View File

@ -7,9 +7,12 @@ submodule(homogenization:thermal) thermal_pass
contains
module subroutine pass_init()
print'(/,1x,a)', '<<<+- homogenization:thermal:pass init -+>>>'
if (homogenization_Nconstituents(1) /= 1) &
call IO_error(211,ext_msg='N_constituents (pass)')
end subroutine pass_init
end submodule thermal_pass

View File

@ -91,6 +91,11 @@ module phase
integer, intent(in) :: ph
end subroutine damage_results
module subroutine thermal_results(group,ph)
character(len=*), intent(in) :: group
integer, intent(in) :: ph
end subroutine thermal_results
module subroutine mechanical_forward()
end subroutine mechanical_forward
@ -487,6 +492,7 @@ subroutine phase_results()
call mechanical_results(group,ph)
call damage_results(group,ph)
call thermal_results(group,ph)
end do

View File

@ -336,7 +336,7 @@ end subroutine damage_results
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!> @brief Constitutive equation for calculating the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
function phase_damage_collectDotState(ph,en) result(broken)

View File

@ -180,11 +180,12 @@ submodule(phase) mechanical
end function elastic_nu
end interface
type :: tOutput !< new requested output (per phase)
type :: tOutput !< requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:) :: output_constituent
type(tOutput), allocatable, dimension(:) :: output_mechanical
procedure(integrateStateFPI), pointer :: integrateState
@ -216,7 +217,7 @@ module subroutine mechanical_init(phases)
print'(/,1x,a)', '<<<+- phase:mechanical init -+>>>'
!-------------------------------------------------------------------------------------------------
allocate(output_constituent(phases%length))
allocate(output_mechanical(phases%length))
allocate(phase_mechanical_Fe(phases%length))
allocate(phase_mechanical_Fi(phases%length))
@ -251,12 +252,12 @@ module subroutine mechanical_init(phases)
allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
phase => phases%get(ph)
mech => phase%get('mechanical')
phase => phases%get(ph)
mech => phase%get('mechanical')
#if defined(__GFORTRAN__)
output_constituent(ph)%label = output_as1dString(mech)
output_mechanical(ph)%label = output_as1dString(mech)
#else
output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
#endif
enddo
@ -330,7 +331,7 @@ module subroutine mechanical_results(group,ph)
integer, intent(in) :: ph
call crystallite_results(group,ph)
call results(group,ph)
select case(phase_plasticity(ph))
@ -879,9 +880,9 @@ end function integrateStateRK
!--------------------------------------------------------------------------------------------------
!> @brief writes crystallite results to HDF5 output file
!> @brief Write mechanical results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results(group,ph)
subroutine results(group,ph)
character(len=*), intent(in) :: group
integer, intent(in) :: ph
@ -891,9 +892,9 @@ subroutine crystallite_results(group,ph)
call results_closeGroup(results_addGroup(group//'/mechanical'))
do ou = 1, size(output_constituent(ph)%label)
do ou = 1, size(output_mechanical(ph)%label)
select case (output_constituent(ph)%label(ou))
select case (output_mechanical(ph)%label(ou))
case('F')
call results_writeDataset(phase_mechanical_F(ph)%data,group//'/mechanical/','F',&
'deformation gradient','1')
@ -919,13 +920,13 @@ subroutine crystallite_results(group,ph)
call results_writeDataset(phase_mechanical_S(ph)%data,group//'/mechanical/','S', &
'second Piola-Kirchhoff stress','Pa')
case('O')
call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_constituent(ph)%label(ou),&
call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_mechanical(ph)%label(ou),&
'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)')
call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_mechanical(ph)%label(ou))
if (any(phase_lattice(ph) == ['hP', 'tI'])) &
call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_mechanical(ph)%label(ou))
end select
enddo
end do
contains
@ -947,7 +948,7 @@ subroutine crystallite_results(group,ph)
end function to_quaternion
end subroutine crystallite_results
end subroutine results
!--------------------------------------------------------------------------------------------------
@ -1335,5 +1336,4 @@ module subroutine phase_set_F(F,co,ce)
end subroutine phase_set_F
end submodule mechanical

View File

@ -6,6 +6,7 @@ submodule(phase) thermal
type :: tThermalParameters
real(pReal) :: C_p = 0.0_pReal !< heat capacity
real(pReal), dimension(3,3) :: K = 0.0_pReal !< thermal conductivity
character(len=pStringLen), allocatable, dimension(:) :: output
end type tThermalParameters
integer, dimension(:), allocatable :: &
@ -108,6 +109,11 @@ module subroutine thermal_init(phases)
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asFloat('K_33')
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
#if defined(__GFORTRAN__)
param(ph)%output = output_as1dString(thermal)
#else
param(ph)%output = thermal%get_as1dString('output',defaultVal=emptyStringArray)
#endif
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
else
@ -381,4 +387,35 @@ function thermal_active(source_label,src_length) result(active_source)
end function thermal_active
!----------------------------------------------------------------------------------------------
!< @brief writes damage sources results to HDF5 output file
!----------------------------------------------------------------------------------------------
module subroutine thermal_results(group,ph)
character(len=*), intent(in) :: group
integer, intent(in) :: ph
integer :: ou
if (allocated(param(ph)%output)) then
call results_closeGroup(results_addGroup(group//'thermal'))
else
return
endif
do ou = 1, size(param(ph)%output)
select case(trim(param(ph)%output(ou)))
case ('T')
call results_writeDataset(current(ph)%T,group//'thermal','T', 'temperature','T')
end select
end do
end subroutine thermal_results
end submodule thermal