Temperature Module (USE_TEMP)

ExaGOOP includes an optional coupled heat-transfer module that solves a Lagrangian energy equation alongside the momentum equation. When enabled, each material point carries thermal state (temperature, specific heat, thermal conductivity, heat flux, and internal heat source), and the nodal grid carries five additional fields that are deposited and interpolated every time step.

The module is compiled in or out at build time through the USE_TEMP preprocessor flag. It is enabled by default in both the CMake and GNUmake build systems.

Enabling the module at build time

CMake

cmake -DEXAGOOP_USE_TEMP=ON  ..   # ON  → compiles temperature module (default)
cmake -DEXAGOOP_USE_TEMP=OFF ..   # OFF → omits temperature module

The option defaults to ON. When ON, the preprocessor symbol USE_TEMP=1 is injected; when OFF, USE_TEMP=0 is injected.

GNUmake

In the relevant GNUmakefile (e.g., Build_Gnumake/GNUmakefile or a test-specific makefile), set the flag before including Make.defs:

USE_TEMP = TRUE   # enables temperature module  (default in build system)
USE_TEMP = FALSE  # disables temperature module

The macro -DUSE_TEMP=$(USE_TEMP) is passed to the compiler automatically by Build_Gnumake/GNUmakefile.

Note

TRUE and FALSE are defined by the build system as 1 and 0 respectively, so #if USE_TEMP guards in the source work correctly for both build paths.

Autogeneration parameters

When mpm.use_autogen = 1 and the temperature module is compiled in, four additional material properties must be supplied for the autogenerated particle block. All four are optional and fall back to the defaults listed below.

Parameter

Type

Default

Description

mpm.T_autogen

Real

0.001

Initial temperature \(T_p^0\) of all autogenerated particles.

mpm.cp_autogen

Real

0.001

Specific heat capacity \(c_p\) (energy per unit mass per degree).

mpm.thermcond_autogen

Real

0.001

Thermal conductivity \(\kappa\) (treated as a scalar; units consistent with the rest of the input).

mpm.heatsrc_autogen

Real

0.001

Volumetric internal heat source \(\dot{q}\).

When particles are loaded from an external file (mpm.particle_file), each particle line must include the four thermal fields in the order: temperature, specific_heat, thermal_conductivity, heat_source. The heat-flux vector components are initialised to zero by the reader regardless of what the file contains.

Temperature boundary conditions

Temperature BCs follow the same per-face naming convention as the momentum BCs. Each domain face is identified by an axis prefix (x, y, z) and a side suffix (lo, hi), giving keys of the form mpm.bc_xlo_temp, mpm.bc_xhi_temp, mpm.bc_ylo_temp, etc.

Parameter

Type

Default

Description

mpm.bc_<face>_temp

string

dirichlet (non-periodic faces); periodic (periodic faces)

BC type for the named face (see keyword table below).

mpm.bc_<face>_temp.T_wall

Real

0.0

Wall temperature \(T_w\) for dirichlet faces.

mpm.bc_<face>_temp.flux

Real

0.0

Prescribed temperature gradient \(\partial T / \partial n\) for heatflux faces. For a dimensional heat flux \(q\), set flux = \(-q / k\) where \(k\) is the material thermal conductivity.

mpm.bc_<face>_temp.h

Real

0.0

Heat transfer coefficient \(h\) for convective faces.

mpm.bc_<face>_temp.T_inf

Real

0.0

Ambient temperature \(T_\infty\) for convective faces.

The five recognised BC type keywords are:

Keyword

Meaning

periodic

Periodic face. Automatically set for directions where mpm.is_it_periodic is 1; no sub-parameters needed.

dirichlet

Fixed-temperature wall. The nodal temperature is overwritten each step with mpm.bc_<face>_temp.T_wall.

adiabatic

Zero heat-flux (insulated) wall. The boundary node temperature is set equal to the first interior neighbour, enforcing \(\partial T / \partial n = 0\). No sub-parameters needed.

