Skip to content

Implicit time integration — phased plan

Status: HISTORICAL / COMPLETE (kept as a design record). The implicit stack this plan scoped has shipped: rk = "implicit_euler" / "bdf2" / "radau_iia" with block-tridiagonal Newton, adaptive Δt, and analytic closure Jacobians are all on main. The optional Phase I-8 JFNK path shipped, benchmarked slower than the direct block-Thomas solve at every tested N, and was removed (see docs/numerics.md). For current usage and every NewtonConfig knob read implicit-solver-reference.md, not this plan.

User request 2026-04-29: add a fully-implicit Euler scheme that allows CFL > 1, plus a BDF2 scheme on top. This document scopes the work into reviewable phases and pins the validation gates each phase must pass before the next starts.

The motivation is real: the convective-only analyzer (commit 8fb4d70a) shows the explicit scheme pays an order-of-magnitude acoustic penalty for low-Mach engine flows (Mach ~0.05 → |u|+a)/|u| ~20×). A working implicit / low-Mach scheme would cut wall-clock by close to that factor for those cases.

This is not on the active critical path. It is a multi-week solver-architecture project. The active near-term work continues to be Cruise M interop + chemistry / catalyst polish; pull from this plan when the user explicitly asks.

What "fully implicit" means here

For implicit Euler the per-step nonlinear system is

U^{n+1} - U^n + Δt · R(U^{n+1}) = 0                 (IE)

where R(U) = (1/Δx) · (F_{i+1/2}(U) - F_{i-1/2}(U)) - S(U) is the spatial residual (HLLC fluxes + closure source terms). For BDF2 with constant Δt:

(3/2) U^{n+1} - 2 U^n + (1/2) U^{n-1} + Δt · R(U^{n+1}) = 0   (BDF2)

Both are solved per-step by Newton iteration:

J(U^{(k)}) · δU = -F(U^{(k)})
U^{(k+1)} = U^{(k)} + δU

with J = ∂F/∂U. The block sparsity of J is what determines whether the linear solve is cheap (block-tridiagonal Thomas, O(N)) or hard (Krylov on the full network).

Validation gates per phase

Every phase must clear:

  • Sanity: at CFL = 0.5, the implicit run matches the explicit RK2 run pointwise to ~10⁻⁶ on the Sod problem at t = 0.2 s. If they don't agree at small CFL the implementation has a bug the gate at CFL > 1 will hide.
  • Stability: at CFL = 10, 100 the run completes without NaN on the same Sod problem. Solution accuracy degrades — that is expected and is the headroom payoff.
  • Order of accuracy: dt halving at fixed CFL > 1 gives the expected convergence rate (1 for IE, 2 for BDF2) on the manufactured smooth problem (Phase 4 test fixture).
  • Closure consistency: the porous-plug bench (commits 5d13e69e + 8fb4d70a regime) reaches the same steady-state Δp at CFL = 50 as it does at CFL = 0.5.
  • Newton off-ramp: when Newton fails to converge (bad initial guess or genuinely stiff transient), Δt is halved and retried; the run must not silently produce garbage.

Phasing

I-1  Implicit Euler scaffolding (single pipe, 1st-order, no closures)
        │
        ▼
I-2  Closures into the residual (friction / heat / porous)
        │
        ▼
I-3  MUSCL with frozen reconstruction
        │
        ▼
I-4  Network coupling — junctions + plena
        │
        ▼
I-5  BDF2 on top of I-1..I-4 infrastructure
        │
        ▼
I-6  Adaptive Δt control (Newton-iteration heuristic)
        │
        ▼
I-7  CLI + GUI surface (scheme selector, Newton diagnostics)
        │
        ▼
I-8  (optional) Matrix-free JFNK for very large networks

Phase I-1 — Implicit Euler scaffolding

Scope. Single pipe, 1st-order spatial (no MUSCL), HLLC flux, ghost-cell BCs. No friction / heat / porous / catalyst yet — those land in I-2. State is just [ρ, ρu, ρE] per cell.

Residual. F(U) = U^{n+1} - U^n + Δt · R(U^{n+1}), R(U)_i = (1/Δx)(HLLC(U_i, U_{i+1}) - HLLC(U_{i-1}, U_i)).

Jacobian. Analytical ∂HLLC/∂U_L and ∂HLLC/∂U_R — each is a 3×3 block. Resulting J is block-tridiagonal with 3×3 blocks. N cells → 3N unknowns. Hand-coded analytic Jacobian is faster + more stable than finite-difference; see Toro §11 for the closed-form expressions.

