Skip to content

What's new in Mallard's numerics

Most of Mallard's numerics are published methods, implemented as their papers describe them and cited on the page that describes each one (References). This page says which parts are not: what Mallard adds to the literature, what it combines in a way we did not find published, and what it only implements. It is meant to be conservative: where we were not sure, an item is classed lower.

Classes.

  • Original: we found no published method that does this (the searches are listed in the audit behind this page).
  • Original extension: a published method, extended in a specific way that we did not find published; the row says exactly what is new.
  • Adapted: a published method, with engineering changes that are Mallard's own but not a new method.
  • Standard: a published method, implemented as published.

"Not found" means not found by the searches made (Google Scholar, arXiv, Crossref and journal pages, October 2026); it is not proof that no one has done it. Corrections are welcome as issues.

The items classed original or original extension were developed for Mallard; their class says only whether we found them published.

Summary

Feature Class Closest prior work What Mallard adds Evidence
Central-stencil size chosen by a Lebesgue-constant bound Original Lebesgue constants rating spectral-volume partitions (Wang 2002; Chen 2005); eigenvalue-based reconstruction stabilization (Zangeneh & Ollivier-Gooch 2019) A per-cell, a priori rule: grow the least-squares stencil until the reconstruction's Lebesgue constant at the face points is at most 10 (2D) or 4 (3D) Lebesgue bound
Stencil ranking in a mesh-spacing metric on thin 3D cells Original extension Local tangential-normal frames for least squares on anisotropic meshes (Jalali & Ollivier-Gooch 2013, 2017) Candidates ranked by the unit-determinant inverse of the second moment of the vertex-neighbor offsets, only where that metric's ratio exceeds 2.5 Ranking metric
TENO-E on 3D hexahedra, prisms, pyramids and mixed meshes, orders 3 to 6 Original extension TENO on unstructured meshes (Ji, Liang & Fu 2022: 2D, one 3D tetrahedral test); TENO-E (Liang, Shyy & Fu 2025); WENO on 3D mixed elements (Tsoutsanis et al. 2011) The TENO-E algorithm on every 3D cell type without splitting cells, with the two items above, mirror stencils and curved walls; stability measured on the linearized operator TENO-E in 3D
Species reconstructed with the flow's contact-field stencil choice Original extension Mass-flux upwinding of species (Larrouturou 1991); primitive-variable reconstruction at material interfaces (Johnsen & Colonius 2006) Every scalar takes, per face, the stencil TENO chose for the entropy/contact field (one byte), so all species share linear weights and one bound-preserving factor Species stencils
Double flux with TENO on unstructured meshes Original extension (narrow) Double flux (Abgrall & Karni 2001; Billet & Abgrall 2003); with WENO on structured grids (Houim & Kuo 2011); in DG on unstructured meshes (Lv & Ihme 2014; Johnson & Kercher 2020) The combination with high-order TENO finite volumes on mixed meshes, and wave speeds from each side's own frozen pair so both fluxes carry the same mass flux Double flux
SIMPLER splitting with an explicit integrator: stability polynomial, and the instability of non-dissipative fluxes Original (analysis) SIMPLER (Wu, Ma & Ihme 2019) and balanced splitting (Speth et al. 2013), analysed with exact substeps The composite polynomial with SSPRK3, which exceeds 1 everywhere on the imaginary axis, so a central flux is unstable at any step SIMPLER stability
RHLL blended toward HLL by a pressure-jump weight Original extension (minor) Rotated hybrid (Nishikawa & Kitamura 2008); the weight of AUSMPW+ (Kim, Kim & Rho 2001); pressure-based shock sensing in RoeM (Kim, Kim & Rho 2003); Rotated-RoeM with a Mach-number weight (Choi et al. 2024) The blend of the rotated flux with face-normal HLL by min(p_L/p_R, p_R/p_L)^3, which removes the stagnation-line carbuncle on tetrahedra RHLL blend
Mirror-image stencil entries at straight boundaries Adapted Mirrored ghost cells in high-order unstructured stencils (Yang, Pan & Xu 2023) Coplanar faces as one plane with several conditions, moments reflected degree by degree, images not counted toward the stencil size TENO-E details
Curved boundaries for k-exact reconstruction Adapted Ollivier-Gooch & Van Altena 2002; Krivodonova & Berger 2006 Quadrature area vectors shifted to sum exactly to the face's; degree cap 3 next to curved walls Curved boundaries
Hybrid KEEP / Riemann flux with a Ducros sensor Standard Kuya, Totani & Kawai 2018; Ducros et al. 1999; Bhagatwala & Lele 2009; Ma, Lv & Ihme 2017 Upwind floor and the rule for wall faces LES design
Low-Mach correction Adapted Thornber et al. 2008; cutoff as Weiss & Smith 1995 Cutoff Mach number Overview
Dynamic thickened-flame efficiency Adapted Charlette et al. 2002 (II); Wang, Boileau & Veynante 2011; Volpiani et al. 2016 Test filter as repeated vertex-neighbor averages sized by their covariance; exact global sums LES design
Conservative remap on an exact supermesh Adapted Farrell & Maddison 2011; Menon & Schmidt 2011 Intersections as lists of clipped simplices; results bitwise independent of both partitions Remap design
Results bitwise independent of rank and thread counts Adapted Order-invariant integer sums (Hallberg & Adcroft 2014); reproducible summation (Demmel & Nguyen 2015) Applied throughout the solver: per-face fluxes gathered in a fixed order, fixed-point global sums, partition-independent setup and restarts Performance
Wall-model exchange point located by a walk Standard Lawson 1977; Devillers, Pion & Teillaud 2002; exchange height of Kawai & Larsson 2012 Wall models

