Lattice Boltzmann methods, automatic code generation, symbolic computation, method of manufactured solutions, high-performance computing

1 Introduction↩︎

Traditional frameworks for approximating continuous systems of conservation laws on unstructured meshes, most prominently finite volume and finite element methods, evaluate spatial derivatives and interfacial flux tensors by gathering extended, irregularly addressed node neighborhoods. While robust and geometrically flexible, these indirect, non-local memory access patterns map poorly onto coalesced, throughput-oriented memory systems on contemporary High-Performance Computing (HPC) architectures. As modern HPC nodes commonly lean on thousands of concurrent execution threads within heterogeneous GPU accelerators, the performance of Partial Differential Equation (PDE) solvers is increasingly bound by memory bandwidth and synchronization limits rather than raw floating-point throughput.

The classical Lattice Boltzmann Method (LBM) [1] circumvents this bottleneck by considering the continuous problem at the mesoscopic level. By tracking particle distribution functions across a discrete velocity space, the algorithm decomposes into a perfectly parallel collision step and a neighborhood-local streaming step along discrete characteristics. This locality yields excellent scalability on heterogeneous supercomputers [2]–[6]. However, standard scalar LBM formulations are structurally rigid. Constructing a scheme for a new target is fundamentally a moment-matching problem: the discrete velocity moments of the equilibrium distribution must reproduce, order by order, the macroscopic densities and fluxes of the target equations. This matching is classically realized through high-order Hermite expansions of the equilibrium, custom-tailored to a limited set of target equations [7], [8], or through entropic equilibria on enlarged velocity sets [9]. Because the number of independent moments a lattice can represent isotropically is fixed by its stencil, supporting the moments a given physics demands typically requires case-specific multi-speed lattices.

To overcome these limitations, a class of LBM schemes relying on a point-wise, state-decoupled kinetic formulation was recently advanced in [10]. By assigning an independent set of distribution functions to each element of the macroscopic state vector, this approach embeds nonlinear physical fluxes directly in the first-order discrete moments, so a single compact lattice serves diverse conservation laws [11]. Such point-wise kinetic schemes have historically remained case-specific, requiring tedious manual derivation and custom code bases for each targeted physical system [10], [12]–[14]. Yet this case-by-case effort is a single mechanical procedure, much as automatic differentiation reduced the hand calculus of derivatives to one graph transformation. The enabling observation is that the flux-in-first-moment construction is an instance of the discrete-kinetic relaxation approximation of conservation laws [15], [16] (Section 2.2.2): a single map from a system of conservation laws to a lattice scheme, whose domain of validity is fixed a priori by a sub-characteristic condition (Section 2.2.4) rather than rediscovered case by case.

The present work goes beyond these specialized formulations by establishing PDE2LBM, a generalized and fully automated symbolic compiler for this class of LBM schemes. Automated code generation for LBM is itself well established: templated GPU solvers such as Sailfish [17] and TCLB [18], the lbmpy/waLBerla toolchain [4], [19], [20], and OpenLB’s expression-level kernel optimization and generated adjoints [6], [21]. All of these optimize or transform the executable form of a lattice scheme the user must still supply: collision model, equilibrium, stencil, and grid scaling enter as givens. Our contribution sits one step upstream, deriving the equilibrium, gradient-tracking cascade, and grid scaling automatically from the continuous PDE, and kernel-level generation can consume the schemes we emit. Driven by a coordinate-free, declarative Domain-Specific Language (DSL) built on top of SymPy [22], our framework automatically deconstructs systems of mixed hyperbolic and parabolic PDEs into point-wise advection-relaxation loops. When higher-order spatial derivatives or nonlinear coupling tensors are present, the compiler recursively extends the global state vector with auxiliary local gradient-tracking states. This reduces multi-node differential operators to nested, first-order kinetic interactions, bypassing macroscopic difference stencils to preserve local concurrency across the calculation.

In this paper we present the theoretical foundation, validation, and hardware execution of this multi-physics compiler, which targets the platform-transparent framework OpenLB [6], [23]: asymptotic convergence across twelve PDE systems, and near-roofline performance of the generated GPU kernels.

The remainder of this manuscript is organized as follows: Section 2 analyzes the kinetic formulation, its macroscopic recovery via Chapman–Enskog expansion, and the automated gradient-tracking cascade. Section 3 states the continuous formulation and compiled configuration of the test cases. Section 4 reports the recursive Method of Manufactured Solutions (MMS) convergence study, and Section 5 evaluates the roofline efficiency of the generated kernels.

Listing [lst:dsl_euler] shows the complete user-facing specification of the compressible Euler equations in our DSL: from these few coordinate-free lines the compiler generates a validated OpenLB solver, with no hand-derived stencils or equilibria.

from pde2lbm import *

eqs = ConservationLaws(dim=2)

# 1. Define Variables
rho, E = eqs.scalars("rho", "E")
rhoU = eqs.vector("rhoU")

gamma = eqs.register_parameter("gamma")

# 2. Dependent Variables
u = rhoU / rho
p = (gamma - 1) * (E - 0.5 * rho * u.dot(u))

# 3. Equation System
eqs.add([
    Eq(dt(rho)  + div(rhoU)                       , 0),
    Eq(dt(rhoU) + div(outer(rhoU, u) + p * eye(2)), 0),
    Eq(dt(E)    + div((E + p) * u)                , 0)
])

# 4. Track Statistics
eqs.track("rho",  rho)
eqs.track("uSqr", u.dot(u))

# 5. Generate OpenLB solver
eqs.compile(class_name="EulerDynamics")

2 Methodology↩︎

2.1 General Conservation Laws↩︎

We consider a system of \(K\) continuous PDEs specified on a \(D\)-dimensional spatial domain \(\Omega \subset \mathbb{R}^D\). The target problem class is cast into the general conservative form: \[\begin{align} \label{eq:pde95target} \partial_t \mathbf{Q} + \nabla \cdot \mathbf{\Phi} = \mathbf{S}, \end{align}\tag{1}\] where \(\mathbf{Q}=[Q_1, Q_2, \dots, Q_K]^T \in \mathbb{R}^K\) denotes the primary macroscopic state vector field. The physical flux tensor \(\mathbf{\Phi} \in \mathbb{R}^{K \times D}\) and the localized source vector \(\mathbf{S} \in \mathbb{R}^K\) capture linear or nonlinear algebraic couplings. These transport terms can depend not only on the primary states \(\mathbf{Q}\) but also on their spatial gradients \(\mathbf{G}=\nabla \mathbf{Q}\) or higher-order derivative tensors, which are resolved through nested gradient tracking hierarchies within a local, component-wise kinetic formulation.

2.2 A Discrete Kinetic Relaxation Scheme↩︎

2.2.1 Vector Populations and the Exact-Flux Equilibrium↩︎

To solve 1 while maintaining spatial locality, the continuous system is mapped into a discrete kinetic phase space by allocating an independent set of particle distribution functions \(f_{k,i}(\mathbf{x},t)\) for each component \(k \in \{1,\dots,K\}\) of the macroscopic state vector. The formulation is dimension- and stencil-agnostic, accommodating any symmetric lattice whose discrete velocities \(\mathbf{c}_i\) and weights \(w_i\) satisfy the tensor isotropy conditions: \[\begin{align} \sum_i w_i = 1, \quad \sum_i w_i \mathbf{c}_i = \mathbf{0}, \quad \sum_i w_i \mathbf{c}_i \otimes \mathbf{c}_i = c_s^2 \mathbf{I}, \end{align}\] where \(\mathbf{I}\) is the identity tensor and \(c_s\) is the dimensionless lattice speed of sound.

The primary macroscopic variables are recovered point-wise via the zeroth moment: \[\begin{align} \label{eq:moments} Q_k(\mathbf{x}, t) = \sum_{i=0}^{I-1} f_{k,i}(\mathbf{x}, t). \end{align}\tag{2}\] Classical LBM equilibria truncate a Hermite expansion of the Maxwell–Boltzmann distribution. We instead embed the analytical physical flux \(\mathbf{\Phi}_k\) directly and exactly into the first-order discrete moment of the equilibrium distribution function \(f_{k,i}^{\text{eq}}\), which is therefore linear in the macroscopic data: \[\begin{align} \label{eq:feq} f_{k,i}^{\text{eq}}(\mathbf{x}, t) = w_i \left( Q_k(\mathbf{x}, t) + \frac{\mathbf{c}_i \cdot \left[ \frac{\Delta t}{\Delta x} \mathbf{\Phi}_k(\mathbf{x}, t) \right]}{c_s^2} \right). \end{align}\tag{3}\] Because no series in \(\mathbf{u}/c_s\) is formed, no low-Mach truncation ceiling arises from this construction: the minimal lattice serves at high Mach number without enlarging the velocity set (wave speeds are instead limited by the sub-characteristic condition of Section 2.2.4), whereas moment-matched Maxwell–Boltzmann constructions reach comparable compressibility through enlarged velocity sets [9] or reformulations that abandon the fixed lattice frame [24]. We adopt this per-field, flux-in-first-moment equilibrium in the form of [10], where it is exercised in the strongly compressible regime. Taking the exact discrete moments verifies the validity of the linear projection: \[\begin{align} \label{eq:eq95moments} \sum_i f_{k,i}^{\text{eq}} = Q_k, \quad \sum_i f_{k,i}^{\text{eq}} \mathbf{c}_i = \frac{\Delta t}{\Delta x} \mathbf{\Phi}_k. \end{align}\tag{4}\] The discrete time-marching follows the single-relaxation-time Bhatnagar–Gross–Krook (BGK) equation with an explicit source term: \[\begin{align} \label{eq:lbe95update} f_{k,i}(\mathbf{x} + \mathbf{c}_i \Delta x, t + \Delta t) = f_{k,i}(\mathbf{x}, t) - \frac{1}{\tau_{\text{LB}}} \left( f_{k,i}(\mathbf{x}, t) - f_{k,i}^{\text{eq}}(\mathbf{x}, t) \right) + \Delta t w_i S_k^{\text{total}}(\mathbf{x}, t), \end{align}\tag{5}\] where \(\tau_{\text{LB}}\) is the relaxation time parameter, and \(S_k^{\text{total}}\) incorporates the target PDE sources. The formulation extends to other collision models such as RLB [25].

2.2.2 Relation to Relaxation Schemes↩︎

The construction of is not merely analogous to a relaxation method but an instance of one. Introducing the physical lattice velocities \(\mathbf{\xi}_i = \mathbf{c}_i \frac{\Delta x}{\Delta t}\), the equilibrium moment constraints 4 read \[\begin{align} \label{eq:moment95constraints} \sum_i f_{k,i}^{\text{eq}} = Q_k, \qquad \sum_i \mathbf{\xi}_i\, f_{k,i}^{\text{eq}} = \mathbf{\Phi}_k, \end{align}\tag{6}\] which are precisely the conditions defining the multidimensional discrete-velocity kinetic approximation of conservation laws analysed by Aregba-Driollet and Natalini [16], a class that includes the Jin–Xin relaxation schemes [15]. In this correspondence the BGK collision is the discretization of a local relaxation toward equilibrium with relaxation time \(\varepsilon = (\tau_{\text{LB}} - \tfrac{1}{2})\Delta t\), while the discrete velocities \(\mathbf{\xi}_i\) act as the constant characteristic speeds of the underlying linear transport operator. All nonlinearity resides in the local equilibrium. The streaming is linear and constant-coefficient.

The canonical scalar instance makes the identification concrete: For a one-dimensional law \(\partial_t u + \partial_x f(u) = 0\), the two-velocity (\(\xi_{1,2} = \pm a\)) scheme has equilibria \[\begin{align} \label{eq:maxwellians} M_{\pm}(u) = \tfrac{1}{2}u \pm \tfrac{1}{2a} f(u), \end{align}\tag{7}\] whose zeroth and first moments return \(u\) and \(f(u)\) exactly, the one-dimensional form of 3 . At full relaxation, \(\tau_{\text{LB}} = 1\), the populations are replaced by their equilibria in every step and this two-velocity scheme collapses to the classical Lax–Friedrichs update, the relaxation time thus tuning the stabilizing diffusion of Section 2.2.3 continuously below the Lax–Friedrichs level as \(\tau_{\text{LB}} \to \tfrac{1}{2}\). Recognizing the scheme as a discrete kinetic relaxation scheme is what grounds the consistency and stability statements that follow: the asymptotic recovery of 1 (Section 2.2.3) and the admissibility boundary (Section 2.2.4).

The correspondence between LBM and relaxation or finite-difference schemes is well established [26]–[29], and vectorial schemes placing physical fluxes in the first moments have been constructed by hand for particular settings [14], [28]. Here the correspondence is used generatively: the compiler derives the scheme directly from the supplied conservation law, whose flux Jacobians fix the admissibility condition a priori.

2.2.3 Asymptotic Consistency↩︎

