pipe-system-simulation-test/docs/superpowers/plans/2026-04-14-0d-1d-tank-pipe-blowdown-mvp.md

0D–1D Tank-Pipe Blowdown MVP Implementation Plan

For agentic workers: REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (- [ ]) syntax for tracking.

Goal: Build a transient simulation of a high-pressure tank blowing down through a 1 m pipe (D=5 mm) into a low-pressure tank, using a 0D lumped-parameter model for each tank coupled to a 1D finite-volume Euler solver on the pipe via ghost-cell + HLL Riemann fluxes.

Architecture: 7 Python source modules (config, riemann, tank, pipe, solver, output, main) under src/, with tests/ for pytest-based unit + integration tests. Test-driven workflow: each module is built by writing failing tests first, then the minimal implementation that passes them, then committing. The solver enforces mass/energy conservation via "flux doubling" — the same boundary flux vector updates both the pipe boundary cell and the adjacent tank.

Tech Stack: Python 3, NumPy, Matplotlib (with PillowWriter for GIF animation, avoiding ffmpeg dependency), pytest.

Spec reference: docs/superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md

---

File Layout (Target State)

pipe-system-simulation-test/
├── src/
│   ├── config.py       # constants (gas props, geometry, ICs, numerics)
│   ├── riemann.py      # hll_flux(W_L, W_R, gamma) pure function
│   ├── tank.py         # Tank class (0D, state = mass + total internal energy U)
│   ├── pipe.py         # Pipe class (1D, conservative W field, step())
│   ├── solver.py       # run(tank1, tank2, pipe, t_end, cfl, ...) time loop
│   ├── output.py       # save_history / plot_* / make_pipe_animation
│   └── main.py         # entry: assemble → run → sanity check → save
├── tests/
│   ├── test_riemann.py         # 3 tests
│   ├── test_tank.py            # 3 tests
│   ├── test_pipe.py            # 2 tests
│   └── test_integration.py     # 2 tests (total mass / energy conservation)
├── conftest.py         # sys.path hack so tests/ can import from src/
├── results/            # runtime outputs (.npz, .png, .gif) — created at runtime
└── docs/superpowers/
    ├── specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md
    └── plans/2026-04-14-0d-1d-tank-pipe-blowdown-mvp.md  (this file)

The src/README.md, tests/, results/README.md, cases/README.md, scripts/README.md existing scaffolding files are untouched; only src/README.md and friends are left as-is.

---

Pre-flight: Dependencies

Verify Python 3 + required packages before starting any task. Run this once at the beginning.

python3 -c "import numpy, matplotlib, pytest; print(numpy.__version__, matplotlib.__version__, pytest.__version__)"

Expected: three version strings printed with no import error.

If matplotlib or pytest is missing:

pip install --user numpy matplotlib pytest

---

Task 1: Project Skeleton

Create the directory structure, conftest.py path shim, and results/ directory. No Python code yet — just scaffolding so subsequent tasks can import cleanly.

Files:

• Create: conftest.py

• Create: tests/ (directory only)

• Create: results/ (directory only; add a .gitkeep so git tracks the folder)

• Step 1: Create conftest.py at project root

Write this file exactly (path shim so tests/test_*.py can do from riemann import ...):

# conftest.py
import os
import sys

_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.join(_HERE, "src"))

• Step 2: Create empty tests/ directory

mkdir -p tests

• Step 3: Create results/ directory with .gitkeep

mkdir -p results
touch results/.gitkeep

• Step 4: Verify pytest can discover the (empty) tests directory

pytest tests/ -v

Expected output: no tests ran in X.XXs (exit code 5 is fine — it means "no tests collected", which is what we want at this stage).

• Step 5: Commit

git add conftest.py tests results/.gitkeep
git commit -m "scaffold: add conftest, tests/, results/ for MVP"

---

Task 2: src/config.py — Constants & Scenario

A pure data module with all physical/numerical constants. Includes assert statements for parameter validation (per spec §6.1). No functions, no classes, no side effects.

Files:

• Create: src/config.py

• Step 1: Write src/config.py

# src/config.py
"""
Physical, geometric, and numerical constants for the 0D-1D tank-pipe
blowdown MVP. Pure data module — no functions, no side effects.
"""

# ---------- Gas properties (ideal air-like) ----------
GAMMA = 1.4
R_GAS = 287.0       # J / (kg K)

# ---------- High-pressure tank (upstream, Tank 1) ----------
V1 = 5.0            # m^3
P1_INIT = 10e6      # Pa  (10 MPa)
T1_INIT = 300.0     # K

# ---------- Low-pressure tank (downstream, Tank 2) ----------
V2 = 10.0           # m^3
P2_INIT = 101325.0  # Pa  (1 atm)
T2_INIT = 300.0     # K

# ---------- Pipe geometry ----------
L = 1.0             # m
D = 5e-3            # m  (5 mm)
N_CELLS = 20        # number of finite-volume cells

# ---------- Simulation control ----------
T_END = 0.1         # s
CFL = 0.5

# ---------- Output & animation ----------
ANIMATION_STRIDE = 10     # keep every Nth frame in the GIF
OUTPUT_DIR = "results"

# ---------- Parameter validation (per spec §6.1) ----------
assert GAMMA > 1, "GAMMA must be > 1"
assert R_GAS > 0, "R_GAS must be > 0"
assert V1 > 0 and V2 > 0, "tank volumes must be > 0"
assert L > 0, "L must be > 0"
assert D > 0, "D must be > 0"
assert N_CELLS >= 2, "N_CELLS must be >= 2"
assert P1_INIT > 0 and P2_INIT > 0, "initial pressures must be > 0"
assert T1_INIT > 0 and T2_INIT > 0, "initial temperatures must be > 0"
assert 0 < CFL <= 1, "CFL must be in (0, 1]"
assert T_END > 0, "T_END must be > 0"
assert ANIMATION_STRIDE >= 1, "ANIMATION_STRIDE must be >= 1"

• Step 2: Verify the module imports cleanly and assertions pass