The TENO-E setup on the GPU, its incremental rebuild and its cache are implementation, not numerics, and are not counted here (performance).

Original and extended items

Lebesgue bound

Mallard fits each cell's central polynomial by least squares over a stencil of nearby cells (Barth & Frederickson 1990), as TENO on unstructured meshes does with a fixed size of about twice the number of coefficients (Ji, Liang & Fu 2022). A fixed size is not enough: a stencil can have full rank yet resolve one direction only through small centroid offsets, and then amplify the cell averages, and the truncation error with them. Mallard grows the stencil one entry at a time and takes the first size whose Lebesgue constant at the cell's face quadrature points,

\[ \begin{equation} \Lambda = \max_q\left(\Bigl\lvert 1 - \sum_s c_{qs}\Bigr\rvert + \sum_s\lvert c_{qs}\rvert\right), \label{eq:lebesgue} \end{equation} \]

is at most 10 (2D) or 4 (3D), where \(c_{qs}\) is the weight of cell \(s\)'s average in the face value at point \(q\) (TENO-E details). In 3D, if no size reaches the bound, it takes the smallest within 1.25 times the best found.

Lebesgue constants have rated reconstruction operators before, to choose the control-volume partitions of the spectral volume method (Wang 2002; Chen 2005), and eigenvalue analysis has been used to modify reconstructions until the discrete operator is stable (Zangeneh & Ollivier-Gooch 2019). We did not find the Lebesgue constant used to choose the size of a reconstruction stencil cell by cell. The bounds themselves are empirical.

Evidence. The largest real part of the eigenvalues of the linearized operator (acoustics at rest, HLLC), with the bound of 10 and of 4: 1.9 to 3.5 against below \(10^{-7}\) on prism meshes with walls, seams or jitter, and 1.1 against \(8\times 10^{-8}\) on 25:1 Kuhn tetrahedra at order 4. Bounds of 5 to 7 do not remove the growing modes. Lowering the bound kept the measured convergence rates of every cell type and order (3D stability). Known exception: pure pyramid lattices at order 4 are less stable with the bound of 4.

Ranking metric

On thin cells, ranking candidates by physical distance takes the whole wall-normal column first. Mallard ranks them instead in the metric of the local mesh spacing,

\[ \begin{equation} \lVert\mathbf{d}\rVert_\mathbf{M}^2 = \mathbf{d}^T\mathbf{M}\,\mathbf{d}, \qquad \mathbf{M} = (\det\mathbf{C})^{1/3}\,\mathbf{C}^{-1}, \qquad \mathbf{C} = \sum_j \mathbf{d}_j\mathbf{d}_j^T, \label{eq:metric} \end{equation} \]

with \(\mathbf{d}_j\) the offsets of the cell's vertex neighbors (mirrored at boundary faces), on cells where \(\sqrt{\lambda_{\max}/\lambda_{\min}}\) of \(\mathbf{C}\) exceeds 2.5, a value above what cell shapes alone give on regular tilings (ranking metric). Least squares on anisotropic meshes has been improved by fitting in a local tangential-normal frame (Jalali & Ollivier-Gooch 2013, 2017); we did not find a metric of this kind used to rank stencil candidates.

Evidence. Without the metric, thin boxes have eigenvalues of the linearized operator up to \(+13\); with it, none above zero. On the examples/sphere_re300 mesh (prism layers of aspect ratio about 7), TENO3 once reached Mach 0.9 at \(t = 1\); it now stays at or below 0.30 to \(t = 10\).

TENO-E in 3D

