Solver method tuning

Aragog dispatches stiff time integration through _EnergyParameters.solver_method (see src/aragog/parser.py). Three values are accepted: "cvode" (default), "radau", and "bdf". This page is the practical recipe for picking one and tuning the associated tolerances.

Available integrators

If solver_method = "cvode" is selected and scikits_odes is missing, the solver logs a warning and falls back to Radau. See Installation: production solver path.

Picking tolerances

The relevant keys are solver.atol (absolute) and solver.rtol (relative); atol is floored at \(10^{-8}\) inside the solver, and rtol is used as given. Suggested starting points:

When a CVODE call stalls or runs far slower than expected, check whether CVODE is locked at the convective switch (the convective mask and the \(\kappa_h\) floor). In a lock the step stays at one small value and the entropy difference across one basic node oscillates with period 2 (around zero in the measured locks, all at \(D = 0\)); in the measured locks CVODE used about 1.4 RHS evaluations per step and evaluated its Jacobian about every 50 to 60 steps (with INFO logging, the CVODE stats line shows rhs_wrap/nst and a small last_step). Two of the runs counted below were run again with every accepted step recorded: both kept a constant step with a period-2 oscillation of the entropy difference at mid-mantle nodes, at 1.41 RHS evaluations per step on the energy_balance core (the difference stayed just below zero) and at 1.14 on the quasi_steady core (around zero); the other counted runs were not checked for the pattern. A lock can run for a long time without an error, or end in CV_TOO_MUCH_WORK once max_steps is used up in one output interval. Locks were seen at rtol 1e-6 and not at 1e-8; the remedy is rtol and atol at 1e-8, and a CVODE run logs one warning when solver.rtol is above 1e-8. The counts below were measured at the default convective_threshold = 0, where the switch is at \(dS/dr = 0\); the probe did not cover \(D > 0\), where CVODE also stalled at 1e-8 (see Heat transport). They come from a probe variant of the energy_balance isentropic start of tests/test_phase_boundary_cap.py (24 nodes, 200 yr, fixed cap unless stated) that draws the initial entropy uniformly between 8000 and 8400 J/kg/K for each run and adds a relative noise of \(10^{-15}\); the probe script is not in the repository, and the runs were on macOS arm64. A run counts as locked when it has not reached 200 yr after 600 s of wall time: every run was screened at 60 s, and each run that missed that was run again alone to 600 s, by which time a locked run had used about \(1.7 \times 10^5\) RHS evaluations. With rtol = atol = 1e-6, CVODE locked in 12 of 1133 runs, in 4 of 474 on the quasi_steady core and in 2 of 286 with the rate cap, where 2 more runs finished only after 233 s and 346 s. With both at 1e-5, all 80 runs were run to 600 s and none locked (too few runs to show it safe; the slowest finished after 100.6 s and \(6.5 \times 10^4\) RHS evaluations), but the median run needed about \(1.6 \times 10^4\) RHS evaluations, against about \(1.3 \times 10^3\) at 1e-8. With both at 1e-8, CVODE locked in none of 1026, 627 and 1224 runs for three setups (fixed and rate cap on the energy_balance core, fixed cap on quasi_steady). With atol at 1e-8 and rtol at 2e-8, 5e-8 or 1e-7, the values of the PROTEUS retry ladder, it locked in none of the runs: 1308 to 1836 for each value and core with the fixed cap, and 620 and 636 with the rate cap at 5e-8 and 1e-7. These counts hold for the 24-node standalone case, not for a coupled solve. At \(D = 0\) the interior flux is continuous at the switch and only its slope jumps, but with core_bc energy_balance or gradient the convective CMB flux jumps by \(\rho T \kappa_h^\mathrm{floor}\,|\partial S/\partial r|_\mathrm{cmb}\) when the node above the CMB crosses 0; at \(D > 0\) the convective flux stays continuous at 0, but the mixing flux of a mush and the CMB flux jump at 0 and the floor flux jumps at \(D\) (see Heat transport), so the \(D = 0\) counts and the 1e-8 remedy do not carry over to \(D > 0\). On CI the AVX-512 numpy path made the test's own case, at its single fixed initial entropy, stall at 1e-6, judged by the CI test timeout rather than by the 600 s criterion. For a stall without this pattern, raise max_steps; the PROTEUS retry ladder raises max_steps and relaxes rtol, and the relaxed rtol can itself be above 1e-8 (see PROTEUS coupling). If the stall persists with an effective rtol of 1e-8, the issue is usually a sharp solidus/liquidus crossing, not the integrator class; consider enabling phi_step_cap. With convective_threshold above 0 it can also be a node whose gradient stays near \(dS/dr = 0\); the default 0 does not have that jump.

