Running neoVULCAN¶

There are two ways to drive the code: the standalone command-line entry point vulcan.py, and the library API in vulcan_api.py for embedding the chemistry in a larger atmospheric model.

Command-line use¶

The script vulcan.py performs a complete simulation in the following stages:

Parse command-line flags (-c, -n), chdir into the vulcan.py directory, and prepend src/ to sys.path.
Load the TOML configuration through neovulcan_config.VulcanConfig.from_toml and install it as the process-wide singleton via neovulcan_runtime.set_cfg().
Regenerate src/chemistry_jax.py from the network file referenced in [network] unless -n was passed.
Build the atmosphere (build_atm.Atm): pressure grid, temperature and \(\Kzz\) profile, mean molecular weight, layer thickness, molecular diffusion coefficients, boundary-condition fluxes.
Parse the reaction network and evaluate all thermal rate coefficients on the grid (rates.ReadRate).
Set the initial number densities (build_atm.InitialAbun) from FastChem equilibrium, a constant mixing ratio, a saved .vul file, or a tabulated initial condition.
Instantiate the chosen ODE solver (Ros2 or Rodas3).
If photochemistry is enabled, bin the stellar flux, load the photolysis cross-sections, construct the radiative-transfer object selected by photochemistry.rt_scheme (two-stream or DisORT++), and call it once.
Run the integration loop (integration.Integration) until either the steady-state criteria are met or one of the hard limits (count_max, runtime) is reached.
Pickle the final state to <output_dir>/<out_name>.

The driver itself is a thin orchestration layer; almost all of the work is done inside the modules in src/ documented in Code architecture.

Command-line options¶

Flag	Purpose
`-c PATH`,	Path to the TOML configuration file. Default
`--config PATH`	`vulcan_cfg.toml` next to `vulcan.py`.
`-n`,	Skip the regeneration of
`--no-remake-chemistry`	`src/chemistry_jax.py`. Use this when the network file has not changed since the last run; it noticeably reduces start-up time.

Library API¶

For coupling neoVULCAN to a three-dimensional general-circulation model (GCM) or another driver, use vulcan_api.VulcanChemistry. A typical GCM time step looks like

import sys
sys.path.insert(0, '/path/to/VULCAN')
from neoVULCAN.vulcan_api import VulcanChemistry

BASE = '/path/to/neoVULCAN'
chem = VulcanChemistry(
    BASE,
    config_path=f'{BASE}/cfg_examples/HD189.cfg',
    cfg_overrides={
        'photochemistry': {'use_photo': True, 'rt_scheme': 'disort'},
        'solver':         {'rtol': 0.5},
    },
)
chem.initialize(regenerate_chemistry=True)

for step in range(n_steps):
    T_new, P_new, Kzz_new = gcm.get_profiles(col)
    chem.set_atmosphere(T=T_new, P=P_new, Kzz=Kzz_new)
    chem.run_to_convergence()

    ymix = chem.get_mixing_ratios()         # shape (nz, ni)
    info = chem.get_convergence_info()
    gcm.update_chemistry(col, ymix, chem.species)

cfg_overrides is a nested dict that follows the TOML structure exactly: top-level keys are section names (solver, atmosphere, …), each holding a sub-dict of field overrides. It is deep-merged on top of the loaded TOML and then re-validated by Pydantic, so the same strictness applies to programmatic overrides as to file values.

Because src/chemistry_jax.py bakes in the species list and the number of layers atmosphere.nz, only one VulcanChemistry instance with a given network and grid can exist in the same Python process at a time.

The configuration singleton¶

Internally, both entry points install the loaded neovulcan_config.VulcanConfig into a process-wide singleton in neovulcan_runtime. Every module under src/ reads its parameters through

from neovulcan_runtime import get_cfg
cfg = get_cfg()
nz   = cfg.atmosphere.nz
rtol = cfg.solver.rtol

so the parameters are typed (Pydantic objects), tab-completable, and shared between modules without import-order surprises.

Stand-alone helper scripts under plot_py/, atm/ and tools/ that are not launched through vulcan.py use neovulcan_runtime.get_cfg_or_load(), which lazily loads vulcan_cfg.toml (or a path you pass) the first time it is called.

Choosing the integrator¶

Two integrators are exposed via the solver.ode_solver parameter:

Ros2: Second-order Rosenbrock W-method. A-stable, two stages, two right-hand-side evaluations per step, one LU factorisation re-used across both stages. Robust default. See Numerical methods.
Rodas3: Third-order Rosenbrock–Wanner method. L-stable and stiffly accurate, with an embedded second-order error estimate. Four stages, three RHS evaluations, four banded back-substitutions. Costs roughly twice as much per step as Ros2 but often takes fewer steps to converge. Requires atmosphere.use_moldiff = true (the Pydantic validator rejects the combination otherwise).

Empirically, Ros2 is preferred for the canonical hot-Jupiter and terrestrial set-ups; Rodas3 becomes attractive when very tight convergence is required or when the chemistry is unusually stiff.

An optional Newton finisher can be enabled with solver.use_newton_finisher = true; it switches from Rosenbrock to a short damped-Newton tail once the per-step change drops below solver.newton_switch_dy to polish the residual without growing the step count. See Configuration reference for the associated tuning knobs.

Choosing the radiative-transfer scheme¶

The photochemistry.rt_scheme parameter selects between two RT backends:

"two-stream" (default): Delta-Eddington two-stream solver implemented in pure NumPy (radiative_transfer.TwoStreamRT). Fast, robust, and good enough for the vast majority of exoplanet runs. The Eddington coefficient is set by photochemistry.edd (default 0.5).
"disort": Per-bin discrete-ordinates solver based on the DisORT++ Python bindings (radiative_transfer.DisortRT). The number of streams is set by photochemistry.disort_nstr (default 8); photochemistry.surface_albedo is honoured for the lower boundary. Requires the disortpp package; the import is deferred to DisortRT.__init__ so two-stream runs work without it.

Both schemes share the optical-depth assembly and the J / J_ion spectral integrals — only the flux step is replaced — so output arrays have the same shape and the rest of the pipeline (convergence checks, plotting) is unchanged. See Mathematical background for the physics and Numerical methods for the implementation.

Performance notes¶

The JAX kernels (chemistry_jax.py and jacobian_jax.py) are jit-compiled on the first call. The first time step of a run therefore pays a one-off cost of a few seconds; subsequent steps are fast.
The Jacobian is stored in LAPACK banded format and factorised once per time step with dgbtrf / dgbtrs, which is markedly faster than the general scipy.linalg.solve_banded.
Radiative transfer is the next biggest cost after the linear solves. The update frequency switches automatically from photochemistry.ini_update_photo_frq to photochemistry.final_update_photo_frq once the chemistry is close to steady state. The DisORT++ backend is several times more expensive per call than the two-stream solver but parallelises across wavelength bins via OpenMP inside disortpp.
atmosphere.update_frq controls how often the layer thickness dz is recomputed from the mean molecular weight; for nearly hydrostatic configurations this can be increased without affecting the solution.