Energy-conservation diagnostics

Aragog emits a two-track energy-conservation diagnostic on every solver call:

1. A physical conservation budget that compares the integrated mantle enthalpy against the cumulative time integral of the boundary fluxes and volumetric sources.
2. A solver-correctness check that verifies the entropy ODE itself was integrated to its stated tolerance, independent of any thermodynamic-frame considerations.

The first track answers did the simulation conserve energy as a physical principle, the second answers did CVODE solve its ODE accurately. Both are exposed as columns in the PROTEUS runtime_helpfile.csv and as fields on Aragog's SolverOutput dataclass.

Per-call integrals (SolverOutput)

Every call to EntropySolver.solve() returns a SolverOutput with the energy-related integrals computed over the actual CVODE sub-step trajectory.

All integrals use the trapezoidal rule on the integrator's accepted sub-steps. Integrating over the actual CVODE trajectory rather than between end-of-step snapshots is required because a single CVODE phase-boundary spike on a coarse end-of-step interpolation can shift the integral by orders of magnitude.

The two surface-flux integrals (step_dE_F_int_J and step_dE_F_cmb_J) are area-weighted only and have no mass-weighting variant. The volumetric source integrals come in two flavours, distinguished by the mass weighting on \(Q_\text{radio} = \rho\,h_\text{radio}\,V\) and \(Q_\text{tidal} = \rho\,h_\text{tidal}\,V\):

• State-mass: \(\rho_\text{state} = \rho(P_\text{stag},\,S(t))\), evolves with the solver state. Used for instantaneous power diagnostics.
• Frozen mass: \(\rho_\text{struct}\,V\) from mesh.staggered_effective_density × mesh.basic.volume, fixed at IC. Used in the conservation-grade integrated enthalpy E_state_cons (see below).

State-heat and compression terms

Two further SolverOutput fields describe where the mantle's internal energy sits and how it moves. They are easy to conflate because both carry units of joules and both go through the equation of state, but only one enters the conservation budget.

step_dE_state_heat_J is the change in the mantle's Lagrangian heat content over the call,

evaluated by direct EOS quadrature along each staggered cell's entropy path with its pressure \(P_i\) held fixed over the call. This is the state side of the conservation budget: the capacitance weighting \(\rho T\) is exactly what the entropy ODE transports, so the coupler differences its cumulative sum against the integrated boundary fluxes and volumetric sources. It is deliberately not an enthalpy change. The two-phase EOS specific enthalpy does not satisfy \(dh = T\,dS\) across the mushy zone, so an enthalpy-snapshot difference carries a spurious flow-work term and cannot close the budget. Sign: negative when the mantle cools.

step_dE_compression_J is the adiabatic compression work from the structure re-solve that precedes the call,

taken at fixed entropy \(S\) and frozen mass: the enthalpy change each mass element sees when the static pressure at its position shifts between the old and new equilibrium structures. It is informational only and is deliberately not added to the predicted energy, for two reasons.

• The conservation residual already closes against the state-heat term \(\sum \rho T\,dS\), which carries the same thermodynamic content as the \(P\,dV\) work; adding compression to the predicted side would double-count it.
• The two-phase specific enthalpy does not satisfy \(dh = T\,dS\) across the mushy zone, so an enthalpy-snapshot difference cannot close the budget on its own in any case.

step_dE_compression_J is zero for a static structure, where the pressures are unchanged, and for the first solve, where there is no prior mesh to difference against.

Cumulative conservation residuals (PROTEUS helpfile)

The PROTEUS coupling layer (proteus.utils.coupler._populate_energy_residual) accumulates the per-call integrals into running cumulative columns on runtime_helpfile.csv:

E_residual_cons_J is the physical conservation residual: the difference between the actual change in mantle enthalpy and the integrated boundary-and-source budget. solver_residual_J is the solver-correctness residual: by construction the entropy ODE satisfies \(\sum_\text{cells}\rho T\,(\partial S/\partial t)\,V = -F_\text{int}\,A_\text{int} + F_\text{cmb}\,A_\text{cmb} + Q_\text{radio} + Q_\text{tidal}\) at every accepted CVODE sub-step, so any drift in solver_residual_J flags real CVODE step rejection or atol/rtol issues.

Closed-mantle energy balance

Aragog's closed-mantle energy balance is

A volumetric-work source \(\Phi_\text{vol}\) does not appear in this budget: the volumetric work done when a melt of different density is transported across a pressure gradient is already implicit in the divergence of the \(\Delta h\)-weighted mass-flux contributions to _heat_flux. By definition \(\Delta h = \Delta u + P\,\Delta v\), and on a hydrostatic column \(\partial \Delta h/\partial r \supset \Delta v\,\partial P/\partial r = -\rho g\,\Delta v\), so \(-\partial/\partial r(j\,\Delta h)\) already carries the volumetric-work term. A separate \(\Phi_\text{vol}\) source would double-count (Bower et al. 2018¹ §3, SPIDER energy.c); a negative regression test in tests/test_jax_no_phi_vol_source.py enforces the absence.

Why frozen-mass weighting

E_state_cons_J integrates the EOS-consistent specific enthalpy \(h(P,S)\) table against a frozen mass-per-shell profile \(\rho_\text{struct}\,V\) taken once from the equilibrium structure at IC, not against the time-varying state-dependent mass \(\rho(P,S(t))\,V\).

