Skip to content

Phase I-Y — Picard outer iteration on the implicit network solver

User-facing guide to the Picard outer loop, the adaptive controller, the diagnostic surfaces, and when to touch each knob. Code reference is aeolus1d.solver.implicit_network; bundled example examples/54_adaptive_picard_turbo.toml / examples/plots/plot_54_adaptive_picard.py.

What the loop is for

The implicit network driver wraps the pipe Newton solve in an outer Picard iteration that refreshes pipe-interior-coupled ghosts between iterations:

for picard_iter in range(picard_max_iter):
    snapshot pipe_Un  (only on iter 0; warm-start otherwise)
    for pipe in network:
        Newton-solve to convergence
    apply_port_bcs over every junction → refresh ghosts from
        the converged pipe interiors
    if max-relative ghost change < picard_tol: break

Components whose ghost state depends on the pipe interior — Compressor, Turbine, LODI ports, py_junction_bc callables, neck pipes — were the reason for adding the loop. Phase I-X (the predecessor) held those ghosts fixed across the Newton solve, giving an O(Δt) splitting error per step. Phase I-Y folds the coupling into the implicit step.

Default behaviour — what you get for free

NewtonConfig.picard_max_iter = 3, picard_tol = 1e-6 are the out-of-the-box defaults. For a typical engine config:

Network type Picard iters per step Cost vs Phase I-X
Plenum-only (no pipe coupling) 2 (iter 0 + 1 check) ~1 % overhead
Compressor / Turbine / LODI 2-3 ~5-15 % overhead
Tight turbo + small inertia saturates at 3 ~25 % + accuracy

The loop's "cost" is the extra Newton solve(s) on Picard re-iterations. Phase I-Y-warm makes those cheaper by warm- starting Newton from the previous iter's converged state (Newton needs ~2 iters per re-iter instead of ~4).

picard_max_iter == 1 reproduces the Phase I-X block-Gauss-Seidel-without-iteration behaviour exactly. Use this if you need to reproduce a Phase I-X result.

TOML knobs

The [case] block exposes:

[case]
# ... existing fields (cfl, t_end, rk, ...) ...

# Phase I-Y outer iteration knobs.
picard_max_iter = 3       # initial / static budget
picard_tol      = 1.0e-6  # max-relative ghost change to exit early

# Phase I-Y-adaptive controller (off by default).
adaptive_picard          = false
adaptive_picard_min      = 1
adaptive_picard_max      = 8
adaptive_picard_n_easy   = 5

When adaptive_picard = true, the driver manages picard_max_iter per step:

  • saturated (used >= max AND ghost_change > picard_tol) → grow by 1, reset easy-streak.
  • easy with headroom (used + 1 <= max for adaptive_picard_n_easy consecutive steps) → shrink by 1, reset streak.
  • otherwise → hold, reset streak unless still in the easy regime.

Diagnostic fields on LiveControl

The implicit driver writes four fields after every completed step. Outer callers (CLI, GUI, custom drivers) can poll them between steps without instrumenting the solver:

Field Meaning
last_picard_iters Number of Picard iters used (1 ≤ N ≤ picard_max_iter_managed).
last_picard_max_iter_managed Controller-managed budget. = static picard_max_iter when adaptive=False.
last_picard_total_newton_iters Σ Newton iters across pipes × Picard iters. Drops under warm-start.
last_picard_ghost_change Max-relative ghost change at exit (the < picard_tol test).

Reading the CLI / GUI surface

When --info-every-n is on (CLI flag, or "Info every N steps" in the runner-dock toolbar), each per-step line carries a trailing diagnostic block:

  step    37230   t = 0.05 s   dt = 1.34e-04 s   picard = 2/3   newton = 8
  • picard = 2/3 — the network used 2 Picard iters out of the managed budget of 3.
  • newton = 8 — total Newton iters summed across pipes × Picard iters this step.

The GUI runner-dock parses both tokens into header-bar badges:

  • Picard: 2/3 — blue when N < M, red when saturated (N == M AND ghost_change > picard_tol).
  • Σ N: 8 — the warm-start payoff signal.

A red Picard badge with elevated Σ N is the textbook signal to raise adaptive_picard_max.

Decision tree — when to touch each knob