CVODE specifics

CVODE is selected by default and is the only path that supports the JAX-traced analytic Jacobian (use_jax_jacobian = true). With JAX absent, CVODE falls back to a finite-difference Jacobian: correct, but \(O(N)\) RHS evaluations per Jacobian build and noisier on stiff profiles. See the CVODE and JAX explainer.

The tolerance warning also names a configured atol above 1e-8, but only on the energy_balance core, the one setup whose lock data show a loose atol adding to the lock. No lock data exist for the other core boundary conditions, core_module included, so the warning makes no atol claim for them.

CVODE also accepts a SUNDIALS root function for melt-fraction step capping; see phi_step_cap how-to. The scipy fallbacks register the equivalent capping logic as a solve_ivp event with terminal=True, so the same TOML config behaves consistently.

Radau and BDF specifics

Both scipy paths use solve_ivp with dense_output=False and a 1 yr max_step near phase boundaries. They do not support the JAX-traced Jacobian; a Jacobian is computed by scipy via finite differences when needed.

With phase_boundary_cap = "fixed" the 1 yr max_step activates at the same trigger on all three integrators. The default "rate" runs only on CVODE: it arms within the larger of phase_boundary_entropy_margin and the stiff-zone half-width and integrates the call in segments; the scipy paths fall back to the 1 yr step.

EOS floor

The entropy state of CVODE has no lower bound, but the EOS tables stop at a minimum entropy \(S_\mathrm{min}\); below it every lookup is clamped to the table edge, so a cell keeps losing entropy at a fixed temperature. This happens when a prescribed surface flux, fixed over a coupling step, drains the top cell after it crystallises. The EOS floor (energy.eos_floor_margin \(m\)) is off by default (\(m = 0\)), so a drain past \(S_\mathrm{min}\) continues at the clamped table-edge temperature. A positive margin below the entropy span of the tables turns it on; 100 J kg\(^{-1}\) K\(^{-1}\) is the tested value. PROTEUS has no key for it, so PROTEUS runs have the floor off. When on, the floor is one more component of the step-cap root function. It watches the smaller of the top basic-node entropy, which sets the surface temperature, and the smallest cell entropy (the CMB basic node is not watched), and it ends the call when that value reaches \(S_\mathrm{min} + m\), or \(S_\mathrm{min} + m/2\) when the call starts below \(S_\mathrm{min} + m\). A start counts as above the level \(S_\mathrm{min} + m\) or \(S_\mathrm{min} + m/2\) only when it lies more than \(10^{-6} m\) above it; the table boundary below uses \(2 \times 10^{-6} m\). CVODE then returns status 0 with a shorter dt_actual, scipy status 1 at the event, and the stop is logged as a WARNING with its level. A call that starts at or below \(S_\mathrm{min} + m/2\) arms a level at half its distance to \(S_\mathrm{min}\); it does not advance (status 1, dt_actual = 0) when the rate of that entropy at the call start, taken over every component no more than \(10^{-6} m\) above that entropy, reaches the level before the call ends. Above \(S_\mathrm{min} + 2 \times 10^{-6} m\) the rate of the top basic node is the larger of its own rate and the rate of the top cell: the node is a linear extrapolation of the top two cells and can fall when only the gradient below it steepens. Inside that band its own rate counts. A start at or below \(S_\mathrm{min} + 2 \times 10^{-6} m\) also arms half its distance to \(S_\mathrm{min}\), but does not advance on any fall of its lowest component, however slow, at every span, nor when no level fits strictly between \(S_\mathrm{min}\) and the start. Every start at or below \(S_\mathrm{min}\) is refused, also a top basic node extrapolated below \(S_\mathrm{min}\) from a steep top gradient while every cell is inside the table and the node or the top cell rises, so no call advances from a top node outside the table. A top node that falls only by extrapolation while the top cell warms therefore advances above the band in calls that each stop at half its distance, at most \(\lceil \log_2 (d_0 / 2 \times 10^{-6} m) \rceil\) calls that stop at the floor from a distance \(d_0 \le m/2\) (a call too short to reach its level adds one that does not stop), and is refused once inside the band. A call with no span (start equal to end) skips the floor. Between \(S_\mathrm{min} + 2 \times 10^{-6} m\) and \(S_\mathrm{min} + m/2\) the refusal depends on the call span, so a caller that retries with a shorter span can then advance. A refused or discarded call returns status 1, dt_actual = 0 and its start state, with a WARNING. PROTEUS runs have the floor off; with the floor on, PROTEUS would treat such a call as a failed attempt: its retry ladder retries with half the span (and up to 5 times the atol), which can then advance in that range, so under a drain that does not stop the surface entropy can come closer to \(S_\mathrm{min}\) than \(m/2\); at the table minimum a shorter span does not advance a falling start. When the ladder runs out, PROTEUS would keep the interior state of the previous step while the other modules advance, and stop the run at the third failure in a row. Above \(S_\mathrm{min} + 2 \times 10^{-6} m\), a drain too slow to reach half the distance in the call, a pause or a warming surface continues. The refusal is logged as a WARNING that says the entropy is at or below the table minimum, leaves no level between it and the table minimum, falls at the table minimum, or falls too fast between the band and \(S_\mathrm{min} + m/2\), also when the start rate is NaN. A floor stop at the call start, where the level lies within the root location error of the start, is a refusal with its own WARNING.

