ADR-0002 — Direct block-Thomas implicit solve; matrix-free JFNK removed

Status: Accepted

Context

The implicit time integrators (Implicit Euler / BDF2 / Radau IIA) solve F(U) = U − U_hist − Δt·R(U) = 0 per step with Newton's method. For the 1-D first-order/frozen-MUSCL spatial residual, the Jacobian J = I − Δt·∂R/∂U is block-tridiagonal with 3×3 blocks (each cell couples only to i−1, i, i+1 through one HLLC face per side).

A matrix-free Newton–Krylov (JFNK / GMRES with a finite-difference J·v) path was prototyped as a potentially more scalable alternative, with both block-Jacobi and ILU(0) preconditioners.

Decision

Use a single direct path: a block-Thomas sweep that exploits the block-tridiagonal structure exactly (no fill-in, O(N) work, captures the off-diagonal acoustic + convective coupling). The JFNK path was removed (v0.9); _jfnk_step and the block-Jacobi apply kernel no longer exist in the tree.

Consequences

• The implicit solve is simple, fast, and exact for the systems aeolus1d targets (≤ a few-thousand cells × 3 unknowns); there is one Newton-solve path to maintain, not two.
• A chord (lazy-Jacobian) heuristic sits on top for the per-iteration cost (NewtonConfig.chord_max_age), reproducing full Newton bit-exactly at chord_max_age=0.
• A future 2-D/3-D extension would reopen the preconditioner question (line preconditioning no longer collapses to a 1-D sweep); the block-Thomas design is forward-compatible but not a 2-D answer.

Alternatives considered

• JFNK + block-Jacobi — rejected: catastrophically worse (hundreds of ×). Dropping the off-diagonal coupling makes the preconditioned operator less diagonally dominant at large Δt, not more — the opposite of what helps.
• JFNK + ILU(0) — rejected: ~1.8–1.9× slower at every N tested, because the only competitive preconditioner (ILU(0)) collapses back to the very block-Thomas factorisation the direct path already performs.

The generalisable lesson: 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.

References

• docs/numerics.md — "Why the direct block-tridiag solve is the right shape" (the full benchmark write-up).
• Related: ADR-0001 (keeping the stiff reaction block out of this Jacobian is what preserves the tridiagonal shape).