Writing and inspecting NetCDF output

Aragog writes its own NetCDF4 snapshots in standalone use, and PROTEUS-coupled runs additionally produce a per-iteration interior snapshot via the framework wrapper. The two paths are independent on purpose: the standalone schema captures the full SolverOutput dataclass (every scalar diagnostic and every profile), while the PROTEUS schema is a curated subset of that, written under a shorter variable-name convention that the PROTEUS post-processors and read_last_Sfield resume helper expect.

This page covers both: how to write a standalone snapshot, and how to inspect either.

Standalone NetCDF write

SolverOutput.to_netcdf writes a self-contained file capturing every field on the dataclass. For convenience EntropySolver.write_netcdf calls self.get_state().to_netcdf(...) so a typical script reads:

from pathlib import Path
from aragog.solver import EntropySolver

solver = EntropySolver.from_file(
    filename="src/aragog/cfg/abe_solid.toml",
    eos_dir="/path/to/eos/tables",
)
solver.initialize()
solver.set_initial_dSdr_cmb(0.0)        # only when core_bc='energy_balance'
solver.set_initial_entropy(2900.0)
solver.solve()

# Single call: dumps the final state to a NetCDF4 file.
solver.write_netcdf(Path("output/snapshot.nc"), time=0.0,
                    description="Final state of abe_solid smoke run")

The writer mkdirs any missing parent directories, opens the file in mode='w' (overwriting any pre-existing file), and stamps the dataset with aragog_version, an ISO-8601 UTC created_utc, the user-supplied description, and the CF-1.8 conventions tag. Every variable carries a units and a long_name attribute so the file is interpretable without consulting the source.

To dump an intermediate state (rather than the final), call to_netcdf directly on a SolverOutput:

out = solver.get_state()
out.to_netcdf("output/snapshot.nc", time=t_yr)

When the time argument is omitted the writer falls back to dt_actual (the per-call integration duration); pass an explicit time whenever you have an absolute simulation clock.

Variables in the standalone schema

Variables are named verbatim after the SolverOutput field names, with their units in SI by default.

Staggered-node profiles (length \(N\))

Variable	Source field	Unit	Description
`r_stag`	`r_stag`	m	Staggered-node radii
`P_stag`	`P_stag`	Pa	Pressure at staggered nodes
`S_final`	`S_final`	J kg⁻¹ K⁻¹	Specific entropy
`T_stag`	`T_stag`	K	Temperature
`phi_stag`	`phi_stag`	--	Melt mass fraction
`rho_stag`	`rho_stag`	kg m⁻³	Density
`visc_stag`	`visc_stag`	Pa s	Dynamic viscosity
`vol`	`vol`	m³	Per-shell volume
`mass_stag`	`mass_stag`	kg	Per-shell mass
`heating`	`heating`	W kg⁻¹	Internal heating (radiogenic + tidal)
`cap_stag`	`cap_stag`	kg K m⁻³	Capacitance \(\rho T\)

Basic-node profiles (length \(N+1\))

Variable	Source field	Unit	Description
`r_basic`	`r_basic`	m	Basic-node radii
`heat_flux`	`heat_flux`	W m⁻²	Total radial heat flux
`eddy_diff`	`eddy_diff`	m² s⁻¹	Thermal eddy diffusivity (\(\kappa_h\))
`jcond_b`	`jcond_b`	W m⁻²	Conductive flux
`jconv_b`	`jconv_b`	W m⁻²	Convective (MLT) flux
`jgrav_b`	`jgrav_b`	W m⁻²	Gravitational-separation flux
`jmix_b`	`jmix_b`	W m⁻²	SPIDER-parity phase-mixing flux
`dSdr_b`	`dSdr_b`	J kg⁻¹ K⁻¹ m⁻¹	Radial entropy gradient
`phi_basic`	`phi_basic`	--	Melt fraction at basic nodes
`T_basic`	`T_basic`	K	Temperature at basic nodes
`cp_basic`	`cp_basic`	J kg⁻¹ K⁻¹	Heat capacity at basic nodes
`rho_basic`	`rho_basic`	kg m⁻³	Density at basic nodes
`porosity_b`	`porosity_b`	--	Density-derived porosity at basic nodes
`rho_solid_b`	`rho_solid_b`	kg m⁻³	Solid density at the phase boundary, at the node pressure
`rho_melt_b`	`rho_melt_b`	kg m⁻³	Melt density at the phase boundary, at the node pressure
`g_b`	`g_b`	m s⁻²	Gravitational acceleration at basic nodes: the magnitude of the mesh gravity, or an external `eos_gravity` profile as given

