Gravity Source Term Discretization

• Gravity source term discretization is a set of numerical techniques that integrate gravitational forces as source terms into PDE models, preserving equilibrium states.
• Techniques recast source terms as discrete flux divergences using methods like finite volume, DG, and lattice Boltzmann to achieve well-balanced and conservative schemes.
• Advanced methods ensure entropy stability, positivity preservation, and precise momentum and energy conservation across complex meshes and varying equations of state.

Gravity source term discretization encompasses a set of numerical techniques designed to accurately and robustly represent gravitational forces as source terms in PDE-based models of fluid and magnetohydrodynamic flows. The goal is to construct discrete approximations that precisely balance gravitational sources with fluxes in the underlying conservation laws, thereby ensuring exact preservation of hydrostatic or more general equilibrium states, and maintaining critical invariants such as total energy, momentum, and entropy. Specialized discretizations have been developed for finite volume, finite element, discontinuous Galerkin, and lattice Boltzmann frameworks, tailored for various equations of state, mesh geometries, and application domains. Rigorous results establish criteria for well-balancedness, conservation, and entropy-stability, with verified numerical tests on highly stratified, turbulent, or astrophysical problems.

1. Conservative Reformulations and Flux-Splitting

Modern schemes often recast the gravitational source terms as discrete divergences of associated fluxes or stress tensors. For example, Mullen et al. (Mullen et al., 2020) rewrite the momentum and energy Euler equations with gravity via:

• Momentum equation: $∂(ρv)/∂t + ∇·(ρvv + PI) = ρg$
• Energy equation: $∂E/∂t + ∇·[(E+P)v] = ρv·g$, where $g = -∇φ$

They define additional fluxes:

• Gravitational stress tensor $T_g = (1/4πG)[∇φ∇φ - ½|∇φ|^2I]$
• Gravitational energy flux $F_g = ρφv$

Such reformulations allow the source terms $ρg$ and $ρv·g$ to be recast as $-∇·T_g$ and $-∇·F_g′+∂_tE_g′$, and thus incorporated directly into the fluxes in finite volume schemes, enabling strict discrete conservation of momentum and total energy.

IMEX relaxation schemes for compressible flows with gravity achieve the same effect by pressure splitting and relaxation, so that the gravitational acceleration enters the explicit source part and is folded into the flux, facilitating exact well-balanced dynamics (see (Thomann et al., 2019)).

2. Well-Balanced Schemes for Hydrostatic and General Equilibria

A defining feature of gravitational source discretization is the construction of well-balanced schemes that exactly preserve discrete hydrostatic (or more general steady) solutions. Various papers (Varma et al., 2018, Berberich et al., 2018, Chandrashekar et al., 2015, Chertock et al., 2017, Kanbar et al., 2022, Liu et al., 29 Jul 2025) present algebraic construction of source terms and reconstructions guaranteeing this property.

• In second-order FV schemes (Varma et al., 2018), the gravitational source in momentum is approximated by the difference of hydrostatically reconstructed pressures at cell faces: $s_i = [p_i^L-p_i^R]/Δx$, with face values constructed so that $w=[ρe^{-ψ},u,pe^{-ψ}]$ is constant for equilibrium states.
• Discontinuous Galerkin methods (Chandrashekar et al., 2015, Liu et al., 29 Jul 2025) project the source onto GLL nodes, with special rewritings (isothermal, polytropic) ensuring the discrete derivative exactly matches the flux derivative at interpolation points.

Table: Key well-balanced discretization components

Method	Source Approximation	Reconstruction/Quadrature
FV (Godunov)	Face pressure difference	Hydrostatic variable
DG (nodal)	GLL projection of rewritten source	Lagrange/Legendre basis
Central schemes	Subtraction of steady state	Staggered grid averages

These constructions yield preservation of isothermal/polytropic equilibria and general equilibrium states to machine precision, even on curvilinear or unstructured meshes.

3. Entropy Stability and Positivity Preservation

Advanced gravitational source term discretizations also integrate mechanisms for fully discrete entropy-stability and positivity preservation. The entropy-stable nodal DG methods (Liu et al., 29 Jul 2025) and FV schemes (Berthon et al., 2024) implement source-flux coupling and entropy correction terms to rigorously enforce discrete entropy inequalities (à la Tadmor).