The reason is dimensional consistency with the entropy ODE's frame. The entropy ODE evolves in fixed-volume Eulerian coordinates with state-dependent \(\rho(P,S(t))\). Its per-cell discretisation gives

and summing over cells yields the surface telescope \(\sum \rho_i T_i (\partial S_i/\partial t)\,V_i = -F_\text{int}\,A_\text{int} + F_\text{cmb}\,A_\text{cmb} + Q_\text{radio} + Q_\text{tidal}\). The mass weighting in this identity is \(\rho(t)\,V\) at every instant.

If the conservation diagnostic uses state-dependent mass too, the actual integrated enthalpy is

whose time derivative picks up an \(h\,\partial\rho/\partial t\) cross term that does not appear in the divergence-of-flux + sources budget on the right-hand side:

The frame artefact grows with mantle cooling because \(\rho(P,S(t))\) increases as the mantle solidifies. On a static-Zalmoxis dry \(a_1\) benchmark, a state-mass residual reaches 51 % of total cooling at 100 kyr, a non-physical signal that masks any real conservation drift.

E_state_cons_J uses the frozen mass-per-shell, eliminating \(h\,\partial\rho/\partial t\) from the diagnostic side. The same trajectory closes to 5.3 % of total cooling at 100 kyr, an order of magnitude smaller and physically meaningful. E_residual_cons_J and E_residual_cons_frac are therefore the canonical conservation columns.

Why two diagnostics

The frozen-mass residual E_residual_cons_frac answers a slightly different question than the solver residual solver_residual_J. Both are useful:

• E_residual_cons_frac quantifies the gap between the physical enthalpy change and the integrated divergence-of-flux + sources budget. Part of it is a frame-correction component: because \(\rho_\text{struct}\,V \ne \rho(t)\,V\), the frozen-mass enthalpy time-derivative is not exactly the entropy-budget LHS. That component is expected accounting from the frozen-mass convention, not lost energy, and it grows as the mantle solidifies. While cooling is fast the frame correction stays small and the fraction dilutes as cumulative cooling accumulates; a fraction that grows over that fast-cooling span signals a real conservation regression.

• solver_residual_J is the cumulative integral of the per-substep mismatch between \(\sum \rho T\,(\partial S/\partial t)\,V\) and the entropy-budget RHS, evaluated at the same state-dependent \(\rho\). By construction the entropy ODE satisfies this identity at every accepted CVODE step; any drift here is purely the gap between the trapezoidal accumulation Aragog uses for diagnostics and CVODE's own internal step accuracy. Closes to ~10⁻⁷ of total cooling on the verified benchmark, at the machine-precision noise floor. A non-trivial drift here flags real solver issues (step rejection, atol/rtol problems, hidden NaN propagation).

In production CHILI runs you should monitor both columns. E_residual_cons_frac is the headline "is the simulation conserving energy?" check; solver_residual_J is the gate that catches numerical pathologies before they show up as a misleading conservation drift.

Verification

Reference smoke runs on a static-Zalmoxis dry \(a_1\) (\(1\,M_\oplus\), no atmosphere, no radio, no tidal) trajectory starting from a fully molten mantle at \(T_\text{magma} \approx 4244\,\text{K}\):

The frozen-mass residual is largest early and dilutes as cumulative cooling accumulates; the solver residual sits at the floating-point noise floor across the same span. The frame correction starts near zero at the molten IC and grows only as the mantle solidifies, so at \(\Phi = 0.97\) it is far too small to set the 11.5 % residual seen there: the large early fraction is the residual divided by a still-small cumulative cooling, not a frame effect. Two points constrain the dilution over this fast-cooling span but not how the residual splits between the frame correction and the non-frame error that dominates it; the trend past 100 kyr is not captured here.

Worked example

For a single 100 kyr coupled step with surface flux \(F_\text{int} \sim 10^3\,\text{W/m}^2\) on an Earth-sized planet:

• \(A_\text{int} \approx 5 \times 10^{14}\,\text{m}^2\), so \(\sim 10^3 \cdot 5 \times 10^{14} \cdot 3.16 \times 10^{12} = 1.6 \times 10^{30}\,\text{J}\).
• step_dE_F_int_J should land near \(-1.6 \times 10^{30}\,\text{J}\) for that step.
• step_dE_F_cmb_J is order \(10^{27}\text{-}10^{28}\,\text{J}\) depending on core thermal state.
• step_dE_Q_radio_J and step_dE_Q_tidal_J are typically order \(10^{25}\text{-}10^{27}\,\text{J}\) for primordial Earth.

While cooling is fast, an E_residual_cons_frac that grows instead of diluting, with solver_residual_J / |dE_actual| still at its ~10⁻⁷ floor, points to a missing energy contribution. As the mantle solidifies the frame correction grows, so a later plateau or slow rise in the fraction can be that expected accounting rather than a leak; the solver residual carries no frame correction and is the reliable gate there. A solver_residual_J / |dE_actual| above ~10⁻⁵ points to an integrator-tolerance problem. Plot both columns over time using tools/plot_energy_balance.py (PROTEUS-side):

python tools/plot_energy_balance.py output/<run_dir> --label <run_label>
# writes output_files/energy_balance/<label>_balance.pdf and _powers.pdf