Skip to content

Pipe time step

One pipe's explicit step: the CFL time step, the residual dU/dt (primitive recovery, reconstruction, Riemann flux, source terms) and the Runge-Kutta stages that combine it. The single-pipe driver advance_single_pipe is drawn here; the network driver builds every pipe's residual through the same make_pipe_residual_fn and sizes its step with the same pipe_min_dt_dispatch (see Explicit network step). The contracts of every function below are in their docstrings (the standard is Function contracts and flow charts).

The step loop

What enters: a built Pipe with its initial state, t_end, the CFL, the scheme and limiter, the RK key, and optionally BCs, dt_fixed, a CFL ramp and live control. What leaves: pipe.U advanced in place and the final time.

flowchart TD
    IN["pipe, t_end, cfl, scheme, limiter, rk, BCs, dt_fixed"] --> GH{"_single_pipe_ghosts"}
    GH -- "one of left_bc / right_bc" --> X1["ValueError"]
    GH --> RF["make_pipe_residual_fn (second chart)"]
    RF --> IM{"rk implicit?"}
    IM -- yes --> IMP{"_advance_single_pipe_implicit"}
    IMP -- "CFL ramp, or scheme mp5" --> X2["NotImplementedError"]
    IMP --> DRV["implicit_single_pipe_driver: implicit Euler, BDF2 or Radau IIA"]
    IM -- no --> RK{"ssp_rk2, ssp_rk3, rk4 or rk5?"}
    RK -- no --> X3["ValueError: unknown rk"]
    RK --> T{"t < t_end?"}
    T -- no --> OUT["t_end reached: pipe.U advanced, clock.t = t"]
    T -- "max_steps spent short of t_end" --> X4["RuntimeError: max_steps exhausted"]
    T -- yes --> LC{"poll_live_control, then the CFL ramp"}
    LC -- stop --> OUT
    LC --> G["ghosts: BCs fill the Ng ghost rows"]
    G --> DT{"_single_pipe_dt: CFL dt (fourth chart), dt_fixed overrides it, clipped to t_end"}
    DT -- "dead cell, or non-finite / 1e300 dt, even under dt_fixed" --> X5["RuntimeError naming the pipe and cells"]
    DT -- "dt <= 0" --> X6["assert_advancing_dt: NonAdvancingStepError"]
    DT --> ST{"_explicit_rk_step"}
    ST -- "ssp_rk2 / ssp_rk3" --> SSP["Shu-Osher stages; ghosts refilled at t + dt (and t + dt/2)"]
    ST -- "rk4 / rk5" --> BT["_butcher_step: stage i at t + c_i dt"]
    SSP --> NX["t += dt; on_step(t)"]
    BT --> NX
    NX --> T

Every stage's residual goes through _diagnosed_residual: a ZeroDivisionError from a kernel (a cell with rho == 0, or p == 0 on a diverging face) becomes a RuntimeError naming the pipe and its bad cells.

Choosing the residual kernel

make_pipe_residual_fn walks the scheme × EOS × geometry × γ tree once and returns a closure bound to one @njit kernel. The kernels hard-code their ghost width, so a wider pad is cut to the inner window (slice_pad_to_ng, slice_pad_to_ng1) rather than handed over whole.

flowchart TD
    P["pipe, scheme, limiter"] --> NG{"check_scheme_ng: Ng >= 1 / 2 / 3?"}
    NG -- no --> X1["ValueError: scheme needs Ng >= n"]
    NG --> S{"scheme"}
    S -- other --> X2["ValueError: unknown scheme"]
    S -- mp5 --> M5{"_residual_fn_mp5"}
    M5 -- "real gas" --> M5E["_residual_fn_mp5_eos: EOS MP5, quasi-1D if needs_quasi1d"]
    M5 -- ideal --> M5I["MP5, quasi-1D and / or per-cell gamma kernel"]
    S -- muscl --> L{"Limiter(code) in 0, 1, 2?"}
    L -- no --> X3["ValueError: unknown limiter"]
    L --> VG{"_residual_fn_muscl_var_gamma: ideal, per-cell gamma, gamma_pad?"}
    VG -- yes --> MVG["var-gamma MUSCL, quasi-1D if needs_quasi1d"]
    VG -- no --> MB{"_residual_fn_muscl"}
    MB -- "real gas" --> MEOS["_residual_fn_muscl_eos: EOS MUSCL on (rho, u, T)"]
    MB -- "ideal, scalar gamma" --> MPL["MUSCL, quasi-1D if needs_quasi1d"]
    MB -- "per-cell gamma, no gamma_pad" --> FO
    S -- first_order --> FO{"_residual_fn_first_order"}
    FO -- "real gas" --> FOE["_residual_fn_first_order_eos"]
    FO -- "needs_quasi1d" --> FOQ["quasi-1D first order"]
    FO -- "hllc_low_mach" --> FOL["pipe_residual_first_order_lm"]
    FO -- "per-cell gamma" --> FOV["first order, per-cell gamma"]
    FO -- else --> FOS["pipe_residual_first_order"]
    M5E --> W["wrap_with_py_overlay: add py_friction / py_heat / py_porous (_overlay_py_sources)"]
    M5I --> W
    MVG --> W
    MEOS --> W
    MPL --> W
    FOE --> W
    FOQ --> W
    FOL --> W
    FOV --> W
    FOS --> W
    W --> OUT["no-arg closure: dU/dt, shape (N, 3)"]