heatflux

Prescribed heat-flux wall. The boundary node temperature is set using a ghost-point formula: \(T_b = T_\text{int} + \texttt{flux} \cdot \Delta x\), where \(T_\text{int}\) is the first interior neighbour and \(\Delta x\) is the cell width in the face-normal direction.

convective

Convective (Robin) wall. The boundary temperature is set via a Biot-number weighting: \(T_b = (T_\text{int} + \text{Bi} \cdot T_\infty) / (1 + \text{Bi})\), where \(\text{Bi} = h \cdot \Delta x\).

Nodal fields added by the temperature module

When USE_TEMP=1, the nodal MultiFab grows from 20 to 25 components (NUM_STATES = 25). The five additional components, their integer indices, and their physical meaning are:

Index

Name (constants.H macro)

Physical quantity

20

MASS_SPHEAT

\(\sum_p m_p c_p N_I(\mathbf{x}_p)\) — thermal mass at node \(I\), used as the denominator when normalising the nodal temperature.

21

MASS_SPHEAT_TEMP

\(\sum_p m_p c_p T_p N_I(\mathbf{x}_p)\) — mass-weighted temperature sum; divided by MASS_SPHEAT to obtain the nodal temperature after deposition.

22

TEMPERATURE

Nodal temperature \(T_I\) after normalisation and boundary condition enforcement.

23

SOURCE_TEMP_INDEX

Net nodal heat source \(S_I = \sum_p (-\nabla N_I \cdot \mathbf{q}_p + \dot{q}_p N_I) V_p\), accumulated during the particle-to-grid pass.

24

DELTA_TEMPERATURE

\(\Delta T_I = T_I^{t+\Delta t} - T_I^t\), computed after the nodal temperature advance and used during the grid-to-particle interpolation step (analogous to DELTA_VEL in the momentum update).

These field names are visible in ParaView when loading grid output files.

Per-particle thermal fields

The following real-data slots are appended to the particle struct when USE_TEMP=1 (indices continue from the non-thermal fields; they are not contiguous with the non-thermal fields in memory — see realData enum in mpm_specs.H for absolute indices):

Field (realData name)

Description

temperature

Current particle temperature \(T_p\).

specific_heat

Specific heat \(c_p\); constant per particle throughout the simulation.

thermal_conductivity

Scalar thermal conductivity \(\kappa\); constant per particle.

heat_flux

Heat-flux vector \(\mathbf{q}_p\) (SPACEDIM components); updated during the grid-to-particle pass by interpolate_from_grid_temperature.

heat_source

Volumetric internal heat source \(\dot{q}_p\); constant per particle.

Time-step sequence for the thermal solve

Each time step, the thermal module executes the following operations inside the main integration loop (after the momentum update):

  1. Particle-to-grid (P2G) — mass and source deposition deposit_onto_grid_temperature() scatters MASS_SPHEAT, MASS_SPHEAT_TEMP, and SOURCE_TEMP_INDEX from particles to grid nodes.

  2. Nodal temperature normalisation For nodes where MASS_SPHEAT ≥ mass_tolerance, the nodal temperature is computed as:

    \[T_I = \frac{\sum_p m_p c_p T_p N_I(\mathbf{x}_p)} {\sum_p m_p c_p N_I(\mathbf{x}_p)}\]

  3. Backup nodal temperature — backup_current_temperature() copies TEMPERATURE into DELTA_TEMPERATURE for use in step 6.

  4. Boundary condition enforcement — nodal_bcs_temperature() applies Dirichlet values from bc_lower_tempval / bc_upper_tempval at the domain boundary nodes.

  5. Nodal temperature advance — advance_nodal_temperature() updates each node:

    \[T_I^{t+\Delta t} = T_I^t + \frac{S_I}{\sum_p m_p c_p N_I} \, \Delta t\]

  6. Store \(\Delta T\) — store_delta_temperature() computes DELTA_TEMPERATURE = T_I^{new} − T_I^{old}.

  7. Grid-to-particle (G2P) — temperature and heat-flux interpolation interpolate_from_grid_temperature() interpolates DELTA_TEMPERATURE back to each particle to update \(T_p\), and computes the new particle heat-flux vector \(\mathbf{q}_p = -\kappa_p \nabla T(\mathbf{x}_p)\) by interpolating the gradient of TEMPERATURE.

