Gravity#
Overview#
Athelas provides Newtonian gravity as an explicit physics package. The default coupled formulation is designed for spherical self-gravity: it holds hydrostatic balance to the order of the scheme and conserves fluid plus gravitational energy in the semi-discrete limit. A constant acceleration model is also available for test problems.
Enable gravity with physics.gravity = true. Coupled gravity is the
recommended mode. gravity.operator_split = true is a deliberately simpler
compatibility mode; it uses a direct velocity kick and gives up well-balancing
and total-energy conservation. The available options are summarized in
Configuration.
Physical model#
Gravity enters through the Newtonian potential \(\varphi\). In spherical symmetry,
where \(M(r)\) is the mass enclosed within \(r\); the constant
model uses \(\varphi(r) = g_{\rm val}\,r\). The acceleration is minus its
gradient,
recovering \(g = -GM/r^2\) and \(g = -g_{\rm val}\) respectively.
Gravity contributes acceleration and work to the Lagrangian fluid equations:
Here \(v\) is velocity and \(E\) is specific total fluid energy.
Numerical model#
Coupled gravity#
A pointwise discretization of these source terms would neither preserve hydrostatic balance nor conserve the discrete potential energy. We seek a representation of gravitational sources that is consistent with the operators used for the hydrodynamic flux divergence and geometric source. Moreover, in the Lagrangian picture, the gravitational potential evolves through mesh advection, so any treatment of gravity must be consistent with mesh motion. The coupled package builds a gravity pressure \(\tilde\Pi\), defined in the continuum by
This is precisely the equation of hydrostatic equilibrium: \(\tilde\Pi\) is the pressure profile that balances the potential \(\varphi\). Equilibrium is therefore the statement \(P = \tilde\Pi + \text{constant}\). In the mass coordinate, dropping factors of \(4\pi\) so that \(\mathrm{d}m = \rho\sqrt{\gamma}\,\mathrm{d}r\) and \(\sqrt{\gamma}=r^2\) in spherical geometry, this reads
Discretely, \(\tilde\Pi\) is built as the gravitational analogue of a pressure. A gas pressure is how the internal energy resists compression; \(\tilde\Pi\) is how the gravitational energy resists it. Let \(W_h=\sum_q M_q\,\varphi(r_q)\) be the discrete gravitational energy (nodal mass \(M_q=w_q\mu_q\) times the potential at each node), and define
Nudging the specific volume \(\tau\) in a cell repacks the gas, which moves the nodes and so shifts where mass sits in the potential; \(\tilde\Pi\) is precisely how much \(W_h\) responds. That is exactly the sense in which the gas and gravity pressures balance in equilibrium.
Computing that response is a single inward sweep. The change in \(W_h\) at a cell accumulates the gravitational pull of everything exterior to it, which the package assembles as one surface-to-centre running sum. Formally that sum is the transpose of the mesh reconstruction \(r(\tau)\) – the map that integrates \(\mu\tau\) to the volume coordinate \(X=r^3/3\). Writing \(J_X=\partial X/\partial\tau\) for its tangent (the same object that advances the grid, \(\dot r=J_X\dot\tau/r^2\)) and \(f=\rho\,\partial_m\varphi=\partial\varphi/\partial X\) for the force per unit volume coordinate,
which is the adjoint of the reconstruction.
Building \(\tilde\Pi\) this way is what makes the scheme conservative. Any \(\tilde\Pi\) obeying the continuum relation above would balance the momentum equation for a hydrostatic state, but only the exact adjoint satisfies \(\dot W_h=\tilde\Pi^{\mathsf T} M\dot\tau\): the gravity work returned to the fluid then equals the rate of change of the discrete potential energy, term for term. That identity is the engine of the energy source below.
Well-balancing#
The momentum source applies the same weak divergence operator to \(\tilde\Pi\) that hydro applies to fluid pressure \(P\). For a discrete state satisfying
the interior (volume) and curvature parts of the two operators cancel term for term to roundoff – the payoff of building \(\tilde\Pi\) as the exact adjoint. What survives is a residual confined to the cell faces, because the two packages compute the face pressure by different rules: gravity uses a single continuous scan value, while hydro uses a Riemann solve of two independently reconstructed traces. These agree at the boundaries – the inner face carries \(\sqrt\gamma = 0\), and the outer face is exact under the surface gauge with a matching surface pressure – but only to truncation order at interior faces.
The balance is therefore high order in the mesh spacing, not exact to roundoff: the spurious equilibrium acceleration converges away under refinement at the order of the scheme. That convergence, rather than a machine-zero residual, is the operative well-balancing criterion. Making the balance exact would additionally require the fluid flux to reuse \(\tilde\Pi\) at the face – a hydrostatic-reconstruction-style well-balanced flux – which is a possible extension, not the current scheme.
Gauge choice#
\(\tilde\Pi\) is defined up to a constant. Both the momentum and energy operators are invariant under \(\tilde\Pi\to\tilde\Pi+c\): the constant cancels in the weak momentum operator, and the energy contribution reduces to the discrete specific-volume identity. We take \(\tilde\Pi=0\) at the stellar surface, where \(P\to0\); this avoids subtractive cancellation in \(P-\tilde\Pi\).
Energy conservation#
The coupled energy source uses a weak product rule,
assembled as (face - volume - dilatation). The face term telescopes across
cells, the summed volume term vanishes, and the dilatation term is
\(\tilde\Pi^{\mathsf T}M\dot\tau\). Therefore the total gravity work is
exactly \(-\mathrm{d}W_h/\mathrm{d}t\) in the semi-discrete sense.
The hydro or IMEX rad-hydro package publishes the complete mass-matrix-scaled
specific-volume RHS in dtau_dt earlier in the same Runge–Kutta stage.
Hydro must be registered before coupled gravity. This same-stage dependency is
why the weak source cannot be used as an operator-split package.
Operator-split gravity#
Split gravity holds the mesh fixed and applies a direct kick after the hydro step. With \(g\) evaluated on that post-hydro mesh, it uses
This leaves specific internal energy unchanged to roundoff. It does not move the mesh or update \(W_h\), so it is not well-balanced and does not conserve fluid plus gravitational energy. The limiter-energy correction is unavailable in split mode.
Note
Operator splitting is not recommended for gravity. The gravity update is not a large computational expense, so there is no real gain and a demonstrable loss of accuracy.
Diagnostics and conservation#
The history file reports Total Gravitational Energy [erg] as the discrete
form of \(W = \int\varphi\,\mathrm{d}m\), using the same potential
\(\varphi\) as the source. Together with total fluid and radiation energy,
this is the conservation audit for coupled gravity.
Limiter mesh work#
Limiters preserve each cell’s mass-weighted average, so they do not move cell
faces. They can nevertheless reshape \(\tau\) within a cell, relocate its
interior nodes, and change \(W_h\) without a physical energy source.
Cumulative Limiter Mesh Work [erg] measures this change for the
end-of-step limiter block on the committed state.
Limiter energy correction#
gravity.limiter_energy_correction returns the committed limiter’s potential
change to the fluid. For each cell it computes
and adds \(-\delta W_i/m_i\) uniformly to the cell’s specific energy. The
correction is clamped at the EOS internal-energy floor. The associated history
columns are Cumulative Limiter Energy Correction [erg] and
Cumulative Limiter Energy Clamp Residual [erg]. They satisfy
Configuration#
Enable the package and select a model:
config.physics = { gravity = true }
config.gravity = {
model = "spherical", -- or "constant"
gval = 1000.0, -- acceleration [cm s^-2], for "constant"
operator_split = false, -- coupled gravity is recommended
limiter_energy_correction = false, -- opt-in
}
See Input Deck Reference for the complete option table.