To evaluate asymptotic consistency, we perform a Chapman–Enskog expansion, perturbing the distribution around its equilibrium state: \(f_{k,i} = f_{k,i}^{(0)} + \Delta t f_{k,i}^{(1)} + \mathcal{O}(\Delta t^2)\), where \(f_{k,i}^{(0)} \equiv f_{k,i}^{\text{eq}}\) and \(\sum_i f_{k,i}^{(1)} = 0\). A Taylor series expansion of the left-hand side of 5 yields: \[\begin{align} \label{eq:taylor} \begin{aligned} f_{k,i}(\mathbf{x} + \mathbf{c}_i \Delta x, t + \Delta t) ={}& f_{k,i} + \Delta t \left( \partial_t + \frac{\Delta x}{\Delta t} \mathbf{c}_i \cdot \nabla \right) f_{k,i} \\ &\hphantom{f_{k,i}} + \frac{\Delta t^2}{2} \left( \partial_t + \frac{\Delta x}{\Delta t} \mathbf{c}_i \cdot \nabla \right)^2 f_{k,i} + \mathcal{O}(\Delta t^3). \end{aligned} \end{align}\tag{8}\] Isolating the first-order non-equilibrium perturbation yields: \[\begin{align} f_{k,i}^{(1)} = -\tau_{\text{LB}} \left( \partial_t f_{k,i}^{(0)} + \frac{\Delta x}{\Delta t} \mathbf{c}_i \cdot \nabla f_{k,i}^{(0)} - w_i S_k^{\text{total}} \right). \end{align}\] Taking the first moment of \(f_{k,i}^{(1)}\) with respect to \(\mathbf{\xi}_i = \mathbf{c}_i \frac{\Delta x}{\Delta t}\) defines the non-equilibrium deviation flux \(\mathbf{\Pi}_k^{(1)} = \sum_i f_{k,i}^{(1)} \mathbf{\xi}_i\). Using the isotropy conditions (odd moments vanish on the symmetric stencil), the second equilibrium moment is \(\sum_i f_{k,i}^{(0)} \mathbf{\xi}_i \otimes \mathbf{\xi}_i = a^2 Q_k \mathbf{I}\) with the squared lattice speed \(a^2 = c_s^2 (\Delta x/\Delta t)^2\), giving: \[\begin{align} \label{eq:deviation95flux} \mathbf{\Pi}_k^{(1)} = -\tau_{\text{LB}} \left( \partial_t \mathbf{\Phi}_k + a^2 \nabla Q_k \right). \end{align}\tag{9}\] The standard second-order accounting reduces the prefactor \(\tau_{\text{LB}}\) to \(\tau_{\text{LB}} - \tfrac{1}{2}\) [1], so the deviation enters the macroscopic balance through the relaxation time \(\varepsilon = (\tau_{\text{LB}} - \tfrac{1}{2})\Delta t\) of Section 2.2.2: \[\begin{align} \label{eq:recovered95raw} \partial_t Q_k + \nabla \cdot \mathbf{\Phi}_k = S_k^{\text{total}} + \nabla \cdot \left( \varepsilon\, \partial_t \mathbf{\Phi}_k + \varepsilon\, a^2\, \nabla Q_k \right) + \mathcal{O}(\varepsilon^2), \end{align}\tag{10}\] The correction is not yet closed: it still contains the time derivative \(\partial_t \mathbf{\Phi}_k\). We take the flux to be an algebraic function of the state, \(\mathbf{\Phi} = \mathbf{\Phi}(\mathbf{Q})\). Gradient-dependent fluxes are reduced to this form by the tracker construction of Section 2.3, and prescribed coefficient fields (the homogenized compressible-flow case of Section 3) add an \(\mathcal{O}(\varepsilon)\) closure term that vanishes under refinement. To leading order the conserved fields then obey \(\partial_t \mathbf{Q} = -\nabla\cdot\mathbf{\Phi} + \mathbf{S} + \mathcal{O}(\varepsilon)\), so with the flux Jacobians \(\mathbf{A}_\alpha = \partial \mathbf{\Phi}_\alpha / \partial \mathbf{Q}\), \[\begin{align} \label{eq:flux95time95closure} \partial_t \mathbf{\Phi}_k = \mathbf{A}_k\, \partial_t \mathbf{Q} = -\,\mathbf{A}_k \mathbf{A}_\ell\, \partial_\ell \mathbf{Q} + \mathbf{A}_k \mathbf{S} + \mathcal{O}(\varepsilon). \end{align}\tag{11}\] Substituting 11 fuses the temporal and spatial corrections into a single second-order operator, the relaxation diffusion, \[\begin{align} \label{eq:relax95viscosity} \partial_t Q_k + \nabla \cdot \mathbf{\Phi}_k &= S_k^{\text{total}} + \nabla \cdot \big( \boldsymbol{\nu}^{\text{eff}}\, \nabla \mathbf{Q} \big)_k + \mathcal{O}(\varepsilon^2), \\ \boldsymbol{\nu}^{\text{eff}}(\mathbf{n}) &= \varepsilon\big( a^2 \mathbf{I} - \mathbf{A}_\mathbf{n}^2 \big), \qquad \mathbf{A}_\mathbf{n} = \textstyle\sum_\alpha n_\alpha \mathbf{A}_\alpha, \end{align}\tag{12}\] anisotropic through the directional flux Jacobian \(\mathbf{A}_\mathbf{n}\), the combination that governs a plane wave propagating along the unit direction \(\mathbf{n}\) (the source coupling \(\mathbf{A}_k \mathbf{S}\) of 11 modifies only the effective source, not the principal symbol, and is absorbed into \(S_k^{\text{total}}\)). For a scalar law this collapses to \(\nu^{\text{eff}} = \varepsilon\,(a^2 - f'(u)^2)\), anticipated from the Jin–Xin equilibria of 7 . Whether \(a^2\mathbf{I} - \mathbf{A}_\mathbf{n}^2\) keeps the net operator dissipative is the question taken up in Section 2.2.4.

Under refinement the correction is controlled by \(\varepsilon = (\tau_{\text{LB}} - \tfrac{1}{2})\Delta t\) together with the lattice speed \(a\). Throughout this work the relaxation time is parametrized by a resolution-independent rate \(\tau_R > 0\) as \(\tau_{\text{LB}} = \tfrac{1}{2} + \tau_R\,\Delta t\), so \(\varepsilon a^2 = c_s^2\,\tau_R\,\Delta x^2\) contracts under any refinement law: the relaxation diffusion is a second-order-consistent numerical stabilizer, never a modelled transport coefficient (at fixed \(\tau_{\text{LB}}\) it would instead converge to a finite limit, the classical route by which lattice Boltzmann schemes model physical diffusion [27]). Physical diffusion, where the target carries it, enters through the tracked-gradient fluxes of Section 2.3.

The two refinement laws differ in which lattice number they hold fixed. Diffusive scaling (\(\Delta t \propto \Delta x^2\)) fixes the lattice diffusion number \(D\,\Delta t/\Delta x^2\) of the tracked diffusive flux, the explicit-stability requirement of a parabolic principal part, while the lattice speed grows as \(a \propto \Delta x^{-1}\) and the sub-characteristic condition of Section 2.2.4 relaxes under refinement. Acoustic scaling (\(\Delta t \propto \Delta x\)) instead fixes the wave Courant number \(\lambda\,\Delta t/\Delta x\) of the physical signal speeds \(\lambda\) while \(\boldsymbol{\nu}^{\text{eff}} = \mathcal{O}(\Delta x^2) \to 0\), the natural regime for hyperbolic targets, which carry no physical diffusion to sustain. In either regime the truncation error is \(\mathcal{O}(\Delta x^2)\). No single refinement holds both numbers fixed as \(\Delta x \to 0\), so each target is refined by whichever timescale binds. The source enters once per fused step, evaluated at the temporal midpoint \(t + \Delta t/2\), which keeps the source coupling second-order accurate.

2.2.4 A Priori Stability: The Sub-Characteristic Condition↩︎

The relaxation structure that supplies the consistency above also delimits the class of systems the framework can stably represent. The leading correction \(\boldsymbol{\nu}^{\text{eff}}\) of 12 is dissipative only where it is positive semi-definite. Where it is indefinite, the leading-order correction is anti-diffusive. For a scalar conservation law with lattice speed \(a\) this is the sign of \(a^2 - f'(u)^2\), i.e. the classical sub-characteristic condition [30], [31] \[\begin{align} \label{eq:subchar95scalar} a \ge |f'(u)| \qquad \text{for all attained } u, \end{align}\tag{13}\] i.e. the lattice must propagate faster than the fastest physical wave.

For a system the requirement is positive semi-definiteness of \(\boldsymbol{\nu}^{\text{eff}}(\mathbf{n}) = \varepsilon(a^2\mathbf{I} - \mathbf{A}_\mathbf{n}^2)\) in every direction. For systems admitting a symmetrizer (equivalently, a convex entropy) [32], [33], \(\mathbf{A}_\mathbf{n} = \sum_\alpha n_\alpha \mathbf{A}_\alpha(\mathbf{Q})\) is diagonalizable with real spectrum in the symmetrized inner product, and positive semi-definiteness of \(a^2\mathbf{I} - \mathbf{A}_\mathbf{n}^2\) is equivalent to the lattice dominating the physical wave speeds in every direction, so that the scalar bound becomes the spectral condition \[\begin{align} \label{eq:subchar95system} \max \left| \operatorname{eig} \sum_\alpha n_\alpha \mathbf{A}_\alpha(\mathbf{Q}) \right| \le a, \qquad \forall\, \mathbf{n}, \;\|\mathbf{n}\| = 1. \end{align}\tag{14}\] Both inequalities are checkable symbolically from the Jacobians of the fluxes the compiler already assembles, yielding an a priori admissibility criterion and, where it fails, a principled rescaling of the lattice speed or time step. The criterion identifies where the leading-order operator is dissipative.

Systems whose principal symbol violates hyperbolicity or whose signal speeds no finite lattice can dominate (elliptic constraints, higher-order dispersive operators) fall outside the admissible class a priori, delimiting the targets the framework represents by construction.

2.3 Recursive Gradient Tracking↩︎

When physical fluxes or sources contain spatial derivatives, traditional methods typically rely on multi-node difference stencils. To preserve locality, the component-wise kinetic formulation maps any spatial gradient requirement onto an independent auxiliary tensor \(\mathbf{G}\) that is appended directly to the global state vector \(\mathbf{Q}\). Following this approach [10], the spatial gradient of a continuous parent field \(V\) is resolved via a dedicated linear advection-relaxation equation: \[\begin{align} \label{eq:gradient95advection} \partial_{t}\mathbf{G} + \nabla \cdot \mathbf{\Phi}_{\mathbf{G}} = \mathbf{S}_{\mathbf{G}}, \end{align}\tag{15}\] where the auxiliary flux tensor and relaxation source vector are defined as: \[\begin{align} \mathbf{\Phi}_{\mathbf{G}} = \mathbf{G} \otimes \mathbf{u} - \mathbf{I} \otimes \frac{V}{\tau_{\nabla}}, \quad \mathbf{S}_{\mathbf{G}} = -\frac{\mathbf{G}}{\tau_{\nabla}}. \end{align}\] Here, \(\mathbf{u}\) represents a background advection velocity vector, \(\mathbf{I}\) is the Kronecker identity tensor, and \(\tau_{\nabla}\) is a small gradient-tracker relaxation time. Pre-scaling the injected field by \(1/\tau_{\nabla}\) makes the tracker carry the physical gradient directly.

In the limit \(\tau_{\nabla} \to 0\) the algebraic relaxation penalty dominates, forcing \(\partial_{t}\mathbf{G} \to 0\). Since \(\nabla \cdot \big(\mathbf{I} \otimes \tfrac{V}{\tau_{\nabla}}\big) \equiv \tfrac{\nabla V}{\tau_{\nabla}}\), the steady-state balance recovers the spatial derivative locally: \[\begin{align} \label{eq:gradient95recovery} \mathbf{G} = \nabla V + \mathcal{O}(\tau_{\nabla}). \end{align}\tag{16}\] Each tracker is a separate linear advection-relaxation system whose equilibrium encodes \(\nabla V\). Unlike the flux-in-first-moment construction of Section 2.2.2, it relaxes toward the gradient of an external field rather than an algebraic flux of a conserved state, so its \(\mathcal{O}(\tau_{\nabla})\) consistency follows from 16 directly rather than from the relaxation-to-conservation-law analysis, with a relaxation lag whenever the parent field is unsteady. The tracker time is fixed at \(\tau_{\nabla} = g_\tau\,\Delta t\) with a resolution-independent step count \(g_\tau\) (supplementary materials), so offset and lag are both \(\mathcal{O}(\Delta t)\): \(\mathcal{O}(\Delta x^2)\) under diffusive scaling and \(\mathcal{O}(\Delta x)\) under acoustic scaling, where it enters only through dissipative fluxes whose small diffusivities keep it below the \(\mathcal{O}(\Delta x^2)\) remainder at the reported resolutions (cf.Table ¿tbl:tab:pde95taxonomy?). First derivatives are alternatively available from the non-equilibrium populations, but we adopt prognostic trackers uniformly here and note that route as future work. For physical systems demanding higher-order derivatives, nested gradient tracking hierarchies are constructed by recursively applying this tracking mechanism to previously generated auxiliary states (e.g., \(\mathbf{G}_{2} = \nabla \mathbf{G}_{1} = \nabla(\nabla V)\)). This reduces multi-node macroscopic differential operators to chains of nested, first-order kinetic interactions, maintaining stencil locality throughout the calculation.

2.4 Declarative DSL and Automated Compiler↩︎

To bridge the gap between abstract formulations and optimized execution on real-world hardware, we implement a symbolic compiler: rather than requiring manual, error-prone decomposition of physical fluxes, spatial derivatives, and gradient-tracking variables, our framework exposes a coordinate-free declarative DSL on top of the SymPy Computer Algebra System (CAS) [22].

The continuous model of any target system is declared by applying coordinate-free operators directly to symbolic state variables: the temporal derivative (\(\mathrm{dt}\)), spatial gradient (\(\mathrm{grad}\)), divergence (\(\mathrm{div}\)), Laplacian (\(\mathrm{laplacian}\)), outer tensor product (\(\mathrm{outer}\)), and identity tensor (\(\mathrm{eye}\)). For instance, the nonlinear Euler equations of compressible fluid dynamics are declared in a few lines (Listing [lst:dsl_euler] in the Introduction). The compiler then processes these declared equations through a sequence of symbolic transformation passes prior to C++ code generation, as detailed in Algorithm 2.

The compiler recursively resolves gradient cascades, which we illustrate on the compressible Navier–Stokes–Fourier system that serves as the running example through the rest of this section. Its energy flux carries a Fourier heat-conduction term \(-K\nabla T\) with temperature \(T = p/(\rho R)\), and its viscous stress \(\mathbf{\tau}\) depends on the velocity gradient \(\nabla\mathbf{u}\). Neither \(\nabla T\) nor \(\nabla\mathbf{u}\) is the gradient of a conserved variable, so the Cascade step of 2 expands them by the chain rule over the conserved fields the scheme can transport (\(\mathbf{G}_\rho = \nabla\rho\), \(\mathbf{G}_{\rho\mathbf{u}} = \nabla(\rho\mathbf{u})\), \(\mathbf{G}_E = \nabla E\)), giving \(\nabla\mathbf{u} = \rho^{-1}(\mathbf{G}_{\rho\mathbf{u}} - \mathbf{u}\otimes\mathbf{G}_\rho)\) for the stress and \(\nabla T = \rho^{-1}(R^{-1}\mathbf{G}_p - T\,\mathbf{G}_\rho)\) for the heat flux (with the pressure gradient \(\mathbf{G}_p=\nabla p\) itself reconstructed from \(\mathbf{G}_\rho,\mathbf{G}_{\rho\mathbf{u}},\mathbf{G}_E\)), then resolves the resulting conserved-field gradients into trackers. A single sweep therefore resolves this chain of nested first-order derivatives without ever forming a multi-node stencil. Systems whose fluxes contain still higher derivatives trigger additional passes of the same loop, each registering the gradient of a previously generated auxiliary state, so the cascade terminates at the highest derivative order present in \(\mathbf{\Phi}\) and \(\mathbf{S}\).

2.4.1 Automated Physical-Unit Scaling↩︎

In traditional hand-coded LBMs, converting physical variables \(Q_{\text{phys}}\) into dimensionless lattice variables \(Q_{\text{LBM}}\) requires solving by hand for the exponents \(p\) and \(q\) in the grid-dependent conversion \(Q_{\text{phys}} = Q_{\text{LBM}} \cdot dx^p \, dt^q\), with \(dx\) and \(dt\) the spatial and temporal lattice spacings. To eliminate this, our symbolic DSL lets the user declare the governing PDE directly in physical (SI) units: each state variable and parameter is registered together with its physical dimension (a velocity as \(\mathrm{length}/\mathrm{time}\), a kinematic transport coefficient as \(\mathrm{length}^2/\mathrm{time}\), and so on). The compiler then reduces every quantity to the lattice scales by a Buckingham-\(\Pi\) projection (the Dedimensionalize step of 2). Each declared quantity \(X\) carries an SI dimension whose exponents over the \(b\) base dimensions (length and time first) form a rational vector \(\mathbf{D}(X)\in\mathbb{Q}^b\). The lattice carries only two independent scales, the spacings \(dx\) and \(dt\), so reducing \(\mathbf{D}(X)\) to them is a linear map \(\pi:\mathbf{D}(X)\mapsto(p,q)\) fixed by three requirements: length and time map to themselves, \(\pi(\mathrm{length})=(1,0)\) and \(\pi(\mathrm{time})=(0,1)\), and each remaining base dimension is eliminated through a characteristic reference quantity \(R_3,\dots,R_b\) already present in the model and held lattice-invariant, \(\pi(\mathbf{D}(R_k))=(0,0)\). These conditions make \(\{\mathrm{length},\mathrm{time},R_3,\dots,R_b\}\) an alternative basis of dimension space, whose exponent vectors form the columns of the invertible \(b\times b\) matrix \(B = [\, \mathbf{e}_1,\;\mathbf{e}_2,\;\mathbf{D}(R_3),\;\dots,\;\mathbf{D}(R_b) \,]\), so that the projection becomes the explicit linear map \[\label{eq:pq95projection} \big(p(X),\,q(X)\big)^{\!\top} = \Pi\,\mathbf{D}(X), \qquad \Pi = [\,I_2 \mid \mathbf{0}\,]\,B^{-1},\tag{17}\] in which \(B^{-1}\) re-expresses any SI dimension in the \(\{\mathrm{length},\mathrm{time},\text{reference}\}\) basis and \(\Pi\) reads off its length and time content. The references need only be dimensionally independent. The exponents are in general non-zero and field-specific (for example \((p, q) = (1, -1)\) for a velocity), degenerating to \((p, q) = (0, 0)\) for already lattice-dimensionless quantities, including each characteristic reference itself: the density, once mass is eliminated through \(\rho_0\), carries no grid-dependent conversion at all.

Fluxes and sources carry their own conversions: a flux is the conserved quantity carried at a velocity, so it inherits the field’s own scaling \(dx^{-p}\,dt^{-q}\) together with one velocity conversion factor \(dt/dx\): \[\mathbf{F}_{\text{LBM}} = \mathbf{F}_{\text{phys}} \cdot dx^{-p}\,dt^{-q}\,\frac{dt}{dx} = \mathbf{F}_{\text{phys}} \cdot dx^{-(p+1)}\,dt^{-(q-1)},\] while a source, being a rate, picks up one time step instead, \(\mathbf{S}_{\text{LBM}} = \mathbf{S}_{\text{phys}} \cdot dx^{-p}\,dt^{1-q}\). Each gradient tracker inherits the field exponents \((p-1, q)\) of its parent, so the gradient cascade is scaled consistently. The compiler injects these transformations automatically: each PDE is written once, in SI form, with no hand rescaling.

Continuing the running example, this machinery is exercised on the same system declared in SI units, the dynamic viscosity (\(\unit{\pascal\second}\)), thermal conductivity (\(\unit[per-mode=power]{\watt\per\metre\per\kelvin}\)), and specific gas constant each carrying named SI dimensions (the complete DSL specification is given in the supplementary material). Eliminating the mass and temperature dimensions through the density and the gas constant and projecting via 17 , both the dynamic viscosity and the thermal conductivity collapse to the identical \((p,q)=(2,-1)\) scaling of a \(\mathrm{length}^2/\mathrm{time}\) diffusivity, recovering the momentum and thermal diffusion scalings with no hand input. The result is the twelve-field compressible Navier–Stokes–Fourier kernel verified at second order in Section 3, generated and run with no hand-derived equilibrium, stencil, or unit conversion.

2.4.2 Automated Reference and Equilibrium Shifting↩︎

Many target systems carry a small fluctuation \(\delta \mathbf{Q}(\mathbf{x}, t)\) on a large background \(\mathbf{Q}_0\), \(|\delta \mathbf{Q}| \ll |\mathbf{Q}_0|\). Populations formed from absolute state values then push the active variations into the low mantissa bits, and spatial differences suffer catastrophic cancellation, especially in single precision (float).

Each component \(Q_k\) is carried by its own substencil whose resting equilibrium \(f_{k, i}^{\text{eq}} = w_i Q_k\) distributes \(Q_k\) over the lattice weights, so near a background reference state \(\mathbf{Q}_0\) even the quiescent populations carry a large static offset \(w_i Q_{k, 0} \neq 0\) [23].

The compiler removes this bias at code-generation time through two coordinated shifts. First, a reference-state shift moves the macroscopic origin off the large background: the compiler evolves only the fluctuation \(\delta\mathbf{Q} = \mathbf{Q} - \mathbf{Q}_0\) about a declared background state \(\mathbf{Q}_0\). Second, an equilibrium-population shift subtracts the matching background equilibrium from every population: evaluating \(\mathbf{Q}_0\) and its flux \(\mathbf{F}_{k, 0}\), the compiler stores only the deviation \[f_{k, i}^* = f_{k, i} - f^{\text{eq}}_{k, i}(\mathbf{Q}_0, \mathbf{F}_0),\] subtracting the background equilibrium of 3 , so the zeroth-moment state offset \(w_i Q_{k, 0}\) and the first-moment flux offset are removed together. The kernels then run entirely on these shifted populations \(f_{k, i}^* = \mathcal{O}(\delta Q)\), which vanish at background equilibrium, so the full mantissa is spent on the active fluctuation. This reduces roundoff and improves stability, decisively so under single precision.

3 Applications↩︎

To validate the generic mapping layer of the compiler framework, a suite of twelve continuous target systems has been compiled, taxonomized in Table ¿tbl:tab:pde95taxonomy?. The systems range over fluid, electromagnetic, and solid mechanics. All are evaluated on a two-dimensional periodic domain \(\Omega=[0,1]^2\), except the three-dimensional Homogenized Compressible Navier–Stokes–Fourier (HCNSF) case on \(\Omega=[0,1]^3\). Each system is refined under the lattice scaling matched to its dominant transport mechanism (Section 2.2.3): diffusive (\(\Delta t \propto \Delta x^2\)) for the diffusion-resolving parabolic systems and acoustic (\(\Delta t \propto \Delta x\)) for the hyperbolic and convection-dominated ones.

3.1 Inviscid Burgers↩︎

The inviscid Burgers equation is the minimal scalar conservation law with a state-dependent characteristic speed \(f'(u) = u\). It isolates the core mechanism of the framework (the exact nonlinear flux \(\tfrac{1}{2}u^2\) embedded directly in the equilibrium first moment) with no gradient trackers, no coupling, and a single field. The scalar law reads \[\begin{align} \partial_t u + \nabla \cdot \begin{pmatrix*}[l] \tfrac{1}{2}u^2 \\ \tfrac{1}{2}u^2 \end{pmatrix*} = 0, \end{align}\] so that \(\nabla \cdot \mathbf{\Phi} = u\left( \partial_{x_0} u + \partial_{x_1} u \right)\), where \(u\) is a velocity-like scalar field. This continuous system maps directly onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] u \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] \tfrac{1}{2}u^2 \\ \tfrac{1}{2}u^2 \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \end{pmatrix*}. \end{align}\] Because \(\max|f'(u)| = \max|u|\), the manufactured amplitude directly sets the physical wave speeds the lattice must dominate, making this the cleanest setting in which to probe the sub-characteristic bound empirically.

3.1.0.1 Empirical Sub-Characteristic Boundary

Figure 3 reports a brute-force stability sweep of the generated scheme over the manufactured Mach number \(\mathrm{Ma} = A/c_s\) and the collision rate \(\omega\). For each pair the simulation is run to \(t = 1.0\) and classified as stable or divergent, while the background shading records the relative \(L^2\) error at \(N = 512\). The empirically resolved boundary tightens under grid refinement toward \(\mathrm{Ma} = 1/\sqrt{3}\) at \(\omega \to 2\): the bound \(\max|f'(u)| \le c_s^2\) at which the exact-flux equilibria 3 cease to be monotone in \(u\), the monotone-equilibria condition of the relaxation framework [16]. It sharpens the sub-characteristic condition of 13 by the factor \(c_s\) in the over-relaxation limit, while under heavier relaxation damping the boundary rises toward the per-axis wave-speed ceiling. The verification working point (\(A = 1.3\), starred) lies safely within the stable interior, consistent with the a priori criterion of Section 2.2.4.

3.2 Compressible Euler↩︎

The Euler equations govern the conservation of mass, momentum, and energy in inviscid compressible fluids, where thermal and mechanical energy modes are coupled dynamically through an equation of state. The system reads \[\begin{align} \begin{aligned} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) &= 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p\mathbf{I}) &= \mathbf{0}, \\ \partial_t E + \nabla \cdot ((E+p)\mathbf{u}) &= 0, \end{aligned} \end{align}\] where \(\rho\) represents the fluid mass density, \(\mathbf{u}\) is the physical velocity vector, \(E\) is the total energy density, and \(p\) is the static thermal pressure. LBM formulations for compressible Euler equations have been developed using multi-speed lattices to decouple fluid velocity from temperature [34]. This system maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] \rho \\ \rho \mathbf{u} \\ E \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p\mathbf{I} \\ (E+p)\mathbf{u} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \\ \mathbf{0} \\ 0 \end{pmatrix*}. \end{align}\] The pressure is dynamically closed via the ideal gas equation of state \(p = (\gamma - 1)(E - \frac{1}{2}\rho |\mathbf{u}|^2)\), where \(\gamma = 1.4\) is the adiabatic index.

