Skip to content

Implicit network step

The implicit network driver (advance_network_implicit, rk = implicit_euler, bdf2 or radau_iia): every pipe is advanced by a Newton-converged implicit step with its port ghosts held fixed, inside a Picard loop over the pipe-to-component coupling; the components, closures, species and chemistry then follow as explicit first-order splits. The contracts of every function below are in their docstrings (the standard is Function contracts and flow charts).

One run

What enters: the built pipes, junctions and far boundary conditions, t_end, the CFL, rk / scheme, the Newton, adaptive-Δt and adaptive-Picard configurations, and the in-Newton / split flags. The junctions' neck pipes are advanced with the user pipes.

flowchart TD
    IN["pipes, junctions, far BCs, t_end, cfl, rk, scheme, configs, flags"] --> R1{"chem_every_n >= 1; _refuse_orifice_cylinders"}
    R1 -- "chem_every_n < 1, an orifice cylinder" --> X1["ValueError"]
    R1 --> F["_stamp_newton_flags on the user pipes; _warn_if_no_substrate_wall"]
    F --> R2{"_check_implicit_rk_and_scheme; _check_phase_i_constraints per pipe"}
    R2 -- "unknown rk / scheme, unsupported pipe" --> X2["ValueError / NotImplementedError"]
    R2 --> ST{"_ImplicitRun.start: ends classified, ghost closures, ops, Picard budget >= 1"}
    ST -- "adaptive_dt with dt_fixed" --> X3["ValueError"]
    ST -- "a port on a pipe outside the network" --> X4["NotImplementedError naming both"]
    ST --> T{"t < t_end?"}
    T -- no --> OUT["t reached; network advanced in place; clock.t = t"]
    T -- "budget spent short of t_end" --> X5["RuntimeError"]
    T -- yes --> LC{"poll_live_control (CFL and adaptive_cfl_max overrides)"}
    LC -- stop --> FL["flush_chemistry: pending super-cycle window reacted"]
    FL --> OUT
    LC --> DT{"refresh gamma; step_dt: ghosts filled, dt fixed / adaptive / CFL, HRR cap"}
    DT -- "non-finite pipe step, dt <= 0" --> X6["RuntimeError / NonAdvancingStepError"]
    DT --> PS["picard_solve (second chart)"]
    PS -- "a pipe fails after every halving" --> X7["RuntimeError: the Newton-failure diagnostic"]
    PS --> SPL["advance_junctions; catalyst_split; dpf_split; substrate_wall_step"]
    SPL --> SPC["species_transport: rebase Y onto the new rho, then sub-steps or the implicit-upwind march"]
    SPC --> CH["chemistry_supercycle; urea spray; advance_thermal_couplings"]
    CH --> UC["update_controllers (adaptive dt, adaptive Picard, dt_prev); HRR cap; t += dt"]
    UC --> T

The Picard loop

picard_solve iterates pipe Newton solves against component ghosts that depend on the pipe interiors. A Newton failure halves Δt for the rest of the step; a pass that halved never counts as converged, and a halving on the last budgeted pass extends the budget, up to a cap.

flowchart TD
    A["dt, start-of-step U_n per pipe"] --> SEED["_seed_jac_chain_cache per pipe (once per step)"]
    SEED --> NP["_newton_pass"]
    NP --> PAR{"parallel_pipes?"}
    PAR -- yes --> TP["_try_parallel_first_pass"]
    TP -- "a pipe failed" --> RS["every pipe restored to U_n"]
    TP -- "all converged" --> GH
    PAR -- no --> SP
    RS --> SP["_solve_pipe, one pipe at a time"]
    SP --> NW{"_newton_step (IE, variable-step BDF2 via _bdf2_var_step_coeffs) or _newton_step_radau"}
    NW -- "not converged" --> HV["pipe restored to U_n; dt halved"]
    HV --> NW
    HV -- "halvings exhausted" --> X["_format_newton_failure_diagnostic -> RuntimeError"]
    NW -- converged --> GH["_fill_component_ghosts and far BCs from the new interiors"]
    GH --> C{"ghost change (_ghost_snapshot_max_rel_change) < picard_tol and no halving?"}
    C -- yes --> DONE["commit dt"]
    C -- "no, budget left (extended by one after a last-pass halving)" --> NP
    C -- "budget spent" --> DONE
    DONE -- "the committed pass halved (only at the cap)" --> W["RuntimeWarning: pipes ran at different dt"]