"My turbo case has noticeably different mdot vs the explicit run"

The Phase I-X frozen-Dirichlet split has O(Δt) error on Compressor/Turbine ports. This is exactly what Phase I-Y was introduced to remove. Confirm with a CFL ramp:

  • picard_max_iter = 1 (Phase I-X equivalent).
  • picard_max_iter = 3 (Phase I-Y default).
  • picard_max_iter = 5, picard_tol = 1e-8 (tight Picard).

The mdot error should drop monotonically. If picard_max_iter = 5 still differs from the explicit run substantially, something else is going on (CFL too high, BC mismatch).

"The Picard badge keeps showing red (saturated)"

The adaptive controller hit its cap without converging the ghost coupling. Two paths:

  1. Raise the cap: adaptive_picard_max = 12. Costs more Newton solves per step (~max × n_pipes × per-Newton).
  2. Loosen the tolerance: picard_tol = 1e-4. The implicit constraint converges to a looser fixed point — fine when the network's natural noise floor is itself ~1e-4 (near surge / choke).

Watch the Σ N badge: if it's also climbing, the work is real. If Σ N stays low while the Picard badge is red, you've got Newton instability inside the outer loop — usually a CFL or initial-condition issue, not a Picard one.

"Picard always exits at 2 iters; can I make it cheaper?"