Linear solve. Block-tridiagonal Thomas algorithm. O(N) flops, O(N) memory. Direct, no preconditioning. Reuse the existing pipe struct-of-arrays for the Jacobian storage: J_lower[N, 3, 3], J_diag[N, 3, 3], J_upper[N, 3, 3].

Newton iteration. Up to 8 inner iterations per step. Convergence test: ||δU|| / ||U|| < 1e-8. Backtracking line- search if the residual norm grows: halve the step until it shrinks. Off-ramp on Newton failure: halve Δt and re-enter the step from the same U^n.

Test fixtures.

  • tests/unit/test_implicit_euler_jacobian.py — verify the analytic Jacobian matches a finite-difference Jacobian on a smooth state to relative 1e-6.
  • tests/integration/test_implicit_euler_sod.py — Sod problem at CFL = 0.5 + CFL = 10. Sanity vs. explicit at CFL = 0.5; stability at CFL = 10.
  • tests/integration/test_implicit_euler_smooth.py — manufactured-solution smooth problem, dt-halving convergence, expect order 1.

Effort. ~5-7 days. The block-tridiagonal Thomas + analytic HLLC Jacobian are both literature-standard; the engineering work is the integration with aeolus1d.solver.timeloop and the test fixture build-out.

Decision points the user will care about.

  • Where to put it: separate aeolus1d.solver.implicit_euler module, dispatched from advance_single_pipe via a new scheme=SCHEME_IMPLICIT_EULER arm? Or a parallel advance_single_pipe_implicit entry point?
  • Whether to expose the Newton tolerance / max-iteration knobs in the case TOML or hard-code them. Recommendation: hard-code in I-1, expose in I-7 once they've earned it.

Phase I-2 — Closures into the implicit residual

Scope. Friction (darcy_const, colebrook), heat (constant_twall), porous (darcy_forchheimer) source terms enter the implicit residual + Jacobian. Catalyst chemistry stays explicit (operator-split) — see "Out of scope" below.

Why it matters. Friction at small Δx + low f is mildly stiff; porous at small K is genuinely stiff (the μ/K term scales as 1/K which can be 10⁹+). Treating these explicitly cancels most of the CFL > 1 headroom: even though the convective dt allows CFL = 100, the porous closure forces Δt < 1/(μ_ref/K) which can be smaller than the original explicit acoustic limit. They have to be implicit to be useful.

Implementation. Each closure exposes a residual_contrib and a jacobian_contrib callback:

def residual_contrib(U: np.ndarray) -> np.ndarray:
    """Return -S(U) added to the cell residual."""

def jacobian_contrib(U: np.ndarray) -> np.ndarray:
    """Return -∂S/∂U as a (N, 3, 3) block array, added to the
    diagonal of J."""

Source terms are local (no flux coupling), so they only touch J_diag. Wire through the closure registry (Phase A1 — commit 86f0be96); friction / heat / porous already live there, just need new methods.

Test fixtures. Extend the I-1 Sod test with friction + porous variants. The porous-plug bench (tests/integration/test_porous_plug.py) is the natural home for the high-stiffness gate.

Effort. ~3-4 days for the three Tier-1 closures.

Phase I-3 — MUSCL with frozen reconstruction

Scope. Add 2nd-order spatial reconstruction (MUSCL or MP5) to the implicit residual with frozen reconstruction — the slope limiter is evaluated once per Newton iteration on the current Newton iterate, not differentiated through.

Why frozen. Differentiating MUSCL through the limiter (e.g. minmod, van Leer) gives a 5-block-wide Jacobian — workable but ~3× the storage + flops of block-tridiagonal. Frozen reconstruction keeps the Jacobian block-tridiagonal at the cost of one extra Newton iteration per step (the slope re-converges as the iterate settles). For most cases the extra Newton iteration is far cheaper than tripling the linear solve.

Implementation. Reconstruct slopes σ_i^{(k)} from U^{(k)} at each Newton iterate; treat them as constants in ∂R/∂U so the Jacobian sparsity stays block-tridiagonal.

Validation. Manufactured smooth problem must show order-2 spatial convergence at fixed CFL > 1 (BDF2 in Phase I-5 will also need this — borrow the fixture).

Effort. ~3-4 days. The frozen-reconstruction pattern is well-understood; the work is bookkeeping + a careful test fixture.

Phase I-4 — Network coupling

Scope. Multiple pipes coupled at junctions / plena. The plenum's ODE in [m, m_species, E] enters the implicit residual; the junction Riemann at each pipe end couples that pipe's last cell to the plenum DOFs.