Functions

  • src/aeolus1d/solver/implicit_network.py::advance_network_implicit
  • src/aeolus1d/solver/implicit_network.py::_refuse_orifice_cylinders
  • src/aeolus1d/solver/implicit_network.py::_stamp_newton_flags
  • src/aeolus1d/solver/implicit_network.py::_warn_if_no_substrate_wall
  • src/aeolus1d/solver/solver_ops.py::network_pipes
  • src/aeolus1d/solver/implicit_network.py::_check_implicit_rk_and_scheme
  • src/aeolus1d/solver/implicit_network.py::_check_network_constraints
  • src/aeolus1d/solver/implicit_kernels.py::_check_phase_i_constraints
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.start
  • src/aeolus1d/solver/implicit_network.py::_classify_pipe_ends
  • src/aeolus1d/solver/implicit_network.py::_port_pipe
  • src/aeolus1d/solver/implicit_network.py::_make_network_apply_ghosts
  • src/aeolus1d/solver/implicit_network.py::_ImplicitOps.resolve
  • src/aeolus1d/solver/implicit_network.py::_get_pipe_pool
  • src/aeolus1d/solver/solver_ops.py::poll_live_control@3
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.flush_chemistry
  • src/aeolus1d/solver/solver_ops.py::refresh_network_gamma
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.step_dt
  • src/aeolus1d/solver/timeloop.py::pipe_min_dt_dispatch
  • src/aeolus1d/solver/solver_ops.py::assert_advancing_dt
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.picard_solve
  • src/aeolus1d/solver/implicit_kernels.py::_seed_jac_chain_cache
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun._newton_pass
  • src/aeolus1d/solver/implicit_network.py::_try_parallel_first_pass
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun._solve_pipe
  • src/aeolus1d/solver/implicit_kernels.py::_newton_step
  • src/aeolus1d/solver/implicit_kernels.py::_bdf2_var_step_coeffs
  • src/aeolus1d/solver/implicit_radau.py::_newton_step_radau
  • src/aeolus1d/solver/implicit_network.py::_format_newton_failure_diagnostic
  • src/aeolus1d/solver/implicit_network.py::_fill_component_ghosts
  • src/aeolus1d/solver/implicit_network.py::_capture_ghost_snapshot
  • src/aeolus1d/solver/implicit_network.py::_ghost_snapshot_max_rel_change
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.report_picard
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.advance_junctions
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun._advance_one_junction
  • src/aeolus1d/solver/implicit_network.py::_advance_junction_volume
  • src/aeolus1d/solver/network_timeloop.py::advance_shaft_housing
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.catalyst_split
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.dpf_split
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun._dpf_regen_step
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.substrate_wall_step
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.species_transport
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun._species_substeps
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun._species_implicit_upwind
  • src/aeolus1d/solver/implicit_network.py::_publish_conserved_face_flux
  • src/aeolus1d/solver/quasi_steady_network.py::_qs_implicit_species
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.chemistry_supercycle
  • src/aeolus1d/solver/chemistry_ops.py::chemistry_step_on_pipes
  • src/aeolus1d/solver/solver_ops.py::apply_urea_spray_to_pipe
  • src/aeolus1d/solver/solver_ops.py::advance_thermal_couplings
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.update_controllers
  • src/aeolus1d/solver/implicit_kernels.py::_adaptive_dt_update
  • src/aeolus1d/solver/implicit_kernels.py::_adaptive_picard_update
  • src/aeolus1d/solver/solver_ops.py::hrr_next_dt_limit

Notes from the kernel contracts

Found writing the contracts of the kernels this chart calls (2026-09-29):

  • The Jacobian chain cache (_seed_jac_chain_cache) is seeded once per step, before the Picard loop. The numerical chain is frozen at the start-of-step state for every Picard and Newton iteration, and once seeded it takes precedence over bc_kind / apply_ghosts until the next seeding.
  • When a Newton step fails, pipe.U holds the last accepted iterate; the caller restores its own base state.
  • The Picard budget is clamped at the start of a step to min(picard_max, max(picard_min, x)). An enabled adaptive band with min < 1 or min > max is refused when the case is built (AdaptivePicardConfig).
  • Refusals inside the kernels: _check_phase_i_constraints and the ValueErrors in _split_bc_kind.