python3 -c "import sys; sys.path.insert(0,'src'); import config; print(config.GAMMA, config.N_CELLS, config.T_END)"

Expected output: 1.4 20 0.1

• Step 3: Commit

git add src/config.py
git commit -m "feat(config): add physical and numerical constants"

---

Task 3: src/riemann.py — HLL Flux (TDD)

HLL Riemann solver as a pure function. Written TDD: 3 tests first, then implementation.

Files:

• Create: tests/test_riemann.py

• Create: src/riemann.py

• Step 1: Write the failing tests

Create tests/test_riemann.py:

# tests/test_riemann.py
import numpy as np
from riemann import hll_flux


GAMMA = 1.4


def _to_conservative(rho, u, P, gamma):
    """(rho, u, P) -> [rho, rho*u, rho*E] where E = e + u^2/2."""
    return np.array([
        rho,
        rho * u,
        P / (gamma - 1) + 0.5 * rho * u ** 2,
    ])


def _physical_flux(W, gamma):
    """F(W) = [rho*u, rho*u^2 + P, u*(rho*E + P)]."""
    rho = W[0]
    u = W[1] / rho
    E = W[2]
    P = (gamma - 1) * (E - 0.5 * rho * u ** 2)
    return np.array([rho * u, rho * u ** 2 + P, u * (E + P)])


def test_hll_identical_states_returns_physical_flux():
    """
    When W_L == W_R, HLL must return the exact physical flux F(W)
    with zero numerical dissipation (the (W_R - W_L) term vanishes).
    """
    W = _to_conservative(rho=1.2, u=50.0, P=2.5e5, gamma=GAMMA)
    F = hll_flux(W, W, GAMMA)
    expected = _physical_flux(W, GAMMA)
    assert np.allclose(F, expected, rtol=1e-12), f"F={F}, expected={expected}"


def test_hll_equal_pressure_equal_energy_gives_exact_pressure_flux():
    """
    Two stationary states (u=0) with same P but different rho:
    - Both have the same energy density E = P/(gamma-1), so the HLL
      (W_R - W_L)[2] term vanishes -> exact zero energy flux.
    - F_L[1] = F_R[1] = P, and W_L[1] = W_R[1] = 0, so the momentum
      flux is exactly P.
    - The mass flux is NOT exactly zero for HLL (the density jump
      triggers the (W_R - W_L)[0] dissipation term) — this is a known
      HLL limitation for stationary contact discontinuities. We do
      not assert on F[0] here.
    """
    W_L = _to_conservative(rho=10.0, u=0.0, P=1e5, gamma=GAMMA)
    W_R = _to_conservative(rho=1.0,  u=0.0, P=1e5, gamma=GAMMA)
    F = hll_flux(W_L, W_R, GAMMA)
    assert abs(F[1] - 1e5) < 1e-6, f"momentum flux should equal P=1e5, got {F[1]}"
    assert abs(F[2]) < 1e-8, f"energy flux should be exactly 0, got {F[2]}"


def test_hll_sod_shock_tube_directional_fluxes_all_positive():
    """
    Classical Sod initial values:
      left:  (rho, u, P) = (1.0,   0, 1.0)
      right: (rho, u, P) = (0.125, 0, 0.1)
    The pressure gradient drives flow from left to right, so the
    HLL flux at the interface should have all three components
    strictly positive:
      F[0] > 0 : mass flux rightward
      F[1] > 0 : momentum flux rightward
      F[2] > 0 : energy flux rightward
    """
    W_L = _to_conservative(rho=1.0,   u=0.0, P=1.0, gamma=GAMMA)
    W_R = _to_conservative(rho=0.125, u=0.0, P=0.1, gamma=GAMMA)
    F = hll_flux(W_L, W_R, GAMMA)
    assert F[0] > 0, f"expected positive mass flux, got {F[0]}"
    assert F[1] > 0, f"expected positive momentum flux, got {F[1]}"
    assert F[2] > 0, f"expected positive energy flux, got {F[2]}"

• Step 2: Run tests to verify they fail

pytest tests/test_riemann.py -v

Expected: 3 failures with ModuleNotFoundError: No module named 'riemann'.

• Step 3: Write the minimal implementation

Create src/riemann.py:

# src/riemann.py
"""
HLL Riemann flux for the 1D compressible Euler equations.

Conservative variable vector: W = [rho, rho*u, rho*E]
  where E = e + u^2/2 is specific total energy,
        e = P / (rho * (gamma - 1)) is specific internal energy.

Physical flux: F(W) = [rho*u, rho*u^2 + P, u*(rho*E + P)]
"""

import numpy as np


def hll_flux(W_L, W_R, gamma):
    """
    Compute the HLL numerical flux at the interface between two states.

    Parameters
    ----------
    W_L, W_R : array-like of shape (3,)
        Left and right conservative state vectors.
    gamma : float
        Ratio of specific heats.

    Returns
    -------
    np.ndarray of shape (3,)
        HLL numerical flux vector.

    Raises
    ------
    ValueError
        If either state has non-positive density or pressure.
    """
    # --- Recover primitives from left state ---
    rho_L = W_L[0]
    if rho_L <= 0:
        raise ValueError(f"Non-positive density in W_L: rho={rho_L}, W_L={W_L}")
    u_L = W_L[1] / rho_L
    E_L = W_L[2]  # total energy density (= rho * E)
    p_L = (gamma - 1) * (E_L - 0.5 * rho_L * u_L ** 2)
    if p_L <= 0:
        raise ValueError(f"Non-positive pressure in W_L: p={p_L}, W_L={W_L}")
    a_L = np.sqrt(gamma * p_L / rho_L)

    # --- Recover primitives from right state ---
    rho_R = W_R[0]
    if rho_R <= 0:
        raise ValueError(f"Non-positive density in W_R: rho={rho_R}, W_R={W_R}")
    u_R = W_R[1] / rho_R
    E_R = W_R[2]
    p_R = (gamma - 1) * (E_R - 0.5 * rho_R * u_R ** 2)
    if p_R <= 0:
        raise ValueError(f"Non-positive pressure in W_R: p={p_R}, W_R={W_R}")
    a_R = np.sqrt(gamma * p_R / rho_R)

    # --- Physical fluxes at left/right states ---
    F_L = np.array([
        rho_L * u_L,
        rho_L * u_L ** 2 + p_L,
        u_L * (E_L + p_L),
    ])
    F_R = np.array([
        rho_R * u_R,
        rho_R * u_R ** 2 + p_R,
        u_R * (E_R + p_R),
    ])

    # --- HLL wave-speed estimates ---
    S_L = min(u_L - a_L, u_R - a_R)
    S_R = max(u_L + a_L, u_R + a_R)

    # --- HLL flux, piecewise on wave configuration ---
    if S_L >= 0:
        return F_L
    if S_R <= 0:
        return F_R
    W_L_arr = np.asarray(W_L, dtype=float)
    W_R_arr = np.asarray(W_R, dtype=float)
    return (S_R * F_L - S_L * F_R + S_L * S_R * (W_R_arr - W_L_arr)) / (S_R - S_L)

