.. highlight:: rst

.. _temperature_module:

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

.. code-block:: bash

   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``:

.. code-block:: makefile

   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.

.. list-table::
   :header-rows: 1
   :widths: 32 12 14 42

   * - Parameter
     - Type
     - Default
     - Description
   * - ``mpm.T_autogen``
     - ``Real``
     - ``0.001``
     - Initial temperature :math:`T_p^0` of all autogenerated particles.
   * - ``mpm.cp_autogen``
     - ``Real``
     - ``0.001``
     - Specific heat capacity :math:`c_p` (energy per unit mass per
       degree).
   * - ``mpm.thermcond_autogen``
     - ``Real``
     - ``0.001``
     - Thermal conductivity :math:`\kappa` (treated as a scalar; units
       consistent with the rest of the input).
   * - ``mpm.heatsrc_autogen``
     - ``Real``
     - ``0.001``
     - Volumetric internal heat source :math:`\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.

.. list-table::
   :header-rows: 1
   :widths: 38 12 14 36

   * - 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 :math:`T_w` for ``dirichlet`` faces.
   * - ``mpm.bc_<face>_temp.flux``
     - ``Real``
     - ``0.0``
     - Prescribed temperature gradient :math:`\partial T / \partial n` for
       ``heatflux`` faces.  For a dimensional heat flux :math:`q`, set
       ``flux`` = :math:`-q / k` where :math:`k` is the material thermal
       conductivity.
   * - ``mpm.bc_<face>_temp.h``
     - ``Real``
     - ``0.0``
     - Heat transfer coefficient :math:`h` for ``convective`` faces.
   * - ``mpm.bc_<face>_temp.T_inf``
     - ``Real``
     - ``0.0``
     - Ambient temperature :math:`T_\infty` for ``convective`` faces.

The five recognised BC type keywords are:

.. list-table::
   :header-rows: 1
   :widths: 20 80

   * - 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
       :math:`\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:
       :math:`T_b = T_\text{int} + \texttt{flux} \cdot \Delta x`,
       where :math:`T_\text{int}` is the first interior neighbour and
       :math:`\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:
       :math:`T_b = (T_\text{int} + \text{Bi} \cdot T_\infty) / (1 + \text{Bi})`,
       where :math:`\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:

.. list-table::
   :header-rows: 1
   :widths: 16 32 52

   * - Index
     - Name (``constants.H`` macro)
     - Physical quantity
   * - 20
     - ``MASS_SPHEAT``
     - :math:`\sum_p m_p c_p N_I(\mathbf{x}_p)` — thermal mass at node
       :math:`I`, used as the denominator when normalising the nodal
       temperature.
   * - 21
     - ``MASS_SPHEAT_TEMP``
     - :math:`\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 :math:`T_I` after normalisation and boundary
       condition enforcement.
   * - 23
     - ``SOURCE_TEMP_INDEX``
     - Net nodal heat source :math:`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``
     - :math:`\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):

.. list-table::
   :header-rows: 1
   :widths: 30 70

   * - Field (``realData`` name)
     - Description
   * - ``temperature``
     - Current particle temperature :math:`T_p`.
   * - ``specific_heat``
     - Specific heat :math:`c_p`; constant per particle throughout
       the simulation.
   * - ``thermal_conductivity``
     - Scalar thermal conductivity :math:`\kappa`; constant per particle.
   * - ``heat_flux``
     - Heat-flux vector :math:`\mathbf{q}_p` (``SPACEDIM`` components);
       updated during the grid-to-particle pass by
       ``interpolate_from_grid_temperature``.
   * - ``heat_source``
     - Volumetric internal heat source :math:`\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:

   .. math::

      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:

   .. math::

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

6. **Store** :math:`\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 :math:`T_p`, and
   computes the new particle heat-flux vector
   :math:`\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 :math:`y`-faces:

.. code-block:: bash

   # -------------------------------------------------------------------
   # 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 :math:`x`- and :math:`z`-faces are left at their
defaults (``dirichlet`` with ``T_wall = 0.0`` for non-periodic directions),
while the lower :math:`y`-face is held at 400 K and the upper
:math:`y`-face at 300 K, imposing a temperature gradient in the
:math:`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
   :math:`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
:math:`\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
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. _test-heatflux:

**Physical problem**

A stationary elastic slab occupies :math:`x \in [0, L]` with :math:`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 :math:`y`-faces are periodic and the domain is thin
(:math:`L_y = 0.1`) so there is no variation in :math:`y`.

Governing equation:

.. math::

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

Boundary and initial conditions:

.. math::

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

**Material properties**

.. list-table::
   :header-rows: 1
   :widths: 40 20 40

   * - Property
     - Value
     - Notes
   * - Density :math:`\rho`
     - 1.0
     - uniform
   * - Specific heat :math:`c_p`
     - 1.0
     - uniform
   * - Thermal conductivity :math:`k`
     - 1.0
     - uniform, giving :math:`\alpha = 1`
   * - Internal heat source :math:`\dot{q}`
     - 0.0
     - no volumetric source
   * - Initial temperature :math:`T_0`
     - 0.0
     - uniform

**Domain and discretisation**

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

**Analytical solution**

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

.. math::

   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
:math:`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:

.. code-block:: json

   {
     "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:

.. code-block:: bash

   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:

.. code-block:: bash

   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 :math:`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 :math:`q / k` accordingly.

**Running**

Build the binary with the temperature module enabled:

.. code-block:: bash

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

Then run:

.. code-block:: bash

   ./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 :math:`t = 0.5`:

.. code-block:: bash

   cd Tests/1D_Heat_Conduction_HeatFlux/PostProcess
   python validate.py

A passing run prints the RMS error and ``PASS``.  The acceptance
criterion is RMS :math:`< 10^{-2}`.

**Sample result**

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

.. list-table::
   :header-rows: 1
   :widths: 20 30 30

   * - :math:`x`
     - :math:`T_\text{exact}(x,\,0.5)`
     - Steady state :math:`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 :math:`\approx 3 \times 10^{-3}`.


1D heat conduction with convective boundary condition
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. _test-convective:

**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:

.. math::

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

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
:math:`\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: :math:`\rho = c_p = k = 1`,
:math:`\dot{q} = 0`, :math:`T_0 = 0`.

**Analytical solution**

The steady state satisfying both boundary conditions is:

.. math::

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

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

.. math::

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

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

.. math::

   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:

.. code-block:: bash

   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``:

.. code-block:: bash

   cd Tests/1D_Heat_Conduction_Convective/PreProcess
   python Generate_MPs_Inputfile_Generic.py

**Running**

.. code-block:: bash

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

**Post-processing and validation**

.. code-block:: bash

   cd Tests/1D_Heat_Conduction_Convective/PostProcess
   python validate.py

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

**Sample result**

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

.. list-table::
   :header-rows: 1
   :widths: 20 30 30

   * - :math:`x`
     - :math:`T_\text{exact}(x,\,0.5)`
     - Steady state :math:`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 :math:`\approx 2 \times 10^{-3}`.