refactor: organize src modules by category

src/tank_pipe/__init__.py

@@ -0,0 +1 @@
+"""0D-1D tank-pipe blowdown simulation package."""

src/tank_pipe/config.py

@@ -0,0 +1,52 @@
+"""
+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 = 2e6       # Pa  (2 MPa)
+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
+
+# ---------- Friction ----------
+MU = 1.8e-5         # Pa·s  (dynamic viscosity of air at ~300 K)
+ROUGHNESS = 0.0     # m  (absolute wall roughness; 0 = smooth pipe)
+
+# ---------- Simulation control ----------
+T_END = 0.1         # s
+CFL = 0.5
+RIEMANN_SOLVER = "roe"  # "hll" or "roe"
+
+# ---------- 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 MU >= 0, "MU must be >= 0"
+assert ROUGHNESS >= 0, "ROUGHNESS must be >= 0"
+assert 0 < CFL <= 1, "CFL must be in (0, 1]"
+assert RIEMANN_SOLVER in ("hll", "roe"), "RIEMANN_SOLVER must be 'hll' or 'roe'"
+assert T_END > 0, "T_END must be > 0"
+assert ANIMATION_STRIDE >= 1, "ANIMATION_STRIDE must be >= 1"

src/tank_pipe/friction.py

@@ -0,0 +1,74 @@
+# src/friction.py
+"""
+Darcy-Weisbach friction factor calculation.
+
+Supports:
+  - Laminar:   f = 64 / Re                           (Re < 2300)
+  - Turbulent: Colebrook-White implicit equation      (Re > 4000)
+  - Transition: linear blend between laminar & turbulent (2300 <= Re <= 4000)
+"""
+
+import numpy as np
+
+
+def _colebrook_white(Re, eps_D, n_iter=10):
+    """
+    Solve the Colebrook-White equation for Darcy friction factor f:
+        1/sqrt(f) = -2 log10( eps_D/3.7 + 2.51/(Re*sqrt(f)) )
+
+    Uses fixed-point iteration seeded with the Swamee-Jain approximation.
+    """
+    # Swamee-Jain initial guess (explicit approximation)
+    A = eps_D / 3.7
+    B = 2.51 / Re
+    f = 0.25 / (np.log10(A + B / np.sqrt(0.02))) ** 2
+
+    for _ in range(n_iter):
+        f = 0.25 / (np.log10(A + B / np.sqrt(f))) ** 2
+
+    return f
+
+
+def darcy_friction_factor(Re, eps_D):
+    """
+    Compute Darcy-Weisbach friction factor for a given Reynolds number
+    and relative roughness eps/D.
+
+    Parameters
+    ----------
+    Re : float or ndarray
+        Reynolds number (ρ|u|D/μ). Values <= 0 return 0 (no flow).
+    eps_D : float
+        Relative roughness ε/D (dimensionless).
+
+    Returns
+    -------
+    f : same shape as Re
+        Darcy friction factor.
+    """
+    Re = np.asarray(Re, dtype=float)
+    scalar = Re.ndim == 0
+    Re = np.atleast_1d(Re)
+
+    f = np.zeros_like(Re)
+
+    lam = Re < 2300
+    turb = Re > 4000
+    trans = ~lam & ~turb  # 2300 <= Re <= 4000
+
+    # Laminar: f = 64/Re (avoid division by zero for Re~0)
+    Re_lam = np.where(Re > 1e-12, Re, 1e-12)
+    f[lam] = 64.0 / Re_lam[lam]
+
+    # Turbulent: Colebrook-White
+    if np.any(turb):
+        f[turb] = _colebrook_white(Re[turb], eps_D)
+
+    # Transition: linear blend
+    if np.any(trans):
+        f_lam = 64.0 / Re_lam[trans]
+        f_turb = _colebrook_white(Re[trans], eps_D)
+        alpha = (Re[trans] - 2300.0) / 1700.0  # 0 at Re=2300, 1 at Re=4000
+        f[trans] = (1.0 - alpha) * f_lam + alpha * f_turb
+
+    return float(f[0]) if scalar else f

src/tank_pipe/main.py

@@ -0,0 +1,130 @@
+# src/tank_pipe/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/tank_pipe/main.py
+"""
+import os
+import sys
+
+_HERE = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.dirname(_HERE))
+
+import numpy as np
+
+from tank_pipe.config import (
+    GAMMA, R_GAS,
+    V1, P1_INIT, T1_INIT,
+    V2, P2_INIT, T2_INIT,
+    L, D, N_CELLS,
+    MU, ROUGHNESS,
+    T_END, CFL, RIEMANN_SOLVER,
+    ANIMATION_STRIDE, OUTPUT_DIR,
+)
+from tank_pipe.tank import Tank
+from tank_pipe.pipe import Pipe
+from tank_pipe.solver import run
+from tank_pipe.output import (
+    save_history,
+    plot_tank_pressure,
+    plot_tank_temperature,
+    plot_pipe_final_profiles,
+    make_pipe_animation,
+    write_summary_report,
+)
+
+
+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, mu=MU, roughness=ROUGHNESS,
+                riemann_solver=RIEMANN_SOLVER)
+
+    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"Riemann solver: {RIEMANN_SOLVER.upper()}")
+    if MU > 0:
+        print(f"Friction: mu={MU:.2e} Pa·s, roughness={ROUGHNESS:.2e} m (eps/D={ROUGHNESS/D:.4f})")
+    else:
+        print("Friction: OFF")
+    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"))
+    plot_pipe_final_profiles(history, pipe,
+                             os.path.join(OUTPUT_DIR, "pipe_final_profiles.png"),
+                             GAMMA, R_GAS)
+    make_pipe_animation(history, pipe,
+                        os.path.join(OUTPUT_DIR, "pipe_animation.gif"),
+                        GAMMA, R_GAS, stride=ANIMATION_STRIDE)
+    write_summary_report(
+        history, pipe,
+        os.path.join(OUTPUT_DIR, "summary_report.html"),
+        GAMMA, R_GAS,
+        config={
+            'V1': V1, 'P1_INIT': P1_INIT, 'T1_INIT': T1_INIT,
+            'V2': V2, 'P2_INIT': P2_INIT, 'T2_INIT': T2_INIT,
+            'L': L, 'D': D, 'N_CELLS': N_CELLS,
+            'T_END': T_END, 'CFL': CFL,
+        },
+    )
+
+    print(f"Outputs written to {OUTPUT_DIR}/")
+    print(f"  - history.npz")
+    print(f"  - tank_pressure.png")
+    print(f"  - tank_temperature.png")
+    print(f"  - pipe_final_profiles.png")
+    print(f"  - pipe_animation.gif")
+    print(f"  - summary_report.html")
+
+
+if __name__ == "__main__":
+    main()