With \(m > 0\) the floor is armed in every solve with the entropy EOS, except on the gradient core; a call in which it does not fire takes the same steps and gives the same result as without it. CVODE locates a root to a time tolerance of about \(100\,u\,(|t| + |h|)\), with \(u\) the unit roundoff, \(t\) the absolute call time and \(h\) the step, so a stop can lie below its level by up to about \(|dS/dt| \times 2.2 \times 10^{-14} \times (|t| + |h|)\) (\(t\) and \(h\) in yr), about \(|dS/dt| \times 2.2 \times 10^{-14} \times t\) for steps much shorter than \(t\); scipy event location also loses precision with \(t\), more slowly; for the 11-cell state of the floor tests under \(10^6\) W m\(^{-2}\), a stop at \(10^9\) yr lay \(1.6 \times 10^{-3}\) J kg\(^{-1}\) K\(^{-1}\) below its level on CVODE, \(6 \times 10^{-6}\) below it on Radau and \(9 \times 10^{-6}\) above it on BDF (a scipy event can stop on either side of its level). A call whose returned states or energy trace hold a top basic-node or cell entropy below \(S_\mathrm{min}\) is discarded: it is solved again as a refusal, so its state, time and energy trace stay at the call start, and the WARNING names the time after the call start, and that entropy and the level as offsets from \(S_\mathrm{min}\). So, with the floor on, every state a call returns is inside the table when the call starts inside it; with the floor off (\(m \le 0\), the default) returned states can lie below \(S_\mathrm{min}\). Trial steps and Newton iterates of the integrator can still evaluate the right-hand side below \(S_\mathrm{min}\), with clamped lookups, at any margin; with the floor on, the returned states and the energy trace cannot hold such a state. The floor guards the table range, not its resolution: where a table is flat in temperature above \(S_\mathrm{min}\), a stop at the level still gives the clamped edge temperature. In the solid table of the floor tests, \(T\) is the same over the first four entropy rows (\(S_\mathrm{min}\) to \(S_\mathrm{min} + 322\) J kg\(^{-1}\) K\(^{-1}\)) at every pressure up to 6 GPa, so with \(m = 100\) the floor keeps the state inside the table but the surface temperature at the edge value; in one PALEOS solid table at the lowest pressure, \(T\) is the same over the rows from about \(S_\mathrm{min} + 70\) to \(S_\mathrm{min} + 620\) J kg\(^{-1}\) K\(^{-1}\).

Cross-references

• CVODE and JAX explainer: why the analytic Jacobian helps and how set_jax_cvode_factory is wired.
• phi_step_cap how-to: SUNDIALS-rootfn capping of \(|\Delta\phi|\) across a step.
• Installation: production solver path: runtime requirement for use_jax_jacobian = true.