TENO on unstructured meshes was introduced in 2D, with one 3D test on tetrahedra, other cell types split into triangles or tetrahedra (Ji, Liang & Fu 2022), and TENO-E adds the adaptive dissipation that Mallard uses (Liang, Shyy & Fu 2025). Mallard runs TENO-E on tetrahedra, hexahedra, prisms, pyramids and mixed meshes without splitting the cells, at orders 3 to 6, with cell moments from one tetrahedral decomposition per cell and the stencil rules above. WENO-type schemes on 3D mixed elements are established (Tsoutsanis, Titarev & Drikakis 2011); what is new is TENO-E there, with its stability measured and the stencil rules it needed. We could not check the full text of Liang, Shyy & Fu 2025, which may include 3D tests.

Evidence. Measured orders on symmetry-walled boxes: 3.83 and 4.85 for orders 4 and 5 on hexahedra, 2.92, 3.90 and 4.84 for orders 3 to 5 on Kuhn tetrahedra (measured orders); polynomial exactness on prisms, pyramids and mixed meshes in the test suite.

Species stencils

In a gas mixture, TENO selects stencils for the flow variables only. Each species then takes, on each face, the stencil that TENO chose for the entropy/contact characteristic field there, the wave that carries composition, stored as one byte per face side (chemistry design). All species are then linear in the same weights, so their face values sum to one, and one Zhang–Shu factor per cell keeps them in \([0, 1]\). Species reconstruction costs only dot products with weights already computed. The species flux is Larrouturou's (1991). We did not find the stencil choice of one characteristic field reused for passive scalars.

Double flux

Double flux (Abgrall & Karni 2001; Billet & Abgrall 2003) keeps pressure and velocity uniform across contacts between gases of different \(\gamma\). It has been combined with WENO on structured grids (Houim & Kuo 2011) and with DG on unstructured meshes (Lv & Ihme 2014; Johnson & Kercher 2020). Mallard combines it with TENO-E on unstructured mixed meshes, for all five Riemann solvers. Both of a face's fluxes estimate their wave speeds from each side's own frozen \((\gamma, e_0)\): with each flux's own pair, HLL, Rusanov and RHLL gave the two fluxes different mass fluxes at a contact and perturbed the pressure by 0.1%. Contacts keep \(p\) and \(\mathbf{u}\) uniform to \(10^{-12}\) with MUSCL and TENO5 (chemistry design).

SIMPLER stability

SIMPLER balanced splitting (Wu, Ma & Ihme 2019) gives every cell the flow's tendency during the reaction substep, then integrates the flow over half a step. With an explicit integrator for the flow, the step is a new Runge-Kutta method; with SSPRK3 its stability polynomial is

\[ \begin{equation} P(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{8} + \frac{z^4}{48}, \qquad \lvert P(iy)\rvert^2 = 1 + \frac{y^4}{24} - \frac{y^6}{192} + \frac{y^8}{2304} > 1 \quad (y \ne 0), \label{eq:simpler} \end{equation} \]

so a non-dissipative operator, such as the central flux of the hybrid scheme, is unstable at any time step, while upwind operators gain a larger real-axis limit (chemistry design). The published analyses of balanced splitting (Speth et al. 2013; Wu, Ma & Ihme 2019) treat substeps solved exactly for dissipative operators. The derivation is elementary; we did not find it published.

Evidence. A 2D Taylor-Green vortex grows exponentially with SIMPLER and the hybrid flux (10 times per 50 steps with MUSCL) and stays bounded with Strang splitting or the Riemann flux. Mallard rejects the combination.

RHLL blend

Nishikawa and Kitamura's rotated hybrid applies HLL along the velocity jump and Roe across it. On tetrahedra this turns HLL toward the shock normal on every face that crosses a shock and grows a carbuncle on the stagnation line of blunt bodies. Mallard blends the rotated flux with face-normal HLL,

\[ \begin{equation} \hat{\mathbf{F}} = w\,\hat{\mathbf{F}}_\mathrm{rot} + (1 - w)\,\hat{\mathbf{F}}_\mathrm{HLL}(\mathbf{n}), \qquad w = \min\left(\frac{p^L}{p^R}, \frac{p^R}{p^L}\right)^3, \label{eq:rhll} \end{equation} \]

with the pressure weight of AUSMPW+ (Kim, Kim & Rho 2001); contacts and shear layers keep the rotated flux (theory). Pressure-based shock sensing already controls the dissipation of Roe-type fluxes (RoeM, Kim, Kim & Rho 2003), and the rotated hybrid has been combined with RoeM through a Mach-number weight (Choi et al. 2024); the blend with the rotated hybrid is the only new part, and it is small.

References