3.3 Shallow Water↩︎

The Shallow Water equations describe flow propagation in open channels and coastal zones, modeling depth-averaged fluid thickness and momentum under hydrostatic pressure. LBM formulations for the shallow water equations identify the fluid height with density [35] and the system is also the code-generation showcase of lbmpy [20]. The depth-averaged nonlinear free-surface model simulates channel hydraulics over a static spatial bottom topography field \(z_b(\mathbf{x})\): \[\begin{align} \begin{aligned} \partial_t h + \nabla \cdot (h \mathbf{u}) &= 0, \\ \partial_t (h \mathbf{u}) + \nabla \cdot \left( h \mathbf{u} \otimes \mathbf{u} + \frac{1}{2}g h^2 \mathbf{I} \right) &= -g h \nabla z_b, \end{aligned} \end{align}\] where \(h\) represents the depth-averaged fluid thickness (water column height), \(\mathbf{u}\) is the depth-averaged horizontal velocity vector, \(g\) is the gravitational acceleration constant, and \(z_b\) is the stationary spatial bottom topography profile. It maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] h \\ h \mathbf{u} \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] h \mathbf{u} \\ h \mathbf{u} \otimes \mathbf{u} + \frac{1}{2}g h^2 \mathbf{I} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \\ -g h \nabla z_b \end{pmatrix*}. \end{align}\]

3.4 Ideal Ultrarelativistic Fluid↩︎

The ultrarelativistic fluid is the ideal gas of massless particles (\(c = 1\), equation of state \(p = \frac{e}{3}\)), whose dynamics is closed by energy–momentum conservation alone: the barotropic relativistic Euler system with sound speed \(\frac{1}{\sqrt{3}}\) [36]. LBM formulations have been developed for relativistic [37] and ultrarelativistic [38] hydrodynamics. Tracking the energy density \(E\) and relativistic momentum density vector \(\mathbf{S}\), \[\begin{align} \begin{aligned} \partial_t E + \nabla \cdot \mathbf{S} &= 0, \\ \partial_t \mathbf{S} + \nabla \cdot \left( \mathbf{S} \otimes \mathbf{v} + p \mathbf{I} \right) &= \mathbf{0}, \end{aligned} \end{align}\] the system maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] E \\ \mathbf{S} \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] \mathbf{S} \\ \mathbf{S} \otimes \mathbf{v} + p \mathbf{I} \end{pmatrix*}, \quad \mathbf{S}_{\text{source}} = \begin{pmatrix*}[l] 0 \\ \mathbf{0} \end{pmatrix*}. \end{align}\] For the massless gas the primitive recovery is algebraic: \[\begin{align} p = \frac{\sqrt{4E^2 - 3|\mathbf{S}|^2} - E}{3}, \quad \mathbf{v} = \frac{\mathbf{S}}{E + p}, \end{align}\] with physicality requiring \(E > |\mathbf{S}|\), the null limit \(|\mathbf{v}| \to 1\) attained at equality.

3.5 Maxwell Electromagnetics↩︎

Maxwell’s equations are a first-order linear hyperbolic system of conservation laws whose defining operator is the curl, carried here as a pure flux divergence in the equilibrium first moment. We treat the two-dimensional Transverse-Magnetic (TM) mode in vacuum, with out-of-plane electric field \(E_z\) and in-plane magnetic field \((H_x, H_y)\) propagating at light speed \(c\) [39]: \[\begin{align} \begin{aligned} \partial_t E_z + \nabla \cdot \big( c\,(-H_y,\, H_x) \big) &= 0, \\ \partial_t H_x + \nabla \cdot \big( c\,(0,\, E_z) \big) &= 0, \\ \partial_t H_y + \nabla \cdot \big( c\,(-E_z,\, 0) \big) &= 0, \end{aligned} \end{align}\] which reproduce Ampère’s and Faraday’s laws and, on eliminating \(\mathbf{H}\), the scalar wave equation \(\partial_{tt} E_z = c^2 \nabla^2 E_z\). The system maps onto 1 with \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] E_z \\ H_x \\ H_y \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] c\,(-H_y,\, H_x) \\ c\,(0,\, E_z) \\ c\,(-E_z,\, 0) \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \\ 0 \\ 0 \end{pmatrix*}. \end{align}\]

3.5.0.1 Admissibility at Constant Wave Speed

For Maxwell the admissibility condition of Section 2.2.4 binds identically at every state: the flux Jacobians are constant, with eigenvalues \(\{+c, -c, 0\}\) independent of the field magnitude. With the naive choice \(c = 1\) the wave speed equals the full lattice speed, so the effective diffusion of 12 carries \(c_s^2 - c^2 = \tfrac{1}{3} - 1 < 0\): the recovered operator is anti-diffusive and the scheme is unstable at every \(\tau\). The remedy is to declare \(c\) as a physical quantity of dimension length/time: the dedimensionalizer then scales the flux by the acoustic factor \(\kappa = \Delta t / \Delta x\), as for every hyperbolic case. The lattice wave speed is then \(c\kappa\), and the sub-characteristic bound of 14 becomes the ordinary lattice-Mach constraint \[\begin{align} c\,\kappa \;\le\; c_s = \tfrac{1}{\sqrt{3}}, \qquad \text{i.e.}\qquad \kappa \;\le\; \frac{c_s}{c}. \end{align}\] At the working point \(c = 1\), \(\kappa = 0.2\) the bound holds with nearly threefold margin, and the scheme runs stably at the BGK edge (\(\tau_{\text{LB}} \approx 0.5\)) at second order.

3.6 Nonlinear Finite-Strain Elasticity↩︎

For large-deformation solid mechanics we take the deformation gradient \(\mathbf{F} = \partial \mathbf{x}/\partial \mathbf{X}\) and the velocity \(\mathbf{v}\) as the conserved state, evolved through the kinematic compatibility relation and the momentum balance: \[\begin{align} \begin{aligned} \partial_t \mathbf{F} - \nabla \mathbf{v} &= \mathbf{0}, \\ \partial_t (\rho_0 \mathbf{v}) - \nabla \cdot \mathbf{P} &= \mathbf{0}, \end{aligned} \end{align}\] where \(\rho_0\) is the constant reference-configuration density and \(\mathbf{P}\) is the first Piola–Kirchhoff stress of a compressible neo-Hookean solid, \[\begin{align} \mathbf{P} = \mu_s \left( \mathbf{F} - \mathbf{F}^{-T} \right) + \lambda_s \ln(J)\, \mathbf{F}^{-T}, \qquad J = \det \mathbf{F}, \end{align}\] with shear modulus \(\mu_s\) and first Lamé parameter \(\lambda_s\). This system maps onto 1 with \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] \mathbf{F} \\ \rho_0 \mathbf{v} \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] - \mathbf{v} \otimes \mathbf{I} \\ - \mathbf{P} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] \mathbf{0} \\ \mathbf{0} \end{pmatrix*}. \end{align}\] Prior LBM treatments of elasticity map Cauchy-stress components to single-distribution moments [40] or build vectorial formulations for linear elastic media [12]. The same deformation-gradient formulation was arrived at concurrently by Feng and Chu [41].

3.7 Advection-Diffusion-Reaction↩︎

The Advection-Diffusion-Reaction (ADR) equations model scalar chemical transport under an invariant velocity transport field coupled with linear isotropic diffusion and a nonlinear logistic Fisher–KPP reaction source term. The governing equation reads \[\begin{align} \partial_t c + \nabla \cdot \left( c \mathbf{u}_{\text{adv}} - D \nabla c \right) = R\,c(1-c), \end{align}\] where \(c\) is the scalar concentration field, \(\mathbf{u}_{\text{adv}}\) is the constant advective velocity vector, \(D\) is the isotropic diffusion coefficient, and \(R\) is the reaction rate constant. The LBM has long been applied to such transport, from the original diffusion model [42] to the single-relaxation-time advection-diffusion scheme [43]. The system maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] c \\ \mathbf{G}_c \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] c \mathbf{u}_{\text{adv}} - D \mathbf{G}_c \\ \mathbf{G}_c \otimes \mathbf{u}_{\text{adv}} - \mathbf{I}\,c/\tau_{\nabla} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] R\,c(1-c) \\ - \mathbf{G}_c/\tau_{\nabla} \end{pmatrix*}. \end{align}\]

3.8 Allen–Cahn Phase-Field↩︎

The Allen–Cahn equation models non-conservative interface dynamics, coupling diffusion with a bistable cubic reaction. The governing equation reads \[\begin{align} \partial_t \phi + \nabla \cdot \left( -\gamma \nabla \phi \right) = r\left( \phi - \phi^3 \right), \end{align}\] where \(\phi\) is the dimensionless order parameter, \(\gamma\) is the interface mobility (with units of \(\text{length}^2/\text{time}\)), and \(r\) is the reaction rate. LBM phase-field models have been developed for interface capturing [44]. The system maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] \phi \\ \mathbf{G}_\phi \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] -\gamma \mathbf{G}_\phi \\ - \mathbf{I}\,\phi/\tau_{\nabla} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] r\left( \phi - \phi^3 \right) \\ - \mathbf{G}_\phi/\tau_{\nabla} \end{pmatrix*}. \end{align}\]

3.9 Incompressible Navier–Stokes↩︎

The Incompressible Navier–Stokes Equations (NSE) model the conservation of momentum and mass in viscous, constant-density fluids. The equations read \[\begin{align} \begin{aligned} \nabla \cdot \mathbf{u} &= 0, \\ \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + \frac{1}{\rho_0} \nabla p - \nu \nabla^2 \mathbf{u} &= \mathbf{0}, \end{aligned} \end{align}\] where \(\rho_0\) is the constant mass density, \(\mathbf{u}\) is the fluid velocity vector, \(p\) is the physical pressure, and \(\nu = \mu / \rho_0\) is the kinematic viscosity. Because the divergence-free condition is a kinematic constraint requiring a non-local elliptic Poisson solver, it cannot be directly represented in our explicit local framework (1 ). Instead, isothermal viscous transport is solved using a weakly compressible formulation, where density fluctuations represent a pressure proxy \(p = \rho c_{s,\text{phys}}^2\) in the low-Mach limit, restoring a local, hyperbolic conservation loop: \[\begin{align} \begin{aligned} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) &= 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + \rho c_{s,\text{phys}}^2 \mathbf{I} - \mathbf{\tau} \right) &= \mathbf{0}, \end{aligned} \end{align}\] where \(\mathbf{\tau} = \mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^T)\) is the viscous shear stress tensor with dynamic viscosity coefficient \(\mu\). This system maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] \rho \\ \rho \mathbf{u} \\ \mathbf{G}_{\rho\mathbf{u}} \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + \rho c_{s,\text{phys}}^2 \mathbf{I} - \mathbf{\tau} \\ \mathbf{G}_{\rho\mathbf{u}} \otimes \mathbf{u} - \mathbf{I} \otimes \rho\mathbf{u}/\tau_{\nabla} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \\ \mathbf{0} \\ - \mathbf{G}_{\rho\mathbf{u}}/\tau_{\nabla} \end{pmatrix*}. \end{align}\]

3.10 Compressible Navier–Stokes–Fourier↩︎

The Navier–Stokes–Fourier (NSF) equations, the running example of Section 2.4, govern viscous, heat-conducting gas flows: the full compressible system closing the momentum balance with the Newtonian viscous stress and the energy balance with Fourier heat conduction \(\mathbf{q}=-\kappa\nabla T\). This is the first-order Chapman–Enskog closure, as opposed to extended-moment methods such as Grad-13 [45]. Prior LBM treatments often rely on multi-speed lattices [46]. The governing equations read: \[\begin{align} \begin{aligned} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) &= 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} - \mathbf{\tau} \right) &= \mathbf{0}, \\ \partial_t E + \nabla \cdot \left( (E + p)\mathbf{u} - \mathbf{\tau} \cdot \mathbf{u} - K \nabla T \right) &= 0, \end{aligned} \end{align}\] where \(\rho\) represents the fluid density, \(\mathbf{u}\) is the physical fluid velocity vector, \(E\) is the total energy density, \(p\) is the static thermal pressure, \(T = \frac{p}{\rho R}\) is the temperature field with specific gas constant \(R\), \(K\) is the isotropic thermal conductivity, and \(\mathbf{\tau} = \mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3}(\nabla \cdot \mathbf{u})\mathbf{I} \right)\) is the viscous stress tensor with dynamic viscosity coefficient \(\mu\). The system maps onto 1 with: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] \rho \\ \rho \mathbf{u} \\ E \\ \mathbf{G}_\rho \\ \mathbf{G}_{\rho\mathbf{u}} \\ \mathbf{G}_E \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} - \mathbf{\tau} \\ (E + p)\mathbf{u} - \mathbf{\tau} \cdot \mathbf{u} - K \mathbf{G}_{T}\\ \mathbf{G}_\rho \otimes \mathbf{u} - \mathbf{I}\,\rho/\tau_{\nabla} \\ \mathbf{G}_{\rho\mathbf{u}} \otimes \mathbf{u} - \mathbf{I} \otimes \rho\mathbf{u}/\tau_{\nabla} \\ \mathbf{G}_E \otimes \mathbf{u} - \mathbf{I}\,E/\tau_{\nabla} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \\ \mathbf{0} \\ 0 \\ - \mathbf{G}_\rho/\tau_{\nabla} \\ - \mathbf{G}_{\rho\mathbf{u}}/\tau_{\nabla} \\ - \mathbf{G}_E/\tau_{\nabla} \end{pmatrix*}, \end{align}\] where the trackers and the chain-rule reconstruction of \(\nabla\mathbf{u}\) and \(\mathbf{G}_T = \nabla T\) are those of the running example (Section 2.4).

3.11 Resistive Compressible Magnetohydrodynamics↩︎

The Resistive Magnetohydrodynamics (MHD) equations couple viscous, thermally conducting, electrically active fluids to magnetic fields, capturing shear stress, heat flux, and magnetic dissipation. Existing LBMs represent the magnetic field by vector-valued distribution functions, with the divergence error controlled rather than enforced [47], [48]. The system, augmented by the Godunov–Powell divergence-cleaning that restores the symmetrizability assumed in Section 2.2.4 and controls the growth of \(\nabla\cdot\mathbf{B}\) [49], reads: \[\begin{align} \begin{aligned} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) &= 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + p_{\text{tot}}\mathbf{I} - \mathbf{B} \otimes \mathbf{B} - \mathbf{\tau} \right) &= -(\nabla \cdot \mathbf{B})\,\mathbf{B}, \\ \partial_t E + \nabla \cdot \left( (E + p_{\text{tot}})\mathbf{u} - \mathbf{B}(\mathbf{u} \cdot \mathbf{B}) - \mathbf{\tau} \cdot \mathbf{u} + \mathbf{q} + \eta (\mathbf{j} \times \mathbf{B}) \right) &= -(\nabla \cdot \mathbf{B})(\mathbf{u} \cdot \mathbf{B}), \\ \partial_t \mathbf{B} + \nabla \cdot \left( \mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u} - \eta \nabla \mathbf{B} \right) &= -(\nabla \cdot \mathbf{B})\,\mathbf{u}, \end{aligned} \end{align}\] where \(\rho\) represents the mass density, \(\mathbf{u}\) is the physical velocity vector, \(E\) is the total energy density, \(\mathbf{B}\) is the magnetic field vector, \(\eta\) is the uniform magnetic resistivity, \(\mathbf{\tau} = \mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3}(\nabla \cdot \mathbf{u})\mathbf{I} \right)\) is the viscous shear stress tensor with dynamic viscosity coefficient \(\mu\), \(\mathbf{q} = - K \nabla T\) is the thermal heat flux vector with conductivity \(K\) and temperature \(T = \frac{p_{\text{gas}}}{\rho}\), \(\mathbf{j} = \nabla \times \mathbf{B}\) is the electrical current density, and \(p_{\text{tot}} = p_{\text{gas}} + \frac{1}{2}|\mathbf{B}|^2\) is the total pressure with \(p_{\text{gas}} = (\gamma - 1)(E - \frac{1}{2}\rho |\mathbf{u}|^2 - \frac{1}{2}|\mathbf{B}|^2)\). To evaluate the dissipative gradients without stencils, the system maps onto 1 with: \[\adjustbox{max width=\textwidth}{\displaystyle \mathbf{Q} = \begin{pmatrix*}[l] \rho \\ \mathbf{\rho} \mathbf{u} \\ E \\ \mathbf{B} \\ \mathbf{G}_{\mathbf{B}} \\ \mathbf{G}_\rho \\ \mathbf{G}_{\rho\mathbf{u}} \\ \mathbf{G}_E \end{pmatrix*}, \; \mathbf{\Phi} = \begin{pmatrix*}[l] \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p_{\text{tot}}\mathbf{I} - \mathbf{B} \otimes \mathbf{B} - \mathbf{\tau} \\ (E + p_{\text{tot}})\mathbf{u} - \mathbf{B}(\mathbf{u} \cdot \mathbf{B}) - \mathbf{\tau} \cdot \mathbf{u} + \mathbf{q} + \eta (\mathbf{j} \times \mathbf{B}) \\ \mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u} - \eta \mathbf{G}_{\mathbf{B}} \\ \mathbf{G}_{\mathbf{B}} \otimes \mathbf{u} - \mathbf{I} \otimes \mathbf{B}/\tau_{\nabla} \\ \mathbf{G}_\rho \otimes \mathbf{u} - \mathbf{I}\,\rho/\tau_{\nabla} \\ \mathbf{G}_{\rho\mathbf{u}} \otimes \mathbf{u} - \mathbf{I} \otimes \rho\mathbf{u}/\tau_{\nabla} \\ \mathbf{G}_E \otimes \mathbf{u} - \mathbf{I}\,E/\tau_{\nabla} \end{pmatrix*}, \; \mathbf{S} = \begin{pmatrix*}[l] 0 \\ -(\nabla \cdot \mathbf{B})\,\mathbf{B} \\ -(\nabla \cdot \mathbf{B})(\mathbf{u} \cdot \mathbf{B}) \\ -(\nabla \cdot \mathbf{B})\,\mathbf{u} \\ - \mathbf{G}_{\mathbf{B}}/\tau_{\nabla} \\ - \mathbf{G}_\rho/\tau_{\nabla} \\ - \mathbf{G}_{\rho\mathbf{u}}/\tau_{\nabla} \\ - \mathbf{G}_E/\tau_{\nabla} \end{pmatrix*}. }\] The divergence-cleaning source is itself stencil-free: \(\nabla\cdot\mathbf{B} = \mathrm{tr}\,\mathbf{G}_{\mathbf{B}}\) is read off the same tracker, so the Powell correction is generated automatically alongside the resistive term. This system is treated by hand in the companion work [10]. The compiler derives the same scheme from the conservative form alone, differing only in the dissipative closure: conserved-field gradient trackers with chain-rule reconstruction in place of the transported viscous stress.

3.12 Homogenized Compressible Navier–Stokes–Fourier↩︎