• Discrete corrections: e.g., in (Liu et al., 29 Jul 2025), any entropy violating component of the source is locally cancelled via $$(S^{\rm corr})_\ell^i = \sigma_i(V_\ell^i - \bar{V}^i)$$, with the projection orthogonal to conserved quantities.
• Positivity limiters: FV and DG schemes enforce cell- and node-wise positivity of density and pressure via scaling limiters post-update, with CFL constraints computable from scheme coefficients.

A proof-of-concept is presented in (Berthon et al., 2024), where interface source terms are designed to ensure that in discrete steady states, the numerical flux and the gravitational source cancel identically. The updated solution remains positive definite and satisfies a discrete entropy law for any convex entropy.

4. Extensions to Arbitrary Equations of State and Meshes

Schemes such as those in (Berberich et al., 2018, Chandrashekar et al., 2015) operate without requiring analytical forms of hydrostatic equilibria or restriction to ideal gas. The gravitational source can be constructed as a divergence of metric-weighted hydrostatic profiles for arbitrary EoS and over general (even curvilinear or unstructured) meshes. Hydrostatic variable reconstructions are EoS-agnostic and mesh-agnostic, with the only requirement being the provision of tabulated or analytic profiles.

Finite element approaches (Lundgren et al., 2024) propose a sharp energy-conserving modification for variable density flows under gravity, compatible with fully unstructured meshes. The additional term $-\frac12 ((ρ_h/ρ_0)gy, ∇·v_h)$ is trivially assembled from local cell matrices and guarantees exact potential energy conservation.

5. Lattice Boltzmann and High-Order Moment Approaches

In the LBM framework, gravity is injected as a Hermite polynomial expansion of the body force source (Li et al., 2024). The first- and second-order Hermite moments for body forces set $s^{(1)}_\alpha = F_\alpha$, $$s^{(2)}_{\alpha \beta}=F_\alpha u_\beta+F_\beta u_\alpha$$; all higher moments are zero. The discrete source $S_i$ is then evaluated on each lattice direction via analytical formulae, and the collision-streaming step employs a second-order trapezoidal update. No modification to the collision operator (BGK) is needed, and the resulting scheme exhibits second-order accuracy, as established by Chapman-Enskog analysis. Benchmark validation confirms exact hydrostatic and gravity-driven flow profiles.

6. Conservation, Curl-Free Gravity, and Validation

Rigorous analyses prove that in conservative formulations (Mullen et al., 2020), the discrete divergence of gravitational stress is exactly equivalent to the discrete gravitational source. Summing over all cells (with appropriate boundary conditions) gives total momentum and energy conservation to machine precision. The curl-free property of gravity is maintained by diagonal-component prescriptions for stress tensors, suppressing grid-aligned circulations and anomalous accelerations typical in naive discretizations.

Validation studies (Spitzer-sheet, Jeans instability, polytropes, glacier-fjord flows) confirm second-order accuracy, well-balanced maintenance over large time intervals, reduction in numerical diffusion, and robust resolution of small perturbations. Schemes demonstrate error levels down to machine roundoff, consistent and sharp energy estimates, and stable performance over physically relevant convective and stratified regimes.

7. Algorithmic Implementations and Practical Considerations

Algorithmic frameworks span stage-wise finite volume updates, DG strong/weak forms with quadrature-matched projections, central-upwind and unstaggered central schemes with equilibrium subtraction, and Hermite-moment-driven LBM kernels. For high-order accuracy, polynomial or piecewise hydrostatic reconstructions are used, with positivity limiters and adaptive viscosity near equilibrium. In multidimensional and multiphysics contexts (MHD, AMR grids), special care is taken to maintain divergence-free magnetic fields, well-balanced mass/momentum/energy, and compatibility with mesh adaptivity and arbitrary boundary conditions. Time integration is commonly via SSP RK methods with source-term consistency and stage-averaged gravity fields.

In summary, gravity source term discretization methods comprise a comprehensive set of mathematical, algorithmic, and physical principles enabling exact conservation, well-balancedness for equilibria, entropy stability, and robust numerical performance across the spectrum of computational fluid, atmospheric, astrophysical, and geophysical simulation paradigms (Mullen et al., 2020, Lundgren et al., 2024, Varma et al., 2018, Chertock et al., 2017, Thomann et al., 2019, Berberich et al., 2018, Chandrashekar et al., 2015, Kanbar et al., 2022, Li et al., 2024, Berthon et al., 2024, Liu et al., 29 Jul 2025).