• Step 4: Run tests to verify they pass

pytest tests/test_riemann.py -v

Expected: 3 passed.

• Step 5: Commit

git add src/riemann.py tests/test_riemann.py
git commit -m "feat(riemann): HLL flux with 3 unit tests"

---

Task 4: src/tank.py — 0D Tank (TDD)

Tank class with (mass, U) as primary state and (rho, T, P) as derived properties. Written TDD: 3 tests first, then implementation.

Files:

• Create: tests/test_tank.py

• Create: src/tank.py

• Step 1: Write the failing tests

Create tests/test_tank.py:

# tests/test_tank.py
import numpy as np
import pytest
from tank import Tank


GAMMA = 1.4
R_GAS = 287.0


def test_tank_initial_state_matches_ideal_gas():
    """
    Given P, T, V construct a Tank; mass and U should match ideal gas:
      rho  = P / (R * T)
      mass = rho * V
      U    = P * V / (gamma - 1)     (since u=0 inside tank)
      T    = (U / mass) * (gamma - 1) / R   (round-trip)
    """
    tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS)
    rho_expected = 10e6 / (R_GAS * 300.0)
    assert abs(tank.rho - rho_expected) < 1e-9
    assert abs(tank.mass - rho_expected * 5.0) < 1e-6
    assert abs(tank.U - 10e6 * 5.0 / (GAMMA - 1)) < 1e-3
    assert abs(tank.P - 10e6) < 1e-3
    assert abs(tank.T - 300.0) < 1e-9


def test_tank_ghost_state_is_stagnation_conservative_vector():
    """
    ghost_state() should return [rho, 0, P/(gamma-1)]
    (u_ghost = 0 per spec §2.3, so total energy density equals
    internal energy density = P/(gamma-1)).
    """
    tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS)
    g = tank.ghost_state()
    assert g.shape == (3,)
    assert abs(g[0] - tank.rho) < 1e-12
    assert g[1] == 0.0
    assert abs(g[2] - 10e6 / (GAMMA - 1)) < 1e-3


def test_tank_apply_flux_outflow_reduces_mass_energy_and_pressure():
    """
    apply_flux(mdot, edot, dt, sign=-1) should subtract mdot*dt from mass
    and edot*dt from U. Derived P should decrease correspondingly.
    """
    tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS)
    mass_before = tank.mass
    U_before = tank.U
    P_before = tank.P
    mdot = 1.0      # kg/s
    edot = 5e5      # J/s (enthalpy rate)
    dt = 1e-3
    tank.apply_flux(mdot=mdot, edot=edot, dt=dt, sign=-1)
    assert abs(tank.mass - (mass_before - mdot * dt)) < 1e-12
    assert abs(tank.U    - (U_before    - edot * dt)) < 1e-9
    assert tank.P < P_before

• Step 2: Run tests to verify they fail

pytest tests/test_tank.py -v

Expected: 3 failures with ModuleNotFoundError: No module named 'tank'.

• Step 3: Write the implementation

Create src/tank.py:

# src/tank.py
"""
0D lumped-parameter tank for ideal gas. The tank's *primary* state is
(mass, U) where U is total internal energy in joules. Pressure, temperature,
and density are derived properties computed from (mass, U) on demand, so
they are always consistent with the conservation-law updates.

Conservation laws (u=0 inside tank):
  dm/dt = mdot_in                     (mass)
  dU/dt = Hdot_in = mdot_in * h_t,in  (energy, open-system first law)

where h_t is specific total enthalpy. When the tank couples to a 1D pipe
through the HLL boundary flux, flux[0]*A = mdot and flux[2]*A = Hdot
automatically — see solver.py.
"""

import numpy as np


class Tank:
    def __init__(self, V, P_init, T_init, gamma, R_gas):
        self.V = V
        self.gamma = gamma
        self.R = R_gas
        rho = P_init / (R_gas * T_init)
        self.mass = rho * V
        # For u=0, total internal energy equals rho*e*V = P*V / (gamma-1)
        self.U = P_init * V / (gamma - 1)

    @property
    def rho(self):
        return self.mass / self.V

    @property
    def T(self):
        return (self.U / self.mass) * (self.gamma - 1) / self.R

    @property
    def P(self):
        return self.rho * self.R * self.T

    def ghost_state(self):
        """
        Return the conservative variable vector [rho, rho*u, rho*E] that
        represents this tank as a ghost cell for the 1D pipe solver.
        Since u_ghost = 0, rho*u = 0 and rho*E = P/(gamma-1).
        """
        return np.array([self.rho, 0.0, self.P / (self.gamma - 1)])

    def apply_flux(self, mdot, edot, dt, sign):
        """
        Update (mass, U) from one time step of boundary flux.

        Parameters
        ----------
        mdot : float
            Mass flux across the interface in kg/s (already multiplied
            by pipe cross-sectional area). Sign is the "outward normal"
            convention of the pipe: positive = pipe-rightward.
        edot : float
            Total enthalpy rate in W (= flux[2] * A), same convention.
        dt : float
            Time-step size in seconds.
        sign : int (+1 or -1)
            Orientation for this tank. For an upstream tank whose gas
            flows "out to the right" into the pipe, the HLL left-boundary
            flux has mdot > 0, so sign = -1 (tank loses mass).
            For a downstream tank receiving gas from the right boundary
            with mdot > 0 entering, sign = +1.

        Raises
        ------
        RuntimeError
            If the tank's mass becomes non-positive after the update.
        """
        self.mass += sign * mdot * dt
        self.U    += sign * edot * dt
        if self.mass <= 0:
            raise RuntimeError(
                f"Tank mass non-positive after apply_flux: mass={self.mass}, "
                f"mdot={mdot}, edot={edot}, dt={dt}, sign={sign}"
            )

