Algebraic flux-corrected transport (FCT) for scalar advection–diffusion on the finite-volume (FV) hex shell grid (fusing two wedges into one hex cell). See fct_advection_diffusion.hpp.
We discretise the scalar transport equation in conservative form
\[ \frac{\partial T}{\partial t} + \nabla \cdot (\mathbf{u}\, T) - \nabla \cdot (\kappa\, \nabla T) = f. \]
Here \(T\) is the transported scalar (e.g. temperature), \(\mathbf{u}\) the given velocity field, \(\kappa \geq 0\) the physical diffusivity, and \(f\) an optional volumetric source term [T/time] (e.g. radiogenic heat production).
The conservative form is used throughout — in the PDE, in the derivation, and in the implementation — because it maps directly onto face fluxes and the finite-volume divergence theorem. The discrete analogue of \(\nabla\cdot(\mathbf{u}T)\) is \(\sum_j \beta_{ij}\,T_{ij}^{\mathrm{face}}\), where the face flux \(\beta_{ij} = \int_{\partial K_{ij}} \mathbf{u}\cdot\hat{n}\,\mathrm{d}S\) is evaluated by numerical quadrature.
On divergence-free vs. compressible velocity fields. Expanding the conservative flux: \(\nabla\cdot(\mathbf{u}T) = \mathbf{u}\cdot\nabla T + T\,\nabla\cdot\mathbf{u}\). For a passive scalar such as temperature the physically correct equation is the advective form \(\partial_t T + \mathbf{u}\cdot\nabla T = \ldots\) — a fluid parcel's temperature changes only through diffusion and sources, not compression. The two forms are equivalent only when \(\nabla\cdot\mathbf{u} = 0\) (Boussinesq).
The subtract_divergence flag (available on all step functions and on UnsteadyAdvectionDiffusion) subtracts the discrete divergence \((\sum_j \beta_{ij})/M_{ii}\) from the upwind operator, recovering the advective form. It defaults to true because:
Additional source terms for compressible flow. Subtracting the divergence error is a necessary but not sufficient correction for compressible or anelastic velocity fields. The full temperature (or entropy) equation under compressible flow contains further source terms — most notably adiabatic heating/cooling (proportional to \(Dp/Dt\), i.e. depending on pressure and velocity) and viscous dissipation. These are problem-specific and their derivation is left to the user; they can be supplied through the source parameter as a pre-computed volumetric rate \(f\) [T/time].
Dirichlet boundary conditions \(T = T_{\mathrm{bc}}\) at the CMB ( \(r = r_{\min}\)) and the outer surface ( \(r = r_{\max}\)) are enforced strongly (pointwise) at the end of each time step by overwriting the boundary cell values with the prescribed constants \(T_{\mathrm{cmb}}\) and \(T_{\mathrm{surface}}\) (see DirichletBCs in helpers.hpp). This is first-order in time: the boundary value is fixed after each step and acts as a Dirichlet datum for the stencil in the next step. No lateral in/outflow boundary conditions are needed since flow does not cross the radial shell boundaries by construction.
Integrate the conservative PDE over a single cell \(K_i\):
\[ \frac{\mathrm{d}}{\mathrm{d}t}\int_{K_i} T\,\mathrm{d}x + \int_{K_i} \nabla\cdot(\mathbf{u}T)\,\mathrm{d}x - \int_{K_i} \nabla\cdot(\kappa\nabla T)\,\mathrm{d}x = \int_{K_i} f\,\mathrm{d}x. \]
Step 1 — Cell average. Represent \(T\) as piecewise constant: \(T|_{K_i} \approx T_i\). Then \(\int_{K_i} T\,\mathrm{d}x = |K_i|\,T_i =: M_{ii}\,T_i\) and \(\int_{K_i} f\,\mathrm{d}x \approx M_{ii}\,f_i\).
Step 2 — Divergence theorem on the advective term \(\nabla\cdot(\mathbf{u}T)\). The volume integral becomes a surface integral, split into one contribution per face:
\[ \int_{K_i} \nabla\cdot(\mathbf{u}T)\,\mathrm{d}x = \oint_{\partial K_i} T\,(\mathbf{u}\cdot\hat{n})\,\mathrm{d}S \approx \sum_{j\in\mathcal{N}(i)} T_{ij}^{\mathrm{face}}\,\beta_{ij}, \]
where \(\beta_{ij} = \int_{\partial K_i \cap \partial K_j} \mathbf{u}\cdot\hat{n}\,\mathrm{d}S\) (positive = net outflow from \(K_i\)) is computed by numerical quadrature over the wedge faces (see GeometryHelper::compute_geometry), and \(T_{ij}^{\mathrm{face}}\) is the face value chosen by the reconstruction (Step 5).
Step 3 — Divergence theorem on the diffusive term \(\nabla\cdot(\kappa\nabla T)\).
\[ \int_{K_i} \nabla\cdot(\kappa\nabla T)\,\mathrm{d}x = \oint_{\partial K_i} \kappa\,\nabla T\cdot\hat{n}\,\mathrm{d}S \approx \sum_{j\in\mathcal{N}(i)} d_{ij}\,(T_j - T_i), \]
where the normal gradient across face \((i,j)\) is approximated by a central difference along the cell-to-cell vector. The two-point diffusion coefficient
\[ d_{ij} = \kappa\, \frac{\mathbf{S}_f^{(j)}\cdot\mathbf{S}_f^{(j)}} {(\mathbf{x}_j - \mathbf{x}_i)\cdot\mathbf{S}_f^{(j)}} \]
accounts for mesh non-orthogonality via the area-weighted face normal \(\mathbf{S}_f^{(j)} = \int_{\partial K_{ij}} \hat{n}_{ij}\,\mathrm{d}S\). It is exact on orthogonal meshes and first-order consistent on smooth non-orthogonal ones.
Step 4 — Collecting terms. Substituting Steps 1–3 gives the semi-discrete ODE for cell \(i\):
\[ M_{ii}\,\dot{T}_i + \underbrace{\sum_j \beta_{ij}\,T_{ij}^{\mathrm{face}}}_{\text{discrete } \nabla\cdot(\mathbf{u}T)} - \underbrace{\sum_j d_{ij}(T_j - T_i)}_{\text{discrete } \nabla\cdot(\kappa\nabla T)} = M_{ii}\,f_i. \]
Forward Euler in time gives the explicit low-order predictor (Stage 1 of FCT); backward Euler gives the semi-implicit system solved with FGMRES.
**Divergence correction (subtract_divergence = true).** With the upwind reconstruction, the advective term expands to \(\sum_j \beta_{ij}^+\,T_i + \sum_j \beta_{ij}^-\,T_j\). Subtracting the discrete divergence \(T_i \sum_j \beta_{ij}\) gives
\[ \sum_j \beta_{ij}^+\,T_i + \sum_j \beta_{ij}^-\,T_j - T_i \underbrace{\left(\sum_j \beta_{ij}^+ + \sum_j \beta_{ij}^-\right)}_{\sum_j \beta_{ij}} = \sum_j \beta_{ij}^-\,(T_j - T_i), \]
which is the discrete advective form \(\mathbf{u}\cdot\nabla T\): only inflow faces contribute, and the difference \(T_j - T_i\) is the upwind gradient. Concretely, the correction reduces to subtracting \(\sum_j \beta_{ij}\) from the diagonal of the upwind operator before the time step — one cheap operation with no additional geometry work, since all \(\beta_{ij}\) are already computed.
Step 5 — Face reconstruction. Two choices for \(T_{ij}^{\mathrm{face}}\):
FCT combines both: the upwind step guarantees monotonicity; the antidiffusive flux \(f_{ij} = \frac{|\beta_{ij}|}{2}(T_i - T_j)\) — the difference between the central and upwind face values, multiplied by the face flux magnitude — is added back in a limited fashion to recover accuracy where the solution is smooth.
The fundamental challenge in scalar transport is the tension between two conflicting goals:
FCT resolves this conflict by treating the two schemes as complementary:
The result is a scheme that is as close to second-order as the local solution smoothness allows, while being provably free of spurious oscillations.
Physical diffusion ( \(\kappa > 0\)) is treated differently: it is built into the low-order predictor (Stage 1) and is not FCT-limited. Physical diffusion is already a smoothing process; limiting it would reduce it below the physical value, which is incorrect. Only the purely numerical antidiffusive flux (the gap between central and upwind advection) is subject to FCT limiting.
Source terms ( \(f \neq 0\)) are added directly to the low-order predictor as \(+\Delta t\,f_i\). They are not FCT-limited, consistent with the physical interpretation that the source genuinely changes the local \(T\) value.
Each cell \(K_i\) carries a cell-averaged value \(T_i\). The semi-discrete system is
\[ M_{ii}\,\dot{T}_i + \sum_j \beta_{ij}\, T_j^{\mathrm{upw}} - \sum_j d_{ij}\,(T_j - T_i) = M_{ii}\,f_i, \]
where
All geometric quantities are computed by GeometryHelper::compute_geometry (see geometry_helper.hpp).
Each explicit time step consists of four stages.
Default (subtract_divergence = true): advective upwind form.
\[ T_i^L = T_i^n - \frac{\Delta t}{M_{ii}} \Bigl[ \underbrace{\sum_j \beta_{ij}^-\,(T_j^n - T_i^n)}_{\text{upwind, advective}} + \underbrace{\sum_j d_{ij}\,(T_i^n - T_j^n)}_{\text{physical diffusion}} \Bigr] + \Delta t\, f_i. \]
With subtract_divergence = false: conservative form (only needed if \(\nabla\cdot(\mathbf{u}T)\) is genuinely the intended equation rather than \(\mathbf{u}\cdot\nabla T\)):
\[ T_i^L = T_i^n - \frac{\Delta t}{M_{ii}} \Bigl[ \underbrace{\beta_{ii}^+\, T_i^n + \sum_{j \neq i} \beta_{ij}^-\, T_j^n}_{\text{upwind, conservative}} + \underbrace{\sum_j d_{ij}\,(T_i^n - T_j^n)}_{\text{physical diffusion}} \Bigr] + \Delta t\, f_i. \]
Both forms use \(\beta^+ = \max(\beta, 0)\), \(\beta^- = \min(\beta, 0)\). Both are monotone (LED) under the CFL condition below. When no source is provided the last term vanishes.
The antidiffusive flux on face \((i,j)\) is the difference between what a central scheme and the upwind scheme would have transported:
\[ f_{ij} = \frac{|\beta_{ij}|}{2}\,(T_i^n - T_j^n). \]
Fluxes are stored pre-scaled as \(\tilde{f}_{ij} = \frac{\Delta t}{M_{ii}}\,f_{ij}\). A positive \(f_{ij}\) means the central scheme would have pushed more of \(T_i^n\) into the face than upwind did (or vice versa for negative).
For each cell \(i\), sum the positive and negative antidiffusive contributions it would receive, and find the local allowable range from the low-order neighbours:
\[ P_i^{\pm} = \sum_{j:\,\pm f_{ij} > 0} f_{ij}, \qquad Q_i^+ = \max_j(T_j^L) - T_i^L \geq 0, \qquad Q_i^- = \min_j(T_j^L) - T_i^L \leq 0. \]
The nodal correction factors cap how much flux can be received before the cell would violate its neighbours' range:
\[ R_i^+ = \min\!\left(1,\,\frac{Q_i^+}{P_i^+}\right), \qquad R_i^- = \min\!\left(1,\,\frac{Q_i^-}{P_i^-}\right) \quad (\text{set to } 1 \text{ if denominator is zero}). \]
For each face, the limiter takes the minimum of what the two cells on either side can handle:
\[ \alpha_{ij} = \begin{cases} \min(R_i^+,\, R_j^-) & \text{if } f_{ij} > 0, \\ \min(R_i^-,\, R_j^+) & \text{if } f_{ij} < 0. \end{cases} \]
The factor \(\alpha_{ij} = 1\) means full second-order correction is safe; \(\alpha_{ij} = 0\) means the flux is entirely suppressed (the solution is locally at an extremum). The final solution
\[ T_i^{n+1} = T_i^L + \sum_j \alpha_{ij}\,\tilde{f}_{ij} \]
satisfies \(\min_j T_j^L \leq T_i^{n+1} \leq \max_j T_j^L\) by construction (LED).
After the FCT correction, Dirichlet boundary conditions are applied by overwriting the boundary radial cells with \(T_{\mathrm{cmb}}\) and \(T_{\mathrm{surface}}\) (see apply_dirichlet_bcs in helpers.hpp).
For large time steps (CFL > 1) the explicit predictor becomes unstable. In the semi-implicit variant the low-order predictor is obtained instead by solving the backward-Euler system
\[ (M + \Delta t\,A_{\mathrm{upw}} + \Delta t\,D)\,T^L = M\,T^n + \Delta t\,M\,f, \]
where \(A_{\mathrm{upw}}\) is the upwind advection–divergence matrix and \(D\) is the diffusion matrix assembled from the \(d_{ij}\) coefficients. When no source is present the RHS reduces to \(M\,T^n\). The solve is performed with FGMRES (see UnsteadyAdvectionDiffusion::compute_rhs and UnsteadyAdvectionDiffusion::apply_impl in advection_diffusion.hpp). The antidiffusive fluxes and Zalesak limiter (Stages 2–4) are identical to the explicit variant.
When subtract_divergence = true is passed to UnsteadyAdvectionDiffusion, the diagonal of \(A_{\mathrm{upw}}\) is reduced by \(\sum_j \beta_{ij}\) per cell, so the implicit system likewise solves the advective form \(\mathbf{u}\cdot\nabla T\) rather than \(\nabla\cdot(\mathbf{u}T)\).
Note on Dirichlet BCs in the implicit solve. The linear system does not enforce boundary conditions; boundary cells are free unknowns during the FGMRES solve and are overwritten only after the FCT correction. This is first-order accurate at the boundary. A more accurate treatment (row replacement in the linear system) is not yet implemented.
The explicit variants require a CFL condition
\[ \mathrm{CFL} = \Delta t\, \frac{\max_i \sum_j \beta_{ij}^+}{M_{ii}} < 1. \]
The semi-implicit variant is unconditionally stable for the low-order predictor; the Zalesak correction cannot introduce new extrema.
Practical time step estimate. The flux-based CFL above requires evaluating all face integrals, which is not convenient at setup time. A safe and practical substitute is
\[ \Delta t < C_{\mathrm{CFL}} \,\frac{h_{\min}}{|\mathbf{u}|_{\max}}, \]
where \(h_{\min}\) is the smallest cell edge length in the mesh and \(|\mathbf{u}|_{\max}\) is the global maximum velocity magnitude. A conservative choice is \(C_{\mathrm{CFL}} \approx 0.1\text{–}0.5\); start at the lower end and increase once the solution behaviour is understood. grid::shell::min_radial_h provides the minimum radial layer thickness, which is typically the smallest length scale on shell meshes. The lateral mesh size at refinement level \(\ell\) scales as \(\pi/(2 \cdot 2^\ell)\) and can be smaller at coarse levels, so verify both directions.
Each hexahedral cell has 6 neighbours, indexed as:
| index | direction |
|---|---|
| 0 | \(x-1\) |
| 1 | \(x+1\) |
| 2 | \(y-1\) |
| 3 | \(y+1\) |
| 4 | \(r-1\) (inner shell) |
| 5 | \(r+1\) (outer shell) |
Use fct_explicit_step for production runs. Pure upwind (upwind_explicit_step) is provided as a debugging baseline and for verifying the geometry and time-step choice, but it is very diffusive: sharp temperature gradients smear out quickly and the thermal structure is lost long before physical equilibrium is reached. FCT recovers near-second-order accuracy in smooth regions while remaining oscillation-free — use it.
The semi-implicit variant (fct_semiimplicit_step + FGMRES) removes the CFL restriction but requires solving a linear system each step; at typical mantle convection resolutions the explicit CFL limit is not prohibitive and the explicit FCT step is simpler, cheaper per step, and free of solver convergence issues. Prefer explicit unless very large time steps are needed.
With a volumetric heat source q [T/time]:
The result will be noticeably more diffusive than FCT — use it to verify that the geometry and time-step logic are correct, then switch to fct_explicit_step.
1.9.8 