This case showcases two capabilities at once: a three-dimensional lattice (D3Q7) and prescribed external coefficient fields that enter the dynamics without being functions of the evolved state. It models compressible viscous flow interacting with a moving solid through a homogenized formulation: a prescribed porosity field \(\alpha(\mathbf{x},t)\in[0,1]\) (\(\alpha=1\) pure fluid, \(\alpha\to 0\) solid) volume-averages the dissipative fluxes, while a velocity coupling drives the momentum toward a prescribed solid velocity \(\mathbf{u}_s(\mathbf{x},t)\) in the solid region. The porosity acts as a numerical volume penalization rather than a physical homogenization, a homogenized form of Brinkman penalization [50], [51] in the tradition of HLBM fluid–structure coupling [52], [53]: \[\begin{align} \begin{aligned} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) &= 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} - \alpha\,\mathbf{\tau} \right) &= \frac{1-\alpha}{\eta}\big( \rho \mathbf{u}_s(\mathbf{x},t) - \rho \mathbf{u} \big), \\ \partial_t E + \nabla \cdot \left( (E + p)\mathbf{u} - \alpha\,\mathbf{\tau} \cdot \mathbf{u} + \alpha\,\mathbf{q} \right) &= \frac{1-\alpha}{\eta}\big( \rho \mathbf{u}_s - \rho \mathbf{u} \big) \cdot \mathbf{u}, \end{aligned} \end{align}\] with viscous stress \(\mathbf{\tau} = \mu\big( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \tfrac{2}{3}(\nabla \cdot \mathbf{u})\mathbf{I} \big)\), Fourier flux \(\mathbf{q} = -\kappa \nabla T\), temperature \(T = p/\rho\), and pressure \(p = (\gamma - 1)\big( E - \tfrac{1}{2}\rho |\mathbf{u}|^2 \big)\), relaxation time \(\eta\), and the limit \(\alpha \to 1\) recovering the plain compressible NSF system. Relative to the two-dimensional compressible NSF case of Section 3, the closure is built from a reduced tracker basis that evolves the symmetric velocity gradient \(\mathbf{G}_{\mathbf{u}} = \operatorname{sym}\nabla\mathbf{u}\) and the temperature gradient \(\mathbf{G}_T = \nabla T\) directly, now in three dimensions and augmented by the homogenization: \[\begin{align} \mathbf{Q} = \begin{pmatrix*}[l] \rho \\ \rho \mathbf{u} \\ E \\ \mathbf{G}_{\mathbf{u}} \\ \mathbf{G}_T \end{pmatrix*}, \quad \mathbf{\Phi} = \begin{pmatrix*}[l] \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} - \alpha\mathbf{\tau} \\ (E + p)\mathbf{u} - \alpha\mathbf{\tau} \cdot \mathbf{u} - \alpha\kappa \mathbf{G}_T \\ \mathbf{G}_{\mathbf{u}} \otimes \mathbf{u} - \operatorname{sym}(\mathbf{I} \otimes \mathbf{u})/\tau_{\nabla} \\ \mathbf{G}_T \otimes \mathbf{u} - \mathbf{I}\,T/\tau_{\nabla} \end{pmatrix*}, \quad \mathbf{S} = \begin{pmatrix*}[l] 0 \\ \tfrac{1-\alpha}{\eta}(\rho\mathbf{u}_s - \rho\mathbf{u}) \\ \tfrac{1-\alpha}{\eta}(\rho\mathbf{u}_s - \rho\mathbf{u})\cdot\mathbf{u} \\ - \mathbf{G}_{\mathbf{u}}/\tau_{\nabla} \\ - \mathbf{G}_T/\tau_{\nabla} \end{pmatrix*}. \end{align}\]

4 Validation↩︎

To verify the convergence and execution precision of the compiled schemes, we track the global \(L_2\) error norm reduction rates across successive grid resolutions (\(N \in \{64, \dots, 512\}\); \(N \in \{32, \dots, 256\}\) for the three-dimensional case). For each system, the compiler evaluates the continuous residual of the analytical profiles to construct a manufactured source correction: \[\begin{align} \mathbf{S}_{\text{MMS}}(\mathbf{x}, t) = \partial_t \mathbf{Q}_{\text{exact}} + \nabla \cdot \mathbf{\Phi}(\mathbf{Q}_{\text{exact}}) - \mathbf{S}(\mathbf{Q}_{\text{exact}}). \end{align}\] This source field is injected directly into the mesoscopic update loop via 5 . At the final time \(t = 1\) we report, per macroscopic component, the discrete relative \(L_2\) error \(\lVert Q_h - Q_{\text{exact}}\rVert_2 / \lVert Q_{\text{exact}}\rVert_2\) with \(\lVert f\rVert_2 = (N_{\text{cells}}^{-1}\sum_{\mathbf{x}} f(\mathbf{x})^2)^{1/2}\), and take the Empirical Order of Convergence (EOC) to be the negative of the least-squares slope of \(\log_2\lVert e\rVert_2\) against \(\log_2 N\) over the full resolution sequence. The manufactured fields are designed to stress the generated operators rather than for analytical convenience, using large-amplitude fluctuations, broken axis symmetry, independent thermodynamic modes, and exactly satisfied involutions. The compressible working points are genuinely supersonic, with per-axis mean-flow Mach numbers up to \(3\) for the Euler, Navier–Stokes–Fourier, and resistive MHD systems on the same minimal five-velocity lattice. The per-case design principles are detailed in the supplementary materials. All cases run on periodic domains, verifying the bulk scheme in isolation. The generation of boundary operators from the same PDE-level representation is in development.

The resulting convergence matrix across IEEE single and double precisions, each system refined under the lattice scaling matched to its dominant transport mechanism (per the Scaling column of Table ¿tbl:tab:pde95taxonomy?), is documented in Table 1. Every case is fully specified by its refinement law and three resolution-independent scalars: the time-step factor \(f\) (\(\Delta t = f\,\Delta x\) acoustic, \(f\,\Delta x^2\) diffusive), the relaxation rate \(\tau_R\) (\(\tau_{\text{LB}} = \tfrac{1}{2} + \tau_R\,\Delta t\), Section 2.2.3), and the tracker step count \(g_\tau\) (\(\tau_\nabla = g_\tau\,\Delta t\), Section 2.3). The refinement law follows the dominant transport mechanism, \(f\) is bounded by the sub-characteristic condition of Section 2.2.4, and \(\tau_R\) trades stability margin against the relaxation-diffusion error constant. No parameter is retuned per grid. Per-case values are tabulated in the supplementary materials.

2pt

4.1 Convergence Characteristics↩︎

4.1.0.1 Second-Order Convergence in Double Precision

Across the twelve physical systems under double-precision arithmetic, the fitted EOC of the macroscopic state variables is at or near second order, with \(L^2\) slopes between \(1.70\) and \(2.51\). The fitted slopes near \(2.5\) of the three diffusively refined light systems (scalar ADR, Allen–Cahn, incompressible NSE) are not asymptotic superconvergence: their pairwise orders decrease monotonically toward two (\(2.8\) down to \(2.2\) pairwise). Physical diffusion is carried by the tracked-gradient fluxes, so under both refinement laws the relaxation diffusion of 12 is a second-order-consistent numerical stabilizer (Section 2.2.3). For the hyperbolic and convection-dominated systems the observed second order is thus a property of the resolved smooth solutions, not a claim of intrinsic second-order accuracy in the inviscid limit. In both regimes the clean rate confirms that the compiler correctly deconstructs the declared PDEs, resolves the automatic gradient cascades, and emits consistent kernels.

4.1.0.2 Single Precision

Comparing the single- and double-precision columns isolates the floating-point floor of the method. Single-precision execution is a standard lever for LBM throughput on GPUs [54], [55], so it matters whether the generated schemes retain their convergence order in float. The reference- and equilibrium-shifted formulation (Section 2.4.2) evolves the fluctuation \(\delta\mathbf{Q}\) about the background state directly, so small fluctuations are never formed as the difference of two large numbers. This removes the catastrophic cancellation that would otherwise destroy the single-precision result, and in some cases is what enables convergence in float at all. Two residual effects remain, measured against double precision at identical configuration. First, the per-step rounding errors of slowly evolving fields accumulate into a relative error floor of order \(10^{-4}\) over the \(\sim\!10^{5}\) steps of the finest diffusive grids. Systems whose discretization error stays above it track their double-precision columns (scalar ADR at \(2.55\), Allen–Cahn at \(2.66\)). The compressible NSF system, whose double error dives below it, tracks double through \(N = 256\) and is floor-limited at \(N = 512\). The incompressible NSE density (\(1.58\)) is floor-limited in the same way. Second, single precision is unstable near the over-relaxation limit \(\omega \to 2\) where double precision is not. The affected systems run a uniformly larger relaxation rate in single precision (supplementary materials) and converge cleanly under that single uniform configuration (incompressible NSE momentum at \(2.23\), resistive MHD at \(2.24\)). The acoustic hyperbolic systems, taking only \(\propto N\) steps, erode where the double-precision error is smallest, most visibly the compressible Euler momentum (\(1.62\)).

5 Performance↩︎

In the final code generation layer (cf. Fig. 1) the compiled schemes are emitted as platform-transparent C++ operators for the open-source multiphysics framework OpenLB [6], [23], whose Common Subexpression Elimination (CSE) machinery compacts them algebraically: OpenLB extracts each operator’s expression tree, reduces it in SymPy, and re-emits it as a per-target template specialization [6]. By abstracting memory layout and execution from the continuous equations, this architecture enables compile-time binding to platform-specific data structures, through which OpenLB targets multi-core CPUs and diverse GPU accelerators [5], [6], [56]. As shown in Figure 4, the generated cell-level kernels approach the memory-bandwidth roofline on an exemplary NVIDIA RTX A5000, reaching up to \(96\%\) of peak in single precision. The per-cell byte traffic reported in Table 2 is obtained by introspecting the compiled OpenLB dynamics, counting the populations and field components actually read and written per cell update. The attained throughput is then measured directly from the kernel runtime on a saturating grid, and the saturation column is its ratio to the measured \(\require{physics} \qty{680}{\giga\byte\per\second}\) peak (BabelStream Triad [57]).

2.5pt

The per-kernel resource columns of Table 2 account for the two single-precision bandwidth outliers. The two largest-state kernels, the 18-field resistive MHD and the 14-field HCNSF system, both reach the 255-register ceiling of the sm_86 architecture and spill to local memory (\(32\) and \(104\) bytes per thread). They are precisely the two whose bandwidth saturation falls below the \(85\)–\(97\%\) band, at \(78\%\) and \(63\%\). Occupancy alone does not explain this: the 12-field compressible NSF kernel runs at comparable occupancy (\(19\%\) against \(17\%\)) but, fitting within the register file without spilling, still sustains \(94\%\) of peak. Of the two, the HCNSF system saturates the least: its 14-field state on the D3Q7 lattice (seven populations per field against five on D2Q5) is the largest live-quantity count in the suite, so at the shared register ceiling it spills the hardest, and the uncounted local-memory traffic of that spill displaces the useful bandwidth. In double precision the three largest-state kernels become compute-bound: compressible NSF, resistive MHD, and HCNSF sustain \(345\), \(361\), and \(\require{physics} \qty{363}{\giga\flop\per\second}\) against the \(\require{physics} \qty{434}{\giga\flop\per\second}\) FP64 peak (\(79\)–\(84\%\)), so their low bandwidth-saturation figures (all under \(38\%\)) reflect the \(1{:}64\) FP64 throughput of the architecture rather than a memory-system inefficiency. Of the remaining double-precision points, Burgers and Maxwell stay at the bandwidth roofline and nonlinear elasticity at the FP64 ceiling, while the mid-intensity kernels sit below both roofs.

6 Conclusion↩︎

We introduced a symbolic compiler that automates the translation of systems of continuous conservation laws into optimized, local lattice Boltzmann schemes, work previously restricted to manual and case-specific derivations. Driven by point-wise advection-relaxation cascades, it maps nonlinear physical fluxes directly onto the first-order discrete moments of independent distribution functions and recursively deconstructs higher-order derivatives into auxiliary states, eliminating non-local finite-difference stencils and preserving local concurrency. Asymptotic spatial convergence at or near second order of the generated models was verified by MMS across twelve PDEs spanning fluid dynamics, resistive MHD, phase-field interface dynamics, electromagnetics, and nonlinear finite-strain elasticity. Compile-time arithmetic minimization and the reference-shifted formulation make single precision usable at close to the double-precision convergence order, where the device kernels generated for the platform-transparent OpenLB framework reach up to 96% of the memory-bandwidth roofline.

PDE2LBM is the first framework in which LBMs themselves are derived automatically from continuous conservation laws rather than supplied as input. Within this class, retargeting LBM to a new system is thereby reduced from a manual derivation to a declarative specification, pairing the rapid PDE retargeting that historically distinguished general-purpose finite-volume and finite-element frameworks with the structural locality and efficiency inherent to LBM.

Building on this, boundary condition generation from the PDE-level representation is in development as the prerequisite for validation against physical, non-manufactured benchmarks, as already demonstrated manually for MHD [10].

Code Availability↩︎

The generated schemes are fully specified in print by Sections 2 and 3 together with the expanded component equations, manufactured solutions, and numerical configurations in the supplementary materials. The PDE2LBM compiler is available from the authors upon reasonable request.

Acknowledgement↩︎

Google Gemini and Anthropic Claude assisted in drafting and revising the manuscript. All content was reviewed, verified, and approved by the authors, who take full responsibility for it.

7 Empirical Convergence Profile Plots↩︎

In this appendix, we present the Empirical Order of Convergence (EOC) profiles for all twelve verified PDE targets. Each plot displays the relative \(L^1\), \(L^2\), and \(L^\infty\) error norms of the macroscopic state components across consecutive grid resolutions, in both double and single precision, together with a second-order reference slope. Each system is refined under the lattice scaling matched to its dominant transport mechanism (diffusive for the diffusion-resolving parabolic systems, acoustic for the hyperbolic and convection-dominated ones). The numerical configuration underlying these runs is tabulated in Appendix 8.

8 Numerical Configuration of the Convergence Study↩︎

This appendix documents the complete numerical configuration behind the convergence results reported in the main article, so that every run is reproducible from the tabulated values. All cases share one study harness parametrized by the refinement scaling and three scalar knobs. The time step follows the per-system scaling law, \(\Delta t = f\,\Delta x\) under acoustic and \(\Delta t = f\,\Delta x^2\) under diffusive refinement, with time-step factor \(f\). The relaxation time of the physical fields is parametrized by a resolution-independent rate \(\tau_R\) as \(\tau_{\text{LB}} = \tfrac{1}{2} + \tau_R\,\Delta t\), so the bare relaxation diffusion contracts as \(\varepsilon a^2 = c_s^2\,\tau_R\,\Delta x^2\) under either scaling. The gradient-tracker fields relax on the tracker time \(\tau_\nabla = g_\tau\,\Delta t\), and their populations collide at the fixed lattice rate \(1/(g_\tau + 1)\). The two-dimensional cases run the D2Q5 stencil and the three-dimensional HCNSF case runs D3Q7. Table 3 lists the discrete velocity sets, weights, and lattice speeds of sound. All populations are initialized at the (shifted) equilibrium of the exact manufactured fields at \(t = 0\).

Table 4 lists the per-case values. Every configuration is uniform across the reported resolutions in both precisions: no parameter is retuned per grid, so each fitted convergence order measures a single fixed scheme family.

4pt

9 Domain-Specific Language Example↩︎

Complementing the compressible Euler specification in the main article, Listing [lst:dsl_cns] gives the complete PDE2LBM specification of the compressible Navier–Stokes–Fourier system declared in SI units. Two non-dimensionalization steps are declared. First, characteristic references (the density and the gas constant) eliminate the mass and temperature base dimensions, from which the compiler derives every lattice grid-scaling exponent automatically, with the dynamic viscosity (\(\unit{\pascal\second}\)) and thermal conductivity (\(\unit[per-mode=power]{\watt\per\metre\per\kelvin}\)) reducing to the diffusive (\(\mathrm{length}^2/\mathrm{time}\)) scaling. Second, the realistic-SI magnitudes (density \(\require{physics} \qty{1.225}{\kilogram\per\cubic\metre}\), sound speed \(\require{physics} \qty{340}{\metre\per\second}\)) are declared as characteristic symbols, and the manufactured solution is then authored as (scale) \(\times\) (dimensionless), so the value non-dimensionalization divides the scales out through the products and places the working point at the \(O(1)\) lattice scale.

from pde2lbm import *
from sympy.physics.units import (length, time, velocity, mass,
        temperature, energy, pressure, power, momentum)

cns = ConservationLaws(dim=2)

# 1. States and parameters with named SI dimensions
rho  = cns.register_state("rho",  dim=mass/length**3)        # density
rhou = cns.register_state("rhou", shape=2, dim=momentum/length**3)
E    = cns.register_state("E",    dim=energy/length**3)      # energy

gamma = cns.register_parameter("gamma")
mu = cns.register_parameter("mu", dim=pressure*time)         # Pa s
K  = cns.register_parameter("K",  dim=power/(length*temperature))  # W/(m K)
R  = cns.register_parameter("R",  dim=energy/(mass*temperature))   # gas const

# 2a. Eliminate mass and temperature base dimensions
cns.set_characteristic(mass, rho)
cns.set_characteristic(temperature, R)

# 2b. Value non-dimensionalization: realistic SI working scales as characteristic symbols
rho0 = cns.characteristic_symbol(mass,     1.225, name="rho0")  # air at ~1 atm
c_s  = cns.characteristic_symbol(velocity, 340.0, name="c_s")   # sound speed

# 3. Dependent variables
u = cns.register_derived("u", rhou / rho)       # registers u_0, u_1
p = (gamma - 1)*(E - 0.5*rho*u.dot(u))
T = p/(rho*R)                                   # temperature

# 4. Newtonian stress + Fourier flux
gu  = cgrad(u)
tau = mu*(gu + gu.T - (2.0/3.0)*trace(gu)*eye(2))
q   = -K*cgrad(T)

# 5. Equation system
cns.add([
  Eq(dt(rho)  + div(rhou), 0),
  Eq(dt(rhou) + div(outer(rhou,u) + p*eye(2) - tau), 0),
  Eq(dt(E)    + div((E+p)*u - tau*u + q), 0),
])

# 6. Generate OpenLB solver
cns.compile(class_name="CompressibleNSF")

10 Fully Expanded Compiled Component Equations↩︎

In this appendix, we present the component-level mathematical expressions of the macroscopic state vectors (\(\mathbf{Q}\)), physical flux tensors (\(\mathbf{\Phi}\)), and local source or sink vectors (\(\mathbf{S}\)) as automatically compiled and executed point-wise on the lattice by our symbolic backend. Each element is listed individually. Compound conserved symbols such as \(\rho u_0\) are atomic state components, so ratios like \(\rho u_0^2/\rho\) are already in final form: no cancellation against the separate state \(\rho\) is available or intended.

10.1 Inviscid Burgers Equation↩︎

The 2D continuous system has a state vector of dimension \(N = 1\).

10.1.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = u \end{align}\]

10.1.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= \frac{u^{2}}{2} \\ \Phi_{0,1} &= \frac{u^{2}}{2} \end{align}\]

10.1.0.3 Source Vector Components (\(S_i\))

All components of the source vector are identically zero: \(\mathbf{S}_{\text{burgers}} = \mathbf{0}\).

10.2 Compressible Euler Equations↩︎

The 2D continuous system has a state vector of dimension \(N = 4\).

10.2.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = \rho , \quad Q_{1} = \rho u_0 , \quad Q_{2} = \rho u_1 \\ Q_{3} = E \end{align}\]

