Troubleshooting

A field guide to the most common failure modes you'll hit running an aeolus1d case, with the diagnostic flags + root causes for each. Organised by what you SEE first.

"My run diverged / blew up"

Symptoms. Solver prints NaNs, exits with RuntimeError: non-finite state, or a --info-every-n line shows p_max = inf. The HDF5 (if written) has NaN values in the final-state pipe arrays.

Most likely cause. CFL violation or unphysical initial conditions.

Fix in order:

1. Check CFL. The default is 0.5 for explicit RK; if you've lowered it, raise back to 0.5. If you've raised it past 0.7, drop it.
2. Drop t_end first. If t_end is 1.0 s, try 0.001 s. A 1-second simulation that diverges at t=10⁻⁵ s reveals the issue much faster.
3. Check the initial state. [pipe.init].T = 0 or p = 0 produces 1/T and 1/p denormals on the first step. Use T ≥ 200 K, p ≥ 100 Pa.
4. Use --info-every-n 1 to see EVERY step's dt. A shrinking dt (e.g. 1e-5 → 1e-7 → 1e-9) signals adaptive-Δt is fighting an unstable closure.
5. Switch to implicit. If explicit RK can't handle your stiffness (typical with porous-Forchheimer or chemistry), set rk = "implicit_euler" or "bdf2" in [case]. The implicit drivers don't have a CFL constraint. One exception: a [[cylinder]] using the default orifice valve coupling (cylinder_orifice_valve_flow, ADR-0008) is rejected at build under an implicit rk — the valve sources are explicit-driver-only — so either keep an explicit rk or set cylinder_orifice_valve_flow = false first.

"Newton didn't converge"

Symptoms. Console shows Newton failed at t = X.X (n_iters = 12, F_norm / U_norm = ...) and the run terminates.

Most likely cause. The Newton residual norm isn't reaching config.tol within config.max_iters, OR the line search exhausts its backtracks.

Fix in order:

1. Loosen the tolerance. Set [case] newton_tol = 1e-6 (default is 1e-8). The default is conservative; 1e-6 is fine for most physics studies.
2. Allow more iterations. [case] newton_max_iters = 50 (default 8). For stiff catalyst kinetics, 50 is often needed.
3. Reduce the Δt. The Newton matrix conditioning degrades with larger Δt × residual; smaller steps stabilise. Three real knobs, in order of preference: lower [case] cfl (implicit Δt = cfl · acoustic dt, so halving cfl halves the step); turn on [case] adaptive_dt = true and lower adaptive_cfl_max so the Newton-difficulty controller caps the step; or pin the step exactly with [case] dt_fixed = 1e-6 (overrides the CFL calculation entirely). Run aeolus1d analyze-cfl case.toml first to see which pipe/cell sets the acoustic dt.
4. Loosen Picard. Network cases use a Picard outer loop; [case] picard_max_iter = 5 (default 3) gives the outer loop more chances before declaring failure.
5. Activate adaptive Picard. [case] adaptive_picard = true automatically grows picard_max_iter when the convergence rate signals saturation. The runner-dock Picard badge surfaces saturation in real time.
6. Check the closure pack. A friction_kind with wrong units (D_h in mm instead of m) produces huge source values that the Newton can't track. Cross-reference closures reference for the expected units.

"QS Picard non-convergence"

Symptoms. solver = "quasi_steady" raises QS junction-pressure Picard failed to converge after qs_picard_max_iter outer iterations. Often the iterate is visibly near the answer — |Δjp| keeps oscillating between two values, or the relative net mass imbalance plateaus at ~0.5–2 % and refuses to drop below qs_picard_tol = 1e-6. Typical on TWC AT systems with a stepwise inlet ṁ transient, or on cases that init at the outlet pressure (zero pressure gradient → blind FD probe at startup).

Most likely cause. A stale-secant 2-cycle in the junction- pressure Newton-damping update jp_new = jp − qs_picard_relax · net / slope. The damped secant can settle into a non-zero fixed point that is converged in the |Δjp| sense but not the net-imbalance sense.

Algorithmic fix (default since MR !131, 2026-05-31). The QS driver now wraps the secant update in a sign-bracket bisection guard: once net flips sign between two iterates, the next step bisects the bracket instead of trusting the stale slope. This kills the 2-cycle by construction and is the right fix for the failure mode the qs_picard_relax knob was originally introduced to paper over. No flag — it is always on.

