First-run tutorial

Welcome to this first-run tutorial. We will load a configuration, integrate a single mantle cooling problem with EntropySolver, and inspect the resulting profiles and scalar diagnostics.

What does Aragog do?

Aragog evolves the specific entropy \(S(r,t)\) at staggered nodes inside a spherical mantle shell. The solver couples conduction, mixing-length convection, gravitational separation of melt, chemical mixing, radiogenic heating, and tidal heating. Temperature, density, and melt fraction are derived from \((P, S)\) via a tabulated equation of state.

Goals

Reach a first successful integration of a coupled-mantle setup and read its results from the SolverOutput dataclass.

Assumptions

• You are using Python 3.12 or 3.13.
• Aragog is installed (see Installation).
• A directory of SPIDER-format pressure-entropy EOS tables is available on disk (see Reference: data for the file format and the canonical PALEOS data sources).

1. Locate or stage the EOS tables

Aragog's entropy solver requires a directory containing the following files:

For a quick standalone walkthrough you can fetch a small SPIDER-format set bundled as a release asset:

mkdir -p /tmp/aragog-test-data
curl -sL -o /tmp/spider_eos.tar.gz \
    https://github.com/FormingWorlds/aragog/releases/download/test-data-v1/spider_eos_test_data.tar.gz
tar xzf /tmp/spider_eos.tar.gz -C /tmp/aragog-test-data/
export ARAGOG_TEST_EOS_DIR=/tmp/aragog-test-data/spider_eos

The same tarball is what the nightly CI populates before the smoke + slow tiers run. In a PROTEUS coupled run these tables are produced from the configured PALEOS or Wolf-Bower P-T file by the PROTEUS wrapper. For standalone work, point eos_dir at any directory containing this set.

2. Scaffold a configuration

The fastest path is the aragog new scaffold:

aragog new first --from abe_solid

This copies the bundled abe_solid template to first.toml in the cwd. Open the file and edit the file paths in [phase_solid], [phase_liquid], and [phase_mixed] so that they point at your local data; nothing else needs to change for a smoke run. A configuration that carries a [scalings] section is rejected at load time; if you copied a legacy config, remove that block.

Confirm the parsed schema is what you expect before solving:

aragog validate first.toml
# OK first.toml: core_bc=energy_balance, IBC=2, OBC=1, eos_method=1, IC_method=1, radionuclides=4

3. Solve from the CLI

# Replace $FWL_DATA/spider/lookup-fs with your local EOS directory if
# the env var is not set.
aragog run first.toml \
    --eos-dir $FWL_DATA/spider/lookup-fs \
    --out first.nc

aragog run mirrors the recipe at the bottom of this page (EntropySolver.from_file → initialize → set_initial_dSdr_cmb → set_initial_entropy → solve → to_netcdf); see Reference: CLI for the full option list. A successful run leaves aragog.log plus first.nc in the working directory.

The CLI derives the initial entropy \(S_0\) from [initial_condition].surface_temperature in the scaffolded config (abe_solid sets it to 3600 K); pass --initial-entropy <S0> to use an explicit value instead.

Inspect the result:

aragog inspect first.nc

The default human format prints status, the headline scalars (T_magma, T_core, Phi_global, F_heat_total, ...), the mesh dimensions, and a S_final / T_basic profile-range block. Add --json for jq-friendly machine-readable output.

For one-off tweaks during debugging, aragog run accepts repeatable --set <key.path>=<value> flags (see the CLI reference for type coercion rules):

aragog run first.toml --eos-dir ... --initial-entropy 2900.0 \
    --set energy.kappah_floor=20.0 --set solver.atol=1e-11

4. Solve from Python

For scripted workflows, parameter sweeps, or richer plotting you can call the same path the CLI uses directly. Create first.py:

from pathlib import Path
from aragog import aragog_file_logger
from aragog.solver import EntropySolver

# Set up combined console + file logging in the current directory.
aragog_file_logger(log_dir=str(Path.cwd()))

# Build the solver from a TOML config and an EOS-table directory.
solver = EntropySolver.from_file(
    filename="first.toml",
    eos_dir="path/to/eos/tables",
)

solver.initialize()

# Set the initial entropy at staggered nodes (J/kg/K). A scalar
# value produces a uniform isentropic profile.
solver.set_initial_entropy(2900.0)

solver.solve()

out = solver.get_state()

print("=== Aragog run summary ===")
print(f"Status:       {out.status}")
print(f"T_magma:      {out.T_magma:.0f} K")
print(f"T_core:       {out.T_core:.0f} K")
print(f"Phi_global:   {out.Phi_global:.4f}")
print(f"M_mantle:     {out.M_mantle:.3e} kg")
print(f"E_th:         {out.E_th:.3e} J")
print(f"F_heat_total: {out.F_heat_total:.3e} W/m^2")
print(f"Cp_eff:       {out.Cp_eff:.0f} J/kg/K")

Run it:

python first.py

out.S_final, out.T_stag, out.phi_stag, and the basic-node fluxes (out.jcond_b, out.jconv_b, out.jgrav_b, out.jmix_b) are NumPy arrays you can plot directly.

5. Plot the entropy and melt-fraction profiles

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True)

axes[0].plot(out.S_final, out.r_stag * 1e-3)
axes[0].set_xlabel("Specific entropy [J/kg/K]")
axes[0].set_ylabel("Radius [km]")

axes[1].plot(out.phi_stag, out.r_stag * 1e-3)
axes[1].set_xlabel("Melt fraction")

fig.tight_layout()
fig.savefig("first_profiles.pdf")

6. Use Aragog inside PROTEUS

For a coupled atmosphere-interior simulation the configuration is built programmatically by the PROTEUS wrapper at src/proteus/interior_energetics/aragog.py. The wrapper drives EntropySolver.set_initial_entropy() from the previous step's profile, supplies the four-column external mesh file from Zalmoxis when eos_method = 2, and reads SolverOutput back into the PROTEUS Interior_t state. See Standalone vs PROTEUS-integrated usage for the path comparison.

For a coupled walkthrough (atmosphere + interior + outgassing), see the PROTEUS usage guide.

7. Troubleshooting

status = -1 from solve(). The solver hit an integration failure: typically the integrator collapsed its step size at a phase boundary or the EOS table was queried outside its \((P, S)\) domain. PROTEUS handles these via a retry ladder that calls set_initial_dSdr_cmb and a tolerance-relaxation setting. In standalone use, inspect aragog.log for the warning trail and consider enabling the SUNDIALS CVODE path with solver_method = "cvode" if scikits.odes is installed. aragog run writes the snapshot and exits 1 on this status; the Python API only sets it, and SolverOutput.failed is true.

Slow integration. Reduce number_of_nodes to 40-80 for first runs, or shorten the time window. Loosening rtol above 1e-8 makes Aragog log one warning that a loose rtol can lock CVODE at the convective switch (the convective mask and the \(\kappa_h\) floor). 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). 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. Above 1e-7 the rate phase-boundary cap also logs a warning, because its accuracy is measured at 1e-8. The phase-aware max_step cap activates automatically near solidus and liquidus crossings.

Entropy out of table range. Check that the EOS table covers the expected \((P, S)\) envelope. For PALEOS tables sampled at 150 points per decade, the entropy axis typically spans roughly [-100, 5500] J/kg/K for MgSiO₃; values outside that range usually indicate an incorrect initial entropy or an oversharp boundary flux.