10.2.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= \rho u_0 \\ \Phi_{0,1} &= \rho u_1 \\ \Phi_{1,0} &= E \gamma - E - \frac{\gamma \rho u_0^{2}}{2 \rho} - \frac{\gamma \rho u_1^{2}}{2 \rho} + \frac{3 \rho u_0^{2}}{2 \rho} \\ &\quad + \frac{\rho u_1^{2}}{2 \rho} \\ \Phi_{1,1} &= \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,0} &= \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,1} &= E \gamma - E - \frac{\gamma \rho u_0^{2}}{2 \rho} - \frac{\gamma \rho u_1^{2}}{2 \rho} + \frac{\rho u_0^{2}}{2 \rho} \\ &\quad + \frac{3 \rho u_1^{2}}{2 \rho} \\ \Phi_{3,0} &= \frac{E \gamma \rho u_0}{\rho} - \frac{\gamma \rho u_0^{3}}{2 \rho^{2}} - \frac{\gamma \rho u_0 \rho u_1^{2}}{2 \rho^{2}} \\ &\quad + \frac{\rho u_0^{3}}{2 \rho^{2}} + \frac{\rho u_0 \rho u_1^{2}}{2 \rho^{2}} \\ \Phi_{3,1} &= \frac{E \gamma \rho u_1}{\rho} - \frac{\gamma \rho u_0^{2} \rho u_1}{2 \rho^{2}} - \frac{\gamma \rho u_1^{3}}{2 \rho^{2}} \\ &\quad + \frac{\rho u_0^{2} \rho u_1}{2 \rho^{2}} + \frac{\rho u_1^{3}}{2 \rho^{2}} \end{align}\]

10.2.0.3 Source Vector Components (\(S_i\))

All components of the source vector are identically zero: \(\mathbf{S}_{\text{compressible\_euler}} = \mathbf{0}\).

10.3 Shallow Water Equations↩︎

The 2D continuous system has a state vector of dimension \(N = 3\).

10.3.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = h , \quad Q_{1} = h u_0 , \quad Q_{2} = h u_1 \end{align}\]

10.3.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= h u_0 \\ \Phi_{0,1} &= h u_1 \\ \Phi_{1,0} &= \frac{g h^{2}}{2} + \frac{h u_0^{2}}{h} \\ \Phi_{1,1} &= \frac{h u_0 h u_1}{h} \\ \Phi_{2,0} &= \frac{h u_0 h u_1}{h} \\ \Phi_{2,1} &= \frac{g h^{2}}{2} + \frac{h u_1^{2}}{h} \end{align}\]

10.3.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= 0 \\ S_{1} &= - \frac{\pi g h \cos{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{50} \\ S_{2} &= \frac{\pi g h \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{50} \end{align}\]

10.4 Ideal Ultrarelativistic Fluid↩︎

The 2D continuous system has a state vector of dimension \(N = 3\).

10.4.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = E , \quad Q_{1} = S_0 , \quad Q_{2} = S_1 \end{align}\]

10.4.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= S_0 \\ \Phi_{0,1} &= S_1 \\ \Phi_{1,0} &= - \frac{E}{3} + \frac{3 S_0^{2}}{2 E + \sqrt{4 E^{2} - 3 S_0^{2} - 3 S_1^{2}}} + \frac{\sqrt{4 E^{2} - 3 S_0^{2} - 3 S_1^{2}}}{3} \\ \Phi_{1,1} &= \frac{3 S_0 S_1}{2 E + \sqrt{4 E^{2} - 3 S_0^{2} - 3 S_1^{2}}} \\ \Phi_{2,0} &= \frac{3 S_0 S_1}{2 E + \sqrt{4 E^{2} - 3 S_0^{2} - 3 S_1^{2}}} \\ \Phi_{2,1} &= - \frac{E}{3} + \frac{3 S_1^{2}}{2 E + \sqrt{4 E^{2} - 3 S_0^{2} - 3 S_1^{2}}} + \frac{\sqrt{4 E^{2} - 3 S_0^{2} - 3 S_1^{2}}}{3} \end{align}\]

10.4.0.3 Source Vector Components (\(S_i\))

All components of the source vector are identically zero: \(\mathbf{S}_{\text{relativistic\_hydro}} = \mathbf{0}\).

10.5 Maxwell Electromagnetics (TM Mode)↩︎

The 2D continuous system has a state vector of dimension \(N = 3\).

10.5.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = E_z , \quad Q_{1} = H_x , \quad Q_{2} = H_y \end{align}\]

10.5.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= - H_y c \\ \Phi_{0,1} &= H_x c \\ \Phi_{1,0} &= 0 \\ \Phi_{1,1} &= E_z c \\ \Phi_{2,0} &= - E_z c \\ \Phi_{2,1} &= 0 \end{align}\]

10.5.0.3 Source Vector Components (\(S_i\))

All components of the source vector are identically zero: \(\mathbf{S}_{\text{maxwell}} = \mathbf{0}\).

10.6 Nonlinear Finite-Strain Elasticity↩︎

The 2D continuous system has a state vector of dimension \(N = 6\).

10.6.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = F_{00} , \quad Q_{1} = F_{01} , \quad Q_{2} = F_{10} \\ Q_{3} = F_{11} , \quad Q_{4} = v_0 , \quad Q_{5} = v_1 \end{align}\]

10.6.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

with the shared subexpressions \[\begin{align} \chi_{0} &= - v_0 \\ \chi_{1} &= - v_1 \\ \chi_{2} &= \frac{\mu_s}{\rho_{0}} \\ \chi_{3} &= F_{00} F_{11} \\ \chi_{4} &= F_{01} F_{10} \\ \chi_{5} &= \frac{1}{\chi_{3} \rho_{0} - \chi_{4} \rho_{0}} \\ \chi_{6} &= \log{(\chi_{3} - \chi_{4} )} \\ \chi_{7} &= \chi_{5} \chi_{6} \lambda_s \\ \chi_{8} &= \chi_{5} \mu_s \end{align}\] the flux components evaluate to \[\begin{align} \Phi_{0,0} &= \chi_{0} \\ \Phi_{0,1} &= 0 \\ \Phi_{1,0} &= 0 \\ \Phi_{1,1} &= \chi_{0} \\ \Phi_{2,0} &= \chi_{1} \\ \Phi_{2,1} &= 0 \\ \Phi_{3,0} &= 0 \\ \Phi_{3,1} &= \chi_{1} \\ \Phi_{4,0} &= - F_{00} \chi_{2} + F_{11} \chi_{5} \mu_s - F_{11} \chi_{7} \\ \Phi_{4,1} &= - F_{01} \chi_{2} + F_{10} \chi_{5} \chi_{6} \lambda_s - F_{10} \chi_{8} \\ \Phi_{5,0} &= F_{01} \chi_{5} \chi_{6} \lambda_s - F_{01} \chi_{8} - F_{10} \chi_{2} \\ \Phi_{5,1} &= F_{00} \chi_{5} \mu_s - F_{00} \chi_{7} - F_{11} \chi_{2} \end{align}\]

10.6.0.3 Source Vector Components (\(S_i\))

All components of the source vector are identically zero: \(\mathbf{S}_{\text{nonlinear\_elasticity}} = \mathbf{0}\).

10.7 Scalar Advection-Diffusion-Reaction↩︎

The 2D continuous system has a state vector of dimension \(N = 3\).

10.7.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = c , \quad Q_{1} = G_{c,0} , \quad Q_{2} = G_{c,1} \end{align}\]

10.7.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= - D G_{c,0} + c u_0 \\ \Phi_{0,1} &= - D G_{c,1} + c u_1 \\ \Phi_{1,0} &= G_{c,0} u_0 - \frac{c}{\tau_{\nabla}} \\ \Phi_{1,1} &= G_{c,0} u_1 \\ \Phi_{2,0} &= G_{c,1} u_0 \\ \Phi_{2,1} &= G_{c,1} u_1 - \frac{c}{\tau_{\nabla}} \end{align}\]

10.7.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= - R c^{2} + R c \\ S_{1} &= - \frac{G_{c,0}}{\tau_{\nabla}} \\ S_{2} &= - \frac{G_{c,1}}{\tau_{\nabla}} \end{align}\]

10.8 Allen–Cahn Phase-Field Equation↩︎

The 2D continuous system has a state vector of dimension \(N = 3\).

10.8.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = \phi , \quad Q_{1} = G_{\phi,0} , \quad Q_{2} = G_{\phi,1} \end{align}\]

10.8.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= - G_{\phi,0} \gamma \\ \Phi_{0,1} &= - G_{\phi,1} \gamma \\ \Phi_{1,0} &= - \frac{\phi}{\tau_{\nabla}} \\ \Phi_{1,1} &= 0 \\ \Phi_{2,0} &= 0 \\ \Phi_{2,1} &= - \frac{\phi}{\tau_{\nabla}} \end{align}\]

10.8.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= - \phi^{3} r + \phi r \\ S_{1} &= - \frac{G_{\phi,0}}{\tau_{\nabla}} \\ S_{2} &= - \frac{G_{\phi,1}}{\tau_{\nabla}} \end{align}\]

10.9 Weakly Compressible Navier–Stokes Equations (Incompressible Limit)↩︎

The 2D continuous system has a state vector of dimension \(N = 7\).

10.9.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = \rho , \quad Q_{1} = \rho u_0 , \quad Q_{2} = \rho u_1 \\ Q_{3} = G_{\rho u, 0, 0} , \quad Q_{4} = G_{\rho u, 0, 1} , \quad Q_{5} = G_{\rho u, 1, 0} \\ Q_{6} = G_{\rho u, 1, 1} \end{align}\]

10.9.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= \rho u_0 \\ \Phi_{0,1} &= \rho u_1 \\ \Phi_{1,0} &= - \frac{2 G_{\rho u, 0, 0} \mu}{\rho} + c_{s,\text{phys}}^{2} \rho + \frac{\rho u_0^{2}}{\rho} \\ \Phi_{1,1} &= - \frac{G_{\rho u, 0, 1} \mu}{\rho} - \frac{G_{\rho u, 1, 0} \mu}{\rho} + \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,0} &= - \frac{G_{\rho u, 0, 1} \mu}{\rho} - \frac{G_{\rho u, 1, 0} \mu}{\rho} + \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,1} &= - \frac{2 G_{\rho u, 1, 1} \mu}{\rho} + c_{s,\text{phys}}^{2} \rho + \frac{\rho u_1^{2}}{\rho} \\ \Phi_{3,0} &= \frac{G_{\rho u, 0, 0} \rho u_0}{\rho} - \frac{\rho u_0}{\tau_{\nabla}} \\ \Phi_{3,1} &= \frac{G_{\rho u, 0, 0} \rho u_1}{\rho} \\ \Phi_{4,0} &= \frac{G_{\rho u, 0, 1} \rho u_0}{\rho} \\ \Phi_{4,1} &= \frac{G_{\rho u, 0, 1} \rho u_1}{\rho} - \frac{\rho u_0}{\tau_{\nabla}} \\ \Phi_{5,0} &= \frac{G_{\rho u, 1, 0} \rho u_0}{\rho} - \frac{\rho u_1}{\tau_{\nabla}} \\ \Phi_{5,1} &= \frac{G_{\rho u, 1, 0} \rho u_1}{\rho} \\ \Phi_{6,0} &= \frac{G_{\rho u, 1, 1} \rho u_0}{\rho} \\ \Phi_{6,1} &= \frac{G_{\rho u, 1, 1} \rho u_1}{\rho} - \frac{\rho u_1}{\tau_{\nabla}} \end{align}\]

10.9.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= 0 \\ S_{1} &= 0 \\ S_{2} &= 0 \\ S_{3} &= - \frac{G_{\rho u, 0, 0}}{\tau_{\nabla}} \\ S_{4} &= - \frac{G_{\rho u, 0, 1}}{\tau_{\nabla}} \\ S_{5} &= - \frac{G_{\rho u, 1, 0}}{\tau_{\nabla}} \\ S_{6} &= - \frac{G_{\rho u, 1, 1}}{\tau_{\nabla}} \end{align}\]

10.10 Compressible Navier–Stokes–Fourier Equations↩︎

The 2D continuous system has a state vector of dimension \(N = 12\).

10.10.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = \rho , \quad Q_{1} = \rho u_0 , \quad Q_{2} = \rho u_1 \\ Q_{3} = E , \quad Q_{4} = G_{\rho u, 1, 0} , \quad Q_{5} = G_{\rho u, 1, 1} \\ Q_{6} = G_{E,0} , \quad Q_{7} = G_{E,1} , \quad Q_{8} = G_{\rho,0} \\ Q_{9} = G_{\rho,1} , \quad Q_{10} = G_{\rho u, 0, 0} , \quad Q_{11} = G_{\rho u, 0, 1} \end{align}\]

10.10.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= \rho u_0 \\ \Phi_{0,1} &= \rho u_1 \\ \Phi_{1,0} &= E \gamma - E + \frac{4 G_{\rho,0} \mu \rho u_0}{3 \rho^{2}} - \frac{2 G_{\rho,1} \mu \rho u_1}{3 \rho^{2}} \\ &\quad - \frac{4 G_{\rho u, 0, 0} \mu}{3 \rho} + \frac{2 G_{\rho u, 1, 1} \mu}{3 \rho} - \frac{\gamma \rho u_0^{2}}{2 \rho} \\ &\quad - \frac{\gamma \rho u_1^{2}}{2 \rho} + \frac{3 \rho u_0^{2}}{2 \rho} + \frac{\rho u_1^{2}}{2 \rho} \\ \Phi_{1,1} &= \frac{G_{\rho,0} \mu \rho u_1}{\rho^{2}} + \frac{G_{\rho,1} \mu \rho u_0}{\rho^{2}} - \frac{G_{\rho u, 0, 1} \mu}{\rho} \\ &\quad - \frac{G_{\rho u, 1, 0} \mu}{\rho} + \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,0} &= \frac{G_{\rho,0} \mu \rho u_1}{\rho^{2}} + \frac{G_{\rho,1} \mu \rho u_0}{\rho^{2}} - \frac{G_{\rho u, 0, 1} \mu}{\rho} \\ &\quad - \frac{G_{\rho u, 1, 0} \mu}{\rho} + \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,1} &= E \gamma - E - \frac{2 G_{\rho,0} \mu \rho u_0}{3 \rho^{2}} + \frac{4 G_{\rho,1} \mu \rho u_1}{3 \rho^{2}} \\ &\quad + \frac{2 G_{\rho u, 0, 0} \mu}{3 \rho} - \frac{4 G_{\rho u, 1, 1} \mu}{3 \rho} - \frac{\gamma \rho u_0^{2}}{2 \rho} \\ &\quad - \frac{\gamma \rho u_1^{2}}{2 \rho} + \frac{\rho u_0^{2}}{2 \rho} + \frac{3 \rho u_1^{2}}{2 \rho} \\ \Phi_{3,0} &= \frac{E G_{\rho,0} K \gamma}{R \rho^{2}} - \frac{E G_{\rho,0} K}{R \rho^{2}} + \frac{E \gamma \rho u_0}{\rho} \\ &\quad - \frac{G_{E,0} K \gamma}{R \rho} + \frac{G_{E,0} K}{R \rho} - \frac{G_{\rho,0} K \gamma \rho u_0^{2}}{R \rho^{3}} \\ &\quad - \frac{G_{\rho,0} K \gamma \rho u_1^{2}}{R \rho^{3}} + \frac{G_{\rho,0} K \rho u_0^{2}}{R \rho^{3}} \\ &\quad + \frac{G_{\rho,0} K \rho u_1^{2}}{R \rho^{3}} + \frac{4 G_{\rho,0} \mu \rho u_0^{2}}{3 \rho^{3}} \\ &\quad + \frac{G_{\rho,0} \mu \rho u_1^{2}}{\rho^{3}} + \frac{G_{\rho,1} \mu \rho u_0 \rho u_1}{3 \rho^{3}} \\ &\quad + \frac{G_{\rho u, 0, 0} K \gamma \rho u_0}{R \rho^{2}} - \frac{G_{\rho u, 0, 0} K \rho u_0}{R \rho^{2}} \\ &\quad - \frac{4 G_{\rho u, 0, 0} \mu \rho u_0}{3 \rho^{2}} - \frac{G_{\rho u, 0, 1} \mu \rho u_1}{\rho^{2}} \\ &\quad + \frac{G_{\rho u, 1, 0} K \gamma \rho u_1}{R \rho^{2}} - \frac{G_{\rho u, 1, 0} K \rho u_1}{R \rho^{2}} \\ &\quad - \frac{G_{\rho u, 1, 0} \mu \rho u_1}{\rho^{2}} + \frac{2 G_{\rho u, 1, 1} \mu \rho u_0}{3 \rho^{2}} \\ &\quad - \frac{\gamma \rho u_0^{3}}{2 \rho^{2}} - \frac{\gamma \rho u_0 \rho u_1^{2}}{2 \rho^{2}} + \frac{\rho u_0^{3}}{2 \rho^{2}} \\ &\quad + \frac{\rho u_0 \rho u_1^{2}}{2 \rho^{2}} \\ \Phi_{3,1} &= \frac{E G_{\rho,1} K \gamma}{R \rho^{2}} - \frac{E G_{\rho,1} K}{R \rho^{2}} + \frac{E \gamma \rho u_1}{\rho} \\ &\quad - \frac{G_{E,1} K \gamma}{R \rho} + \frac{G_{E,1} K}{R \rho} + \frac{G_{\rho,0} \mu \rho u_0 \rho u_1}{3 \rho^{3}} \\ &\quad - \frac{G_{\rho,1} K \gamma \rho u_0^{2}}{R \rho^{3}} - \frac{G_{\rho,1} K \gamma \rho u_1^{2}}{R \rho^{3}} \\ &\quad + \frac{G_{\rho,1} K \rho u_0^{2}}{R \rho^{3}} + \frac{G_{\rho,1} K \rho u_1^{2}}{R \rho^{3}} \\ &\quad + \frac{G_{\rho,1} \mu \rho u_0^{2}}{\rho^{3}} + \frac{4 G_{\rho,1} \mu \rho u_1^{2}}{3 \rho^{3}} \\ &\quad + \frac{2 G_{\rho u, 0, 0} \mu \rho u_1}{3 \rho^{2}} + \frac{G_{\rho u, 0, 1} K \gamma \rho u_0}{R \rho^{2}} \\ &\quad - \frac{G_{\rho u, 0, 1} K \rho u_0}{R \rho^{2}} - \frac{G_{\rho u, 0, 1} \mu \rho u_0}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 1, 0} \mu \rho u_0}{\rho^{2}} + \frac{G_{\rho u, 1, 1} K \gamma \rho u_1}{R \rho^{2}} \\ &\quad - \frac{G_{\rho u, 1, 1} K \rho u_1}{R \rho^{2}} - \frac{4 G_{\rho u, 1, 1} \mu \rho u_1}{3 \rho^{2}} \\ &\quad - \frac{\gamma \rho u_0^{2} \rho u_1}{2 \rho^{2}} - \frac{\gamma \rho u_1^{3}}{2 \rho^{2}} \\ &\quad + \frac{\rho u_0^{2} \rho u_1}{2 \rho^{2}} + \frac{\rho u_1^{3}}{2 \rho^{2}} \\ \Phi_{4,0} &= \frac{G_{\rho u, 1, 0} \rho u_0}{\rho} - \frac{\rho u_1}{\tau_{\nabla}} \\ \Phi_{4,1} &= \frac{G_{\rho u, 1, 0} \rho u_1}{\rho} \\ \Phi_{5,0} &= \frac{G_{\rho u, 1, 1} \rho u_0}{\rho} \\ \Phi_{5,1} &= \frac{G_{\rho u, 1, 1} \rho u_1}{\rho} - \frac{\rho u_1}{\tau_{\nabla}} \\ \Phi_{6,0} &= - \frac{E}{\tau_{\nabla}} + \frac{G_{E,0} \rho u_0}{\rho} \\ \Phi_{6,1} &= \frac{G_{E,0} \rho u_1}{\rho} \\ \Phi_{7,0} &= \frac{G_{E,1} \rho u_0}{\rho} \\ \Phi_{7,1} &= - \frac{E}{\tau_{\nabla}} + \frac{G_{E,1} \rho u_1}{\rho} \\ \Phi_{8,0} &= \frac{G_{\rho,0} \rho u_0}{\rho} - \frac{\rho}{\tau_{\nabla}} \\ \Phi_{8,1} &= \frac{G_{\rho,0} \rho u_1}{\rho} \\ \Phi_{9,0} &= \frac{G_{\rho,1} \rho u_0}{\rho} \\ \Phi_{9,1} &= \frac{G_{\rho,1} \rho u_1}{\rho} - \frac{\rho}{\tau_{\nabla}} \\ \Phi_{10,0} &= \frac{G_{\rho u, 0, 0} \rho u_0}{\rho} - \frac{\rho u_0}{\tau_{\nabla}} \\ \Phi_{10,1} &= \frac{G_{\rho u, 0, 0} \rho u_1}{\rho} \\ \Phi_{11,0} &= \frac{G_{\rho u, 0, 1} \rho u_0}{\rho} \\ \Phi_{11,1} &= \frac{G_{\rho u, 0, 1} \rho u_1}{\rho} - \frac{\rho u_0}{\tau_{\nabla}} \end{align}\]