Diagnostic scripts. Two helpers ship under scripts/ for when a stuck case still falls through the guard:

• scripts/diagnose_qs_picard_2cycle.py — runs the case with the Picard iterate trace enabled and flags any junction that visits the same (jp, net) pair twice within qs_picard_max_iter. Useful as a first triage: a 2-cycle signal means the bracket guard never engaged (likely a zero-gradient bootstrap blind spot), a monotone-but-slow trace means the relaxation is too aggressive for the transient.
• scripts/analyze_qs_picard_trace.py — post-processes the trace HDF5 / JSON and produces a per-junction convergence rate, the final |net|, and the iterate at which the bracket first formed. Plot this when triaging a case that fails late (well after the inlet transient settled).

qs_picard_relax tuning ladder. When the bracket guard gets you converged but slowly, or for legacy cases that need the old behaviour preserved, the relaxation knob is still exposed and follows a simple ladder (src/aeolus1d/io/schema.py:247, default 0.7):

1. qs_picard_relax = 0.7 — default. Works on near-steady or smoothly-evolving cases; matches the historical QSNetworkConfig design and is bit-equivalent for every validated v1 case.
2. qs_picard_relax = 0.3 — sharp inlet-condition transients (e.g. a stepwise ṁ drop/rise over ~1 s on a TWC AT system). Cuts the secant amplitude so the bracket guard sees the sign flip earlier and lands on the true zero-imbalance solution at a modest convergence cost.
3. qs_picard_relax = 0.1 — very stiff regimes (cold-start light-off ramp, large porous-α + transient ṁ at the same time). Treat as the last-resort knob; convergence cost scales roughly inversely with the relaxation.

In all three cases the sign-bracket guard is doing the algorithmic heavy lifting — the knob just trims the secant's overshoot amplitude when the bracket forms late.

"My case used to load and now it is rejected"

Symptoms. CaseValidationError naming an unknown key, with a did you mean …? suggestion and the file:line — on a case that ran before.

Most likely cause. It was never doing what you thought. Unknown keys in [[bc]], [[junction]], [[turbo]] and the nested sub-blocks used to be a warning: the parser dropped the key and left the field at its default. [case] and [[pipe]] always raised; the split was the dangerous half, because every BC field is read with a default, so a mistyped pO silently loaded p0 = 0.0 and ran a near-vacuum boundary.

Fix. Take the suggestion — it now normalises confusable glyphs, so pO → p0 and TO → T0 are matched even though a two-character key scores below a plain fuzzy-match cutoff. If a key is genuinely obsolete, delete it; the error lists the valid keys for that block.

Related refusal: p0 must be positive / T0 must be positive on a BC. That is the consequence the unknown-key error prevents, caught directly — a stagnation inlet computes a0 = sqrt(γ·R·T0) and divides by it, so a zero used to kill the run at step 0 with a bare ZeroDivisionError naming nothing.

"Everything is at near-vacuum / density ~1e-5"

Symptoms. The run completes, no warning, and the final state has rho around 1e-5 kg/m³ with pressures that look four or five orders of magnitude too small.

Most likely cause. Pressures written in bar. p = 1.0 means one pascal. This is the single most common unit slip for anyone arriving from GT-POWER or AVL, and nothing in the schema catches it — 1.0 is a perfectly valid float for a pressure.

Fix. Multiply by 1e5. aeolus1d is strict SI throughout: Pa, K, m, s, kg. A load-time warning now fires for this ("p = 1 Pa with T = 300 K — pressures are in PASCAL"), but only when a temperature in the Kelvin band shows the case is in physical units at all: several bundled examples (01_single_pipe_sod among them) are the classic non-dimensional Riemann problems with p = 1.0 and no temperature, and warning on those would train you to ignore the message.

"negative density (or e_internal ≤ 0)"

Symptoms. _line_search_trial_jit rejects every trial, exits with Newton step produced non-physical state on cell N. Or HDF5 has rho < 0 in a final cell.

Most likely cause. Either an extreme expansion wave reached near-vacuum conditions, OR a momentum-conservation BC mismatch is forcing mass out of a closed pipe.

Fix in order:

1. Check far BCs. Both ends of every pipe need something filling ghost cells. A wall on a pipe that should have a pressure-outlet creates a closed reservoir; small numerical error grows over many steps.
2. Refine cells near shocks. A first-order scheme on a sharp shock can produce rho ≈ 1e-3 and trip the positivity guard. Either use scheme = "muscl" or double n_cells.
3. Check the slope limiter. The limiter is a case-level field ([case] limiter, applied to every MUSCL pipe) and already defaults to "minmod" — the safest choice. If your case sets limiter = "van_leer" (sharper but less robust near genuine vacuum) or "none", set it back to "minmod".
4. Reduce the step near the expansion. There is no density floor knob — clamping ρ would hide conservation bugs. Lower [case] cfl (or, on the implicit drivers, enable adaptive_dt = true so Newton difficulty shrinks Δt automatically) so the expansion is resolved in time; aeolus1d analyze-cfl case.toml --cells shows which cells are closest to the positivity limit.

"CFL exploded"

Symptoms. --info-every-n shows dt < 1e-15. The adaptive-Δt path is stuck shrinking forever.

Most likely cause. A high-frequency oscillation has set in, typically from a discontinuity-amplifying scheme on near-discontinuous initial data.

Fix in order:

1. Run aeolus1d analyze-cfl case.toml. It reports the acoustic-CFL dt per pipe before you run anything, and --cells pinpoints the constraining cells — answering "which pipe is driving dt to zero" without a solver run.
2. Switch to MP5 or MUSCL. First-order is dissipative but stable; MP5 + a limiter handles sharp gradients cleanly.
3. Check the initial state for steps. A Riemann problem needs [pipe.init] kind = "riemann" with explicit left/right states, not a hand-built step in a uniform init.
4. Add a small artificial viscosity. Set [pipe.friction] kind = "darcy_const" f = 0.01 D_h = (pipe D) — tiny friction smooths high-frequency artifacts at the cost of some physical damping.

"Solver hangs"

Symptoms. No output for minutes; the runner dock's progress bar stops advancing.

Most likely cause. The Numba JIT is compiling, OR a Python callback is blocking.

Fix in order:

1. First run is slow. Numba compiles every kernel on first call — a one-time ~19 s (measured 2026-08-08) for the full quasi-steady + chemistry kernel set on a cold cache (~8 s warm). Run aeolus1d warmup once to do the compile out loud and prime the on-disk cache; aeolus1d doctor reports whether the cache is warm or cold. Subsequent runs reuse the cache.
2. Check py_friction / py_heat / py_porous. A user callable that calls back into matplotlib, GUI, or a slow external library can block. Test the callable in isolation first.
3. GIL contention. With parallel_pipes = true and a Python override on any pipe, the threads collide on the GIL. Use native closure kinds where possible.

"GUI dependencies missing (No module named '_lzma')"

That message was wrong, and is fixed as of 2026-08-17. If you are on an older build and see it, do not run the pip install ... [gui] it suggests: it cannot help.

_lzma is part of the Python standard library, not the [gui] extra. It is a CPython extension module compiled only when liblzma headers are present at interpreter build time, so it is absent from an interpreter built without them — most often a pyenv / python-build compile on a machine that did not have liblzma-dev installed. Nothing you pip install can supply it.

Confirm the diagnosis:

python -c "import lzma"     # ModuleNotFoundError: No module named '_lzma'
aeolus1d doctor             # reports INCOMPLETE PYTHON BUILD and what to do

Fix it by giving the build the headers and rebuilding:

sudo apt install liblzma-dev      # Debian / Ubuntu
sudo dnf install xz-devel         # Fedora / RHEL
brew install xz                   # macOS / Homebrew

pyenv install --force 3.11.9      # then rebuild the interpreter

A distribution Python (apt install python3) or a python.org installer already has it, so switching interpreters also works. The same applies to _bz2, _ssl, _sqlite3, _ctypes, _tkinter, readline and zlib — aeolus1d doctor lists every one it finds missing, with the package that supplies it.

Note that aeolus1d's solver does not need _lzma: a case loads and runs on an interpreter without it. It surfaces through a GUI dependency's import chain, which is why it looked like a GUI problem.

"GUI runs the case but I see no plot"

Symptoms. The runner shows "Run completed"; the network view loads colours, but the Plot dock shows "select a component".

Most likely cause. The case has no [[tracer]] or --trace-dt wasn't set, so only the final state was written. The Plot dock falls back to "click a component" prompt rather than auto-selecting.

