Skip to content

ADR-0001 — Operator-split stiff chemistry (LSODA + Bulirsch–Stoer), not folded into the advection Newton matrix

Status: Accepted

Context

Surface and gas-phase chemistry (TWC/SCR/DOC Brinkmeier kinetics, etc.) is stiff: reaction timescales span many orders of magnitude and are far shorter than the acoustic/advection timestep. There were two structural choices for how chemistry interacts with the compressible-flow solve:

  1. Fold the reaction source into the advection Newton matrix so one implicit solve advances flow + chemistry together.
  2. Operator-split chemistry from transport (Lie/Strang split) and integrate the per-cell reaction ODE with a dedicated stiff integrator.

Option 1 is attractive on paper (one coupled solve, no splitting error) but in practice it makes the Jacobian dense and badly conditioned: the stiff reaction block destroys the clean block-tridiagonal structure the transport solve relies on (see ADR-0002), and the linear solve becomes the bottleneck. Multiple QS-chemistry integrator candidates were compared on the stiff Brinkmeier cases.

Decision

Chemistry stays operator-split from transport. The per-cell stiff ODE is integrated by:

  • scipy.integrate.solve_ivp(method="LSODA") on the explicit / quasi-steady path (LSODA auto-switches between non-stiff and stiff BDF internally) — it won on the stiff QS chemistry across the integrators compared; and
  • a semi-implicit Bulirsch–Stoer (Bader–Deuflhard) @njit kernel (src/aeolus1d/solver/chem_bs_jit.py, after Bader & Deuflhard 1983 / NR3 §17.5) for the implicit-Euler / BDF2 split, where calling back into SciPy per cell is too costly.

The reaction source is never assembled into the transport Newton matrix.

Consequences

  • The transport Jacobian keeps its block-tridiagonal shape and O(N) direct solve.
  • A first-order (Lie) or second-order (Strang) splitting error is accepted; it is controlled by the macro-step size and the chem_every_n super-cycle knob, and is negligible when the chemistry timescale comfortably exceeds the macro step.
  • LSODA is a hard dependency of the stiff path; the learned-closure / linear-op accelerators (chem_linop) must keep an LSODA fallback for out-of-distribution states (src/aeolus1d/closures/chem_linop.py).

Alternatives considered

  • Fold chemistry into the advection Newton matrix — rejected: dense, stiff, ill-conditioned Jacobian; defeats the block-tridiagonal direct solve.
  • A single fixed-step implicit integrator for all chemistry — rejected for stiff cases (accuracy/robustness); a fixed-step backward-Euler path is shipped only for explicitly non-stiff chemistry.

References

  • src/aeolus1d/solver/chem_bs_jit.py — the Bader–Deuflhard kernel.
  • src/aeolus1d/closures/chem_linop.py — the LSODA-fallback contract.
  • Related: ADR-0002.