10.10.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= 0 \\ S_{1} &= 0 \\ S_{2} &= 0 \\ S_{3} &= 0 \\ S_{4} &= - \frac{G_{\rho u, 1, 0}}{\tau_{\nabla}} \\ S_{5} &= - \frac{G_{\rho u, 1, 1}}{\tau_{\nabla}} \\ S_{6} &= - \frac{G_{E,0}}{\tau_{\nabla}} \\ S_{7} &= - \frac{G_{E,1}}{\tau_{\nabla}} \\ S_{8} &= - \frac{G_{\rho,0}}{\tau_{\nabla}} \\ S_{9} &= - \frac{G_{\rho,1}}{\tau_{\nabla}} \\ S_{10} &= - \frac{G_{\rho u, 0, 0}}{\tau_{\nabla}} \\ S_{11} &= - \frac{G_{\rho u, 0, 1}}{\tau_{\nabla}} \end{align}\]

10.11 Resistive Compressible Magnetohydrodynamics↩︎

The 2D continuous system has a state vector of dimension \(N = 18\).

10.11.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = \rho , \quad Q_{1} = \rho u_0 , \quad Q_{2} = \rho u_1 \\ Q_{3} = E , \quad Q_{4} = B_0 , \quad Q_{5} = B_1 \\ Q_{6} = G_{\rho,0} , \quad Q_{7} = G_{\rho,1} , \quad Q_{8} = G_{\rho u, 1, 0} \\ Q_{9} = G_{\rho u, 1, 1} , \quad Q_{10} = G_{E,0} , \quad Q_{11} = G_{E,1} \\ Q_{12} = G_{B, 1, 0} , \quad Q_{13} = G_{B, 1, 1} , \quad Q_{14} = G_{\rho u, 0, 0} \\ Q_{15} = G_{\rho u, 0, 1} , \quad Q_{16} = G_{B, 0, 0} , \quad Q_{17} = G_{B, 0, 1} \end{align}\]

10.11.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

\[\begin{align} \Phi_{0,0} &= \rho u_0 \\ \Phi_{0,1} &= \rho u_1 \\ \Phi_{1,0} &= - \frac{B_0^{2} \gamma}{2} - \frac{B_1^{2} \gamma}{2} + B_1^{2} + E \gamma - E + \frac{4 G_{\rho,0} \mu \rho u_0}{3 \rho^{2}} \\ &\quad - \frac{2 G_{\rho,1} \mu \rho u_1}{3 \rho^{2}} - \frac{4 G_{\rho u, 0, 0} \mu}{3 \rho} + \frac{2 G_{\rho u, 1, 1} \mu}{3 \rho} \\ &\quad - \frac{\gamma \rho u_0^{2}}{2 \rho} - \frac{\gamma \rho u_1^{2}}{2 \rho} + \frac{3 \rho u_0^{2}}{2 \rho} \\ &\quad + \frac{\rho u_1^{2}}{2 \rho} \\ \Phi_{1,1} &= - B_0 B_1 + \frac{G_{\rho,0} \mu \rho u_1}{\rho^{2}} + \frac{G_{\rho,1} \mu \rho u_0}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 0, 1} \mu}{\rho} - \frac{G_{\rho u, 1, 0} \mu}{\rho} + \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,0} &= - B_0 B_1 + \frac{G_{\rho,0} \mu \rho u_1}{\rho^{2}} + \frac{G_{\rho,1} \mu \rho u_0}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 0, 1} \mu}{\rho} - \frac{G_{\rho u, 1, 0} \mu}{\rho} + \frac{\rho u_0 \rho u_1}{\rho} \\ \Phi_{2,1} &= - \frac{B_0^{2} \gamma}{2} + B_0^{2} - \frac{B_1^{2} \gamma}{2} + E \gamma - E - \frac{2 G_{\rho,0} \mu \rho u_0}{3 \rho^{2}} \\ &\quad + \frac{4 G_{\rho,1} \mu \rho u_1}{3 \rho^{2}} + \frac{2 G_{\rho u, 0, 0} \mu}{3 \rho} - \frac{4 G_{\rho u, 1, 1} \mu}{3 \rho} \\ &\quad - \frac{\gamma \rho u_0^{2}}{2 \rho} - \frac{\gamma \rho u_1^{2}}{2 \rho} + \frac{\rho u_0^{2}}{2 \rho} \\ &\quad + \frac{3 \rho u_1^{2}}{2 \rho} \\ \Phi_{3,0} &= - \frac{B_0^{2} G_{\rho,0} K \gamma}{2 \rho^{2}} + \frac{B_0^{2} G_{\rho,0} K}{2 \rho^{2}} \\ &\quad - \frac{B_0^{2} \gamma \rho u_0}{2 \rho} - \frac{B_0 B_1 \rho u_1}{\rho} + \frac{B_0 G_{B, 0, 0} K \gamma}{\rho} \\ &\quad - \frac{B_0 G_{B, 0, 0} K}{\rho} - \frac{B_1^{2} G_{\rho,0} K \gamma}{2 \rho^{2}} + \frac{B_1^{2} G_{\rho,0} K}{2 \rho^{2}} \\ &\quad - \frac{B_1^{2} \gamma \rho u_0}{2 \rho} + \frac{B_1^{2} \rho u_0}{\rho} + B_1 G_{B, 0, 1} \eta \\ &\quad + \frac{B_1 G_{B, 1, 0} K \gamma}{\rho} - \frac{B_1 G_{B, 1, 0} K}{\rho} - B_1 G_{B, 1, 0} \eta \\ &\quad + \frac{E G_{\rho,0} K \gamma}{\rho^{2}} - \frac{E G_{\rho,0} K}{\rho^{2}} + \frac{E \gamma \rho u_0}{\rho} \\ &\quad - \frac{G_{E,0} K \gamma}{\rho} + \frac{G_{E,0} K}{\rho} - \frac{G_{\rho,0} K \gamma \rho u_0^{2}}{\rho^{3}} \\ &\quad - \frac{G_{\rho,0} K \gamma \rho u_1^{2}}{\rho^{3}} + \frac{G_{\rho,0} K \rho u_0^{2}}{\rho^{3}} \\ &\quad + \frac{G_{\rho,0} K \rho u_1^{2}}{\rho^{3}} + \frac{4 G_{\rho,0} \mu \rho u_0^{2}}{3 \rho^{3}} \\ &\quad + \frac{G_{\rho,0} \mu \rho u_1^{2}}{\rho^{3}} + \frac{G_{\rho,1} \mu \rho u_0 \rho u_1}{3 \rho^{3}} \\ &\quad + \frac{G_{\rho u, 0, 0} K \gamma \rho u_0}{\rho^{2}} - \frac{G_{\rho u, 0, 0} K \rho u_0}{\rho^{2}} \\ &\quad - \frac{4 G_{\rho u, 0, 0} \mu \rho u_0}{3 \rho^{2}} - \frac{G_{\rho u, 0, 1} \mu \rho u_1}{\rho^{2}} \\ &\quad + \frac{G_{\rho u, 1, 0} K \gamma \rho u_1}{\rho^{2}} - \frac{G_{\rho u, 1, 0} K \rho u_1}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 1, 0} \mu \rho u_1}{\rho^{2}} + \frac{2 G_{\rho u, 1, 1} \mu \rho u_0}{3 \rho^{2}} \\ &\quad - \frac{\gamma \rho u_0^{3}}{2 \rho^{2}} - \frac{\gamma \rho u_0 \rho u_1^{2}}{2 \rho^{2}} + \frac{\rho u_0^{3}}{2 \rho^{2}} \\ &\quad + \frac{\rho u_0 \rho u_1^{2}}{2 \rho^{2}} \\ \Phi_{3,1} &= - \frac{B_0^{2} G_{\rho,1} K \gamma}{2 \rho^{2}} + \frac{B_0^{2} G_{\rho,1} K}{2 \rho^{2}} \\ &\quad - \frac{B_0^{2} \gamma \rho u_1}{2 \rho} + \frac{B_0^{2} \rho u_1}{\rho} - \frac{B_0 B_1 \rho u_0}{\rho} \\ &\quad + \frac{B_0 G_{B, 0, 1} K \gamma}{\rho} - \frac{B_0 G_{B, 0, 1} K}{\rho} - B_0 G_{B, 0, 1} \eta + B_0 G_{B, 1, 0} \eta \\ &\quad - \frac{B_1^{2} G_{\rho,1} K \gamma}{2 \rho^{2}} + \frac{B_1^{2} G_{\rho,1} K}{2 \rho^{2}} \\ &\quad - \frac{B_1^{2} \gamma \rho u_1}{2 \rho} + \frac{B_1 G_{B, 1, 1} K \gamma}{\rho} - \frac{B_1 G_{B, 1, 1} K}{\rho} \\ &\quad + \frac{E G_{\rho,1} K \gamma}{\rho^{2}} - \frac{E G_{\rho,1} K}{\rho^{2}} + \frac{E \gamma \rho u_1}{\rho} \\ &\quad - \frac{G_{E,1} K \gamma}{\rho} + \frac{G_{E,1} K}{\rho} + \frac{G_{\rho,0} \mu \rho u_0 \rho u_1}{3 \rho^{3}} \\ &\quad - \frac{G_{\rho,1} K \gamma \rho u_0^{2}}{\rho^{3}} - \frac{G_{\rho,1} K \gamma \rho u_1^{2}}{\rho^{3}} \\ &\quad + \frac{G_{\rho,1} K \rho u_0^{2}}{\rho^{3}} + \frac{G_{\rho,1} K \rho u_1^{2}}{\rho^{3}} \\ &\quad + \frac{G_{\rho,1} \mu \rho u_0^{2}}{\rho^{3}} + \frac{4 G_{\rho,1} \mu \rho u_1^{2}}{3 \rho^{3}} \\ &\quad + \frac{2 G_{\rho u, 0, 0} \mu \rho u_1}{3 \rho^{2}} + \frac{G_{\rho u, 0, 1} K \gamma \rho u_0}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 0, 1} K \rho u_0}{\rho^{2}} - \frac{G_{\rho u, 0, 1} \mu \rho u_0}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 1, 0} \mu \rho u_0}{\rho^{2}} + \frac{G_{\rho u, 1, 1} K \gamma \rho u_1}{\rho^{2}} \\ &\quad - \frac{G_{\rho u, 1, 1} K \rho u_1}{\rho^{2}} - \frac{4 G_{\rho u, 1, 1} \mu \rho u_1}{3 \rho^{2}} \\ &\quad - \frac{\gamma \rho u_0^{2} \rho u_1}{2 \rho^{2}} - \frac{\gamma \rho u_1^{3}}{2 \rho^{2}} \\ &\quad + \frac{\rho u_0^{2} \rho u_1}{2 \rho^{2}} + \frac{\rho u_1^{3}}{2 \rho^{2}} \\ \Phi_{4,0} &= - G_{B, 0, 0} \eta \\ \Phi_{4,1} &= \frac{B_0 \rho u_1}{\rho} - \frac{B_1 \rho u_0}{\rho} - G_{B, 0, 1} \eta \\ \Phi_{5,0} &= - \frac{B_0 \rho u_1}{\rho} + \frac{B_1 \rho u_0}{\rho} - G_{B, 1, 0} \eta \\ \Phi_{5,1} &= - G_{B, 1, 1} \eta \\ \Phi_{6,0} &= \frac{G_{\rho,0} \rho u_0}{\rho} - \frac{\rho}{\tau_{\nabla}} \\ \Phi_{6,1} &= \frac{G_{\rho,0} \rho u_1}{\rho} \\ \Phi_{7,0} &= \frac{G_{\rho,1} \rho u_0}{\rho} \\ \Phi_{7,1} &= \frac{G_{\rho,1} \rho u_1}{\rho} - \frac{\rho}{\tau_{\nabla}} \\ \Phi_{8,0} &= \frac{G_{\rho u, 1, 0} \rho u_0}{\rho} - \frac{\rho u_1}{\tau_{\nabla}} \\ \Phi_{8,1} &= \frac{G_{\rho u, 1, 0} \rho u_1}{\rho} \\ \Phi_{9,0} &= \frac{G_{\rho u, 1, 1} \rho u_0}{\rho} \\ \Phi_{9,1} &= \frac{G_{\rho u, 1, 1} \rho u_1}{\rho} - \frac{\rho u_1}{\tau_{\nabla}} \\ \Phi_{10,0} &= - \frac{E}{\tau_{\nabla}} + \frac{G_{E,0} \rho u_0}{\rho} \\ \Phi_{10,1} &= \frac{G_{E,0} \rho u_1}{\rho} \\ \Phi_{11,0} &= \frac{G_{E,1} \rho u_0}{\rho} \\ \Phi_{11,1} &= - \frac{E}{\tau_{\nabla}} + \frac{G_{E,1} \rho u_1}{\rho} \\ \Phi_{12,0} &= - \frac{B_1}{\tau_{\nabla}} + \frac{G_{B, 1, 0} \rho u_0}{\rho} \\ \Phi_{12,1} &= \frac{G_{B, 1, 0} \rho u_1}{\rho} \\ \Phi_{13,0} &= \frac{G_{B, 1, 1} \rho u_0}{\rho} \\ \Phi_{13,1} &= - \frac{B_1}{\tau_{\nabla}} + \frac{G_{B, 1, 1} \rho u_1}{\rho} \\ \Phi_{14,0} &= \frac{G_{\rho u, 0, 0} \rho u_0}{\rho} - \frac{\rho u_0}{\tau_{\nabla}} \\ \Phi_{14,1} &= \frac{G_{\rho u, 0, 0} \rho u_1}{\rho} \\ \Phi_{15,0} &= \frac{G_{\rho u, 0, 1} \rho u_0}{\rho} \\ \Phi_{15,1} &= \frac{G_{\rho u, 0, 1} \rho u_1}{\rho} - \frac{\rho u_0}{\tau_{\nabla}} \\ \Phi_{16,0} &= - \frac{B_0}{\tau_{\nabla}} + \frac{G_{B, 0, 0} \rho u_0}{\rho} \\ \Phi_{16,1} &= \frac{G_{B, 0, 0} \rho u_1}{\rho} \\ \Phi_{17,0} &= \frac{G_{B, 0, 1} \rho u_0}{\rho} \\ \Phi_{17,1} &= - \frac{B_0}{\tau_{\nabla}} + \frac{G_{B, 0, 1} \rho u_1}{\rho} \end{align}\]

10.11.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= 0 \\ S_{1} &= - B_0 G_{B, 0, 0} - B_0 G_{B, 1, 1} \\ S_{2} &= - B_1 G_{B, 0, 0} - B_1 G_{B, 1, 1} \\ S_{3} &= - \frac{B_0 G_{B, 0, 0} \rho u_0}{\rho} - \frac{B_0 G_{B, 1, 1} \rho u_0}{\rho} - \frac{B_1 G_{B, 0, 0} \rho u_1}{\rho} \\ &\quad - \frac{B_1 G_{B, 1, 1} \rho u_1}{\rho} \\ S_{4} &= - \frac{G_{B, 0, 0} \rho u_0}{\rho} - \frac{G_{B, 1, 1} \rho u_0}{\rho} \\ S_{5} &= - \frac{G_{B, 0, 0} \rho u_1}{\rho} - \frac{G_{B, 1, 1} \rho u_1}{\rho} \\ S_{6} &= - \frac{G_{\rho,0}}{\tau_{\nabla}} \\ S_{7} &= - \frac{G_{\rho,1}}{\tau_{\nabla}} \\ S_{8} &= - \frac{G_{\rho u, 1, 0}}{\tau_{\nabla}} \\ S_{9} &= - \frac{G_{\rho u, 1, 1}}{\tau_{\nabla}} \\ S_{10} &= - \frac{G_{E,0}}{\tau_{\nabla}} \\ S_{11} &= - \frac{G_{E,1}}{\tau_{\nabla}} \\ S_{12} &= - \frac{G_{B, 1, 0}}{\tau_{\nabla}} \\ S_{13} &= - \frac{G_{B, 1, 1}}{\tau_{\nabla}} \\ S_{14} &= - \frac{G_{\rho u, 0, 0}}{\tau_{\nabla}} \\ S_{15} &= - \frac{G_{\rho u, 0, 1}}{\tau_{\nabla}} \\ S_{16} &= - \frac{G_{B, 0, 0}}{\tau_{\nabla}} \\ S_{17} &= - \frac{G_{B, 0, 1}}{\tau_{\nabla}} \end{align}\]

10.12 3D Homogenized Compressible Navier–Stokes–Fourier Equations↩︎

The 3D continuous system has a state vector of dimension \(N = 14\).

10.12.0.1 Macroscopic State Components (\(Q_i\))

\[\begin{align} Q_{0} = \rho , \quad Q_{1} = \rho u_0 , \quad Q_{2} = \rho u_1 \\ Q_{3} = \rho u_2 , \quad Q_{4} = E , \quad Q_{5} = G_{u,1,1} \\ Q_{6} = G_{u,2,2} , \quad Q_{7} = G^{\text{aux}}_{0} , \quad Q_{8} = G^{\text{aux}}_{1} \\ Q_{9} = G^{\text{aux}}_{2} , \quad Q_{10} = G_{u,0,0} , \quad Q_{11} = S_{u,01} \\ Q_{12} = S_{u,02} , \quad Q_{13} = S_{u,12} \end{align}\]

10.12.0.2 Flux Tensor Components (\(\Phi_{i,j}\))