Fix. In the GUI runner dock, tick "Trace every Δt:" and set the cadence (defaults to 1e-4 s, capped at 1e6 s), then re-run. On the command line pass --trace-dt 1e-4 to aeolus1d run. There is no [case] TOML key for this — tracing is a run-time option (CLI flag / GUI control), not a case field, and adding trace_dt under [case] is rejected as an unknown key. Default is no trace, so only the final state is written.

"Cylinder pressure looks wrong"

Symptoms. The cylinder p(θ) trace shows compression but combustion adds way too much pressure (or too little).

Most likely cause. Wiebe heat-release parameters off, or the combustion m_fuel is misset.

Fix in order:

1. Check m_fuel in [cylinder.combustion]. Units are kg per cycle per cylinder. A 500 cc cylinder at λ = 1 gets ~5e-5 kg/cycle of fuel; an order-of-magnitude error is obvious in p_max.
2. Two-zone vs single-zone. Bulk T (single-zone) is inflated near peak combustion because constant γ doesn't capture the burnt-mixture composition. Set [cylinder.combustion] two_zone = true and read last_T_burnt (the GUI plot dock surfaces it as T_burnt when traced).
3. Burn duration. Diesel: 40–60° CA. Naturally aspirated SI: 30–50° CA. Turbocharged SI: 25–35° CA. If you're off-window, the Wiebe shape produces wrong p_max.
4. Spark / injection timing. theta_start_deg in [cylinder.combustion] is the crank angle (degrees) at which the heat release ramp starts. TDC is 360 deg by convention; advance is < 360, retard is > 360.

"Everything conserves, but the answer is still wrong"

Symptoms. Mass and energy balance to machine precision, the golden snapshots pass, nothing warns — and a temperature, a residual mass or an emissions number is plainly unphysical.

Most likely cause. An attribution error rather than a conservation one. A transfer can move the wrong amount and still balance perfectly, because the same quantity is credited to one side and debited from the other. Two shipped examples, both fixed: the exhaust valve shed mass at cp(gamma_advection) while the cylinder held cp_eff, retaining 26 % of the heat release; and volume.x_burnt carried the burn progress, so gamma_eff reinterpreted an unchanged (m, E, V) and pressure jumped 29.6 % at every crank wrap.

No conservation check can see either. Tightening one will not help — it is the wrong instrument.

Fix in order:

1. Turn on the attribution checks and re-run:

   AEOLUS_PHYSICS_ASSERTS=1 aeolus1d run case.toml
   

   This asserts two things a conservation check cannot: that every unit of mass a cylinder sheds carries the cylinder's own specific enthalpy, and that no derived quantity (p, T) jumps while its arguments (m, E, V) stand still. Both are off by default and cost nothing when off.

   Use AEOLUS_PHYSICS_ASSERTS=warn to print to stderr and keep going, which is the useful mode when bringing up a case that trips the check immediately.

2. Read the ratio it reports. A shed-enthalpy ratio near cp(gamma_advection)/cp_eff (≈ 0.74 for 1.4 against 1.27) means something is rebuilding enthalpy from the advection gamma. A pressure jump at fixed (m, E, V) means something changed the mapping rather than the state — gamma_eff and its burnt fraction, the EOS branch, or a resume path applying fields in the wrong order.

3. Suspect any cp = gamma*R/(gamma-1) near a component boundary. Crossing between a 0-D volume and a 1-D pipe must be written in a form that presumes neither side's gamma — h = (E + pV)/m names no thermodynamic coefficient at all and is exact under gamma_eff, under caloric_thermo, and under a real-gas EOS alike. See the "Two γ's, and the seam between them" section of cylinder-reference.md.

"Chemistry too slow on AVL-style runs"

Symptoms. A QS driver case with a Brinkmeier TWC (AVL Cruise-M MFR_transition_SC630, MFR EAT bench variants, anything that bolts the oxygen-storage / NOx mechanism onto a catalyst pipe) walls at ≫ 1× realtime even after JIT compile has settled. Profiler shows _qs_surface_chem_step dominating, with the per-substep LSODA RHS closure showing up as a hot Python frame.

Most likely cause. Either the fused @njit Brinkmeier RHS isn't engaging (so LSODA is paying the per-substep Python closure tax), or the chemistry sub-step is firing every macro step through a regime where it could safely super-cycle.

Fix in order:

1. Verify qs_chem_fused_rhs = True. This is the default since MR !134 (2026-05-31, src/aeolus1d/io/schema.py:359). It swaps the per-substep Python RHS closure in _qs_surface_chem_step (dict comprehension + reshape + extras-reorder + ascontiguousarray casts) for a single @njit call into _qs_brinkmeier_rhs_fused. Measured 2.08× speedup on AVL Cruise-M MFR_transition_SC630 (97.7 s → 47.0 s on a 12 s sim, 8.15× → 3.92× realtime), bit-equivalent to the Python RHS — LSODA takes the identical step sequence. Anchored by tests/integration/test_qs_chem_fused_rhs.py. Non-Brinkmeier closures (generic builder, NH3, py_reaction kinds) lack _jit_params and silently fall back to the Python RHS — no per-case opt-out needed. If you have manually set this to False in a case TOML, undo that unless you are diagnosing a Brinkmeier-specific regression. scripts/profile_post_fused.py is the reference profile to confirm the fused path is engaged.

2. Enable the chemistry super-cycle with a settle gate. chem_every_n = N (default 1) fires the chemistry sub-step once every N macro steps over the accumulated Σ Δt (transport + steady-energy still run every step). The same knob drives the Strang and implicit-network Lie super-cycle, not just QS. It is an N× lever on chemistry-dominated walls, but it is catastrophic during light-off — at the transient a chem_every_n = 4 run took 2.25 M steps. The fix is to combine it with a settle gate so the super-cycle self-suspends during transients (the two settle gates below are consulted only when chem_every_n > 1):

3. qs_chem_settle_tol (default 0.0 = disabled): fire chemistry on EVERY macro step while the lagged HRR fractional change exceeds the tol; only honour the every-N skip once HRR settles below it. A small positive value like 1e-2 (1 %/step) is a safe start.

4. qs_chem_coverage_settle_tol (default 0.0): companion gate, coverage-aware. The HRR gate is energy-based and blind to a still-relaxing oxygen-storage θ on a hot catalyst — set this to a small absolute per-fire Δθ (e.g. 1e-3) to keep NOx / θ fidelity while still super-cycling once both gates settle.
5. chem_interval (default 0.0, seconds): simulated-time companion to chem_every_n, shared by the Strang, implicit-network, AND quasi-steady drivers. Forces a chemistry fire whenever the accumulated Σ Δt since the last chem call exceeds the interval, no matter how cheap the macro step is. Typical 1e-3–1e-1 s for transient AT — track the chemistry timescale, not the acoustic / advective one. Fires on EITHER the step counter or the interval.

The gates are independent and fire on EITHER — set the ones that match the regime you're super-cycling through, leave the rest at 0.0.

1. Pick chem_solver for the mechanism. The default chem_solver = "auto" (src/aeolus1d/io/schema.py:375) routes to "bs_jit" only when the pipe matches the JIT kernel's Arrhenius-only, single-species contract (py_reaction_params set, no multi / surface / extras); everything else (including Brinkmeier surface chem) falls back to LSODA. For an Arrhenius-only mechanism that auto is not routing to bs_jit (e.g. you've added extras for diagnostics that block the contract), set chem_solver = "bs_jit" explicitly to force the @njit Bader-Deuflhard semi-implicit Bulirsch-Stoer kernel. Do not set it on stiff oxygen-storage TWC (Brinkmeier) — LSODA's variable-order BDF is the right integrator there, and the fused-RHS path already delivers the 2.08× win without changing the integrator choice. chem_solver = "fixed_step" is the backward-Euler path and is opt-in for non-stiff mechanisms only (overshoots the Brinkmeier Ce-OSC transient).

Got an issue that isn't here?

• Check the validator output. Phase B10 emits a list of every spec problem before the build path runs — fixing those usually fixes the symptom.
• Run aeolus1d analyze-cfl case.toml for a pre-run picture of the acoustic time-step budget (per pipe; --cells for the constraining cells).
• Run with --info-every-n 100 to see the time-stepping cadence. A near-zero dt or growing Newton iteration count is the strongest signal of a physics or closure problem.
• File an issue with the case TOML + run command + the last few lines of solver output, and paste the aeolus1d doctor block (Python / package versions, env vars, numba-cache state) — it is the environment dump SUPPORT.md asks for. The bundled examples are the reference for "what's expected to work."