Linear solver options.

  • Schur complement on plenum DOFs. Eliminate each pipe's interior DOFs against its block-tridiagonal piece, leaving a small (~plenum-count × N_species sized) reduced system on the plenum DOFs. Solve that, back-substitute. Direct, exact, scales linearly with pipe count. Best fit when plena are "few + small" — typical engine cases.
  • Block-Jacobi GMRES. Each pipe's tridiagonal as a block preconditioner; GMRES handles the coupling iteratively. Better fit for dense networks (catalyst with hundreds of parallel channels); the Schur block gets large.

Recommendation. Schur complement first — it's exact and matches typical engine network shapes. Defer GMRES to I-8.

Effort. ~7-10 days. The plenum residual + Jacobian is straightforward; the Schur reduction is the new piece.

Phase I-5 — BDF2

Scope. BDF2 time discretization on top of the I-1..I-4 infrastructure. Reuse the same residual (modulo the time- discretization terms) and Jacobian.

Implementation.

  • Residual: replace U^{n+1} - U^n + Δt · R(U^{n+1}) with (3/2)U^{n+1} - 2 U^n + (1/2) U^{n-1} + Δt · R(U^{n+1}). Jacobian: J = (3/2) I + Δt · ∂R/∂U instead of I + Δt · ∂R/∂U.
  • Startup: first step uses implicit Euler (so U^{n-1} is available for step 2). Crank-Nicolson would also work but is not L-stable; IE is the safer startup.
  • Variable-Δt BDF2 coefficients: if the Phase I-6 controller is in place, use the variable-step BDF2 form (Hairer-Wanner II.5.20).

Test fixtures. Order-2 dt-halving convergence on the manufactured smooth problem at fixed CFL > 1.

Effort. ~2-3 days on top of I-1..I-4 — the heavy infrastructure is already there.

Phase I-6 — Adaptive Δt

Scope. Auto-adjust Δt based on Newton iteration count. PI controller:

Δt_{n+1} = Δt_n · (target_iters / actual_iters)^p