The last four rows are for consumers that compute melt drainage. porosity_b is the porosity the gravitational-separation velocity uses, \(\phi = (\rho_s - \rho)/(\rho_s - \rho_l)\), smoothly clipped to \([0, 1]\) with a width of \(10^{-3}\). Its range is \([-2.5 \times 10^{-7}, 1]\) and it is \(5.0 \times 10^{-4}\) where \(\rho = \rho_s\). It can be negative, so floor it before a logarithm or a power. It is not the melt fraction: in a solid node it is near zero where the solid density is at least \(\rho_s\), and tables whose solid density rises with entropy towards the solidus give larger values there, up to 1. Use phi_basic to decide whether a node holds melt. Where rho_solid_b <= rho_melt_b a mush is denser than rho_solid_b, so porosity_b is about 0 although the node holds melt: do not use it for melt drainage there, and test the sign of the contrast first. Phi_global_vol counts that melt, because it uses the unfloored ratio, hard-clipped, at the staggered nodes, with \(\rho\) built from the melt fraction. The permeability over porosity that goes with it is aragog.eos.mobility_function(porosity, grain_size), which floors it internally. The denominator of \(\phi\) is floored at about 1 kg m⁻³, which acts where rho_solid_b - rho_melt_b is below that, including where the tables give a melt denser than the solid (the difference is negative). There porosity_b is \((\rho_s - \rho)/(1\ \mathrm{kg\,m^{-3}})\) clipped to \([0, 1]\): about 0 where rho_basic is at or above rho_solid_b, 1 where it is more than 1 kg m⁻³ below, and intermediate in between. With const_properties there is no phase contrast: porosity_b is 1 and both densities equal the constant density.

Scalars

Time, surface and CMB temperatures, both averaging conventions of the global melt fraction, the energy-balance integrals, the solver status code, the raw CVODE return flag, and the per-solve core-temperature excursion:

time, T_magma, T_core, Phi_global, Phi_global_vol, M_mantle, M_mantle_liquid, M_mantle_solid, RF_depth, E_th, E_state, E_state_cons, Cp_eff, F_heat_total, F_cmb, Q_radio_total, Q_tidal_total, step_dE_F_int_J, step_dE_F_cmb_J, step_dE_core_J, step_dE_Q_radio_J, step_dE_Q_tidal_J, step_dE_Q_radio_cons_J, step_dE_Q_tidal_cons_J, step_solver_residual_J, step_dE_compression_J, step_dE_state_heat_J, dt_actual, status (i4), cvode_flag (i4), tcore_change_max (f8), tcore_change_exceeded (i4).

step_dE_core_J is the core's heat change over the call for core_bc = "core_module" (the effective-capacity integral over the T_core path) and "bower2018" (the reservoir capacity times the T_core change); the other modes write 0.0, since their core is not integrated as a heat content.

status is the normalised solve outcome (0 success, negative failure). cvode_flag is the raw integer flag CVODE returned for the same solve, kept distinct from status so that a CV_TOO_MUCH_WORK (-1) step-budget stop is told apart from a CV_CONV_FAILURE (-4) convergence failure; the scipy solve_ivp path reports 0. Its name is written as the string global attribute cvode_flag_name. tcore_change_max is the largest change of the core temperature from the solve-entry value over the returned grid [K]; tcore_change_exceeded is 1 when tcore_change_limit is set and that change exceeds it, or when any sampled core temperature is non-finite (a corrupted-solve signal, independent of the limit), else 0.

Inspecting a snapshot

Either via the netCDF4 CLI, the Python netCDF4 library, or xarray:

ncdump -h output/snapshot.nc                     # header + variable list
ncdump output/snapshot.nc | less                 # full contents

import xarray as xr

ds = xr.open_dataset("output/snapshot.nc")
print(ds)                                          # one-line summary
print(ds["Phi_global"].item())                     # scalar access
print(ds["S_final"].values.shape, ds["r_stag"].values.shape)

ncdump is bundled with most netCDF4 distributions; on macOS brew install netcdf provides it, on Linux apt install netcdf-bin or dnf install netcdf. xarray is the recommended tool for analysis-style work because it keeps coordinates and units alongside the data.

Plotting from a snapshot

A minimal plot of temperature, melt fraction, and total heat flux versus pressure:

import matplotlib.pyplot as plt
import xarray as xr

ds = xr.open_dataset("output/snapshot.nc")

fig, axes = plt.subplots(1, 3, figsize=(11, 4), sharey=True)
P_GPa = ds["P_stag"].values / 1e9