src/tank_pipe/output.py

@@ -0,0 +1,300 @@
+# src/tank_pipe/output.py
+"""
+Output helpers: persistence (.npz), static plots (.png/.html), 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 html
+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 _pipe_primitives(history, gamma, R_gas):
+    W_hist = history['W_hist']
+    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)
+    a = np.sqrt(gamma * P / rho)
+    Ma = u / a
+    return rho, u, P, T, a, Ma
+
+
+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']
+    x = pipe.x_centers
+    n_steps = history['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
+    _, u, P, T, _, Ma = _pipe_primitives(history, gamma, R_gas)
+
+    fig, axes = plt.subplots(4, 1, figsize=(10, 12), 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_ylim(T.min() * 0.95, T.max() * 1.05)
+    axes[2].grid(True)
+
+    # Mach number subplot
+    line_Ma, = axes[3].plot(x, Ma[0])
+    axes[3].set_ylabel("Mach [-]")
+    axes[3].set_xlabel("x [m]")
+    Ma_min, Ma_max = float(Ma.min()), float(Ma.max())
+    pad_Ma = max(0.05, 0.05 * (Ma_max - Ma_min) if Ma_max > Ma_min else 0.05)
+    axes[3].set_ylim(Ma_min - pad_Ma, Ma_max + pad_Ma)
+    axes[3].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])
+        line_Ma.set_ydata(Ma[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, line_Ma, 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)
+
+
+def plot_pipe_final_profiles(history, pipe, path, gamma, R_gas):
+    """Plot final pipe P(x), u(x), T(x), and Mach(x) as a static PNG."""
+    x = pipe.x_centers
+    _, u, P, T, _, Ma = _pipe_primitives(history, gamma, R_gas)
+
+    fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True)
+    axes = axes.ravel()
+
+    axes[0].plot(x, P[-1] / 1e6)
+    axes[0].set_ylabel("P [MPa]")
+    axes[0].set_title("Final pressure profile")
+    axes[0].grid(True)
+
+    axes[1].plot(x, u[-1])
+    axes[1].set_ylabel("u [m/s]")
+    axes[1].set_title("Final velocity profile")
+    axes[1].grid(True)
+
+    axes[2].plot(x, T[-1])
+    axes[2].set_xlabel("x [m]")
+    axes[2].set_ylabel("T [K]")
+    axes[2].set_title("Final temperature profile")
+    axes[2].grid(True)
+
+    axes[3].plot(x, Ma[-1])
+    axes[3].axhline(1.0, color="r", linestyle="--", linewidth=1, label="Mach 1")
+    axes[3].set_xlabel("x [m]")
+    axes[3].set_ylabel("Mach [-]")
+    axes[3].set_title("Final Mach profile")
+    axes[3].grid(True)
+    axes[3].legend()
+
+    fig.tight_layout()
+    fig.savefig(path, dpi=120)
+    plt.close(fig)
+
+
+def write_summary_report(history, pipe, path, gamma, R_gas, config):
+    """Write an HTML summary report for the latest simulation outputs."""
+    t = history['t']
+    P1 = history['P1']
+    T1 = history['T1']
+    P2 = history['P2']
+    T2 = history['T2']
+    W_hist = history['W_hist']
+    x = pipe.x_centers
+
+    _, u, P, T, _, Ma = _pipe_primitives(history, gamma, R_gas)
+
+    V1 = config['V1']
+    V2 = config['V2']
+
+    m1_init = P1[0] * V1 / (R_gas * T1[0])
+    m1_final = P1[-1] * V1 / (R_gas * T1[-1])
+    m2_init = P2[0] * V2 / (R_gas * T2[0])
+    m2_final = P2[-1] * V2 / (R_gas * T2[-1])
+    mpipe_init = float(np.sum(W_hist[0, 0, :] * pipe.area * pipe.dx))
+    mpipe_final = float(np.sum(W_hist[-1, 0, :] * pipe.area * pipe.dx))
+    m_total_init = m1_init + m2_init + mpipe_init
+    m_total_final = m1_final + m2_final + mpipe_final
+    rel_err_m = abs(m_total_final - m_total_init) / m_total_init
+
+    U1_init = P1[0] * V1 / (gamma - 1)
+    U1_final = P1[-1] * V1 / (gamma - 1)
+    U2_init = P2[0] * V2 / (gamma - 1)
+    U2_final = P2[-1] * V2 / (gamma - 1)
+    Upipe_init = float(np.sum(W_hist[0, 2, :] * pipe.area * pipe.dx))
+    Upipe_final = float(np.sum(W_hist[-1, 2, :] * pipe.area * pipe.dx))
+    U_total_init = U1_init + U2_init + Upipe_init
+    U_total_final = U1_final + U2_final + Upipe_final
+    rel_err_U = abs(U_total_final - U_total_init) / U_total_init
+
+    idx_ma = np.unravel_index(np.argmax(Ma), Ma.shape)
+    idx_p = np.unravel_index(np.argmax(P), P.shape)
+
+    sections = [
+        ("Simulation setup", [
+            ("Gamma", f"{gamma:.3f}"),
+            ("R_gas", f"{R_gas:.3f} J/(kg K)"),
+            ("Tank 1", f"V={V1:.3f} m^3, P0={config['P1_INIT']/1e6:.6f} MPa, T0={config['T1_INIT']:.3f} K"),
+            ("Tank 2", f"V={V2:.3f} m^3, P0={config['P2_INIT']/1e6:.6f} MPa, T0={config['T2_INIT']:.3f} K"),
+            ("Pipe", f"L={config['L']:.3f} m, D={config['D']*1e3:.3f} mm, N={config['N_CELLS']}"),
+            ("Run control", f"t_end={config['T_END']:.6f} s, CFL={config['CFL']:.3f}, steps={len(t)}"),
+        ]),
+        ("Tank states", [
+            ("Tank 1 pressure", f"{P1[0]/1e6:.6f} -> {P1[-1]/1e6:.6f} MPa"),
+            ("Tank 2 pressure", f"{P2[0]/1e6:.6f} -> {P2[-1]/1e6:.6f} MPa"),
+            ("Tank 1 temperature", f"{T1[0]:.6f} -> {T1[-1]:.6f} K"),
+            ("Tank 2 temperature", f"{T2[0]:.6f} -> {T2[-1]:.6f} K"),
+        ]),
+        ("Pipe extrema", [
+            ("Max pressure", f"{P[idx_p]/1e6:.6f} MPa at t={t[idx_p[0]]:.6e} s, x={x[idx_p[1]]:.6f} m"),
+            ("Max velocity", f"{u.max():.6f} m/s"),
+            ("Min / max temperature", f"{T.min():.6f} / {T.max():.6f} K"),
+            ("Max Mach", f"{Ma[idx_ma]:.6f} at t={t[idx_ma[0]]:.6e} s, x={x[idx_ma[1]]:.6f} m"),
+            ("Final Mach range", f"{Ma[-1].min():.6f} -> {Ma[-1].max():.6f}"),
+            ("Final supersonic cells", f"{int(np.sum(Ma[-1] > 1.0))} / {Ma.shape[1]}"),
+        ]),
+        ("Conservation check", [
+            ("Total mass", f"{m_total_init:.12e} -> {m_total_final:.12e} kg (rel err {rel_err_m:.3e})"),
+            ("Total energy", f"{U_total_init:.12e} -> {U_total_final:.12e} J (rel err {rel_err_U:.3e})"),
+        ]),
+    ]
+
+    parts = [
+        "<!doctype html>",
+        "<html lang='en'>",
+        "<head>",
+        "<meta charset='utf-8'>",
+        "<title>Pipe system simulation summary</title>",
+        "<style>",
+        "body { font-family: Arial, sans-serif; margin: 24px; line-height: 1.45; }",
+        "h1, h2 { margin-bottom: 0.3em; }",
+        "table { border-collapse: collapse; width: 100%; margin: 12px 0 24px; }",
+        "th, td { border: 1px solid #ccc; padding: 8px 10px; text-align: left; vertical-align: top; }",
+        "th { background: #f5f5f5; width: 28%; }",
+        "img { max-width: 100%; height: auto; border: 1px solid #ddd; margin: 8px 0 24px; }",
+        "code { background: #f5f5f5; padding: 1px 4px; }",
+        "</style>",
+        "</head>",
+        "<body>",
+        "<h1>Pipe system simulation summary</h1>",
+        f"<p>Generated from <code>results/history.npz</code>. Final simulation time: {t[-1]:.6f} s.</p>",
+    ]
+
+    for title, rows in sections:
+        parts.append(f"<h2>{html.escape(title)}</h2>")
+        parts.append("<table>")
+        for key, value in rows:
+            parts.append(
+                f"<tr><th>{html.escape(str(key))}</th><td>{html.escape(str(value))}</td></tr>"
+            )
+        parts.append("</table>")
+
+    parts.extend([
+        "<h2>Figures</h2>",
+        "<p><img src='tank_pressure.png' alt='Tank pressure history'></p>",
+        "<p><img src='tank_temperature.png' alt='Tank temperature history'></p>",
+        "<p><img src='pipe_final_profiles.png' alt='Final pipe profiles'></p>",
+        "<p>Animation: <a href='pipe_animation.gif'>pipe_animation.gif</a></p>",
+        "</body>",
+        "</html>",
+    ])
+
+    dirname = os.path.dirname(path)
+    if dirname:
+        os.makedirs(dirname, exist_ok=True)
+    with open(path, "w", encoding="utf-8") as f:
+        f.write("\n".join(parts))

src/tank_pipe/pipe.py

@@ -0,0 +1,126 @@
+# src/tank_pipe/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 tank_pipe.riemann import hll_flux, get_riemann_solver
+from tank_pipe.friction import darcy_friction_factor
+
+
+class Pipe:
+    def __init__(self, L, D, N, P_init, T_init, gamma, R_gas,
+                 mu=0.0, roughness=0.0, riemann_solver="hll"):
+        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.mu = mu                    # dynamic viscosity [Pa·s]
+        self.roughness = roughness      # absolute wall roughness [m]
+        self.eps_D = roughness / D if D > 0 else 0.0  # relative roughness
+        self._flux_fn = get_riemann_solver(riemann_solver)
+        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] = self._flux_fn(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])
+
+        # --- Friction source term (operator splitting, explicit Euler) ---
+        # S = [0, -f/D * rho*u*|u|/2, 0]
+        # Energy source = 0 for adiabatic wall (KE dissipated → internal energy)
+        if self.mu > 0:
+            rho_s = self.W[0, :]
+            u_s = self.W[1, :] / rho_s
+            abs_u = np.abs(u_s)
+            Re = rho_s * abs_u * self.D / self.mu
+            f = darcy_friction_factor(Re, self.eps_D)
+            S_mom = -f / self.D * rho_s * u_s * abs_u / 2.0
+            self.W[1, :] += dt * S_mom
+
+        # 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()}"
+            )

src/tank_pipe/riemann.py

@@ -0,0 +1,206 @@
+# src/tank_pipe/riemann.py
+"""
+Riemann flux solvers 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)]
+
+Available solvers:
+  - hll_flux:  HLL (Harten-Lax-van Leer) two-wave approximate solver
+  - roe_flux:  Roe linearized solver with Harten-Hyman entropy fix
+"""
+
+import numpy as np
+
+
+# ---------------------------------------------------------------------------
+# Helper: recover primitives + physical flux from a conservative state
+# ---------------------------------------------------------------------------
+def _primitives(W, gamma):
+    """Return (rho, u, P, a, H) from conservative W = [rho, rho*u, rho*E]."""
+    rho = W[0]
+    if rho <= 0:
+        raise ValueError(f"Non-positive density: rho={rho}, W={W}")
+    u = W[1] / rho
+    E = W[2]
+    P = (gamma - 1) * (E - 0.5 * rho * u ** 2)
+    if P <= 0:
+        raise ValueError(f"Non-positive pressure: P={P}, W={W}")
+    a = np.sqrt(gamma * P / rho)
+    H = (E + P) / rho  # specific total enthalpy
+    return rho, u, P, a, H
+
+
+def _physical_flux(rho, u, P, E):
+    """Physical Euler flux from primitives + total energy density."""
+    return np.array([
+        rho * u,
+        rho * u ** 2 + P,
+        u * (E + P),
+    ])
+
+
+# ---------------------------------------------------------------------------
+# HLL solver
+# ---------------------------------------------------------------------------
+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.
+    """
+    rho_L, u_L, p_L, a_L, _ = _primitives(W_L, gamma)
+    rho_R, u_R, p_R, a_R, _ = _primitives(W_R, gamma)
+    E_L, E_R = W_L[2], W_R[2]
+
+    F_L = _physical_flux(rho_L, u_L, p_L, E_L)
+    F_R = _physical_flux(rho_R, u_R, p_R, E_R)
+
+    # --- HLL wave-speed estimates (Davis) ---
+    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)
+
+
+# ---------------------------------------------------------------------------
+# Roe solver with Harten-Hyman entropy fix
+# ---------------------------------------------------------------------------
+def roe_flux(W_L, W_R, gamma):
+    """
+    Compute the Roe linearized numerical flux with Harten-Hyman entropy fix.
+
+    The Roe solver resolves all three waves (left acoustic, contact/entropy,
+    right acoustic) and is more accurate than HLL at contact discontinuities.
+
+    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,)
+        Roe numerical flux vector.
+
+    Raises
+    ------
+    ValueError
+        If either state has non-positive density or pressure.
+    """
+    rho_L, u_L, p_L, a_L, H_L = _primitives(W_L, gamma)
+    rho_R, u_R, p_R, a_R, H_R = _primitives(W_R, gamma)
+
+    F_L = _physical_flux(rho_L, u_L, p_L, W_L[2])
+    F_R = _physical_flux(rho_R, u_R, p_R, W_R[2])
+
+    # --- Roe-averaged quantities (density-weighted) ---
+    sqrt_rL = np.sqrt(rho_L)
+    sqrt_rR = np.sqrt(rho_R)
+    denom = sqrt_rL + sqrt_rR
+
+    rho_hat = sqrt_rL * sqrt_rR           # geometric mean density
+    u_hat = (sqrt_rL * u_L + sqrt_rR * u_R) / denom
+    H_hat = (sqrt_rL * H_L + sqrt_rR * H_R) / denom
+    a_hat_sq = (gamma - 1) * (H_hat - 0.5 * u_hat ** 2)
+    if a_hat_sq <= 0:
+        return hll_flux(W_L, W_R, gamma)  # fallback
+    a_hat = np.sqrt(a_hat_sq)
+
+    # --- Eigenvalues of the Roe matrix ---
+    lam1 = u_hat - a_hat   # left acoustic
+    lam2 = u_hat            # entropy / contact
+    lam3 = u_hat + a_hat   # right acoustic
+
+    # --- Wave strengths (jump decomposition onto eigenvectors) ---
+    dp = p_R - p_L
+    du = u_R - u_L
+    drho = rho_R - rho_L
+
+    alpha_1 = (dp - rho_hat * a_hat * du) / (2.0 * a_hat ** 2)
+    alpha_2 = drho - dp / (a_hat ** 2)
+    alpha_3 = (dp + rho_hat * a_hat * du) / (2.0 * a_hat ** 2)
+
+    # --- Right eigenvectors ---
+    r1 = np.array([1.0, u_hat - a_hat, H_hat - u_hat * a_hat])
+    r2 = np.array([1.0, u_hat,          0.5 * u_hat ** 2])
+    r3 = np.array([1.0, u_hat + a_hat, H_hat + u_hat * a_hat])
+
+    # --- Harten-Hyman entropy fix ---
+    # Prevents unphysical expansion shocks at sonic points
+    eps1 = max(0.0, lam1 - (u_L - a_L), (u_R - a_R) - lam1)
+    eps3 = max(0.0, lam3 - (u_L + a_L), (u_R + a_R) - lam3)
+
+    abs_lam1 = abs(lam1)
+    abs_lam2 = abs(lam2)
+    abs_lam3 = abs(lam3)
+
+    if abs_lam1 < eps1:
+        abs_lam1 = (lam1 ** 2 + eps1 ** 2) / (2.0 * eps1)
+    if abs_lam3 < eps3:
+        abs_lam3 = (lam3 ** 2 + eps3 ** 2) / (2.0 * eps3)
+
+    # --- Roe flux: F = 0.5*(F_L + F_R) - 0.5 * sum(alpha_k |lam_k| r_k) ---
+    return 0.5 * (F_L + F_R) - 0.5 * (
+        alpha_1 * abs_lam1 * r1 +
+        alpha_2 * abs_lam2 * r2 +
+        alpha_3 * abs_lam3 * r3
+    )
+
+
+# ---------------------------------------------------------------------------
+# Dispatcher
+# ---------------------------------------------------------------------------
+_SOLVERS = {
+    'hll': hll_flux,
+    'roe': roe_flux,
+}
+
+
+def get_riemann_solver(name):
+    """
+    Return the Riemann flux function for the given solver name.
+
+    Parameters
+    ----------
+    name : str
+        Solver name: 'hll' or 'roe'.
+
+    Returns
+    -------
+    callable
+        A function with signature (W_L, W_R, gamma) -> np.ndarray(3,).
+    """
+    key = name.lower()
+    if key not in _SOLVERS:
+        raise ValueError(
+            f"Unknown Riemann solver '{name}'. Available: {list(_SOLVERS.keys())}"
+        )
+    return _SOLVERS[key]

src/tank_pipe/solver.py

@@ -0,0 +1,115 @@
+"""
+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
+
+
+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, same solver as pipe) ---
+        flux_L = pipe._flux_fn(W_ghost_L,     pipe.W[:, 0],  pipe.gamma)
+        flux_R = pipe._flux_fn(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

src/tank_pipe/tank.py

@@ -0,0 +1,82 @@
+# src/tank_pipe/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}"
+            )