with target_iters = 4, p = 0.4. If Newton fails, halve Δt and retry. If it converges in 1 iteration two steps in a row, double Δt (capped at the user's max).

Optionally: an embedded local-error estimator (compare BDF2 step against an IE step) for proper local-error control. Defer the embedded estimator unless the iteration-count heuristic proves inadequate.

Effort. ~2-3 days for the iteration-count heuristic; ~3-4 days extra for the embedded error estimator if needed.

Phase I-7 — CLI + GUI surface

Scope. Expose the schemes through the user-facing edges.

CLI. New TOML field on [case]:

time_scheme = "rk2_explicit"  # default — current behaviour
# or:
time_scheme = "implicit_euler"
time_scheme = "bdf2"

Plus tuning knobs (newton_tol, newton_max_iters, adaptive_dt).

GUI. Scheme selector on the run-settings strip. Live Newton-iteration count in the progress info. CFL Analyzer dock gains a "with implicit, dt headroom = X×" banner reading the ratio computed by format_report (commit 8fb4d70a).

Effort. ~3-4 days, mostly GUI plumbing.

Phase I-8 — Matrix-free JFNK (optional)

Scope. Replace the analytic Jacobian assembly with finite- difference Jacobian-vector products inside GMRES. Skips the storage cost of the full Jacobian (relevant for networks of hundreds of pipes); pays a ~10× hit per linear iteration due to the FD evaluations.

Pull this only if Phase I-4's Schur-complement direct solver runs out of memory or scales poorly on a real case the user brings. For typical engine networks (10s of pipes, ~10 plena), the direct solver is faster and simpler.

Effort. ~7-10 days. Heavy testing burden — JFNK adds many degrees of freedom (FD step size, preconditioner choice, GMRES restart, etc.) that all interact.

Phase I-9 — Radau IIA(3) (added 2026-04-30)

Scope. A 3-stage 5th-order fully-implicit Runge-Kutta time integrator on top of the I-2/I-3 spatial + closure infrastructure. L-stable; stiffly accurate so Y_3 = U^{n+1} at convergence without a separate weighted-sum update step. Single-pipe at first; network follow-up later.

Why. BDF2 is 2nd-order. For accuracy-critical smooth-flow runs (acoustic-mode validation, smooth wave-propagation benchmarks, model-fitting where the temporal-error floor is the bottleneck) the 2nd-order truncation caps how big a Δt you can take before accuracy degrades visibly. Radau IIA is 5th-order on smooth flow — at the same large CFL it produces a ~1000× smaller error than BDF2.

Implementation. The coupled stage system has 3 stages × N cells × 3 conservatives = 9N unknowns per pipe. The Jacobian in (cell, stage·m) layout is block-tridiagonal in cells with 9×9 blocks per cell; the Kronecker structure J = I − Δt · A ⊗ J_spatial lets us reuse the existing spatial-Jacobian assembly verbatim and just feed bigger blocks to the same block-Thomas solver.

Simplified Newton. Jacobian computed once at the predictor state U^n and reused across all Newton iterations + all stages. Standard for Radau — the linear solve dominates per- step cost and refreshing the Jacobian per iteration buys very little for the identity-dominant systems at moderate Δt.

Validation.

  • Stiffly-accurate Butcher tableau pin: b = A[-1].
  • Sanity at low CFL: matches BDF2 to ~1e-3 on smooth flow.
  • Stability at CFL > 1: Sod-100 at CFL = 2 stays positive.
  • Order benefit: at the same large CFL, Radau err ≪ BDF2 err (the formal 5⁵/2² ratio per Δt halving) — the comparative test is cleaner than measuring 5th-order convergence in isolation since Radau hits the FD/Newton noise floor at small Δt.
  • Closure compatibility: composes with porous Darcy via the shared _residual_full / _augment_jacobian_with_closures path, recovering the analytic exponential decay to ~1% better than BDF2 at the same Δt.

Effort. ~1 day to integrate cleanly; the Kronecker structure + block-Thomas reuse made this much smaller than the I-1 first- cut.

Phase I-9 follow-ups.

  • Network coupling. ✅ Landed alongside the single-pipe driver. advance_network_implicit dispatches to the Radau Newton step when rk = "radau_iia"; per-pipe block-split coupling reuses the same plenum-ghost / Dirichlet / NUMERICAL chain-rule machinery as IE/BDF2. Radau is an RK scheme (no multi-step history) so dt-halvings don't need the BDF2 fallback the IE/ BDF2 path has.
  • Eigen-decomposition simplified Newton. Hairer-Wanner II.5.5: A^{-1} has eigenstructure (one real eigenvalue + one complex pair) so the 9×9 linear solve decouples into a (3×3) real solve + a (3×3) complex solve. ~3× speedup at the cost of complex arithmetic. Not on the critical path until a benchmark shows the 9×9 solve is the bottleneck.

Out of scope (deliberately deferred)

  • Catalyst chemistry inside the implicit residual. Chemistry already runs as a Strang-split sub-step driven by an independent stiff-implicit integrator — aeolus1d.solver.chem_bs_jit is a Bader-Deuflhard semi- implicit Bulirsch-Stoer extrapolation kernel with a per-cell frozen Jacobian and polynomial-extrapolation step control. That sub-step does not change in the I-1..I-5 phases; the new implicit Euler / BDF2 schemes only solve the advection half (conservative variables [ρ, ρu, ρE], 3×3 per-cell Jacobian blocks). The control flow stays:

    for each timestep: advection_step(U) # I-1..I-5: implicit Euler / BDF2 chemistry_step(Y, E) # unchanged: Bader-Deuflhard

Folding chemistry into the same Newton matrix as the conservative variables would inflate the per-cell block to 13×13 or more (5-10 gas species + 3-5 surface coverages + the [ρ, ρu, ρE] triple), and chemistry rate-coefficient stiffness (rate constants spanning 10⁶ orders) mixes badly with the acoustic stiffness inside one matrix. Full coupling is a separate research-scale project — not a blocker for the user's stated low-Mach motivation. If a case ever needs sub-acoustic time scales coupled to fast surface kinetics (e.g. ignition fronts) that's the time to revisit; for steady-state catalyst light-off / conversion sweeps the operator-split path is empirically tight. * Limiter-aware MUSCL Jacobian. See Phase I-3 — the frozen- reconstruction trick avoids it; pursue only if the extra Newton iteration becomes painful. * Compressible-low-Mach asymptotic preserving (AP) schemes. Genuinely low-Mach (M < 0.01) is where the pressure-projection family (Klein, Munz et al.) wins decisively. The implicit- HLLC family in this plan handles M ~ 0.05-0.5 well; AP is a separate ~3-month research project, justified only if a real Cruise M / catalyst case needs M < 0.01 resolution.

Total effort estimate

  • Phases I-1 through I-5 (the user's stated request, end-to-end working): ~3-4 weeks of focused work.
  • Phase I-6 + I-7 (production polish): ~1 week extra.
  • Phase I-8 (optional, only if needed): ~1.5 weeks extra.

Whole stack: ~5-6 weeks. The minimum viable feature (single-pipe implicit Euler with closures, no GUI) is ~2 weeks (I-1 + I-2).