with the shared subexpressions \[\begin{align} \chi_{0} &= \rho u_0^{2} \\ \chi_{1} &= \frac{1}{\rho} \\ \chi_{2} &= \frac{3 \chi_{1}}{2} \\ \chi_{3} &= \alpha \mu \\ \chi_{4} &= \frac{4 \chi_{3}}{3} \\ \chi_{5} &= - E \\ \chi_{6} &= E \gamma \\ \chi_{7} &= \rho u_2^{2} \\ \chi_{8} &= \frac{\chi_{1}}{2} \\ \chi_{9} &= \chi_{7} \chi_{8} \\ \chi_{10} &= \chi_{0} \chi_{8} \\ \chi_{11} &= - \chi_{10} \gamma \\ \chi_{12} &= \rho u_1^{2} \\ \chi_{13} &= \chi_{12} \chi_{8} \\ \chi_{14} &= - \chi_{13} \gamma \\ \chi_{15} &= - \chi_{9} \gamma \\ \chi_{16} &= \frac{2 \chi_{3}}{3} \\ \chi_{17} &= G_{u,2,2} \chi_{16} + \chi_{11} + \chi_{14} + \chi_{15} + \chi_{5} + \chi_{6} + \chi_{9} \\ \chi_{18} &= G_{u,1,1} \chi_{16} + \chi_{13} \\ \chi_{19} &= 2 \chi_{3} \\ \chi_{20} &= - S_{u,01} \chi_{19} + \chi_{1} \rho u_0 \rho u_1 \\ \chi_{21} &= - S_{u,02} \chi_{19} + \chi_{1} \rho u_0 \rho u_2 \\ \chi_{22} &= G_{u,0,0} \chi_{16} + \chi_{10} \\ \chi_{23} &= - S_{u,12} \chi_{19} + \chi_{1} \rho u_1 \rho u_2 \\ \chi_{24} &= \alpha \kappa \\ \chi_{25} &= \frac{1}{\rho^{2}} \\ \chi_{26} &= \rho u_0^{3} \\ \chi_{27} &= \frac{\chi_{25}}{2} \\ \chi_{28} &= \chi_{1} \rho u_0 \\ \chi_{29} &= G_{u,0,0} \chi_{28} \\ \chi_{30} &= \chi_{1} \rho u_1 \\ \chi_{31} &= S_{u,01} \chi_{30} \\ \chi_{32} &= \chi_{1} \rho u_2 \\ \chi_{33} &= S_{u,02} \chi_{32} \\ \chi_{34} &= \chi_{27} \rho u_0 \\ \chi_{35} &= \rho u_1^{3} \\ \chi_{36} &= G_{u,1,1} \chi_{30} \\ \chi_{37} &= S_{u,01} \chi_{28} \\ \chi_{38} &= S_{u,12} \chi_{32} \\ \chi_{39} &= \chi_{27} \rho u_1 \\ \chi_{40} &= \rho u_2^{3} \\ \chi_{41} &= G_{u,2,2} \chi_{32} \\ \chi_{42} &= S_{u,02} \chi_{28} \\ \chi_{43} &= S_{u,12} \chi_{30} \\ \chi_{44} &= \chi_{27} \rho u_2 \\ \chi_{45} &= \frac{1}{\tau_{\nabla}} \\ \chi_{46} &= \chi_{1} \chi_{45} \\ \chi_{47} &= \chi_{27} \chi_{45} \\ \chi_{48} &= \chi_{0} \chi_{47} \\ \chi_{49} &= \chi_{12} \chi_{47} \\ \chi_{50} &= \chi_{47} \chi_{7} \\ \chi_{51} &= E \chi_{46} - \chi_{46} \chi_{6} + \chi_{48} \gamma - \chi_{48} + \chi_{49} \gamma - \chi_{49} + \chi_{50} \gamma - \chi_{50} \\ \chi_{52} &= \chi_{45} \chi_{8} \\ \chi_{53} &= - \chi_{52} \rho u_1 \\ \chi_{54} &= - \chi_{52} \rho u_0 \\ \chi_{55} &= - \chi_{52} \rho u_2 \end{align}\] the flux components evaluate to \[\begin{align} \Phi_{0,0} &= \rho u_0 \\ \Phi_{0,1} &= \rho u_1 \\ \Phi_{0,2} &= \rho u_2 \\ \Phi_{1,0} &= - G_{u,0,0} \chi_{4} + \chi_{0} \chi_{2} + \chi_{17} + \chi_{18} \\ \Phi_{1,1} &= \chi_{20} \\ \Phi_{1,2} &= \chi_{21} \\ \Phi_{2,0} &= \chi_{20} \\ \Phi_{2,1} &= - G_{u,1,1} \chi_{4} + \chi_{12} \chi_{2} + \chi_{17} + \chi_{22} \\ \Phi_{2,2} &= \chi_{23} \\ \Phi_{3,0} &= \chi_{21} \\ \Phi_{3,1} &= \chi_{23} \\ \Phi_{3,2} &= - G_{u,2,2} \chi_{4} + \chi_{11} + \chi_{14} + \chi_{15} + \chi_{18} + \chi_{2} \chi_{7} + \chi_{22} + \chi_{5} + \chi_{6} \\ \Phi_{4,0} &= E \chi_{1} \gamma \rho u_0 - G^{\text{aux}}_{0} \chi_{24} + \frac{2 G_{u,1,1} \alpha \chi_{1} \mu \rho u_0}{3} \\ &\quad + \frac{2 G_{u,2,2} \alpha \chi_{1} \mu \rho u_0}{3} + \frac{\chi_{12} \chi_{25} \rho u_0}{2} - \chi_{12} \chi_{34} \gamma \\ &\quad - \chi_{19} \chi_{31} - \chi_{19} \chi_{33} + \frac{\chi_{25} \chi_{26}}{2} + \frac{\chi_{25} \chi_{7} \rho u_0}{2} \\ &\quad - \chi_{26} \chi_{27} \gamma - \chi_{29} \chi_{4} - \chi_{34} \chi_{7} \gamma \\ \Phi_{4,1} &= E \chi_{1} \gamma \rho u_1 - G^{\text{aux}}_{1} \chi_{24} + \frac{2 G_{u,0,0} \alpha \chi_{1} \mu \rho u_1}{3} \\ &\quad + \frac{2 G_{u,2,2} \alpha \chi_{1} \mu \rho u_1}{3} + \frac{\chi_{0} \chi_{25} \rho u_1}{2} - \chi_{0} \chi_{39} \gamma \\ &\quad - \chi_{19} \chi_{37} - \chi_{19} \chi_{38} + \frac{\chi_{25} \chi_{35}}{2} + \frac{\chi_{25} \chi_{7} \rho u_1}{2} \\ &\quad - \chi_{27} \chi_{35} \gamma - \chi_{36} \chi_{4} - \chi_{39} \chi_{7} \gamma \\ \Phi_{4,2} &= E \chi_{1} \gamma \rho u_2 - G^{\text{aux}}_{2} \chi_{24} + \frac{2 G_{u,0,0} \alpha \chi_{1} \mu \rho u_2}{3} \\ &\quad + \frac{2 G_{u,1,1} \alpha \chi_{1} \mu \rho u_2}{3} + \frac{\chi_{0} \chi_{25} \rho u_2}{2} - \chi_{0} \chi_{44} \gamma \\ &\quad + \frac{\chi_{12} \chi_{25} \rho u_2}{2} - \chi_{12} \chi_{44} \gamma - \chi_{19} \chi_{42} - \chi_{19} \chi_{43} \\ &\quad + \frac{\chi_{25} \chi_{40}}{2} - \chi_{27} \chi_{40} \gamma - \chi_{4} \chi_{41} \\ \Phi_{5,0} &= G_{u,1,1} \chi_{28} \\ \Phi_{5,1} &= - \chi_{30} \chi_{45} + \chi_{36} \\ \Phi_{5,2} &= G_{u,1,1} \chi_{32} \\ \Phi_{6,0} &= G_{u,2,2} \chi_{28} \\ \Phi_{6,1} &= G_{u,2,2} \chi_{30} \\ \Phi_{6,2} &= - \chi_{32} \chi_{45} + \chi_{41} \\ \Phi_{7,0} &= G^{\text{aux}}_{0} \chi_{28} + \chi_{51} \\ \Phi_{7,1} &= G^{\text{aux}}_{0} \chi_{30} \\ \Phi_{7,2} &= G^{\text{aux}}_{0} \chi_{32} \\ \Phi_{8,0} &= G^{\text{aux}}_{1} \chi_{28} \\ \Phi_{8,1} &= G^{\text{aux}}_{1} \chi_{30} + \chi_{51} \\ \Phi_{8,2} &= G^{\text{aux}}_{1} \chi_{32} \\ \Phi_{9,0} &= G^{\text{aux}}_{2} \chi_{28} \\ \Phi_{9,1} &= G^{\text{aux}}_{2} \chi_{30} \\ \Phi_{9,2} &= G^{\text{aux}}_{2} \chi_{32} + \chi_{51} \\ \Phi_{10,0} &= - \chi_{28} \chi_{45} + \chi_{29} \\ \Phi_{10,1} &= G_{u,0,0} \chi_{30} \\ \Phi_{10,2} &= G_{u,0,0} \chi_{32} \\ \Phi_{11,0} &= \chi_{37} + \chi_{53} \\ \Phi_{11,1} &= \chi_{31} + \chi_{54} \\ \Phi_{11,2} &= S_{u,01} \chi_{32} \\ \Phi_{12,0} &= \chi_{42} + \chi_{55} \\ \Phi_{12,1} &= S_{u,02} \chi_{30} \\ \Phi_{12,2} &= \chi_{33} + \chi_{54} \\ \Phi_{13,0} &= S_{u,12} \chi_{28} \\ \Phi_{13,1} &= \chi_{43} + \chi_{55} \\ \Phi_{13,2} &= \chi_{38} + \chi_{53} \end{align}\]

10.12.0.3 Source Vector Components (\(S_i\))

\[\begin{align} S_{0} &= 0 \\ S_{1} &= - \frac{\alpha u_{s,0} \rho}{\eta} + \frac{\alpha \rho u_0}{\eta} + \frac{u_{s,0} \rho}{\eta} - \frac{\rho u_0}{\eta} \\ S_{2} &= - \frac{\alpha u_{s,1} \rho}{\eta} + \frac{\alpha \rho u_1}{\eta} + \frac{u_{s,1} \rho}{\eta} - \frac{\rho u_1}{\eta} \\ S_{3} &= - \frac{\alpha u_{s,2} \rho}{\eta} + \frac{\alpha \rho u_2}{\eta} + \frac{u_{s,2} \rho}{\eta} - \frac{\rho u_2}{\eta} \\ S_{4} &= - \frac{\alpha u_{s,0} \rho u_0}{\eta} - \frac{\alpha u_{s,1} \rho u_1}{\eta} - \frac{\alpha u_{s,2} \rho u_2}{\eta} \\ &\quad + \frac{\alpha \rho u_0^{2}}{\eta \rho} + \frac{\alpha \rho u_1^{2}}{\eta \rho} + \frac{\alpha \rho u_2^{2}}{\eta \rho} \\ &\quad + \frac{u_{s,0} \rho u_0}{\eta} + \frac{u_{s,1} \rho u_1}{\eta} + \frac{u_{s,2} \rho u_2}{\eta} - \frac{\rho u_0^{2}}{\eta \rho} \\ &\quad - \frac{\rho u_1^{2}}{\eta \rho} - \frac{\rho u_2^{2}}{\eta \rho} \\ S_{5} &= - \frac{G_{u,1,1}}{\tau_{\nabla}} \\ S_{6} &= - \frac{G_{u,2,2}}{\tau_{\nabla}} \\ S_{7} &= - \frac{G^{\text{aux}}_{0}}{\tau_{\nabla}} \\ S_{8} &= - \frac{G^{\text{aux}}_{1}}{\tau_{\nabla}} \\ S_{9} &= - \frac{G^{\text{aux}}_{2}}{\tau_{\nabla}} \\ S_{10} &= - \frac{G_{u,0,0}}{\tau_{\nabla}} \\ S_{11} &= - \frac{S_{u,01}}{\tau_{\nabla}} \\ S_{12} &= - \frac{S_{u,02}}{\tau_{\nabla}} \\ S_{13} &= - \frac{S_{u,12}}{\tau_{\nabla}} \end{align}\]

11 Method of Manufactured Solutions (MMS) Exact Fields↩︎

In this appendix, we document the analytical expressions of the Method of Manufactured Solutions (MMS) exact fields (\(\mathbf{Q}^{\text{exact}}\)) for all twelve verified physical PDE targets. The volume forcing terms are derived from these fields symbolically at compile time as the continuous residual \(\mathbf{S}_{\text{mms}} = \partial_t \mathbf{Q}^{\text{exact}} + \nabla \cdot \mathbf{\Phi}(\mathbf{Q}^{\text{exact}}) - \mathbf{S}(\mathbf{Q}^{\text{exact}})\). We list the conserved macroscopic components. The gradient-tracker components carried by some systems are initialized from the spatial derivatives of these fields and are omitted here. Cases declared with SI-valued scales author their fields as (scale) \(\times\) (dimensionless), and the value non-dimensionalization divides the scales out to reach the \(O(1)\) lattice working point tabulated per case.

11.1 Design Principles for the Manufactured Solutions↩︎

1. Smoothness. All fields are \(C^\infty\) trigonometric–exponential profiles, so the measured order reflects the truncation error of the scheme.

2. Large amplitudes. Fluctuations reach a finite fraction of the mean state, so the nonlinearities are exercised.

3. Independent modes. Fields coupled by a constitutive law carry independent fluctuation modes rather than being slaved to one another.

4. Broken axis symmetry. Per-axis wavenumbers, Mach numbers, or phase patterns are distinct.

5. Off-resonance forcing. Temporal frequencies avoid the free eigenfrequencies, and the evaluation instant falls at a generic phase.

6. Involutions by construction. Differential constraints carried by the physics are satisfied analytically by the exact fields.

11.2 Inviscid Burgers Equation MMS↩︎

11.2.0.1 Manufactured solution (physical form)

The exact manufactured profile is a smooth periodic wave of amplitude \(A = \max|u|\): \[\begin{align} u_{\text{exact}}(\mathbf{x}, t) = A\,\sin(2\pi x_0)\sin(2\pi x_1)\cos t. \end{align}\]

11.2.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= A \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(t )} \end{align}\]

11.3 Compressible Euler Equations MMS↩︎

11.3.0.1 Manufactured solution (physical form)

The manufactured solution superimposes smooth fluctuations on a supersonic advection background at per-axis Mach number \(3.0\). The reference state has density \(\rho_0\) and sound speed \(c_s = \sqrt{\gamma p_0/\rho_0}\), with reference pressure \(p_0 = \rho_0 c_s^2/\gamma\) and background velocity \(u_0 = v_0 = 3\,c_s\). With wavenumber \(k = 2\pi\) the fluctuations are \[\begin{align} \begin{aligned} \Delta\rho &= 0.3\,\rho_0\,\sin(k x_0)\cos(k x_1)\cos t, \\ \Delta u &= 0.3\,c_s\,\cos(k x_0)\sin(k x_1)\sin t, \\ \Delta v &= 0.3\,c_s\,\sin(k x_0)\sin(k x_1)\sin t, \\ \Delta p &= c_s^2\,\Delta\rho + 0.2\,p_0\,\cos(k x_0)\sin(k x_1)\sin t, \end{aligned} \end{align}\] so the pressure carries an independent mode in addition to the linear response \(c_s^2\Delta\rho\). The conserved state is assembled as \(\rho_{\text{exact}} = \rho_0 + \Delta\rho\), \((\rho\mathbf{u})_{\text{exact}} = \rho_{\text{exact}}(\mathbf{u}_0 + \Delta\mathbf{u})\), \(p_{\text{exact}} = p_0 + \Delta p\), and \(E_{\text{exact}} = p_{\text{exact}}/(\gamma-1) + \tfrac12\rho_{\text{exact}}|\mathbf{u}_{\text{exact}}|^2\).

11.3.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} + \rho_{0} \\ Q_{1}^{\text{exact}} &= (\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} \\ &\quad + 3 c_{s}) (\frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} + \rho_{0}) \\ Q_{2}^{\text{exact}} &= (\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{10} \\ &\quad + 3 c_{s}) (\frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} + \rho_{0}) \\ Q_{3}^{\text{exact}} &= \frac{5 c_{s}^{2} \rho_{0} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{14} \\ &\quad + \frac{3 c_{s}^{2} \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{4} + \frac{25 c_{s}^{2} \rho_{0}}{14} \\ &\quad + (\frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{20} \\ &\quad + \frac{\rho_{0}}{2}) ((\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{10} + 3 c_{s})^{2} \\ &\quad + (\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} + 3 c_{s})^{2}) \end{align}\]

11.4 Shallow Water Equations MMS↩︎

11.4.0.1 Manufactured solution (physical form)

The manufactured fields evolve the water depth \(h\) and discharge \(h\mathbf{u}\) over a static bathymetry \(z_b(\mathbf{x}) = 0.01\sin(2\pi x_0)\cos(2\pi x_1)\). Depth fluctuations are \(40\%\) of the mean depth, giving a subcritical Froude number \(\approx 0.3\) with \(h \in [0.06,\,0.14]\) so \(h > 0\) everywhere: \[\begin{align} \begin{aligned} h_{\text{exact}} &= 0.04\,\sin(2\pi x_0)\sin(2\pi x_1)\cos t + 0.1, \\ (hu_0)_{\text{exact}} &= 0.02\,\cos(2\pi x_0)\cos(2\pi x_1)\sin t + 0.03, \\ (hu_1)_{\text{exact}} &= 0.02\,\sin(2\pi x_0)\sin(2\pi x_1)\sin t + 0.03. \end{aligned} \end{align}\] The characteristic gravity-wave speed \(c_s = \sqrt{g H}\) at reference depth \(H = 0.1\) sets the convective time scale and reduces the realistic-SI gravity to the working point. The values are listed below.