Minimal input file snippet

The block below shows the additional lines required on top of the standard input file to activate a uniform-temperature elastic-solid simulation with fixed temperatures at the lower and upper \(y\)-faces:

# -------------------------------------------------------------------
# Autogeneration: thermal properties
# (only needed when mpm.use_autogen = 1)
# -------------------------------------------------------------------
mpm.T_autogen        = 300.0     # initial temperature (K)
mpm.cp_autogen       = 500.0     # specific heat  (J / kg K)
mpm.thermcond_autogen = 50.0     # thermal conductivity  (W / m K)
mpm.heatsrc_autogen  = 0.0       # volumetric heat source (W / m^3)

# -------------------------------------------------------------------
# Temperature boundary conditions
# -------------------------------------------------------------------
mpm.bc_ylo_temp          = dirichlet
mpm.bc_ylo_temp.T_wall   = 400.0

mpm.bc_yhi_temp          = dirichlet
mpm.bc_yhi_temp.T_wall   = 300.0

In this example the \(x\)- and \(z\)-faces are left at their defaults (dirichlet with T_wall = 0.0 for non-periodic directions), while the lower \(y\)-face is held at 400 K and the upper \(y\)-face at 300 K, imposing a temperature gradient in the \(y\)-direction.

Warning

The default values for bc_lower_tempval and bc_upper_tempval are 0.0 in all directions. Any face for which these are not explicitly set in the input file will have its boundary nodes forced to \(T = 0\) each step. Always specify all six values (three lower, three upper) when running with USE_TEMP=1.

Verification test cases

The following two test cases verify the heatflux and convective temperature boundary conditions against analytical solutions. Both are located under Tests/ and follow the standard ExaGOOP layout: a PreProcess/ directory containing config.json and Generate_MPs_Inputfile_Generic.py, a top-level input file (Inputs_*.inp), a particle file (mpm_particles.dat), and a PostProcess/ directory containing validate.py.

The thermal diffusivity for both cases is \(\alpha = k / (\rho c_p) = 1 / (1 \times 1) = 1\), which gives unit-diffusivity heat conduction — a clean choice that makes the analytical series converge quickly.

1D heat conduction with prescribed heat flux

Physical problem

A stationary elastic slab occupies \(x \in [0, L]\) with \(L = 1\). The left face is held at a fixed temperature, and the right face has a prescribed outward heat flux. The problem is effectively one-dimensional: the \(y\)-faces are periodic and the domain is thin (\(L_y = 0.1\)) so there is no variation in \(y\).

Governing equation:

\[\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}, \qquad \alpha = \frac{k}{\rho c_p}\]

Boundary and initial conditions:

\[\begin{split}T(0,\, t) &= 0 \quad &&\text{(Dirichlet)} \\ \frac{\partial T}{\partial x}\bigg|_{x=L} &= \frac{q}{k} = 1 \quad &&\text{(prescribed flux, } q = 1,\; k = 1\text{)} \\ T(x,\, 0) &= 0 \quad &&\text{(initial condition)}\end{split}\]

Material properties

Property

Value

Notes

Density \(\rho\)

1.0

uniform

Specific heat \(c_p\)

1.0

uniform

Thermal conductivity \(k\)

1.0

uniform, giving \(\alpha = 1\)

Internal heat source \(\dot{q}\)

0.0

no volumetric source

Initial temperature \(T_0\)

0.0

uniform

Domain and discretisation