• Step 4: Run tests to verify they pass

pytest tests/test_tank.py -v

Expected: 3 passed.

• Step 5: Commit

git add src/tank.py tests/test_tank.py
git commit -m "feat(tank): 0D tank with (mass, U) primary state and 3 tests"

---

Task 5: src/pipe.py — 1D Pipe (TDD)

Pipe class storing the conservative-variable field W of shape (3, N). step() takes left/right boundary fluxes and advances the field by one explicit Euler step.

Files:

• Create: tests/test_pipe.py

• Create: src/pipe.py

• Step 1: Write the failing tests

Create tests/test_pipe.py:

# tests/test_pipe.py
import numpy as np
import pytest
from pipe import Pipe


GAMMA = 1.4
R_GAS = 287.0


def test_pipe_uniform_initialization_all_cells_identical():
    """
    Initial state: uniform P, T, u=0 across all N cells.
    All cells should have identical W. primitives() should round-trip
    back to u=0 and P=P_init.
    """
    pipe = Pipe(L=1.0, D=5e-3, N=20,
                P_init=1e5, T_init=300.0,
                gamma=GAMMA, R_gas=R_GAS)

    # All cells identical
    for i in range(1, pipe.N):
        assert np.allclose(pipe.W[:, i], pipe.W[:, 0], rtol=1e-14)

    rho, u, P, a = pipe.primitives()
    rho_expected = 1e5 / (R_GAS * 300.0)
    assert np.allclose(u, 0.0)
    assert np.allclose(P, 1e5, rtol=1e-10)
    assert np.allclose(rho, rho_expected, rtol=1e-10)

    # Sound speed a = sqrt(gamma * P / rho)
    a_expected = np.sqrt(GAMMA * 1e5 / rho_expected)
    assert np.allclose(a, a_expected, rtol=1e-10)

    # Geometric sanity
    assert pipe.dx == 1.0 / 20
    assert pipe.x_centers.shape == (20,)
    assert abs(pipe.x_centers[0] - 0.025) < 1e-14
    assert abs(pipe.x_centers[-1] - 0.975) < 1e-14


def test_pipe_uniform_state_plus_matching_boundary_flux_is_static():
    """
    If all cells have identical W (so all internal HLL fluxes are
    identical to the physical flux F(W) = [0, P, 0] for u=0 state),
    and we pass in boundary fluxes equal to [0, P, 0] as well, then
    every difference (flux[i+1] - flux[i]) is zero, so W must not
    change after one step. Verify to machine precision.
    """
    pipe = Pipe(L=1.0, D=5e-3, N=20,
                P_init=1e5, T_init=300.0,
                gamma=GAMMA, R_gas=R_GAS)

    W_before = pipe.W.copy()
    # Boundary flux matching the static interior: [rho*u, rho*u^2+P, u*(E+P)]
    # with u=0 -> [0, P, 0]
    flux_boundary = np.array([0.0, 1e5, 0.0])
    pipe.step(flux_boundary, flux_boundary, dt=1e-5)
    assert np.allclose(pipe.W, W_before, atol=1e-6, rtol=1e-12)

• Step 2: Run tests to verify they fail

pytest tests/test_pipe.py -v

Expected: 2 failures with ModuleNotFoundError: No module named 'pipe'.

• Step 3: Write the implementation

Create src/pipe.py:

# src/pipe.py
"""
1D finite-volume pipe for compressible Euler equations:
  dW/dt + dF(W)/dx = 0
  W = [rho, rho*u, rho*E],  F = [rho*u, rho*u^2 + P, u*(rho*E + P)]

Discretization:
  - N uniform cells, cell-averaged piecewise-constant reconstruction
  - HLL numerical flux at all interior interfaces
  - Boundary (tank-side) interface fluxes are provided by the caller
    via step(flux_L, flux_R, dt)
"""

import numpy as np
from riemann import hll_flux


