Gradient Reconstruction in Lattice Boltzmann Methods for Systems of Conservation Laws

1. Introduction.

• The first non-equilibrium moment carries the gradient to leading order, but feedback through the reconstructed flux makes it an operator-transformed gradient.Treating that moment directly as the gradient leaves an error that does not vanish under refinement.Source: §1 · pp. 1–2
• The compiler pipeline replaces gradient tracking with in-collision recovery, so relaxation and lattice stages act on conserved state Q alone.The alternative tracked state is Q+ = (Q, ∇Q).Source: Figure 1 · p. 2
• The reconstructed closure gradients are recovered from conserved fields inside the collision through a single declaration, trading field transport for local computation.The approach was noted in but not pursued there.Source: §1 · pp. 1–2
• The reconstruction is the algebraic inverse of that linear operator, generated in closed form from each declared flux using the compiler’s symbolic coefficients.The compiler already assembles the operator coefficient symbolically with SymPy.Source: §1 · p. 2

2. Methodology.

• The scheme assigns populations to each conserved field, embeds physical flux in the linear equilibrium’s first moment, and recovers the field from its zeroth moment.It advances through a BGK collide-and-stream update with the source projected onto lattice weights.Source: §2.1 · pp. 2–3

2.1. Gradient Reconstruction.

• Subtracting the advective flux yields a direct gradient estimate from the non-equilibrium moment, while the flux time derivative is neglected as a higher-order term under diffusive and acoustic scalings.The paper reports second-order convergence across the tested systems despite this approximation.Source: §2.1 · p. 3

2.2. The Reconstruction Operator.

• The moment contains the diffusive-flux Jacobian’s image of the gradient, so reconstruction applies the algebraic inverse of that linear operator.Near the stability limit, the scaling factor is approximately 0.26 for advection–diffusion–reaction and 0.21 for Allen–Cahn; direct moment readings overshoot and do not converge under refinement.Source: §2.2 · p. 4

2.3. Compilation.

• The compiler generates reconstruction for declared fluxes, forms the required Jacobians symbolically, and rejects requested gradients that cannot be reduced to functions of the conserved state.Constant operators are compile-time constants; state-dependent operators are solved per cell using small block solves and scalar divisions.Source: §2.2 · p. 5, §2.3 · pp. 5–6
• Verification reuses the five systems’ exact solutions, lattices, and refinement laws, replacing only gradient tracking with reconstruction and measuring relative L2 and L∞ errors and empirical convergence order.Refinement follows each system’s declared diffusive or acoustic scaling.Source: §3 · p. 6

3. Verification by Manufactured Solutions.

• Purely hyperbolic systems are outside the operator’s scope because they have no diffusive flux or gradient closure to reconstruct.The operator targets systems with gradient-form constitutive closures.Source: §3 · p. 7
• Every conserved field in every system converges at second order in double precision.Table 1 reports relative errors and observed convergence order at the finest resolution.Source: §3 · p. 7, Table 1 · p. 7
• Reconstruction has about 1.4 times lower relative L2 error on average than gradient tracking, except for HCNSF density.The comparison is at matched scaling and resolution in double precision.Source: §3 · p. 7, Table 2 · p. 8
• At that accuracy, reconstruction carries only conserved fields: 3 versus 7 for Navier–Stokes and 6 versus 18 for magnetohydrodynamics.The comparison is with the gradient-tracking scheme.Source: §3 · p. 7
• In single precision, resistive magnetohydrodynamics and Navier–Stokes retain second-order convergence near double-precision error, while HCNSF has about twice that error.Advection–diffusion–reaction and Allen–Cahn fall below second order on the finest grids; double precision retains second order throughout.Source: §3 · p. 7

4. Performance.

• The reconstruction is faster for every target equation in both precisions, with gains scaling with the reduction in transported fields.The comparison is with the gradient-tracking scheme.Source: §4 · p. 8
• Byte traffic falls by up to a factor of three, setting the speedup in the memory-bound single-precision regime.Fewer populations also reduce register use, increase occupancy, and avoid local-memory spills.Source: §4 · p. 8
• A single declaration generates gradient reconstruction generically from the declared flux, inverting the diffusive-flux image of the non-equilibrium moment without auxiliary state.The reconstructed gradients are recovered inside the collision.Source: §5 · p. 8

• Figure 2 places the single-precision reconstruction kernels on or near the NVIDIA RTX A5000's 680 GB/s memory roof.Source: Figure 2 · p. 9
• Up to 3.7 times faster in single precision on an NVIDIA RTX A5000, with memory-bound kernels reaching up to 97% of peak bandwidth.Source: §5 · p. 10
• Second-order convergence in double precision spans five systems and matches gradient tracking at equal resolution while carrying only the conserved fields.The systems are advection-diffusion-reaction, Allen–Cahn, Navier–Stokes, resistive magnetohydrodynamics, and three-dimensional homogenized compressible Navier–Stokes–Fourier.Source: §5 · pp. 8, 10