Yes — picard_max_iter = 1 skips the convergence check and runs Phase I-X. For a plenum-only network this is exact (the ghosts don't depend on pipe interior). For Compressor/Turbine networks this re-introduces the O(Δt) splitting error.

Alternative: adaptive_picard = true lets the controller pick the budget. On a plenum-only network it walks down to 2 (the floor where used + 1 > current) and holds. On Compressor networks it ratchets to whatever the natural budget is.

"I want to reproduce a Phase I-X result for regression"

[case]
picard_max_iter = 1
adaptive_picard = false

This restores the exact Phase I-X path: one Newton solve per pipe per step, no ghost refresh between, no convergence check.

Bundled example walkthrough

Two bundled cases exercise the adaptive controller from different angles.

Example 54 — quiet ratcheting

examples/54_adaptive_picard_turbo.toml runs a turbine + compressor + boost-PI network at CFL = 2 with the adaptive controller on at picard_tol = 1e-6 (the production default). examples/plots/plot_54_adaptive_picard.py samples LiveControl per step and plots the controller's trajectory.

What you should see:

  1. Initial transient (~0-5 ms): controller saturates at the initial 4 iters, immediately grows to 6 (the cap). Ghost change is ~1e-3 (hot turbo coupling).
  2. Mid-transient (5-15 ms): controller ratchets back: 6 → 5 → 4 → 3 → 2.
  3. Steady state (15+ ms): managed budget settles at 2. Eventually (~35 ms) drops further to 1 as Newton converges in 0 iterations on each step (fully steady).

Example 57 — tuning sweep with stiffness step

examples/57_picard_stiffness_step.toml is the canonical tuning recipe: same turbo network as example 54, but picard_tol = 1e-9 (tight) and the exhaust inlet is stepped mid-run (3.0 bar → 9.0 bar at t = 25 ms, over a 0.1 ms ramp). The plot shows the full controller life cycle and provides two explicit tuning signals.

What you should see in plot_57_picard_stiffness_step.py:

  1. Cold-start (t < 15 ms) — managed budget saturates at adaptive_picard_max = 8 for ~10 ms. Tuning signal #1 — raise the cap. A network whose initial transient pegs the cap for many ms in a row needs more headroom.
  2. Mid-run settle (t = 15-25 ms) — controller walks down to a steady-state floor (4 here). If the floor in your network sits exactly at adaptive_picard_min, lower the floor for less Newton work per step at steady state.
  3. Stiffness step (t ≈ 25-30 ms) — picard_iters_used bumps up briefly (4 → 5), ghost-change residual spikes, controller raises the managed budget. Tuning signal #2 — confirms the controller responds to load disturbances. If your equivalent stiffness event saturates the cap, raise it again.
  4. Re-settle (t > 30 ms) — controller walks back down to the new operating point's floor.

Run either via the CLI to watch the badges live:

aeolus1d run examples/54_adaptive_picard_turbo.toml \
    --info-every-n 200
aeolus1d run examples/57_picard_stiffness_step.toml \
    --info-every-n 200

Failure diagnostics — what the Newton-failure message tells you

When Newton inside a Picard iter fails (after exhausting config.dt_halvings_max), the diagnostic message includes a Picard-state line:

Solver state at failure:
  rk = 'implicit_euler', scheme = 'first_order', cfl = 2
  adaptive_dt = False
  Newton: tol = 1e-08, max_iters = 8, dt_halvings_max = 5
  last Newton iteration count before each halving: 8 (max = 8)
  Picard (Phase I-Y): on iteration 3 of 4 (tol = 1e-06) when Newton failed

This tells you whether the failure was in the FIRST Picard iter (inner step / CFL is the problem) or on a LATER iter (outer coupling didn't converge). Combined with the suggested remedies, particularly the new

[F] **Raise picard_max_iter.** Newton failed on the last
   Picard iteration of a 4-iter budget — the outer loop never
   got to tighten the ghost coupling further.

you usually have enough to pick the right knob without re-running with a debugger.

QS-driver Picard (sign-bracket-bisection guard)

The quasi-steady driver (solver = "quasi_steady") runs its own junction-pressure Picard, separate from the implicit-network loop above: each outer iter evaluates the steady per-pipe SIMPLE solve at the current junction pressures jp, reads the net mass imbalance net(jp) at every junction, refreshes a secant slope from the previous-iter (jp, net) pair, and takes a damped-Newton step

jp_new[jj] = jp[jj] − qs_picard_relax · net[jj] / slope[jj]

with qs_picard_relax = 0.7 (default, src/aeolus1d/io/schema.py:247). The per-junction slope is refreshed from the previous-iter secant pair when |Δjp| > dp_fd and the pair is well-conditioned (negative, |slope| > 1e-20); when |Δjp| ≤ dp_fd the prior slope is simply kept. Whenever the slope is invalid — the first iter (secant zero-initialised), a sign flip / |slope| → 0, or after a bracket-bisection reset — it is (re)seeded by a junction-local FD perturbation of size dp_fd = max(|jp|·1e-6, 1e-3 Pa). Knobs live on [case]:

[case]
qs_picard_max_iter = 100      # default
qs_picard_tol      = 1.0e-6   # relative ṁ imbalance
qs_picard_relax    = 0.7      # damped-Newton relax factor

The stale-secant 2-cycle and the guard

On AVL Cruise-M's MFR_transition_SC630 case the loop was trapped in a stable damped-Newton 2-cycle past the inlet ṁ drop at t ∈ [10, 11] s: net[jj] flipped sign every iter, jp[jj] settled into exactly two alternating values, and the secant slope could not refresh because |Δjp| sat below the FD-refresh gate. The mechanism — a stale-slope underestimate of ~2.86× putting the contraction at λ ≈ -1 (perfectly reflecting) — is documented in the diagnostic-script docstrings linked below. The original workaround (MR !130) exposed qs_picard_relax so users could break the cycle by dropping it to 0.3 (or 0.1 on very stiff transients).

MR !131 ships the algorithmic fix: an iter-local sign-bracket guard built from the previous-iter and current-iter (jp[jj], net[jj]) pair. When those two nets have opposite signs the iteration has straddled the root and the two iterates form a finite bracket spanning one Newton step of off-diagonal drift (small by construction; exact for SC630's single-junction cycle). The guard then fires whenever the damped-Newton step would either

  • land outside the known sign bracket (Newton wandering past the root), or
  • be below the FD-refresh gate while the bracket is still wider than the step (the λ ≈ -1 2-cycle lock — the secant cannot update because |Δjp| ≤ dp_fd, and the two iterates ARE the bracket),

in which case the Newton step is replaced by the bisection midpoint 0.5 · (lo + hi) and slope_valid[jj] is cleared so the secant re-bootstraps from the bracket neighbourhood on the next iter. On smooth convergence consecutive nets share the same sign so the bracket never forms and the guard is silent — bit-equivalent to the legacy path in that regime. Code reference: src/aeolus1d/solver/quasi_steady_network.py (_picard_jp body, the comment block headed by "Per-junction sign-bracket guard (iter-local)").

The cycle's two iterates ARE the bracket and the root is the midpoint by construction, so the fix converges in O(1) extra evals on the trapped case. The iter-local scope keeps the bracket coherent under off-diagonal multi-junction coupling (the bracket spans one Newton step of drift, not many).

qs_picard_relax — still useful, no longer load-bearing

The knob from MR !130 stays on the dataclass as belt-and-braces. With the sign-bracket guard in place the default 0.7 no longer gets trapped on the SC630 transient (the bracket forms on iter 2 and bisection converges immediately). Per-case lowering to 0.3 remains documented in the schema comment for users facing a "failed to converge" with a stuck relative imbalance on networks the guard's bracket somehow misses; in practice the sign-bracket path makes the relax knob a diagnostic-only lever.

Permanent diagnostics — when convergence still looks weird

Two scripts ship under scripts/ and are the durable on-ramp for any future QS-Picard convergence story:

  • scripts/diagnose_qs_picard_2cycle.py — drives the AVL Cruise-M MFR_transition_SC630 case past the trapping transient at multiple qs_picard_relax values and writes a JSONL trace of every Picard iter via the env-gated hook in advance_network_quasi_steady (search the source for _qs_dbg_writer). Outputs qs_picard_trace_relax{0.7,0.3}.jsonl under --out-dir (default /tmp). Prints one line per variant with the path, the RuntimeError (if Picard failed to converge), and the wall-clock time.
  • scripts/analyze_qs_picard_trace.py — reads the JSONL traces and asserts the 2-cycle signature on the trailing --last-k iters of each macro step: net[jj] sign-flip every iter, exactly two alternating jp[jj] values, and the reflection ratio λ = 1 − relax · S_true / σ near -1. Use it to confirm whether a "stuck Picard" on a new case is the same regime-switch 2-cycle or a different failure mode.

Usage::

python scripts/diagnose_qs_picard_2cycle.py \
    --case local_cases/avl_cruise_m/aeolus_optimized.toml \
    --t-end 11.0 --out-dir /tmp
python scripts/analyze_qs_picard_trace.py \
    /tmp/qs_picard_trace_relax0.7.jsonl \
    /tmp/qs_picard_trace_relax0.3.jsonl --last-k 6

If the 2-cycle signature appears AND the sign-bracket guard fails to clear it, capture the trace and file an issue: the bracket-coherence assumption (one Newton step of off-diagonal drift) has hit a multi-junction topology that breaks it.

  • GUI run-settings strip and Inspector form — see the GUI page. The QS-driver knobs (qs_picard_max_iter, qs_picard_tol, qs_picard_relax) live on CaseHeader and sit under the Inspector's "Advanced settings…" expander introduced in MR !132 / !135 / !136 (CaseHeader's 10 BASIC_FIELDS are name, t_end, cfl, scheme, limiter, rk, solver, adaptive_dt, dt_fixed, chem_solver; everything else, including the QS-Picard knobs, is collapsed by default).
  • CLI run controls — see CLI reference. The QS knobs are TOML-only at the [case] level; no aeolus1d run flag overrides them.
  • Phase I-Y (2bce33e) — the outer loop itself.
  • Phase I-Y-diag (afff9ef) — LiveControl surfaces.
  • Phase I-Y-toml (a112d6b) — [case] schema knobs.
  • Phase I-Y-cli (3d4209c) — CLI picard = N token.
  • Phase I-Y-warm (30f776a) — Newton warm-start across iters.
  • Phase I-Y-cache (8e7e468) — BC chain cache hoisted out of the loop.
  • Phase I-Y-adaptive (e6b9552) — controller + 4 TOML knobs.
  • Phase I-Y-adaptive-demo (f950e98) — bundled example 54.
  • Phase I-Y-cli2 (0f54685) — runner-dock Picard badge.
  • Phase I-Y-cli3 (132eca0) — runner-dock Σ Newton badge.
  • MR !130 — qs_picard_relax knob on CaseHeader (initial workaround for the QS-driver stale-secant 2-cycle).
  • MR !131 — sign-bracket-bisection guard inside _picard_jp (algorithmic fix; closes the workaround story). Permanent diagnostics: scripts/diagnose_qs_picard_2cycle.py + scripts/analyze_qs_picard_trace.py.