class Pipe:
    def __init__(self, L, D, N, P_init, T_init, gamma, R_gas):
        self.L = L
        self.D = D
        self.N = N
        self.dx = L / N
        self.area = np.pi * (D / 2) ** 2
        self.gamma = gamma
        self.R = R_gas
        self.x_centers = np.linspace(self.dx / 2, L - self.dx / 2, N)

        # Uniform initial state, u = 0
        rho = P_init / (R_gas * T_init)
        E_density = P_init / (gamma - 1)  # since u=0, total energy density = internal

        self.W = np.zeros((3, N))
        self.W[0, :] = rho
        self.W[1, :] = 0.0
        self.W[2, :] = E_density

    def primitives(self):
        """
        Return (rho, u, P, a) each of shape (N,), computed from W.
        """
        rho = self.W[0, :]
        u = self.W[1, :] / rho
        P = (self.gamma - 1) * (self.W[2, :] - 0.5 * rho * u ** 2)
        a = np.sqrt(self.gamma * P / rho)
        return rho, u, P, a

    def max_wave_speed(self):
        """
        Return max over cells of |u| + a, used for CFL dt calculation.
        """
        _, u, _, a = self.primitives()
        return float(np.max(np.abs(u) + a))

    def step(self, flux_L, flux_R, dt):
        """
        Advance W by one explicit Euler step. The caller provides the
        two boundary interface fluxes (with ghost states already folded
        in); internal interface fluxes are computed here with HLL.

        Parameters
        ----------
        flux_L, flux_R : np.ndarray of shape (3,)
            Numerical fluxes at the leftmost and rightmost interfaces
            (cell -1/2 and cell N-1/2, i.e. the tank-facing boundaries).
        dt : float
            Time-step size.

        Raises
        ------
        RuntimeError
            If the updated state has any non-positive density or pressure.
        """
        N = self.N
        W_snap = self.W.copy()

        # Internal fluxes: flux_int[:, k] is the flux at the interface
        # between cell k and cell k+1, for k = 0 .. N-2  (total N-1 of them)
        flux_int = np.zeros((3, N - 1))
        for k in range(N - 1):
            flux_int[:, k] = hll_flux(W_snap[:, k], W_snap[:, k + 1], self.gamma)

        # First cell: left face = flux_L, right face = flux_int[:, 0]
        self.W[:, 0] = W_snap[:, 0] - (dt / self.dx) * (flux_int[:, 0] - flux_L)

        # Interior cells: left face = flux_int[:, i-1], right face = flux_int[:, i]
        for i in range(1, N - 1):
            self.W[:, i] = W_snap[:, i] - (dt / self.dx) * (flux_int[:, i] - flux_int[:, i - 1])

        # Last cell: left face = flux_int[:, N-2], right face = flux_R
        self.W[:, N - 1] = W_snap[:, N - 1] - (dt / self.dx) * (flux_R - flux_int[:, N - 2])

        # Physical-state sanity check
        rho_new = self.W[0, :]
        if np.any(rho_new <= 0):
            bad = np.where(rho_new <= 0)[0]
            raise RuntimeError(
                f"Non-positive density after pipe step at cells {bad.tolist()}: "
                f"rho={rho_new[bad].tolist()}"
            )
        u_new = self.W[1, :] / rho_new
        P_new = (self.gamma - 1) * (self.W[2, :] - 0.5 * rho_new * u_new ** 2)
        if np.any(P_new <= 0):
            bad = np.where(P_new <= 0)[0]
            raise RuntimeError(
                f"Non-positive pressure after pipe step at cells {bad.tolist()}: "
                f"P={P_new[bad].tolist()}"
            )

• Step 4: Run tests to verify they pass

pytest tests/test_pipe.py -v

Expected: 2 passed.

• Step 5: Commit

git add src/pipe.py tests/test_pipe.py
git commit -m "feat(pipe): 1D FV pipe with HLL internal fluxes and 2 tests"

---

Task 6: src/solver.py + Integration Tests (TDD)

The time-loop driver that orchestrates tank ↔ pipe coupling. The integration tests verify total mass and total energy are conserved to 1e-10 over a short run — the end-to-end correctness check for the "flux doubling" mechanism.

Files:

• Create: tests/test_integration.py

• Create: src/solver.py

• Step 1: Write the failing integration tests

Create tests/test_integration.py:

# tests/test_integration.py
"""
End-to-end conservation tests: after a short run, total system mass
and total system energy must be conserved to machine precision. These
are the load-bearing tests for the "flux doubling" coupling mechanism.
"""
import numpy as np
import pytest
from tank import Tank
from pipe import Pipe
from solver import run


GAMMA = 1.4
R_GAS = 287.0


def _build_scenario():
    """Default blowdown scenario at smaller scale — same physics, short run."""
    tank1 = Tank(V=5.0,  P_init=10e6,    T_init=300.0, gamma=GAMMA, R_gas=R_GAS)
    tank2 = Tank(V=10.0, P_init=101325.0, T_init=300.0, gamma=GAMMA, R_gas=R_GAS)
    pipe = Pipe(L=1.0, D=5e-3, N=20,
                P_init=101325.0, T_init=300.0, gamma=GAMMA, R_gas=R_GAS)
    return tank1, tank2, pipe


def _total_mass(tank1, tank2, pipe):
    return tank1.mass + tank2.mass + float(np.sum(pipe.W[0, :] * pipe.area * pipe.dx))


def _total_energy(tank1, tank2, pipe):
    return tank1.U + tank2.U + float(np.sum(pipe.W[2, :] * pipe.area * pipe.dx))


def test_total_mass_conserved_over_short_run():
    tank1, tank2, pipe = _build_scenario()
    m_init = _total_mass(tank1, tank2, pipe)
    run(tank1, tank2, pipe, t_end=1e-3, cfl=0.5)
    m_final = _total_mass(tank1, tank2, pipe)
    rel_err = abs(m_final - m_init) / m_init
    assert rel_err < 1e-10, (
        f"Total mass not conserved: m_init={m_init:.6e}, "
        f"m_final={m_final:.6e}, rel_err={rel_err:.2e}"
    )


def test_total_energy_conserved_over_short_run():
    tank1, tank2, pipe = _build_scenario()
    U_init = _total_energy(tank1, tank2, pipe)
    run(tank1, tank2, pipe, t_end=1e-3, cfl=0.5)
    U_final = _total_energy(tank1, tank2, pipe)
    rel_err = abs(U_final - U_init) / U_init
    assert rel_err < 1e-10, (
        f"Total energy not conserved: U_init={U_init:.6e}, "
        f"U_final={U_final:.6e}, rel_err={rel_err:.2e}"
    )

• Step 2: Run the tests to verify they fail

pytest tests/test_integration.py -v

Expected: 2 failures with ModuleNotFoundError: No module named 'solver'.

• Step 3: Write the solver implementation

Create src/solver.py:

# src/solver.py
"""
Time-loop driver for the 0D-1D coupled tank-pipe simulation.

Per time step (per spec §3.1):
  1. Compute CFL-limited dt from pipe's max wave speed
  2. Freeze ghost states from current tank states
  3. Compute two boundary HLL fluxes (left and right)
  4. Advance pipe by one step using those two fluxes (pipe.step handles
     the internal fluxes itself)
  5. Advance both tanks using the SAME two boundary fluxes * area
     -> this "flux doubling" is the mechanism that makes system mass
        and energy strictly conserved to machine precision
  6. Advance time
  7. Append snapshot to history
"""
import numpy as np
from riemann import hll_flux


