Numerics — implicit Newton, preconditioners, and the block-tridiag shape

This document records why each Newton-solve path in the implicit network solver is the shape it is. The point isn't to re-derive the math (the references at the end do that better); it's to capture which design decisions came from theory and which came from benchmarking, so the next person reading _newton_step understands why each piece is shaped the way it is.

For the practical what to set — knob defaults, TOML examples, cost-shape table, decision flowchart — see implicit-solver-reference.md. For the audit-driven evidence trail behind the bench numbers, see perf-history.md.

Choosing a time integrator (the [case] rk knob)

aeolus1d ships seven time integrators behind the single rk knob — four explicit Runge-Kutta schemes and three implicit ones. They trade three things against each other: temporal accuracy (formal order), stability (does it unlock CFL > 1, and is it L-stable for stiff terms), and numerical damping (how much spurious high-frequency content it dissipates). GT-POWER frames the same trade in Mechanical §1.7.8–9: a controllable-dissipation spectrum from "fully damped" (implicit-Euler-like) to "sharp" (higher-order, near-energy-conserving). aeolus1d's discrete dial:

rk	Order	Class	CFL	Numerical damping	Pick it for
ssp_rk2	2	Explicit, TVD	≤ 1	only the limiter's	Default. Shocked / discontinuous flow, general use — TVD-safe
ssp_rk3	3	Explicit, TVD	≤ 1	only the limiter's	Smooth flow wanting lower temporal error while staying TVD
rk4	4	Explicit, not TVD	≤ 1	none (anti-dissipative near shocks)	Smooth benchmarks where the temporal-error floor dominates
rk5	5	Explicit, not TVD	≤ 1	none	Smooth, temporal-error-limited benchmarks only
implicit_euler	1	Implicit, L-stable	> 1	maximal — damps HF content	Robust startup, very stiff, low-Mach where you want acoustics damped
bdf2	2	Implicit, A-stable	> 1	light	General implicit default; ~RK2 accuracy at CFL > 1
radau_iia	5 — but 1 on any junction network, see below	Implicit, L-stable + stiffly accurate	> 1	minimal (sharp)	Stiff chemistry-coupled single-pipe runs; a true error controller

How to read it:

• Shocks → stay explicit and TVD. rk4/rk5 are not TVD; they will oscillate at discontinuities (a Sod tube run with rk4 rings). Use them only on provably smooth flow where the temporal-error floor — not the spatial limiter — is the bottleneck. For everything with a shock, contact, or sharp gradient, ssp_rk2 (default) or ssp_rk3 is the right call.
• CFL > 1 → go implicit. All three implicit schemes unlock CFL > 1; the explicit ones are bounded by |u| + a (run aeolus1d analyze-cfl case.toml to see the bound and which pipe pins it). Pick the implicit scheme by the damping you want.
• The damping dial. implicit_euler is the fully-damped end — first-order and L-stable, it scrubs the spurious high-frequency content that a low-Mach or stiff run otherwise rings with, at the cost of accuracy. radau_iia is the sharp end — fifth-order, L-stable and stiffly accurate, it keeps resolved content crisp while still handling stiffness. bdf2 sits in between: second-order and A-stable (lightly damped), the pragmatic implicit default.
• Stiff source terms (chemistry, Forchheimer) → L-stable. Prefer implicit_euler or radau_iia over bdf2 when the stiffness is in the source terms — A-stability alone (bdf2) can let stiff modes ring; L-stability damps them dead.
• Adaptive error control ([case] adaptive_error_control) uses step-doubling on one-step implicit stencils — it pairs with implicit_euler and radau_iia but raises on the multi-step bdf2 (which carries history and can't be Richardson-doubled cleanly). If you want a true accuracy controller, pick radau_iia.

Implicit schemes on a junction network

The Order column above is the scheme's order on a single pipe. On any case that has a junction, every implicit scheme is first-order.

The network coupling is a block-split: each junction's volume is held frozen through the pipe Newton solves and then advanced by forward Euler, and junction ends are Dirichlet. So all three Radau stages see one ghost built from V(t_n), and the split — not the scheme — sets the observed order.

Measured on the two-pipe-plenum network at picard_max_iter = 12:

rk	formal order	observed on a network	interior delta at cfl = 1
implicit_euler	1	1.00	6.50e1
bdf2	2	1.87–1.94 (not capped)	1.51e1
radau_iia	5	0.94–1.00	6.28e1

bdf2's row was 0.98–1.10 until 2026-08-25, and it did not reproduce

That range was prose whose harness was never written down. Promoted to an executable refinement sweep (tests/integration/test_implicit_network_split_order.py — closed two-pipe-plenum, wall far-BCs, N = 20, t_end = 2e-4, same Picard setting), bdf2 rises toward its formal order as the window tightens while radau_iia stays pinned at 1 on the same windows:

refinement window	bdf2	radau_iia
(2, 4, 8)	1.671	—
(4, 8, 16)	1.867	0.975
(8, 16, 32)	1.941	1.001
(4, 8, 16, 32)	1.905	—

Rising under refinement is what a second-order method does; a capped scheme cannot, and radau demonstrably does not. The radau_iia claim — the one this section is about — is unaffected.

The same schemes reach their formal order on a single pipe — radau_iia measures 4.86–4.99, including through advance_network_implicit with junctions=[] — and the tableau itself is verified against Hairer & Wanner II.7 (B(5), C(3), D(2) to ~1e-17, R(∞) = 1.1e-16). Switching the volume update to implicit Euler changes nothing; the deltas are bit-identical.

Practical consequence. On a network, radau_iia delivers implicit-Euler accuracy at roughly 3× the per-step cost. Prefer bdf2 there — it is not capped (observed order ~1.9 against a formal 2) and carries a ~4× smaller error constant at 1× the cost. The advice is unchanged from before 2026-08-25; the reason for it is. It used to read "same order 1, but a ~4× smaller error constant", which was right advice for the wrong reason. validate_case emits this as a warning when it sees radau_iia on a case with junctions.

Radau conserves on a network since 2026-10-08; the order is unchanged

Until then a Radau pipe lost dt·Σ bᵢ·Fᵢ through a junction face (the three stage fluxes) while the junction gained dt·F₃ (the end state), so a closed network created mass and energy at O(dt): 2.6e-5 and 5.1e-5 of the total on the two-pipe-plenum box, 1.1e-4 with a DPF wall exchange in the Newton residual. The junction now evaluates its component at each stored stage state, weighted by the Radau bᵢ, so the same quadrature is applied on both sides of the face: 2.6e-8 and 3.7e-8, the implicit-Euler floor. What is left is the Picard ghost lag. The observed order on a network is still 1, because the junction volume is still held frozen through the stage solve, so the advice above is unchanged.

Lifting the cap needs stage-resolved junction state — fill port ghosts from V_i per stage and close V_i = V^n + dt · Σ_j A_ij G(Y_j, V_j), either as an arrowhead border on the block-tridiagonal solve or by iterating the stage-resolved volume in the existing Picard loop. Tracked in BACKLOG.md.

Practical knob defaults and TOML examples live in implicit-solver-reference.md; this table is the which integrator layer above those.

Verifying the order of accuracy

The Order column above is not a claim on faith — each entry is pinned by a convergence test, and the three discretisation paths (flux, time integrator, source term) are verified separately:

Path	What is measured	Test
Spatial (flux)	mesh refinement at fixed CFL on a smooth density wave → observed spatial order (MUSCL+van Leer ≥ 1.8)	tests/integration/test_smooth_advection.py, test_mp5.py
Temporal (each rk)	Δt refinement on a fixed mesh (self-convergence) → observed temporal order, one row per integrator	tests/integration/test_temporal_order.py
Source term	manufactured-solution truncation error → formal order of a source closure's discretisation	tests/integration/test_mms_source.py

Three things worth knowing when reading those orders:

• Spatial order is a lower bound away from smooth extrema. A TVD limiter must clip the reconstruction at a smooth extremum, so the local order drops toward 1 there even for a formally 2nd-order scheme — minmod clips hardest, van_leer least (hence the ≥ 1.8 gate uses van Leer), and MP5 relaxes monotonicity near extrema. The quoted orders are the formal orders on smooth regions, not a per-point guarantee in a limited run.
• "Radau IIA(3)" is stage-count notation — 3 stages, formal order 2s − 1 = 5, not order 3. The temporal test confirms the 5.
• The temporal test isolates the time error by comparing each Δt to a fine-Δt reference on the same mesh, so the (identical) spatial truncation cancels. The implicit drivers need a tight Newton tolerance there: the default simplified-Newton — a chord-aged Jacobian (chord_max_age = 2) held across inner iterations, terminated at a finite residual tolerance — leaves a ~1e-6 fixed-point bias that is fine for engineering runs but swamps Radau's minute order-5 time error, a property worth remembering when chasing high temporal accuracy implicitly. (The HLLC flux Jacobian itself is analytic on both the uniform-γ and the per-cell-γ ideal-gas paths — Phase B13's _hllc_flux_jac_analytic_impl and the Phase B15 ..._var_gamma_impl variants — and only falls back to a finite-difference Jacobian on the real-gas / EOS-consistent frozen-MUSCL assembler, so the residual bias is a property of the chord/simplified iteration, not of the Jacobian entries.)

