Configuration
Aragog reads its input from a single configuration file. Two formats are accepted: TOML (.toml, recommended for new work) and INI (.cfg).
When Aragog is invoked from PROTEUS the configuration is built programmatically by the PROTEUS wrapper and the file-based path is bypassed; the same key names and value semantics apply in both cases.
TOML configuration
A TOML configuration file specifies all model parameters grouped by section. Below is a representative coupled-mantle setup; the reference tables that follow document every key.
[solver]
start_time = 0
end_time = 1.0e9 # yr
atol = 1e-9
rtol = 1e-9
tsurf_poststep_change = 30
cvode_output_points = 65 # dense-output grid per CVODE macro-step
max_steps = 100000 # internal CVODE steps per output interval
# tcore_change_limit = 3000 # optional; flag a per-solve core-T change above this [K] (off by default)
[boundary_conditions]
outer_boundary_condition = 4 # 4 = prescribed flux (PROTEUS coupling)
outer_boundary_value = 280.0 # W/m^2
inner_boundary_condition = 1 # 1 = core cooling
inner_boundary_value = 0.0 # ignored when inner_bc=1
emissivity = 1.0
equilibrium_temperature = 273.0
core_heat_capacity = 880.0
tfac_core_avg = 1.147
param_utbl = false
param_utbl_const = 1.0e-7
core_bc = "energy_balance" # CMB BC mode; default matches PROTEUS production
[mesh]
outer_radius = 6.371e6 # m
inner_radius = 3.480e6 # m
number_of_nodes = 200
mixing_length_profile = "nearest_boundary"
core_density = 11000.0 # kg/m^3
surface_density = 4078.95095544 # kg/m^3 (PROTEUS production default; SPIDER -adams_williamson_rhos)
gravitational_acceleration = 9.81 # m/s^2 (scalar fallback)
adiabatic_bulk_modulus = 260e9 # Pa
adams_williamson_beta = 0.0
surface_pressure = 0.0 # Pa
mass_coordinates = true # PROTEUS production default; SPIDER mesh layout
eos_method = 1 # 1 = Adams-Williamson, 2 = external file
eos_file = "" # used when eos_method = 2
[energy]
conduction = true
convection = true
gravitational_separation = true
mixing = true
radionuclides = true
tidal = false
eddy_diffusivity_thermal = 1.0
eddy_diffusivity_chemical = 1.0
kappah_floor = 10.0 # PROTEUS production default; 0.0 = textbook MLT
convective_threshold = 0.0 # 0.0 = SPIDER switch; > 0: minimum drive gradient and floor band [J/kg/K/m]
phi_step_cap = 0.0 # 0 = disabled; 0.05 caps per-call melt-fraction change in mushy band
eos_floor_margin = 0.0 # 0 (default) = off; > 0: J/kg/K above the entropy table minimum, 100 tested
phase_boundary_cap = "rate" # "rate" (default) or "fixed" (1 yr step near a phase boundary)
bottom_up_grav_sep = true
phase_smoothing = "tanh" # "tanh" (default, SPIDER-parity) or "cubic_hermite"
solver_method = "cvode" # "cvode" | "radau" | "bdf"
use_jax_jacobian = true
tidal_array = [0.0]
[initial_condition]
initial_condition = 1 # 1=linear T, 2=user file, 3=adiabat
surface_temperature = 3600.0 # K
basal_temperature = 4000.0 # K
init_file = "" # used when initial_condition = 2
[phase_solid]
density = "data/lookup/density_solid.dat"
heat_capacity = "data/lookup/heat_capacity_solid.dat"
melt_fraction = 0.0
thermal_conductivity = 4.0
thermal_expansivity = 1.0e-5
viscosity = 1.0e21
[phase_liquid]
density = "data/lookup/density_melt.dat"
heat_capacity = "data/lookup/heat_capacity_melt.dat"
melt_fraction = 1.0
thermal_conductivity = 4.0
thermal_expansivity = 1.0e-5
viscosity = 1.0e2
[phase_mixed]
latent_heat_of_fusion = 4.0e5
rheological_transition_melt_fraction = 0.4
rheological_transition_width = 0.15
solidus = "data/lookup/solidus.dat"
liquidus = "data/lookup/liquidus.dat"
phase = "composite"
phase_transition_width = 0.01
grain_size = 1.0e-3
matprop_smooth_width = 0.0
separation_viscosity = "melt" # "melt" or "mixture"
cp_blend = "latent" # "latent" or "linear"
const_properties = false
[radionuclide_K40]
name = "K40"
t0_years = 4.55e9
abundance = 1.1668e-4
concentration = 310.0
heat_production = 2.8761e-5
half_life_years = 1.248e9
Any number of [radionuclide_*] sections can be added; each section name must start with radionuclide_.
Each isotope heats at heat_production * abundance * 1e-6 * concentration * 2**((t0_years - t) / half_life_years) per kg, with concentration in ppm; an isotope whose amplitude is zero contributes no heating at any time.
With radionuclides = true, loading raises a ValueError when a half-life is not positive, when heat_production, abundance or concentration is negative, or when an isotope's heating is not finite at start_time.
Configuration sections
A configuration file containing a [scalings] block is rejected at load time with a ValueError.
Aragog applies its internal nondimensionalisation around the integrator only (entropy reference \(S_\mathrm{ref}\) and time reference \(t_\mathrm{ref}\)); it is not user-configurable.
[solver]
Time-integration controls.
| Key | Unit | Description |
|---|---|---|
start_time |
yr | Start of integration |
end_time |
yr | End of integration |
atol |
-- | Absolute tolerance, default \(10^{-8}\) (floored at \(10^{-8}\)) |
rtol |
-- | Relative tolerance, default \(10^{-8}\), the value at which phase_boundary_cap = "rate" is verified. Above \(10^{-8}\) a CVODE run logs one warning: 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. |
tsurf_poststep_change |
K | Maximum allowed surface-temperature change per coupling step (PROTEUS use) |
cvode_output_points |
-- | Number of points on the CVODE dense-output grid returned per macro-step (default 65, minimum 2). The per-call energy integrals take CVODE's accepted internal steps and these output points as nodes, so a coarse grid does not under-resolve them, and their cost follows the number of accepted steps, not this grid. The core-temperature change check (tcore_change_limit) reads only the output points. The grid feeds back into CVODE stepping, so the accepted step count and the final state shift weakly with it. Used only when solver_method = "cvode". |
max_steps |
-- | Maximum number of internal CVODE steps between two consecutive output times (SUNDIALS mxstep; default 100000, minimum 1). CVODE returns CV_TOO_MUCH_WORK and stops once one output interval reaches this count. With phase_boundary_cap = "rate" each CVODE segment applies the same per-interval limit. A stiff phase-change window can need more internal steps than the default budget; raise this value to let such a solve complete. Used only when solver_method = "cvode". |
tcore_change_limit |
K | Optional limit on the per-solve core-temperature change (unset by default; must be positive when set). Aragog always measures the largest change of the core temperature from the solve-entry value over the returned grid and reports it as tcore_change_max. When this limit is set and that change exceeds it, Aragog sets the tcore_change_exceeded flag on the result; the flag is also set whenever any sampled core temperature is non-finite, independent of whether a limit is set, since that signals a corrupted solve. The solve never raises. A caller can use the flag to reject a solve whose core temperature jumps as it crosses a phase boundary in a single accepted step, or whose core temperature is non-finite. |
[boundary_conditions]
Thermal boundary conditions at the surface and CMB.
| Key | Unit | Description |
|---|---|---|
outer_boundary_condition |
int | Surface BC mode. 1 = grey-body atmosphere, 4 = prescribed flux, 5 = prescribed temperature |
outer_boundary_value |
W/m² or K | Surface flux (modes 1, 4) or temperature (mode 5) |
inner_boundary_condition |
int | CMB BC mode. 1 = core cooling, 2 = prescribed flux, 3 = prescribed temperature |
inner_boundary_value |
W/m² or K | CMB flux (modes 1, 2) or temperature (mode 3) |
emissivity |
-- | Surface emissivity (used in mode 1) |
equilibrium_temperature |
K | Radiative equilibrium temperature (mode 1) |
core_heat_capacity |
J/kg/K | Core specific heat capacity |
tfac_core_avg |
-- | Core adiabat correction factor (default 1.147; Bower+2018 Table 2) |
param_utbl |
bool | Enable upper-thermal-boundary-layer parameterisation (default false) |
param_utbl_const |
Kā»Ā² | UTBL constant in \(\Delta T = b\,T^3\) |
core_bc |
str | CMB BC mode. energy_balance (SPIDER bit-parity, length \(N+1\), default and PROTEUS production), quasi_steady (alpha-factor, state length \(N\)), gradient (entropy gradient as state, length \(N+2\), experimental), bower2018 (parity-only, not recommended), or core_module (the core evolution module, length \(N+2\) plus the shell cells with stratification = true, needs inner_boundary_condition = 1). Any other value is refused. |
top_node |
str | Entropy of the top basic node. linear (default) extrapolates the top 2 cells; limited reduces that half step smoothly near an edge of the EOS entropy table, so the node rises with the top cell and, with the top cell in the table, stays in the table (aragog.mesh.top_node_entropy). limited is refused on core_bc = "gradient": when the config is read, on a later assignment of either key on a BoundaryConfig, and when the mesh is built, which also covers a parsed parameter object changed after loading. Any other value is refused. |
[mesh]
Spatial discretisation and pressure-density profile.
| Key | Unit | Description |
|---|---|---|
outer_radius |
m | Planet (mantle top) radius |
inner_radius |
m | CMB radius |
number_of_nodes |
-- | Number of basic-grid nodes (cell faces) |
mixing_length_profile |
str | "nearest_boundary" (distance to nearer mesh boundary) or "constant" (a fixed fraction of mantle thickness; see mixing_length_constant_fraction below) |
mixing_length_constant_fraction |
-- | Fraction of mantle thickness used as the mixing length when mixing_length_profile = "constant". Ignored otherwise. Default 0.25 |
core_density |
kg/m³ | Mean core density |
eos_method |
int | 1 = analytic AdamsāWilliamson; 2 = external file (eos_file) |
surface_density |
kg/m³ | Surface mantle density (Adams-Williamson). Default 4078.95095544 (PROTEUS production; SPIDER -adams_williamson_rhos) |
gravitational_acceleration |
m/s² | Scalar gravity for AdamsāWilliamson; superseded by the per-node profile from eos_file when eos_method = 2. Default 9.81 |
adiabatic_bulk_modulus |
Pa | \(K_S\) in AdamsāWilliamson. Default 260e9 |
adams_williamson_beta |
-- | AāW exponent \(\beta\); 0.0 derives it from \(K_S\). Default 0.0 |
surface_pressure |
Pa | Atmospheric overburden added to the pressure integration. Default 0.0 |
mass_coordinates |
bool | If true, use a SPIDER-parity mass-coordinate grid with Newton-solved spatial radii. Default true (PROTEUS production); set false for uniform spacing in radius |
eos_file |
str | Path to a four-column file (r [m], P [Pa], rho [kg/m³], g [m/s²]) used when eos_method = 2. PROTEUS supplies a Zalmoxis-generated file via this key. A relative path resolves against the directory of the config file first, then the working directory. |
[energy]
Heat-transport switches, transport parameters, and integrator selection.
| Key | Type | Default | Description |
|---|---|---|---|
conduction |
bool | -- | Conductive flux \(-k\,[(T/c_p)\partial S/\partial r + (\partial T/\partial P)_S\,\partial P/\partial r]\) |
convection |
bool | -- | Convective flux \(\rho T\kappa_h(-\partial S/\partial r)\), driven when \(\partial S/\partial r < 0\), with the floor up to convective_threshold |
gravitational_separation |
bool | -- | Stokes-regime separation of melt and solid; phase-boundary smoothed |
mixing |
bool | -- | Chemical mixing flux (SPIDER bracket form involving \(\partial S_\mathrm{liq}/\partial P\) and \(\partial S_\mathrm{sol}/\partial P\)) |
radionuclides |
bool | -- | Internal radiogenic heating from any [radionuclide_*] sections |
tidal |
bool | -- | Tidal heating from tidal_array |
eddy_diffusivity_thermal |
float | 1.0 | Scalar multiplier on \(\kappa_h\). Negative values pin \(\kappa_h\) to the absolute value (SPIDER convention) |
eddy_diffusivity_chemical |
float | 1.0 | Scalar multiplier on \(\kappa_c\). Negative values pin to absolute |
kappah_floor |
m²/s | 10.0 | Phase-modulated lower bound on \(\kappa_h\). Default 10.0 (PROTEUS production); set 0.0 for textbook MLT |
convective_threshold |
J kgā»Ā¹ Kā»Ā¹ mā»Ā¹ | 0 | Entropy gradient \(D\): the mixing-length drive acts only where \(\partial S/\partial r < 0\) and uses at least \(D\) there, and the kappah_floor applies up to \(\partial S/\partial r = D\), so nodes with \(0 \le \partial S/\partial r < D\) get the floor alone. Finite and \(\ge 0\); the default 0.0 gives the \(\partial S/\partial r < 0\) switch of SPIDER; \(D = 10^{-6}\) gives one end state in 3 replayed coupled calls that end in one of 2 states under that switch, but \(D > 0\) can stall CVODE at a node whose gradient stays near 0, also at rtol = \(10^{-8}\). See Heat transport |
phi_step_cap |
-- | 0.0 | Per-call \(\Delta\Phi_\mathrm{global}\) cap. When > 0 and the mantle straddles the rheological transition, a SUNDIALS root function fires at the step where the mass-weighted global melt fraction \(\Phi_\mathrm{global}\) has changed by cap from its value at solve() entry. 0.05 is a useful upper bound for 1 M\(_\oplus\) runs. Default 0.0 (disabled). |
eos_floor_margin |
J kgā»Ā¹ Kā»Ā¹ | 0 | EOS floor margin \(m\); the default 0 leaves the floor off, and a positive margin turns it on (100 is the tested value). With the floor off, a drain past the table minimum continues at the clamped table-edge temperature and the EOS lookup logs an out-of-range WARNING (on the 1st, 10th, 100th, ... occurrence of each lookup); it raises instead only for an EOS that a caller builds with strict_range=True. With the floor on, a root function (CVODE) or a terminal event (scipy) ends the call when the smaller of the top basic-node entropy, which sets the surface temperature, and the smallest cell entropy reaches the table minimum plus \(m\), or plus \(m/2\) when the call starts below the minimum plus \(m\); the call returns a shorter step. A call that starts at or below the minimum plus \(m/2\) arms a level at half its distance to the minimum and does not advance (status 1) when its entropy falls fast enough to reach that level before the call ends, and a caller that retries with a shorter span, such as the PROTEUS retry ladder, can then advance. A start at the table minimum (within \(2 \times 10^{-6} m\)) does not advance on any fall at any span, and a start at or below the minimum, also of the extrapolated top node alone, or so close to it that no level fits between them, does not advance at all. A call whose stop lands below the minimum through the root location error, which grows with the call time, is discarded as a refusal (see the solver tuning page). A missing key keeps the default 0; 0 or a negative value disables it (None too, outside TOML, which has no None), and a margin as large as the entropy span of the tables is refused. PROTEUS has no key for it, so PROTEUS runs have the floor off. See Solver method tuning |
phase_boundary_cap |
str | "rate" |
Step-size cap near or inside the two-phase band. "fixed" uses 1 yr. "rate", the default (also when unset), uses event-driven CVODE segments with max_step set to 0.1 of the shortest time to reach a phase boundary, clipped to [1, 100] yr. Inside the stiff zone (\(\delta = \max(3w \max_j(S_{\mathrm{liq},j}-S_{\mathrm{sol},j}), 10\ \mathrm{J\,kg^{-1}\,K^{-1}})\), with \(w = \mathtt{matprop\_smooth\_width}\)) or below the rate floor, the estimate uses the distance to the nearer boundary. The gradient core and the scipy integrators use 1 yr in both modes. On CVODE the segments of a "rate" call arm within the larger of phase_boundary_entropy_margin and \(\delta\); the step caps, "fixed", the gradient core and the scipy integrators use phase_boundary_entropy_margin. At the default rtol \(= 10^{-8}\), "rate" stays within \(1.3 \times 10^{-4}\) K of a "fixed" run at \(10^{-10}\) in the cases on the energy equation page; an rtol above \(10^{-7}\) makes the rate cap log one warning. |
bottom_up_grav_sep |
bool | true | Apply SPIDER's bottom-up gating on the gravitational-separation flux |
phase_smoothing |
str | "tanh" |
Phase-boundary smoothing. "tanh" (default) is SPIDER's two-branch get_smoothing with width matprop_smooth_width = 0.01; "cubic_hermite" is the fallback \(16\,g\phi^2(1-g\phi)^2\) form. |
solver_method |
str | "cvode" |
ODE integrator. "cvode" selects SUNDIALS CVODE via scikits.odes (default); "radau" and "bdf" use scipy solve_ivp. When scikits.odes is not installed the solver falls back to Radau with a warning. |
use_jax_jacobian |
bool | true | When solver_method = "cvode", install a JAX-traced analytic Jacobian instead of the finite-difference Jacobian. Requires JAX and the PROTEUS-side factory registration; falls back to FD Jacobian otherwise. |
tidal_array |
array[float] | [0.0] |
Per-staggered-node tidal heating in W/kg. Length must be 1 (broadcast scalar) or number_of_nodes - 1 |
[initial_condition]
Initial entropy profile (mapped from the temperature inputs via the EOS).
| Key | Unit | Description |
|---|---|---|
initial_condition |
int | 1 = linear T(r) profile (default), 2 = read from init_file, 3 = adiabatic |
surface_temperature |
K | Top temperature for modes 1 and 3 |
basal_temperature |
K | Bottom temperature for mode 1 |
init_file |
str | Path to a two-column r, T (or r, S) file used when initial_condition = 2. A relative path resolves against the directory of the config file first, then the working directory. |
When the PROTEUS wrapper drives Aragog the initial entropy is supplied directly via EntropySolver.set_initial_entropy(); the section keys are still parsed but the values are overridden by the call.
[phase_liquid] and [phase_solid]
End-member phase properties. Float values mean a constant; string values are file paths to one-column lookup tables.
| Key | Unit | Description |
|---|---|---|
density |
kg/m³ | Reference density (or path to lookup) |
viscosity |
PaĀ·s | Dynamic viscosity (or path) |
heat_capacity |
J/kg/K | Specific heat capacity |
melt_fraction |
-- | Reference melt fraction (0 for solid, 1 for liquid) |
thermal_conductivity |
W/m/K | Thermal conductivity |
thermal_expansivity |
1/K | Thermal expansivity |
entropy |
J/kg/K | Optional reference entropy |
In the production PROTEUS path the per-phase keys are not consumed for material properties; the EOS tables provide \(\rho\), \(c_p\), \(\alpha\), \(k\), and \(T\) as functions of \((P, S)\). The values are kept for the standalone constant-properties path (see const_properties below).
[phase_mixed]
Mixed-phase (mushy) parameters and the optional constant-properties analytical mode.
| Key | Unit | Default | Description |
|---|---|---|---|
latent_heat_of_fusion |
J/kg | -- | Reference latent heat (legacy; production runs use \(L(P)\) from the EOS) |
rheological_transition_melt_fraction |
-- | -- | Critical melt fraction for the viscosity transition |
rheological_transition_width |
-- | -- | Width of the tanh viscosity blend |
solidus |
str | -- | Path to a solidus reference (legacy T-based; not consumed by the entropy solver) |
liquidus |
str | -- | Path to a liquidus reference (legacy T-based) |
phase |
str | -- | Phase label |
phase_transition_width |
-- | -- | Legacy smoothing width |
grain_size |
m | 1.0e-3 |
Grain size \(a\) in the melt-solid mobility \(F(\zeta)=K(\zeta)/\zeta\) |
matprop_smooth_width |
-- | 0.0 | SPIDER-parity blending width across phase boundaries |
separation_viscosity |
str | "melt" |
Drag viscosity in \(v_\mathrm{rel} = \lvert\Delta\rho\rvert g F(\phi)/\eta\): "melt" (fixed single-phase liquid viscosity, SPIDER parity) keeps separation active below phi_rheo, which collapses the CMB temperature in coupled caps-off runs; "mixture" (rheological-transition-blended bulk viscosity) locks separation up below phi_rheo and avoids that collapse, which is why PROTEUS defaults to "mixture" |
cp_blend |
str | "latent" |
"latent" (SPIDER parity) or "linear" for the mushy \(C_p\) rule |
const_properties |
bool | false | Enable the analytic constant-properties path (no EOS tables needed) |
const_rho |
kg/m³ | 4000 | Constant density when const_properties = true |
const_Cp |
J/kg/K | 1000 | Constant heat capacity |
const_alpha |
1/K | 1e-5 | Constant thermal expansivity |
const_cond |
W/m/K | 4.0 | Constant thermal conductivity |
const_log10visc |
-- | 2.0 | Constant logāā viscosity |
const_T_ref |
K | 3500 | Reference temperature for \(T(S) = T_\mathrm{ref}\exp((S - S_\mathrm{ref})/C_p)\) |
const_S_ref |
J/kg/K | 3000 | Reference entropy in the same expression |
Loading configuration in Python
from aragog.parser import Parameters
# From a TOML or INI file (auto-detected by suffix)
params = Parameters.from_file("path/to/config.toml")
# From a dictionary (used by PROTEUS via Config.from_dict)
from aragog.config import Config
params = Config.from_dict(config_dict)
Both calls return a Parameters dataclass; the Config facade exists primarily for the dict-driven path used by the PROTEUS wrapper.
Parameters.from_file dispatches on file suffix:
| Suffix | Loader | Semantics |
|---|---|---|
.toml |
stdlib tomllib |
TOML: quoted strings ("nearest_boundary") and # ... inline comments. Strings round-trip without surrounding quotes; comments are stripped before type coercion. |
.cfg, .ini, others |
typed_configparser |
INI: bare unquoted strings (nearest_boundary). No inline comment support; values are read verbatim, quotes included. |
Use the format matching the file you are writing: a .toml file must use TOML conventions, an .ini / .cfg file must use INI conventions. Mixing the two (e.g. quoted strings in a .cfg) is a common cause of Mixing length profile is unknown errors.
Legacy INI format
The legacy .cfg format is still parsed. Example files live under src/aragog/cfg/. Keys and section names match the TOML schema; only the surface syntax differs (no quotes around strings, no # inline comments).