def run(tank1, tank2, pipe, t_end, cfl, verbose=False, log_every=100):
    """
    Run the coupled tank-pipe simulation from t=0 to t=t_end.

    Parameters
    ----------
    tank1, tank2 : Tank
        Upstream and downstream tanks. tank1 connects to pipe.W[:, 0],
        tank2 connects to pipe.W[:, -1].
    pipe : Pipe
        1D pipe instance with initial state already set.
    t_end : float
        End time in seconds.
    cfl : float
        CFL number in (0, 1].
    verbose : bool, default False
        If True, print step-progress info every `log_every` steps.
    log_every : int, default 100
        Logging interval when verbose=True.

    Returns
    -------
    dict
        History with keys 't', 'P1', 'T1', 'P2', 'T2' (all 1D arrays
        of shape (n_steps,)), and 'W_hist' of shape (n_steps, 3, N).
    """
    history = {
        't': [],
        'P1': [], 'T1': [],
        'P2': [], 'T2': [],
        'W_hist': [],
    }

    t = 0.0
    step = 0

    while t < t_end:
        # --- Phase 1: CFL time step ---
        a_max = pipe.max_wave_speed()
        dt = cfl * pipe.dx / a_max
        dt = min(dt, t_end - t)
        if dt < 1e-12:
            raise RuntimeError(
                f"dt degenerate at step {step}: dt={dt:.3e}, a_max={a_max:.3e}"
            )

        # --- Phase 2: freeze tank ghost states (snapshot for this step) ---
        W_ghost_L = tank1.ghost_state()
        W_ghost_R = tank2.ghost_state()

        # --- Phase 3: two boundary fluxes (solver-level) ---
        flux_L = hll_flux(W_ghost_L,     pipe.W[:, 0],  pipe.gamma)
        flux_R = hll_flux(pipe.W[:, -1], W_ghost_R,     pipe.gamma)

        # --- Phase 4: advance pipe (internal fluxes handled inside) ---
        pipe.step(flux_L, flux_R, dt)

        # --- Phase 5: advance tanks with the SAME boundary fluxes * area ---
        fL_A = flux_L * pipe.area
        fR_A = flux_R * pipe.area
        # Left boundary flux is "rightward positive"; tank1 loses that mass
        tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=dt, sign=-1)
        # Right boundary flux is "rightward positive"; tank2 gains that mass
        tank2.apply_flux(mdot=fR_A[0], edot=fR_A[2], dt=dt, sign=+1)

        # --- Phase 6: advance time ---
        t += dt
        step += 1

        # --- Phase 7: record history ---
        history['t'].append(t)
        history['P1'].append(tank1.P)
        history['T1'].append(tank1.T)
        history['P2'].append(tank2.P)
        history['T2'].append(tank2.T)
        history['W_hist'].append(pipe.W.copy())

        if verbose and step % log_every == 0:
            _, u, _, _ = pipe.primitives()
            print(
                f"step={step:6d}  t={t:.5f}  dt={dt:.2e}  "
                f"P1={tank1.P/1e6:7.4f}MPa  P2={tank2.P/1e6:7.4f}MPa  "
                f"max|u|={float(np.max(np.abs(u))):7.1f}m/s"
            )

    if step == 0:
        raise RuntimeError("solver.run() exited without taking any step")

    # Convert lists to arrays for downstream consumers
    history['t']      = np.asarray(history['t'])
    history['P1']     = np.asarray(history['P1'])
    history['T1']     = np.asarray(history['T1'])
    history['P2']     = np.asarray(history['P2'])
    history['T2']     = np.asarray(history['T2'])
    history['W_hist'] = np.stack(history['W_hist'])   # shape (n_steps, 3, N)

    return history

• Step 4: Run integration tests to verify they pass

pytest tests/test_integration.py -v

Expected: 2 passed.

• Step 5: Run the full test suite to verify nothing regressed

pytest tests/ -v

Expected: 10 passed total (3 riemann + 3 tank + 2 pipe + 2 integration).

• Step 6: Commit

git add src/solver.py tests/test_integration.py
git commit -m "feat(solver): coupled time loop with mass/energy conservation tests"

---

Task 7: src/output.py — Persistence, Plotting, Animation

Utility module for .npz persistence, static time-series plots, and the pipe evolution GIF. No unit tests per spec §7.3 (visual outputs verified by eye).

Files:

• Create: src/output.py

• Step 1: Write src/output.py

# src/output.py
"""
Output helpers: persistence (.npz), static plots (.png), animation (.gif).
Uses matplotlib's Agg backend so it works in headless environments.
The PillowWriter is used for GIF output to avoid an ffmpeg dependency.
"""
import os

import numpy as np
import matplotlib
matplotlib.use("Agg")  # headless-safe; must be set before pyplot import
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter


def save_history(history, pipe, path, gamma, R_gas):
    """
    Persist the full simulation history + pipe geometry to a .npz file.

    Loadable later with:
        d = np.load("results/history.npz")
        rho = d['W_hist'][:, 0, :]
        u   = d['W_hist'][:, 1, :] / rho
        P   = (d['gamma'] - 1) * (d['W_hist'][:, 2, :] - 0.5 * rho * u**2)
    """
    dirname = os.path.dirname(path)
    if dirname:
        os.makedirs(dirname, exist_ok=True)
    np.savez_compressed(
        path,
        t=history['t'],
        P1=history['P1'], T1=history['T1'],
        P2=history['P2'], T2=history['T2'],
        W_hist=history['W_hist'],
        x=pipe.x_centers,
        dx=pipe.dx,
        area=pipe.area,
        gamma=gamma,
        R_gas=R_gas,
    )


def plot_tank_pressure(history, path):
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(history['t'], history['P1'] / 1e6, label="Tank 1 (high pressure)")
    ax.plot(history['t'], history['P2'] / 1e6, label="Tank 2 (low pressure)")
    ax.set_xlabel("Time [s]")
    ax.set_ylabel("Pressure [MPa]")
    ax.set_title("Tank pressures vs time")
    ax.grid(True)
    ax.legend()
    fig.tight_layout()
    fig.savefig(path, dpi=120)
    plt.close(fig)