Inside one residual kernel

The same four steps in every kernel; the scheme decides the second, the gas model the first and third, the geometry the last.

flowchart TD
    U["U_pad: the scheme's ghost width, ghosts filled"] --> W{"primitives"}
    W -- "ideal" --> W1["U_to_W / U_to_W_var_gamma: (rho, u, p)"]
    W -- "real gas" --> W2["_W_eos_from_U_jit: (rho, u, T)"]
    W -- "rho == 0 in any row" --> XZ["ZeroDivisionError, named by the driver"]
    W1 --> R{"reconstruction"}
    W2 --> R
    R -- "first order" --> R0["cell states as they are"]
    R -- muscl --> R1["muscl_slopes with limited_slope (minmod, van_leer, central)"]
    R -- mp5 --> R2["_mp5_face: 5th-order interpolant, median-limited"]
    R1 --> F{"face states"}
    R2 --> F
    F -- ideal --> F1["W_to_U_single: rho, p floored at 1e-12"]
    F -- "real gas" --> F2["_U_eos_from_W_jit: rho, T floored at 1e-12"]
    R0 --> H{"Riemann flux"}
    F1 --> H
    F2 --> H
    H --> H1["hllc_flux / hllc_flux_var_gamma / hllc_flux_eos_jit / hllc_flux_lm"]
    H1 --> HC["_hllc_core_jit: Davis speeds, den floored at 1e-30"]
    HC -- "a == 0 on a diverging face" --> XZ
    HC --> D{"divergence"}
    D -- "uniform area" --> D1["dU/dt = -(F_R - F_L) / dx"]
    D -- "needs_quasi1d" --> D2["-(A_R F_R - A_L F_L)/(A_c dx) + p dA/(A_c dx) + friction, heat, porous"]
    D1 --> OUT["R, shape (N, 3)"]
    D2 --> OUT

The CFL time step

pipe_min_dt_dispatch picks the sound speed the pipe's own flux uses. A cell with rho <= 0 or p <= 0 is not skipped: its sound speed is NaN, and NaN loses every comparison, so a skipped cell would leave the step at the healthy cells' value (ADR-0011). The kernels return the finite 1e300 sentinel instead, and the Python layer turns it into an error.

flowchart TD
    P["pipe, cfl"] --> K{"gas model"}
    K -- "real gas" --> E["pipe_min_dt_eos_jit: |u| + a from the EOS"]
    K -- "per-cell gamma" --> V["pipe_min_dt_checked with pipe_min_dt_var_gamma"]
    K -- "scalar gamma" --> S["pipe_min_dt_checked with pipe_min_dt"]
    E --> Q{"any cell rho <= 0 or p <= 0, N == 0, or max wave speed 0 / NaN / inf?"}
    V --> Q
    S --> Q
    Q -- no --> DT["dt = cfl dx / max(|u| + a)"]
    Q -- yes --> SEN["1e300 sentinel"]
    SEN --> D{"diagnose"}
    D -- "all NaN" --> X1["nonfinite_dt_error"]
    D -- "rho <= 0 or NaN" --> X2["nonpositive_density_error"]
    D -- "rho > 0, p <= 0" --> X3["nonpositive_pressure_error; on a real-gas pipe a RuntimeError naming the cells by EOS pressure"]
    D -- else --> X1

The network driver's packed first-order pipes take the same test in _packed_min_dt, one sentinel per pipe, so one dead pipe cannot hide behind the minimum over the others.