The implicit Newton system

Each implicit time step solves

F(U) = U − U_hist − Δt · R(U) = 0

for the new state U (U_hist is the BDF2 / Implicit-Euler history term and R(U) = -∇·F_flux + S_closure is the spatial residual). Newton's method linearises around the current iterate

J · δ = −F,    J = ∂F/∂U = I − Δt · ∂R/∂U

and updates U ← U + α·δ with a backtracking line search. For the canonical first-order scheme the spatial residual depends on each cell i only through cells i−1, i, i+1 (one HLLC face on each side), so

∂R/∂U[i, j] = 0        for |i − j| > 1

and J is block-tridiagonal with 3×3 blocks. Three arrays of shape (N, 3, 3) carry it through the code:

• L[i] — the lower block (∂R[i] / ∂U[i−1])
• D[i] — the diagonal block (∂R[i] / ∂U[i])
• Up[i] — the upper block (∂R[i] / ∂U[i+1])

After the chord-aged broadcast (Phase Perf-2-quick) the implicit Jacobian becomes L_J = -dt·L, D_J = I - dt·D, Up_J = -dt·Up. The direct path then solves J · δ = -F via block-Thomas (_block_thomas_solve_3x3_jit).

Why the direct block-tridiag solve is the right shape

The direct path solves J · δ = -F with a block-Thomas sweep that exploits the block-tridiagonal structure exactly — no fill-in, O(N) work, and it captures the off-diagonal acoustic + convective coupling (the L / Up blocks) that dominates J at cfl > 1. For the ≤ a few-thousand-cells × 3-unknowns systems aeolus1d solves, that is the correct cost shape and there is nothing cheaper.

A matrix-free Newton-Krylov (JFNK / GMRES with FD J·v) path was prototyped and benchmarked against the direct solve, with both block-Jacobi and ILU(0) preconditioners. It lost at every N tested (~1.8–1.9× slower with ILU(0)/pure; block-Jacobi was catastrophically worse — hundreds of × — because dropping the off-diagonal coupling makes the preconditioned operator less diagonally dominant at large Δt, not more). The lesson that generalises: for a hyperbolic problem at Δt-CFL > 1, a preconditioner that ignores the convective + acoustic off-diagonal coupling is at best useless and often actively harmful, so the only competitive matrix-free preconditioner (ILU(0)) collapses back to the very block-Thomas factorisation the direct path already does. The JFNK path was therefore removed in v0.9 (see the changelog); the direct + chord solve below is the single implicit Newton path.

The chord-Jacobian heuristic (Phase Perf-3)

The full Newton method rebuilds J every iteration. For most steady-state-ish steps the convergence rate after the first Newton iter is already linear with a good constant — the expensive Jacobian assembly is rebuilding almost the same matrix. The chord (lazy-Jacobian) variant holds J fixed across multiple inner iterations and accepts linear-with-good- constant convergence in exchange for fewer assembly calls.

Two NewtonConfig knobs control the trade-off:

• chord_max_age (default 2) — hard cap on iterations between rebuilds. 0 reproduces the legacy full-Newton path bit-exactly.
• chord_rate_threshold (default 0.5) — if F_norm_new / F_norm_old > threshold the chord has stopped matching the local Newton direction; force a rebuild on the next iteration. Line-search rejection on a chord-aged iteration also forces a rebuild.

The cache holds the post-broadcast (L_J, D_J, Up_J) triple, so chord-reused iterations skip both the Jacobian assembly and the closure-Jac augmentation (which Phase Perf-2 made into a JIT'd kernel but still has non-zero cost). In the canonical case 54, the chord saves ~50 % of Jacobian-assembly calls without changing the converged solution.

When each Newton-solve path wins

There is a single implicit Newton-solve path — direct + chord. It was confirmed as the right shape empirically over many bench iterations; the matrix-free Newton-Krylov variants that were prototyped alongside it lost at every N tested and were removed in v0.9 (see the changelog). The cost-shape table is in implicit-solver-reference.md.

What's not yet measured

• Very large N (≫ 1000 cells per pipe). The direct path's block-Thomas solve scales O(N), which is the right shape at the problem sizes aeolus1d targets; no real workload has pushed N far enough to make the per-step Jacobian assembly the bottleneck.
• Stiff non-hyperbolic regimes. When friction or porous-Forchheimer source terms dominate, the diagonal block D becomes more diagonally-dominant relative to the off-diagonals — a regime where the implicit step is easiest, not hardest.

A 2-D / 3-D extension would reopen the preconditioner question: line preconditioning (block-Thomas along each characteristic direction) no longer collapses to a 1-D sweep, so a proper line-Jacobi or red-black-line solver would become the right shape. The current 1-D codebase doesn't need that; the block-Thomas design is forward-compatible.

Quasi-steady DAE robustness — the limit-point / concentration-jump problem

The low-Mach aftertreatment solver (solver = "quasi_steady") drops the gas-phase accumulation terms, turning the channel PDEs into a system of differential-algebraic equations (DAEs) in which the surface coverages and solid temperature are the only states carrying explicit time derivatives. Millo et al. [Ch. 11 §11.2.6, ref. [41] in theory.md] document a generic failure mode of this formulation that aeolus1d's quasi-steady kernel has had to engineer around: at a limit point — where the steady surface-balance Jacobian becomes singular — the surface concentration jumps discontinuously (their Fig. 11.8 shows CO surface concentration stepping from 0.25 % to 0.04 % across one channel location). A naive DAE/Picard solve drives its step size toward zero as it approaches such a jump and often fails the error test there. The literature remedy is to monitor the Jacobian and apply special logic at the limit point rather than let the adaptive controller collapse. aeolus1d's SIMPLE/Picard robustness guards (the zero-response dead-zone bracket+bisect bootstrap, the strong-exotherm guard, and the junction-Picard zero-gradient bootstrap) are concrete instances of exactly this class of mitigation. A related guard hardens the species march: near such a jump the Arrhenius source can transiently overflow, and the convex-combination update is bounded only while its source is finite, so the quasi-steady kernel rejects a non-finite primitive/extra source (substitutes zero) and the implicit species floor rejects +inf rather than letting the positivity clip pass it through — see BACKLOG.md and the related memory notes.

References

• Saad, Y. Iterative Methods for Sparse Linear Systems, 2nd ed. (SIAM, 2003) — chapters 6 (Krylov methods) and 10 (preconditioning) for the general framework.
• Toro, E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed. (Springer, 2009) — chapter 10 for HLLC algebra; chapter 13 for implicit time-stepping.
• Millo, F., Gundlapally, S., Wang, W., Wahiduzzaman, S. 1D Simulation Models for Aftertreatment Components, Ch. 11 in 1D and Multi-D Modeling Techniques for IC Engine Simulation (SAE, 2020) — §11.2.5–11.2.7 for the Method-of-Lines / BDF DAE solution, the limit-point robustness problem, and the inner/outer adaptive mesh.

The source-side documentation tracker is in docs/perf-history.md (the audit-driven evidence trail phase by phase) and BACKLOG.md (open items, including the ones this doc records as deferred or rejected).