def plot_tank_temperature(history, path):
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(history['t'], history['T1'], label="Tank 1 (high pressure)")
    ax.plot(history['t'], history['T2'], label="Tank 2 (low pressure)")
    ax.set_xlabel("Time [s]")
    ax.set_ylabel("Temperature [K]")
    ax.set_title("Tank temperatures vs time")
    ax.grid(True)
    ax.legend()
    fig.tight_layout()
    fig.savefig(path, dpi=120)
    plt.close(fig)


def make_pipe_animation(history, pipe, path, gamma, R_gas, stride=10):
    """
    Render a GIF of the pipe's P(x), u(x), T(x) evolution over time.
    Uses PillowWriter so no ffmpeg is needed.
    """
    t = history['t']
    W_hist = history['W_hist']  # (n_steps, 3, N)
    x = pipe.x_centers
    n_steps = W_hist.shape[0]

    # Frame indices: every `stride`th snapshot, plus the final one
    frames = list(range(0, n_steps, stride))
    if frames[-1] != n_steps - 1:
        frames.append(n_steps - 1)

    # Precompute primitives for all frames in one vectorized pass
    rho = W_hist[:, 0, :]
    u = W_hist[:, 1, :] / rho
    P = (gamma - 1) * (W_hist[:, 2, :] - 0.5 * rho * u ** 2)
    T = P / (rho * R_gas)

    fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)

    # Pressure subplot
    line_P, = axes[0].plot(x, P[0] / 1e6)
    axes[0].set_ylabel("P [MPa]")
    axes[0].set_ylim(P.min() / 1e6 * 0.95, P.max() / 1e6 * 1.05)
    axes[0].grid(True)

    # Velocity subplot
    line_u, = axes[1].plot(x, u[0])
    axes[1].set_ylabel("u [m/s]")
    u_min, u_max = float(u.min()), float(u.max())
    pad = max(1.0, 0.05 * (u_max - u_min) if u_max > u_min else 1.0)
    axes[1].set_ylim(u_min - pad, u_max + pad)
    axes[1].grid(True)

    # Temperature subplot
    line_T, = axes[2].plot(x, T[0])
    axes[2].set_ylabel("T [K]")
    axes[2].set_xlabel("x [m]")
    axes[2].set_ylim(T.min() * 0.95, T.max() * 1.05)
    axes[2].grid(True)

    title = fig.suptitle("")

    def update(frame_idx):
        line_P.set_ydata(P[frame_idx] / 1e6)
        line_u.set_ydata(u[frame_idx])
        line_T.set_ydata(T[frame_idx])
        title.set_text(
            f"t = {t[frame_idx]:.5f} s   (step {frame_idx + 1}/{n_steps})"
        )
        return line_P, line_u, line_T, title

    anim = FuncAnimation(fig, update, frames=frames, interval=50, blit=False)

    dirname = os.path.dirname(path)
    if dirname:
        os.makedirs(dirname, exist_ok=True)
    anim.save(path, writer=PillowWriter(fps=20))
    plt.close(fig)

• Step 2: Smoke-test that the module imports without errors

python3 -c "import sys; sys.path.insert(0,'src'); import output; print('output OK:', dir(output))"

Expected: output OK: [...] 'make_pipe_animation', 'plot_tank_pressure', 'plot_tank_temperature', 'save_history', ...

• Step 3: Commit

git add src/output.py
git commit -m "feat(output): npz persistence, timeseries plots, pipe GIF animation"

---

Task 8: src/main.py — Entry Point + Sanity Check

Wires everything together: build objects from config, call solver.run, assert conservation, write all outputs.

Files:

• Create: src/main.py

• Step 1: Write src/main.py

# src/main.py
"""
Entry point: assemble tanks + pipe from config constants, run the solver,
verify total mass/energy conservation, persist history, and generate
plots + animation.

Run from project root:
    python3 src/main.py
"""
import os
import sys

# Ensure imports work when running from project root
_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, _HERE)

import numpy as np

from config import (
    GAMMA, R_GAS,
    V1, P1_INIT, T1_INIT,
    V2, P2_INIT, T2_INIT,
    L, D, N_CELLS,
    T_END, CFL,
    ANIMATION_STRIDE, OUTPUT_DIR,
)
from tank import Tank
from pipe import Pipe
from solver import run
from output import (
    save_history,
    plot_tank_pressure,
    plot_tank_temperature,
    make_pipe_animation,
)


def _total_mass(tank1, tank2, pipe):
    pipe_mass = float(np.sum(pipe.W[0, :] * pipe.area * pipe.dx))
    return tank1.mass + tank2.mass + pipe_mass


def _total_energy(tank1, tank2, pipe):
    pipe_energy = float(np.sum(pipe.W[2, :] * pipe.area * pipe.dx))
    return tank1.U + tank2.U + pipe_energy