The computational domain is \([0,1] \times [0,0.1]\) with \(40 \times 4\) background grid cells. Two material points per cell in each direction gives \(80 \times 8 = 640\) material points in total. The time step is fixed at \(\Delta t = 10^{-5}\) and the simulation is advanced to \(t_\text{final} = 2.0\).

Analytical solution

Decomposing \(T = T_s(x) + v(x,t)\) where \(T_s = (q/k)\,x\) is the steady state, the transient part satisfies the heat equation with a Dirichlet condition at \(x = 0\) and a zero-flux condition at \(x = L\). The eigenvalues of this Sturm-Liouville problem are \(\lambda_n = (2n-1)\pi / (2L),\; n = 1,2,\ldots\), giving:

\[T(x, t) = \frac{q}{k}\,x \;+\; \sum_{n=1}^{\infty} \frac{2\,(-1)^n}{L\,\lambda_n^2} \sin(\lambda_n x)\, e^{-\lambda_n^2 \alpha t}\]

At long times the solution approaches the linear steady state \(T_s(x) = x\).

Setting up the test case

The particle file and input file are generated from the configuration in PreProcess/config.json. The key editable fields are:

{
  "grid": { "nx": 40, "ny": 4 },
  "ppc": [2, 2],
  "bodies": [{
    "temperature": {
      "T": 0.0,
      "spheat": 1.0,
      "thermcond": 1.0,
      "heatsrc": 0.0
    }
  }]
}

To regenerate the particle and input files, run from the test directory:

cd Tests/1D_Heat_Conduction_HeatFlux/
bash Generate_MPs_and_InputFiles.sh

This writes mpm_particles.dat and Inputs_1DHeatConduction_HeatFlux.inp. The boundary conditions in the generated input file are:

mpm.bc_xlo_temp          = dirichlet
mpm.bc_xlo_temp.T_wall   = 0.0
mpm.bc_xhi_temp          = heatflux
mpm.bc_xhi_temp.flux     = 1.0

To change the imposed flux, set mpm.bc_xhi_temp.flux to the desired value of \(q\). If the material thermal conductivity is also changed (via thermcond in config.json), the flux entry in the input file should be updated to \(q / k\) accordingly.

Running

Build the binary with the temperature module enabled:

cd Tests/1D_Heat_Conduction_HeatFlux
# copy make file here.
make USE_TEMP=TRUE -j4

Then run:

./ExaGOOP2d.<suffix>.ex Inputs_1DHeatConduction_HeatFlux.inp

Output particle snapshots are written to the subdirectory specified by mpm.prefix_asciifilename inside the test directory.

Post-processing and validation

The validation script compares the numerical temperature profile against the analytical Fourier-series solution at \(t = 0.5\):

cd Tests/1D_Heat_Conduction_HeatFlux/PostProcess
python validate.py

A passing run prints the RMS error and PASS. The acceptance criterion is RMS \(< 10^{-2}\).

Sample result

At \(t = 0.5\) the transient has partially decayed. The profile rises from \(T(0) = 0\) and curves upward, approaching the linear steady state \(T = x\). The first-mode relaxation time is \(1/\lambda_1^2 = 4/\pi^2 \approx 0.41\), so at \(t = 0.5\) roughly \(e^{-1.23} \approx 29\%\) of the leading-mode amplitude remains. Representative exact values are:

\(x\)

\(T_\text{exact}(x,\,0.5)\)

Steady state \(x\)

0.00

0.000

0.000

0.25

0.160

0.250

0.50

0.333

0.500

0.75

0.532

0.750

1.00

0.764

1.000

The numerical solution produced by ExaGOOP matches these values to within RMS \(\approx 3 \times 10^{-3}\).

1D heat conduction with convective boundary condition

Physical problem

The same slab geometry and material as the heat-flux case, but with different boundary conditions: the left face has a fixed temperature and the right face is subject to Newton cooling (a Robin condition).

Boundary and initial conditions:

\[\begin{split}T(0,\, t) &= T_w = 1 \quad &&\text{(Dirichlet)} \\ -k\,\frac{\partial T}{\partial x}\bigg|_{x=L} &= h\bigl(T(L,t) - T_\infty\bigr) \quad &&\text{(convective, } h = 2,\; T_\infty = 0\text{)} \\ T(x,\, 0) &= 0 \quad &&\text{(initial condition)}\end{split}\]

This models a slab heated from one end by a fixed wall and cooled from the other by a convective fluid. The Biot number is \(\text{Bi} = h L / k = 2\), placing this case firmly in the mixed-convection regime (neither thermally thin nor thick limit).

Material properties

Identical to the heat-flux case: \(\rho = c_p = k = 1\), \(\dot{q} = 0\), \(T_0 = 0\).

Analytical solution

The steady state satisfying both boundary conditions is:

\[T_s(x) = T_w + \frac{(T_\infty - T_w)\,(h/k)\,x}{1 + h L / k} = 1 - \frac{2x}{3}\]

Decomposing \(T = T_s + w\), the transient part \(w\) satisfies homogeneous boundary conditions of the same type (Dirichlet at \(x = 0\), Robin at \(x = L\)). The eigenfunctions are \(\sin(\lambda_n x)\) where the eigenvalues satisfy the transcendental equation:

\[\lambda_n \cot(\lambda_n L) = -\frac{h}{k}\]

For \(h/k = 2\) the first few roots are approximately \(\lambda_1 \approx 1.077\), \(\lambda_2 \approx 3.644\), \(\lambda_3 \approx 6.578\), … . The full solution is:

\[T(x, t) = T_s(x) \;+\; \sum_{n=1}^{\infty} C_n \sin(\lambda_n x)\,e^{-\lambda_n^2 \alpha t}, \qquad C_n = \frac{\displaystyle\int_0^L -T_s(\xi)\sin(\lambda_n\xi)\,\mathrm{d}\xi} {\displaystyle\int_0^L \sin^2(\lambda_n\xi)\,\mathrm{d}\xi}\]

Setting up the test case

The config.json structure is identical to the heat-flux case. The key boundary condition lines in the generated input file are:

mpm.bc_xlo_temp          = dirichlet
mpm.bc_xlo_temp.T_wall   = 1.0
mpm.bc_xhi_temp          = convective
mpm.bc_xhi_temp.h        = 2.0
mpm.bc_xhi_temp.T_inf    = 0.0

To change the heat transfer coefficient, edit mpm.bc_xhi_temp.h. To change the far-field temperature, edit mpm.bc_xhi_temp.T_inf. Regenerate input and particle files after editing config.json:

cd Tests/1D_Heat_Conduction_Convective/PreProcess
python Generate_MPs_Inputfile_Generic.py

Running

cd Tests/1D_Heat_Conduction_Convective
make USE_TEMP=TRUE -j4
./ExaGOOP2d.gnu.MPI.ex Inputs_1DHeatConduction_Convective.inp

Post-processing and validation

cd Tests/1D_Heat_Conduction_Convective/PostProcess
python validate.py

The script finds eigenvalues numerically (scipy.optimize.brentq), evaluates the Fourier series at \(t = 0.5\), and reports the RMS error. Acceptance criterion: RMS \(< 10^{-2}\).

Sample result

The steady state \(T_s(x) = 1 - 2x/3\) drops from 1 at the heated wall to \(1/3 \approx 0.333\) at the convective boundary. The first eigenvalue \(\lambda_1 \approx 2.289\) gives a relaxation time of \(1/\lambda_1^2 \approx 0.19\), so by \(t = 0.5\) the profile is close to steady state. Representative exact values:

\(x\)

\(T_\text{exact}(x,\,0.5)\)

Steady state \(T_s(x)\)

0.00

1.000

1.000

0.25

0.805

0.833

0.50

0.619

0.667

0.75

0.448

0.500

1.00

0.294

0.333

The numerical solution produced by ExaGOOP matches these values to within RMS \(\approx 2 \times 10^{-3}\).