axes[0].plot(ds["T_stag"], P_GPa)
axes[0].set_xlabel("Temperature [K]")
axes[0].set_ylabel("Pressure [GPa]")
axes[0].invert_yaxis()

axes[1].plot(ds["phi_stag"], P_GPa)
axes[1].set_xlabel("Melt fraction")
axes[1].invert_yaxis()

axes[2].plot(ds["heat_flux"], ds["r_basic"].values / 1e3)
axes[2].set_xlabel(r"Total heat flux [W m$^{-2}$]")
axes[2].set_ylabel("Radius [km]")
axes[2].set_xscale("symlog", linthresh=1e-3)

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

The per-component flux decomposition lives in the same file as jcond_b, jconv_b, jgrav_b, jmix_b; overlay them on the heat-flux panel for the full Figure 2 of Heat transport.

PROTEUS-coupled snapshots

When Aragog is driven by PROTEUS, the framework wrapper (proteus.interior_energetics.aragog.AragogRunner._write_output_ncdf) writes a separate, more compact per-iteration snapshot at output/data/<iter>_int.nc (where <iter> is the integer time in years), using shorter *_s and *_b variable names that the resume helper recognises. This schema is independent of SolverOutput.to_netcdf so that PROTEUS post-processors and the resume path do not break when the standalone schema evolves.

PROTEUS variable name convention

Always written:

Variable	Source field	Unit	Description
`entropy_s`	`S_final`	J kg⁻¹ K⁻¹	Specific entropy
`temp_s`	`T_stag`	K	Temperature derived from \((P, S)\)
`phi_s`	`phi_stag`	--	Melt fraction
`density_s`	`rho_stag`	kg m⁻³	Density
`log10visc_s`	\(\log_{10}(\texttt{visc\_stag})\)	\(\log_{10}\) Pa s	Base-10 dynamic viscosity
`pres_s`	`P_stag / 1e9`	GPa	Pressure (note: GPa, not Pa)
`radius_s`	`r_stag / 1e3`	km	Staggered-node radii (km)
`mass_s`	`mass_stag`	kg	Cell mass
`Htotal_s`	`heating`	W kg⁻¹	Per-node internal heating
`radius_b`	`r_basic / 1e3`	km	Basic-node radii (km)
`Ftotal_b`	`heat_flux`	W m⁻²	Total radial heat flux
`time`	--	yr	Snapshot time
`phi_global`	`Phi_global`	--	Mass-weighted mean melt fraction

Written only when interior_energetics.write_flux_diagnostics = true:

Variable	Source field	Unit	Description
`Jcond_b`	`jcond_b`	W m⁻²	Conductive flux
`Jconv_b`	`jconv_b`	W m⁻²	Convective (MLT) flux
`Jgrav_b`	`jgrav_b`	W m⁻²	Gravitational-separation flux
`Jmix_b`	`jmix_b`	W m⁻²	Chemical-mixing flux
`dSdr_b`	`dSdr_b`	J kg⁻¹ K⁻¹ m⁻¹	Entropy gradient
`eddy_diff_b`	`eddy_diff`	m² s⁻¹	Thermal eddy diffusivity (\(\kappa_h\))
`phi_basic_b`	`phi_basic`	--	Melt fraction at basic nodes
`T_basic_b`	`T_basic`	K	Temperature at basic nodes
`cp_basic_b`	`cp_basic`	J kg⁻¹ K⁻¹	Heat capacity at basic nodes
`rho_basic_b`	`rho_basic`	kg m⁻³	Density at basic nodes

Per-iteration scalar diagnostics (\(T_\mathrm{magma}\), \(T_\mathrm{core}\), \(\Phi_\mathrm{global}\), \(E_\mathrm{th}\), \(C_p^\mathrm{eff}\), \(F_\mathrm{heat}^\mathrm{total}\), the four step_dE_*_J integrals) live in output/runtime_helpfile.csv rather than the snapshot, since PROTEUS aggregates them across all interior steps.

Resuming from a PROTEUS snapshot

PROTEUS itself uses proteus.interior_energetics.aragog.read_last_Sfield to re-load the previous step's entropy field at the start of the next iteration. The same helper is the canonical way to resume a run from a saved snapshot:

from proteus.interior_energetics.aragog import read_last_Sfield

S_prev = read_last_Sfield(output_dir="output", time=t_prev)
solver.set_initial_entropy(S_prev)

The helper reads entropy_s (PROTEUS schema) from the snapshot, validates the array length against the live mesh, and returns the entropy in J/kg/K. For standalone snapshots written via to_netcdf, the equivalent field is S_final and you can re-load it with one line of xarray or netCDF4.