Functions

  • src/aeolus1d/solver/timeloop.py::advance_single_pipe
  • src/aeolus1d/solver/timeloop.py::_single_pipe_ghosts
  • src/aeolus1d/solver/timeloop.py::_advance_single_pipe_implicit
  • src/aeolus1d/solver/timeloop.py::implicit_single_pipe_driver
  • src/aeolus1d/solver/solver_ops.py::poll_live_control@3
  • src/aeolus1d/solver/timeloop.py::_single_pipe_dt
  • src/aeolus1d/solver/solver_ops.py::assert_advancing_dt
  • src/aeolus1d/solver/timeloop.py::_explicit_rk_step
  • src/aeolus1d/solver/timeloop.py::_butcher_step
  • src/aeolus1d/solver/timeloop.py::_diagnosed_residual
  • src/aeolus1d/solver/timeloop.py::make_pipe_residual_fn
  • src/aeolus1d/solver/timeloop.py::check_scheme_ng
  • src/aeolus1d/solver/timeloop.py::slice_pad_to_ng
  • src/aeolus1d/solver/timeloop.py::slice_pad_to_ng1
  • src/aeolus1d/solver/timeloop.py::_residual_fn_mp5
  • src/aeolus1d/solver/timeloop.py::_residual_fn_mp5_eos
  • src/aeolus1d/solver/timeloop.py::_residual_fn_muscl_var_gamma
  • src/aeolus1d/solver/timeloop.py::_residual_fn_muscl
  • src/aeolus1d/solver/timeloop.py::_residual_fn_muscl_eos
  • src/aeolus1d/solver/timeloop.py::_residual_fn_first_order
  • src/aeolus1d/solver/timeloop.py::_residual_fn_first_order_eos
  • src/aeolus1d/solver/timeloop.py::wrap_with_py_overlay
  • src/aeolus1d/solver/timeloop.py::_overlay_py_sources
  • src/aeolus1d/solver/timeloop.py::closure_args
  • src/aeolus1d/solver/reconstruct.py::U_to_W
  • src/aeolus1d/solver/reconstruct.py::U_to_W_var_gamma
  • src/aeolus1d/solver/hllc_eos_njit.py::_W_eos_from_U_jit
  • src/aeolus1d/solver/reconstruct.py::muscl_slopes
  • src/aeolus1d/solver/reconstruct.py::limited_slope
  • src/aeolus1d/solver/kernels_mp5.py::_mp5_face
  • src/aeolus1d/solver/reconstruct.py::W_to_U_single
  • src/aeolus1d/solver/hllc_eos_njit.py::_U_eos_from_W_jit
  • src/aeolus1d/solver/hllc.py::hllc_flux
  • src/aeolus1d/solver/hllc.py::hllc_flux_var_gamma
  • src/aeolus1d/solver/hllc_eos_njit.py::hllc_flux_eos_jit
  • src/aeolus1d/solver/hllc_lm.py::hllc_flux_lm
  • src/aeolus1d/solver/hllc.py::_hllc_core_jit
  • src/aeolus1d/solver/kernels.py::pipe_residual_first_order
  • src/aeolus1d/solver/kernels.py::pipe_residual_first_order_var_gamma
  • src/aeolus1d/solver/kernels.py::pipe_residual_quasi1d_first_order
  • src/aeolus1d/solver/kernels.py::pipe_residual_muscl
  • src/aeolus1d/solver/kernels.py::pipe_residual_quasi1d_muscl
  • src/aeolus1d/solver/kernels_mp5.py::pipe_residual_mp5
  • src/aeolus1d/solver/kernels_mp5.py::pipe_residual_quasi1d_mp5
  • src/aeolus1d/solver/hllc_lm.py::pipe_residual_first_order_lm
  • src/aeolus1d/solver/timeloop.py::pipe_min_dt_dispatch
  • src/aeolus1d/solver/hllc_eos_njit.py::pipe_min_dt_eos_jit
  • src/aeolus1d/solver/kernels.py::pipe_min_dt_checked
  • src/aeolus1d/solver/kernels.py::pipe_min_dt
  • src/aeolus1d/solver/kernels.py::pipe_min_dt_var_gamma
  • src/aeolus1d/solver/kernels.py::nonfinite_dt_error
  • src/aeolus1d/solver/kernels.py::nonpositive_density_error
  • src/aeolus1d/solver/kernels.py::nonpositive_pressure_error
  • src/aeolus1d/solver/packed_pipes.py::_packed_min_dt