def main():
    os.makedirs(OUTPUT_DIR, exist_ok=True)

    # --- Assemble ---
    tank1 = Tank(V=V1, P_init=P1_INIT, T_init=T1_INIT, gamma=GAMMA, R_gas=R_GAS)
    tank2 = Tank(V=V2, P_init=P2_INIT, T_init=T2_INIT, gamma=GAMMA, R_gas=R_GAS)
    pipe = Pipe(L=L, D=D, N=N_CELLS, P_init=P2_INIT, T_init=T2_INIT,
                gamma=GAMMA, R_gas=R_GAS)

    m_init = _total_mass(tank1, tank2, pipe)
    U_init = _total_energy(tank1, tank2, pipe)
    print(f"Initial total mass:   {m_init:.6e} kg")
    print(f"Initial total energy: {U_init:.6e} J")
    print(f"Initial P1 = {tank1.P/1e6:.3f} MPa, P2 = {tank2.P/1e6:.3f} MPa")
    print(f"Pipe: L={L} m, D={D*1e3:.1f} mm, N={N_CELLS} cells, dx={pipe.dx*1e3:.1f} mm")
    print(f"Running to t_end={T_END} s with CFL={CFL}...")
    print()

    # --- Run ---
    history = run(tank1, tank2, pipe,
                  t_end=T_END, cfl=CFL,
                  verbose=True, log_every=200)

    n_steps = len(history['t'])
    print()
    print(f"Simulation complete: {n_steps} steps")

    # --- Conservation sanity check (per spec §6.1, §8) ---
    m_final = _total_mass(tank1, tank2, pipe)
    U_final = _total_energy(tank1, tank2, pipe)
    rel_err_m = abs(m_final - m_init) / m_init
    rel_err_U = abs(U_final - U_init) / U_init
    print(f"Final total mass:     {m_final:.6e} kg   (rel err = {rel_err_m:.2e})")
    print(f"Final total energy:   {U_final:.6e} J    (rel err = {rel_err_U:.2e})")
    print(f"Final P1 = {tank1.P/1e6:.3f} MPa, P2 = {tank2.P/1e6:.3f} MPa")
    assert rel_err_m < 1e-10, f"Total mass not conserved: rel_err={rel_err_m:.2e}"
    assert rel_err_U < 1e-10, f"Total energy not conserved: rel_err={rel_err_U:.2e}"

    # --- Persist + visualize ---
    save_history(history, pipe,
                 os.path.join(OUTPUT_DIR, "history.npz"),
                 GAMMA, R_GAS)
    plot_tank_pressure(history,
                       os.path.join(OUTPUT_DIR, "tank_pressure.png"))
    plot_tank_temperature(history,
                          os.path.join(OUTPUT_DIR, "tank_temperature.png"))
    make_pipe_animation(history, pipe,
                        os.path.join(OUTPUT_DIR, "pipe_animation.gif"),
                        GAMMA, R_GAS, stride=ANIMATION_STRIDE)

    print(f"Outputs written to {OUTPUT_DIR}/")
    print(f"  - history.npz")
    print(f"  - tank_pressure.png")
    print(f"  - tank_temperature.png")
    print(f"  - pipe_animation.gif")


if __name__ == "__main__":
    main()

• Step 2: Verify the module imports cleanly (no run)

python3 -c "import sys; sys.path.insert(0,'src'); import main; print('main module OK')"

Expected: main module OK (importing must not trigger main() because of if __name__ == '__main__').

• Step 3: Commit

git add src/main.py
git commit -m "feat(main): entry point with conservation sanity check"

---

Task 9: End-to-End Verification Run

Run the full program and verify all acceptance criteria from spec §8.

Files:

• No new files; verifies existing ones.

• Step 1: Run the full simulation

python3 src/main.py

Expected behavior:

• Prints "Initial total mass", "Initial total energy", initial pressures

• Prints verbose step logs every 200 steps (roughly 10–15 log lines total)

• Prints "Simulation complete: NNNN steps" where NNNN ≈ 2000–3000

• Prints "Final total mass" and "rel err" < 1e-10

• Prints "Final total energy" and "rel err" < 1e-10

• Prints "Outputs written to results/"

• Exit code 0

• Total wall time < 30 seconds

• Step 2: Verify all 4 output files exist and are non-empty

ls -la results/history.npz results/tank_pressure.png results/tank_temperature.png results/pipe_animation.gif

Expected: all four files present, each with non-zero size (history.npz ≈ 1-3 MB, PNGs tens of KB, GIF 5-15 MB).

• Step 3: Verify the .npz file is loadable and has correct shapes

python3 -c "
import numpy as np
d = np.load('results/history.npz')
print('keys:', sorted(d.files))
print('t shape:', d['t'].shape)
print('W_hist shape:', d['W_hist'].shape)
print('P1 range:', float(d['P1'].min())/1e6, '->', float(d['P1'].max())/1e6, 'MPa')
print('P2 range:', float(d['P2'].min())/1e6, '->', float(d['P2'].max())/1e6, 'MPa')
# Physical monotonicity
assert d['P1'][0] > d['P1'][-1], 'P1 should decrease'
assert d['P2'][0] < d['P2'][-1], 'P2 should increase'
print('monotonicity OK')
"

Expected:

• keys: ['P1', 'P2', 'R_gas', 'T1', 'T2', 'W_hist', 'area', 'dx', 'gamma', 't', 'x']

• t shape: (NNNN,) where NNNN ≈ 2000–3000

• W_hist shape: (NNNN, 3, 20)

• P1 range: starts near 10, decreases toward some intermediate value

• P2 range: starts near 0.1, increases toward some intermediate value

• monotonicity OK

• Step 4: Run the full pytest suite one more time

pytest tests/ -v

Expected: 10 passed.

• Step 5: Commit a .gitignore entry for the bulky result files

Check whether .gitignore already excludes results/:

cat .gitignore

If results/ is not already ignored, update .gitignore to include it (but keep .gitkeep):

Add these two lines to .gitignore:

results/*
!results/.gitkeep

Then:

git add .gitignore
git commit -m "chore: ignore runtime outputs under results/"

If results/ was already handled, skip this step with git status confirming no .gitignore changes are staged.

• Step 6: Final status check

git status
git log --oneline -12

Expected:

• git status shows a clean working tree
• git log shows the sequence of feature commits: scaffold → config → riemann → tank → pipe → solver → output → main → (optionally) .gitignore

---

Acceptance Criteria Recap (from spec §8)

After Task 9 is complete, verify each item is satisfied:

1. python3 src/main.py runs to completion without exceptions → Task 9 Step 1
2. pytest tests/ -v → 10 passed → Task 9 Step 4
3. main.py sanity check: total mass & total energy rel err < 1e-10 → Task 9 Step 1 output
4. P1(t) monotone decreasing, P2(t) monotone increasing → Task 9 Step 3
5. pipe_animation.gif is a valid GIF openable by standard viewers → Task 9 Step 2 + manual inspection
6. history.npz loadable via np.load with all fields accessible → Task 9 Step 3
7. Wall time < 30 seconds → Task 9 Step 1