11.4.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{\sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(t )}}{25} + \frac{1}{10} \\ Q_{1}^{\text{exact}} &= \frac{\sin{(t )} \cos{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{50} + \frac{3}{100} \\ Q_{2}^{\text{exact}} &= \frac{\sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{50} + \frac{3}{100} \end{align}\]

11.5 Ideal Ultrarelativistic Fluid MMS↩︎

11.5.0.1 Manufactured solution (physical form)

The manufactured fields for the ideal ultrarelativistic fluid are smooth fluctuations about a uniform background, kept inside the physicality region \(E > |\mathbf{S}|\) (three-velocities up to \(|\mathbf{v}| \approx 0.79\), Lorentz factor \(\Gamma \approx 1.6\) where \(E\) is low and \(|\mathbf{S}|\) high): \[\begin{align} \begin{aligned} E_{\text{exact}}(\mathbf{x}, t) &= 0.4\,\sin(2\pi x_0)\sin(2\pi x_1)\cos t + 3.0, \\ S_{0,\text{exact}}(\mathbf{x}, t) &= 0.4\,\cos(2\pi x_0)\cos(2\pi x_1)\sin t + 1.2, \\ S_{1,\text{exact}}(\mathbf{x}, t) &= 0.4\,\sin(2\pi x_0)\sin(2\pi x_1)\sin t + 1.2, \end{aligned} \end{align}\] with conserved energy density \(E\) and momentum \(\mathbf{S}\) in natural units (\(c = 1\)).

11.5.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{2 \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(t )}}{5} + 3 \\ Q_{1}^{\text{exact}} &= \frac{2 \sin{(t )} \cos{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{5} + \frac{6}{5} \\ Q_{2}^{\text{exact}} &= \frac{2 \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{5} + \frac{6}{5} \end{align}\]

11.6 Maxwell Electromagnetics (TM Mode) MMS↩︎

11.6.0.1 Manufactured solution (physical form)

The manufactured solution is a TM pattern with distinct wavenumbers \(k_x = 2\pi\) and \(k_y = 4\pi\), amplitude \(A = \tfrac{1}{2}\), and temporal frequency \(\omega = 4\), with the electric and magnetic fields in the physical \(90^\circ\) quadrature (\(E_z \sim \cos\omega t\), \(\mathbf{H} \sim \sin\omega t\)): \[\begin{align} \begin{aligned} E_{z,\text{exact}}(\mathbf{x}, t) &= A \sin(k_x x_0)\sin(k_y x_1)\cos(\omega t), \\ H_{x,\text{exact}}(\mathbf{x}, t) &= A \sin(k_x x_0)\cos(k_y x_1)\sin(\omega t), \\ H_{y,\text{exact}}(\mathbf{x}, t) &= -A \frac{k_x}{k_y} \cos(k_x x_0)\sin(k_y x_1)\sin(\omega t). \end{aligned} \end{align}\] The field satisfies \(\nabla \cdot \mathbf{H} = 0\) analytically.

11.6.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{\sin{(2 \pi x_0 )} \sin{(4 \pi x_1 )} \cos{(4 t )}}{2} \\ Q_{1}^{\text{exact}} &= \frac{\sin{(4 t )} \sin{(2 \pi x_0 )} \cos{(4 \pi x_1 )}}{2} \\ Q_{2}^{\text{exact}} &= - \frac{\sin{(4 t )} \sin{(4 \pi x_1 )} \cos{(2 \pi x_0 )}}{4} \end{align}\]

11.7 Nonlinear Finite-Strain Elasticity MMS↩︎

11.7.0.1 Manufactured solution (physical form)

The case is the deformation-gradient (total-Lagrangian) formulation of compressible neo-Hookean elastodynamics of the main article, with conserved state \((\mathbf{F}, \mathbf{v})\). The manufactured fields derive from a smooth, non-separable displacement potential \(\mathbf{w}\) of amplitude \(A = 0.06\) via \(\mathbf{F} = \mathbf{I} + \nabla\mathbf{w}\) (curl-free, so the compatibility involution holds exactly) and \(\mathbf{v} = \partial_t\mathbf{w}\), reaching a \(37\%\) finite strain with \(\det\mathbf{F}\ge 0.68 > 0\) everywhere: \[\begin{align} \begin{aligned} w_0 &= A\big(\sin(2\pi x_0)\sin(2\pi x_1)\cos t + \tfrac{1}{2}\sin(2\pi x_0)\cos(2\pi x_1)\sin t \\ &\qquad + \tfrac{1}{4}\cos(2\pi(x_0+x_1))\cos t\big), \\ w_1 &= A\big(\cos(2\pi x_0)\cos(2\pi x_1)\sin t + \tfrac{1}{2}\cos(2\pi x_0)\sin(2\pi x_1)\cos t \\ &\qquad + \tfrac{1}{4}\sin(2\pi(x_0-x_1))\sin t\big). \end{aligned} \end{align}\] The material values (aluminium-like) are listed below. The rest state \(\mathbf{F} = \mathbf{I}\), \(\mathbf{v} = 0\) gives \(\mathbf{P}(\mathbf{I}) = 0\), so the reference state is the rest configuration.

11.7.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{3 \pi \sin{(t )} \cos{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{50} \\ &\quad + \frac{3 \pi \sin{(2 \pi x_1 )} \cos{(t )} \cos{(2 \pi x_0 )}}{25} - \frac{3 \pi \sin{(\pi (2 x_0 + 2 x_1) )} \cos{(t )}}{100} \\ &\quad + 1 \\ Q_{1}^{\text{exact}} &= - \frac{3 \pi \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{50} \\ &\quad + \frac{3 \pi \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{25} - \frac{3 \pi \sin{(\pi (2 x_0 + 2 x_1) )} \cos{(t )}}{100} \\ Q_{2}^{\text{exact}} &= - \frac{3 \pi \sin{(t )} \sin{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{25} + \frac{3 \pi \sin{(t )} \cos{(\pi (2 x_0 - 2 x_1) )}}{100} \\ &\quad - \frac{3 \pi \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(t )}}{50} \\ Q_{3}^{\text{exact}} &= - \frac{3 \pi \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{25} - \frac{3 \pi \sin{(t )} \cos{(\pi (2 x_0 - 2 x_1) )}}{100} \\ &\quad + \frac{3 \pi \cos{(t )} \cos{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{50} + 1 \\ Q_{4}^{\text{exact}} &= - \frac{3 \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{50} - \frac{3 \sin{(t )} \cos{(\pi (2 x_0 + 2 x_1) )}}{200} \\ &\quad + \frac{3 \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{100} \\ Q_{5}^{\text{exact}} &= - \frac{3 \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{100} + \frac{3 \sin{(\pi (2 x_0 - 2 x_1) )} \cos{(t )}}{200} \\ &\quad + \frac{3 \cos{(t )} \cos{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{50} \end{align}\]

11.8 Scalar Advection-Diffusion-Reaction MMS↩︎

11.8.0.1 Manufactured solution (physical form)

The exact manufactured verification profile is defined as: \[\begin{align} c_{\text{exact}}(\mathbf{x}, t) = 0.5 \sin(2\pi(x_0 - u_{0}t))\cos(2\pi(x_1 - u_{1}t))\exp(-t) + 0.5. \end{align}\] System parameters are listed below.

11.8.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{1}{2} + \frac{e^{- t} \sin{(\pi (- 2 t u_0 + 2 x_0) )} \cos{(\pi (- 2 t u_1 + 2 x_1) )}}{2} \end{align}\]

11.9 Allen–Cahn Phase-Field Equation MMS↩︎

11.9.0.1 Manufactured solution (physical form)

The exact manufactured verification profile is a smooth periodic order parameter whose amplitude is kept below unity: \[\begin{align} \phi_{\text{exact}}(\mathbf{x}, t) = 0.7 \sin(2\pi x_0)\cos(2\pi x_1)\cos(t). \end{align}\] System parameters are listed below.

11.9.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{7 \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} \end{align}\]

11.10 Weakly Compressible Navier–Stokes Equations (Incompressible Limit) MMS↩︎

11.10.0.1 Manufactured solution (physical form)

The manufactured solution is the two-dimensional Taylor–Green vortex, an exact solution of the incompressible Navier–Stokes equations, used here to verify the weakly compressible (artificial-compressibility) scheme in its incompressible limit. The working point is weakly compressible at Mach number \(\text{Ma} = U/c_s = 0.1\) and Reynolds number \(\text{Re} = U L/\nu = 10\), with reference density \(\rho_0\), artificial sound speed \(c_s\), velocity amplitude \(U\), and kinematic viscosity \(\nu\) (\(\mu = \rho_0\nu\)). With wavenumber \(k = 2\pi\) and the exact viscous decay \(\mathcal{D}(t) = \exp(-2(\nu/c_s)k^2 t)\) in characteristic time, the velocity and pressure are \[\begin{align} \begin{aligned} u &= -U\,\cos(k x_0)\sin(k x_1)\,\mathcal{D}(t), \\ v &= U\,\sin(k x_0)\cos(k x_1)\,\mathcal{D}(t), \\ p &= -\tfrac14\rho_0 U^2\,[\cos(2 k x_0) + \cos(2 k x_1)]\,\mathcal{D}(t)^2. \end{aligned} \end{align}\] The density is slaved to the dynamic pressure, \(\rho_{\text{exact}} = \rho_0 + p/c_s^2\), so the density fluctuation is \(O(\text{Ma}^2)\), and the conserved momentum is \((\rho\mathbf{u})_{\text{exact}} = \rho_{\text{exact}}\mathbf{u}\).

11.10.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= - \frac{\rho_{0} (\cos{(4 \pi x_0 )} + \cos{(4 \pi x_1 )}) e^{- \frac{4 \pi^{2} t}{25}}}{400} + \rho_{0} \\ Q_{1}^{\text{exact}} &= - \frac{c_{s} \rho_{0} e^{- \frac{2 \pi^{2} t}{25}} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} \\ &\quad + \frac{c_{s} \rho_{0} e^{- \frac{6 \pi^{2} t}{25}} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )} \cos{(4 \pi x_0 )}}{4000} \\ &\quad + \frac{c_{s} \rho_{0} e^{- \frac{6 \pi^{2} t}{25}} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )} \cos{(4 \pi x_1 )}}{4000} \\ Q_{2}^{\text{exact}} &= \frac{c_{s} \rho_{0} e^{- \frac{2 \pi^{2} t}{25}} \sin{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{10} \\ &\quad - \frac{c_{s} \rho_{0} e^{- \frac{6 \pi^{2} t}{25}} \sin{(2 \pi x_0 )} \cos{(4 \pi x_0 )} \cos{(2 \pi x_1 )}}{4000} \\ &\quad - \frac{c_{s} \rho_{0} e^{- \frac{6 \pi^{2} t}{25}} \sin{(2 \pi x_0 )} \cos{(2 \pi x_1 )} \cos{(4 \pi x_1 )}}{4000} \end{align}\]

11.11 Compressible Navier–Stokes–Fourier Equations MMS↩︎

11.11.0.1 Manufactured solution (physical form)

The manufactured solution models a supersonic viscous flow with an anisotropic mean velocity at per-axis Mach numbers \(3.0\) and \(2.0\), at Prandtl number \(\text{Pr} = 0.71\) (air). The reference state has density \(\rho_0\), gas constant \(R\), and sound speed \(c_s = \sqrt{\gamma R T_0}\) at temperature \(\require{physics} T_0 = \qty{288}{\kelvin}\), with reference pressure \(p_0 = \rho_0 R T_0\) and background velocity \(u_0 = 3\,c_s\), \(v_0 = 2\,c_s\). With wavenumber \(k = 2\pi\) the fluctuations use a distinct amplitude and spatial structure per field, \[\begin{align} \begin{aligned} \Delta\rho &= 0.20\,\rho_0\,\sin(k x_0)\cos(k x_1)\cos t, \\ \Delta u &= 0.20\,c_s\,\cos(k x_0)\sin(k x_1)\sin t, \\ \Delta v &= 0.15\,c_s\,\sin(k x_0)\sin(k x_1)\cos t, \\ \Delta p &= c_s^2\,\Delta\rho + 0.20\,p_0\,\cos(k x_0)\sin(k x_1)\sin t, \end{aligned} \end{align}\] so the pressure carries an independent mode beyond the linear response. The conserved state is \(\rho_{\text{exact}} = \rho_0 + \Delta\rho\), \((\rho\mathbf{u})_{\text{exact}} = \rho_{\text{exact}}(\mathbf{u}_0 + \Delta\mathbf{u})\), and \(E_{\text{exact}} = p_{\text{exact}}/(\gamma-1) + \tfrac12\rho_{\text{exact}}|\mathbf{u}_{\text{exact}}|^2\) with \(p_{\text{exact}} = p_0 + \Delta p\). The thermal conductivity follows from the Prandtl number, \(K = \mu\,c_p/\text{Pr}\) with \(c_p = \gamma R/(\gamma-1)\), placing the flow in a resolved viscous regime (\(\text{Re} = \rho_0 u_0 L/\mu \approx 3.5\times10^2\) on the unit domain \(L = 1\)).

11.11.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{\rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + \rho_{0} \\ Q_{1}^{\text{exact}} &= (\frac{c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{5} \\ &\quad + 3 c_{s}) (\frac{\rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + \rho_{0}) \\ Q_{2}^{\text{exact}} &= (\frac{3 c_{s} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(t )}}{20} \\ &\quad + 2 c_{s}) (\frac{\rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + \rho_{0}) \\ Q_{3}^{\text{exact}} &= \frac{5 c_{s}^{2} \rho_{0} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{14} \\ &\quad + \frac{c_{s}^{2} \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{2} + \frac{25 c_{s}^{2} \rho_{0}}{14} \\ &\quad + (\frac{\rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} \\ &\quad + \frac{\rho_{0}}{2}) ((\frac{c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{5} + 3 c_{s})^{2} \\ &\quad + (\frac{3 c_{s} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(t )}}{20} + 2 c_{s})^{2}) \end{align}\]

11.12 Resistive Compressible Magnetohydrodynamics MMS↩︎

11.12.0.1 Manufactured solution (physical form)

The manufactured solution is a supersonic magnetized flow at per-axis Mach number \(3.0\) with smooth fluctuations about a uniform background. The magnetic field is carried in Alfven (velocity) units. The reference state has density \(\rho_0\) and sound speed \(c_s = \sqrt{\gamma p_0/\rho_0}\) (reference pressure \(p_0 = \rho_0 c_s^2/\gamma\)), with background velocity \(u_0 = v_0 = 3\,c_s\) and background field \(B_{x,0} = B_{y,0} = 0.5\,c_s\) (plasma \(\beta = 4/\gamma \approx 2.9\)). With wavenumber \(k = 2\pi\) the fluctuations are \[\begin{align} \begin{aligned} \Delta\rho &= 0.3\,\rho_0\,\sin(k x_0)\cos(k x_1)\cos t, \\ \Delta p &= c_s^2\,\Delta\rho + 0.2\,p_0\,\cos(k x_0)\sin(k x_1)\sin t, \\ \Delta u &= 0.3\,c_s\,\cos(k x_0)\sin(k x_1)\sin t, \\ \Delta v &= 0.3\,c_s\,\sin(k x_0)\sin(k x_1)\sin t, \\ \Delta B_x &= -0.3\,c_s\,\cos(k x_0)\sin(k x_1)\sin t, \\ \Delta B_y &= 0.3\,c_s\,\sin(k x_0)\cos(k x_1)\sin t, \end{aligned} \end{align}\] with the magnetic perturbation divergence-free by construction (\(\partial_{x_0}\Delta B_x + \partial_{x_1}\Delta B_y = 0\)), and the pressure carrying an independent mode beyond the linear response \(c_s^2\Delta\rho\). The conserved energy is \(E_{\text{exact}} = p_{\text{exact}}/(\gamma-1) + \tfrac12\rho_{\text{exact}}|\mathbf{u}_{\text{exact}}|^2 + \tfrac12|\mathbf{B}_{\text{exact}}|^2\), and a Powell source proportional to \(\nabla\cdot\mathbf{B}\) controls the solenoidal error. The magnetic diffusivity, dynamic viscosity, and heat-flux coefficient are set to equal diffusivities.

11.12.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} + \rho_{0} \\ Q_{1}^{\text{exact}} &= (\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} \\ &\quad + 3 c_{s}) (\frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} + \rho_{0}) \\ Q_{2}^{\text{exact}} &= (\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{10} \\ &\quad + 3 c_{s}) (\frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{10} + \rho_{0}) \\ Q_{3}^{\text{exact}} &= \frac{5 c_{s}^{2} \rho_{0} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{14} \\ &\quad + \frac{3 c_{s}^{2} \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{4} + \frac{25 c_{s}^{2} \rho_{0}}{14} \\ &\quad + \frac{(\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{10} + \frac{c_{s}}{2})^{2}}{2} \\ &\quad + \frac{(- \frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} + \frac{c_{s}}{2})^{2}}{2} \\ &\quad + (\frac{3 \rho_{0} \sin{(2 \pi x_0 )} \cos{(t )} \cos{(2 \pi x_1 )}}{20} \\ &\quad + \frac{\rho_{0}}{2}) ((\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )}}{10} + 3 c_{s})^{2} \\ &\quad + (\frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} + 3 c_{s})^{2}) \\ Q_{4}^{\text{exact}} &= - \frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_0 )}}{10} + \frac{c_{s}}{2} \\ Q_{5}^{\text{exact}} &= \frac{3 c_{s} \sin{(t )} \sin{(2 \pi x_0 )} \cos{(2 \pi x_1 )}}{10} + \frac{c_{s}}{2} \end{align}\]

11.13 3D Homogenized Compressible Navier–Stokes–Fourier Equations MMS↩︎

11.13.0.1 Manufactured solution (physical form)

The manufactured solution is fully three-dimensional with active compressibility (\(\nabla \rho \neq 0\)), evolved as fluctuations of amplitude \(A = 0.2\) about a subsonic mean state at per-axis Mach number \(\text{Ma} = 0.5\) (air at standard conditions; reference density \(\rho_0\), isothermal sound-speed scale \(c_s = \sqrt{p_0/\rho_0}\), reference pressure \(p_0 = \rho_0 c_s^2\), mean (adiabatic) sound speed \(\bar{c} = \sqrt{\gamma}\,c_s\), drift \(U_0 = \text{Ma}\,\bar{c}\)). The fluid primitives (\(k = 2\pi\)) are \[\begin{align} \begin{aligned} \rho &= \rho_0 + A\,\sin(k x_0)\cos(k x_1)\sin(k x_2)\cos t, \\ u_0 &= U_0 + A\,\cos(k x_0)\sin(k x_1)\sin(k x_2)\sin t, \\ u_1 &= U_0 + A\,\sin(k x_0)\cos(k x_1)\sin(k x_2)\sin t, \\ u_2 &= U_0 + A\,\sin(k x_0)\sin(k x_1)\cos(k x_2)\sin t, \\ p &= p_0 + \bar{c}^2(\rho - \rho_0) + A\,p_0\,\cos(k x_0)\sin(k x_1)\cos(k x_2)\sin t, \end{aligned} \end{align}\] so the pressure combines the linear response \(\bar{c}^2(\rho-\rho_0)\) with an independent mode, and the conserved state is assembled algebraically as \(\mathbf{Q} = (\rho,\;\rho\mathbf{u},\;E)\) with \(E = p/(\gamma-1) + \tfrac{1}{2}\rho|\mathbf{u}|^2\). The prescribed solid-velocity input field superimposes a rotation (amplitude \(0.2\)), a uniform translation \(\mathbf{c} = (0.1,\,-0.1,\,0.05)\), and an independent oscillation (amplitude \(0.1\)), \[\begin{align} \begin{aligned} u_{s,0} &= 0.2\,\sin(k x_0)\cos(k x_1)\cos(k x_2)\cos t + c_0 + 0.1\,\sin t\,\cos(k x_2), \\ u_{s,1} &= -0.2\,\cos(k x_0)\sin(k x_1)\cos(k x_2)\cos t + c_1 + 0.1\,\cos t\,\cos(k x_0), \\ u_{s,2} &= 0.2\,\cos(k x_0)\cos(k x_1)\sin(k x_2)\sin t + c_2 + 0.1\,\sin t\,\cos(k x_1), \end{aligned} \end{align}\] and the porosity is a smooth, fully three-dimensional and non-separable sinusoid advected along the solid displacement \(\boldsymbol{\xi} = \mathbf{x} - \mathbf{d}(\mathbf{x},t)\), \[\begin{align} \alpha &= \frac{1}{2} + \frac{1}{10}\Big[\sin(k\xi_0) + \sin(k\xi_1) + \sin(k\xi_2) \nonumber \\ &\qquad\quad + 2\sin(k\xi_0)\sin(k\xi_1)\sin(k\xi_2)\Big] \in [0,\,1], \end{align}\] with \(d_i = (\text{rotation})\,\sin t + c_i\,t\) a smooth solid displacement, so \(\alpha\) spans the full solid-to-fluid range and is transported.

11.13.0.2 Compiler-expanded fields (\(Q_i^{\text{exact}}\))

\[\begin{align} Q_{0}^{\text{exact}} &= \rho_{0} (\frac{\sin{(2 \pi x_0 )} \sin{(2 \pi x_2 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + 1) \\ Q_{1}^{\text{exact}} &= c_{s} \rho_{0} (\frac{\sin{(t )} \sin{(2 \pi x_1 )} \sin{(2 \pi x_2 )} \cos{(2 \pi x_0 )}}{5} \\ &\quad + 0.591607978309962) (\frac{\sin{(2 \pi x_0 )} \sin{(2 \pi x_2 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + 1) \\ Q_{2}^{\text{exact}} &= c_{s} \rho_{0} (\frac{\sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_2 )} \cos{(2 \pi x_1 )}}{5} \\ &\quad + 0.591607978309962) (\frac{\sin{(2 \pi x_0 )} \sin{(2 \pi x_2 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + 1) \\ Q_{3}^{\text{exact}} &= c_{s} \rho_{0} (\frac{\sin{(t )} \sin{(2 \pi x_0 )} \sin{(2 \pi x_1 )} \cos{(2 \pi x_2 )}}{5} \\ &\quad + 0.591607978309962) (\frac{\sin{(2 \pi x_0 )} \sin{(2 \pi x_2 )} \cos{(t )} \cos{(2 \pi x_1 )}}{5} + 1) \end{align}\]

References↩︎