fix: 11 verified bugs from codebase audit

src/venting/cli.py

@@ -287,13 +287,13 @@ def main() -> None:
         preset = get_default_panel_preset_v9()
         net_cfg = NetworkConfig(
             N_chain=10,
-            N_par=args.n_int,
+            N_par=2,
             V_cell=preset.V_cell,
             V_vest=preset.V_vest,
             A_wall_cell=preset.A_wall_cell,
             A_wall_vest=preset.A_wall_vest,
             d_int_mm=args.d_int,
-            n_int_per_interface=1,
+            n_int_per_interface=args.n_int,
             d_exit_mm=args.d_exit,
             n_exit=args.n_exit,
             Cd_int=args.cd_int,

src/venting/compare.py

@@ -39,10 +39,17 @@ def load_run(run_dir: Path) -> dict:
 
 def compare_runs(run_a: dict, run_b: dict) -> dict:
     """Compare run metrics and return structured differences."""
-    t_final_a = run_a["P"][:, -1]
-    t_final_b = run_b["P"][:, -1]
-    tA = run_a["T"][:, -1]
-    tB = run_b["T"][:, -1]
+    Pa = run_a["P"][:, -1]
+    Pb = run_b["P"][:, -1]
+    Ta = run_a["T"][:, -1]
+    Tb = run_b["T"][:, -1]
+
+    if Pa.shape == Pb.shape:
+        max_final_P_diff = float(np.max(np.abs(Pb - Pa)))
+        max_final_T_diff = float(np.max(np.abs(Tb - Ta)))
+    else:
+        max_final_P_diff = float("nan")
+        max_final_T_diff = float("nan")
 
     delta_summary = {}
     by_edge_a = {row["edge"]: row for row in run_a.get("summary", [])}
@@ -68,8 +75,8 @@ def compare_runs(run_a: dict, run_b: dict) -> dict:
         "run_a": str(run_a["dir"]),
         "run_b": str(run_b["dir"]),
         "edge_deltas": delta_summary,
-        "max_final_P_diff": float(np.max(np.abs(t_final_b - t_final_a))),
-        "max_final_T_diff": float(np.max(np.abs(tB - tA))),
+        "max_final_P_diff": max_final_P_diff,
+        "max_final_T_diff": max_final_T_diff,
         "rel_tau_exit_delta": float(rel_tau),
     }
 

src/venting/flow.py

@@ -65,13 +65,18 @@ def mdot_slot_pos(
 
 
 def friction_factor(Re: float, eps_over_D: float) -> float:
-    """Darcy friction factor (not Fanning)."""
+    """Darcy friction factor (not Fanning), with laminar-turbulent blending."""
     if Re <= 0.0:
         return 0.0
-    if Re < 2300.0:
-        return 64.0 / Re
+    f_lam = 64.0 / Re
+    if Re < 2000.0:
+        return f_lam
     term = eps_over_D / 3.7 + 5.74 / (Re**0.9)
-    return 0.25 / (math.log10(max(term, 1e-20)) ** 2)
+    f_turb = 0.25 / (math.log10(max(term, 1e-20)) ** 2)
+    if Re > 4000.0:
+        return f_turb
+    w = (Re - 2000.0) / 2000.0
+    return (1.0 - w) * f_lam + w * f_turb
 
 
 def mdot_short_tube_pos(

src/venting/io.py

@@ -32,7 +32,7 @@ def package_version() -> str:
     try:
         return version("venting")
     except PackageNotFoundError:
-        return "9.0.0"
+        return "10.0.0"
 
 
 def write_run_json(outdir: Path, params: dict, solver_settings: dict) -> None:

src/venting/profiles.py

@@ -1,4 +1,5 @@
 import math
+import warnings
 from collections.abc import Callable
 from dataclasses import dataclass
 from pathlib import Path
@@ -46,12 +47,16 @@ def p_fn(t: float) -> float:
 def make_profile_exponential(
     p0: float, rate0_mmhg_per_s: float, p_floor: float = 10.0
 ) -> Profile:
+    if rate0_mmhg_per_s <= 0.0:
+        raise ValueError("rate0_mmhg_per_s must be positive")
+    if p_floor <= 0.0:
+        raise ValueError("p_floor must be positive")
     rate0 = rate0_mmhg_per_s * 133.322
     tau = p0 / rate0
     t_floor = tau * math.log(p0 / p_floor)
 
     def p_fn(t: float) -> float:
-        return p0 * math.exp(-t / tau)
+        return max(p0 * math.exp(-t / tau), p_floor)
 
     return Profile(
         "barometric_exp",
@@ -84,9 +89,10 @@ def make_profile_from_table(
 
     p_max = float(np.max(p_tab))
     if pressure_unit.lower() == "pa" and 200.0 <= p_max <= 2000.0:
-        raise ValueError(
+        warnings.warn(
             "Profile table pressure looks like mmHg values provided as Pa. "
-            "Use pressure_unit='mmHg' or convert CSV to Pa."
+            "Use pressure_unit='mmHg' or convert CSV to Pa.",
+            stacklevel=2,
         )
 
     def p_fn(t: float) -> float:

src/venting/solver.py

@@ -72,7 +72,15 @@ def _property_model(thermo: str):
             lambda T: C_P * T,
             lambda T: C_V * T,
         )
-    raise ValueError("Property model valid for 'intermediate' or 'variable'")
+    if thermo == "isothermal":
+        raise ValueError(
+            "_property_model should not be called for 'isothermal' mode; "
+            "isothermal is handled separately in build_rhs"
+        )
+    raise ValueError(
+        f"Unknown thermo mode {thermo!r}; "
+        "expected 'isothermal', 'intermediate', or 'variable'"
+    )
 
 
 def build_rhs(
@@ -415,23 +423,8 @@ def rhs(t: float, y: np.ndarray) -> np.ndarray:
     )
 
 
-def solve_case(
-    nodes: list[GasNode], edges: list, bcs: list[ExternalBC], case: CaseConfig
-):
-    ctx = _prepare_solve(nodes, edges, bcs, case)
-    nodes_local = ctx["nodes_local"]
-    edges_local = ctx["edges_local"]
-    bcs_local = ctx["bcs_local"]
-    ext_idx = ctx["ext_idx"]
-    N = ctx["N"]
-    rhs = ctx["rhs"]
-    p_from_state = ctx["p_from_state"]
-    y0 = ctx["y0"]
-    max_step = ctx["max_step"]
-    t_eval = ctx["t_eval"]
-
-    profile = bcs_local[0].profile if bcs_local else None
-
+def _build_events(N, p_from_state, ext_idx, profile, case):
+    """Build terminal events shared by solve_case and solve_case_stream."""
     events = []
     for i in range(N):
 
@@ -469,6 +462,26 @@ def ev_p_stop(t, y):
     ev_p_stop.terminal = True
     ev_p_stop.direction = -1
     events.append(ev_p_stop)
+    return events
+
+
+def solve_case(
+    nodes: list[GasNode], edges: list, bcs: list[ExternalBC], case: CaseConfig
+):
+    ctx = _prepare_solve(nodes, edges, bcs, case)
+    nodes_local = ctx["nodes_local"]
+    edges_local = ctx["edges_local"]
+    bcs_local = ctx["bcs_local"]
+    ext_idx = ctx["ext_idx"]
+    N = ctx["N"]
+    rhs = ctx["rhs"]
+    p_from_state = ctx["p_from_state"]
+    y0 = ctx["y0"]
+    max_step = ctx["max_step"]
+    t_eval = ctx["t_eval"]
+
+    profile = bcs_local[0].profile if bcs_local else None
+    events = _build_events(N, p_from_state, ext_idx, profile, case)
 
     sol = solve_ivp(
         rhs,
@@ -512,10 +525,14 @@ def solve_case_stream(
     ext_idx = ctx["ext_idx"]
     N = ctx["N"]
     rhs = ctx["rhs"]
+    p_from_state = ctx["p_from_state"]
     t_eval = ctx["t_eval"]
     max_step = ctx["max_step"]
     y = ctx["y0"].copy()
 
+    profile = bcs_local[0].profile if bcs_local else None
+    events = _build_events(N, p_from_state, ext_idx, profile, case)
+
     if n_chunks is not None:
         effective_n_chunks = max(int(n_chunks), 1)
     else:
@@ -548,12 +565,20 @@ def solve_case_stream(
             rtol=1e-7 if case.thermo == "isothermal" else 1e-6,
             atol=1e-10 if case.thermo == "isothermal" else 1e-8,
             max_step=max_step,
+            events=events,
         )
         if not sol_seg.success:
             success = False
             message = str(sol_seg.message)
             break
 
+        terminated = False
+        if hasattr(sol_seg, "t_events"):
+            for ev_t in sol_seg.t_events:
+                if ev_t is not None and len(ev_t) > 0:
+                    terminated = True
+                    break
+
         seg_t = sol_seg.t
         seg_y = sol_seg.y
         if i > 0 and seg_t.size > 0:
@@ -572,7 +597,7 @@ def solve_case_stream(
                     "t": t_cat,
                     "y": y_cat,
                     "progress": float((i + 1) / effective_n_chunks),
-                    "done": bool(i == effective_n_chunks - 1),
+                    "done": bool(i == effective_n_chunks - 1 or terminated),
                     "node_count": N,
                     "n_nodes": N,
                     "thermo": case.thermo,
@@ -581,6 +606,10 @@ def solve_case_stream(
                 }
             )
 
+        if terminated:
+            message = "terminated by event"
+            break
+
     if t_out:
         t_final = np.concatenate(t_out)
         y_final = np.concatenate(y_out, axis=1)

src/venting/thermo.py

@@ -61,7 +61,7 @@ def h_air(T: float) -> float:
 
 def u_air(T: float) -> float:
     t = min(max(float(T), 1.0), T_FIT_HIGH)
-    return h_air(t) - R_GAS * (t - _T_REF)
+    return h_air(t) - R_GAS * t
 
 
 def speed_of_sound(T: float) -> float:

src/venting/validity.py

@@ -4,10 +4,11 @@
 
 import numpy as np
 
-from .constants import D_MOL_AIR, GAMMA, K_BOLTZMANN, R_GAS, T_SAFE
+from .constants import D_MOL_AIR, GAMMA, K_BOLTZMANN, R_GAS, T0, T_SAFE
 from .flow import (
     fanno_choked_state,
     friction_factor,
+    mdot_fanno_tube,
     mdot_short_tube_pos,
     mu_air_sutherland,
 )
@@ -152,18 +153,32 @@ def _short_tube_flag(
                 t_up = float(T[up_idx, k])
 
             t_up = max(t_up, T_SAFE)
-            md = mdot_short_tube_pos(
-                p_up,
-                t_up,
-                p_dn,
-                cd0,
-                e.A_total,
-                e.D,
-                e.L,
-                e.eps,
-                e.K_in,
-                e.K_out,
-            )
+            if e.fanno:
+                md = mdot_fanno_tube(
+                    p_up,
+                    t_up,
+                    p_dn,
+                    cd0,
+                    e.A_total,
+                    e.D,
+                    e.L,
+                    e.eps,
+                    e.K_in,
+                    e.K_out,
+                )
+            else:
+                md = mdot_short_tube_pos(
+                    p_up,
+                    t_up,
+                    p_dn,
+                    cd0,
+                    e.A_total,
+                    e.D,
+                    e.L,
+                    e.eps,
+                    e.K_in,
+                    e.K_out,
+                )
             rho = p_up / (R_GAS * t_up)
             u = md / max(rho * e.A_total, 1e-18)
             a_s = math.sqrt(GAMMA * R_GAS * t_up)
@@ -256,7 +271,11 @@ def _knudsen_flag(
                 t_up = float(T[a, k])
             else:
                 p_up = pb
-                t_up = float(T[a, k] if b == EXT_NODE else T[b, k])
+                if b == EXT_NODE:
+                    # External atmosphere T not in T array; use T0 as estimate
+                    t_up = float(T0)
+                else:
+                    t_up = float(T[b, k])
             p_up = max(p_up, 1e-9)
             t_up = max(t_up, T_SAFE)
             mfp = (

tests/test_solver_stream.py

@@ -77,7 +77,8 @@ def cb(payload):
         seen.append(payload["progress"])
 
     solve_case_stream(nodes, edges, bcs, case, callback=cb, dt_chunk_s=2.0)
-    assert 4 <= len(seen) <= 6
+    # May terminate early due to equilibrium/p_stop events (3-5 chunks typical)
+    assert 1 <= len(seen) <= 6
 
 
 def test_dt_chunk_small_duration():

tests/test_v9_features.py

@@ -23,8 +23,8 @@ def _single_node():
 
 def test_variable_thermo_regression_near_300k():
     nodes, edges, bcs = _single_node()
-    c_int = CaseConfig("intermediate", 0.0, T0, 0.005, 40)
-    c_var = CaseConfig("variable", 0.0, T0, 0.005, 40)
+    c_int = CaseConfig("intermediate", 0.0, T0, 0.3, 100)
+    c_var = CaseConfig("variable", 0.0, T0, 0.3, 100)
     r_int = summarize_result(
         nodes, edges, bcs, c_int, solve_case(nodes, edges, bcs, c_int)
     )
@@ -33,6 +33,10 @@ def test_variable_thermo_regression_near_300k():
     )
     rel = np.max(np.abs(r_var.P[0] - r_int.P[0]) / np.maximum(r_int.P[0], 1.0))
     assert float(rel) <= 0.01
+    rel_T = np.max(
+        np.abs(r_var.T[0] - r_int.T[0]) / np.maximum(np.abs(r_int.T[0]), 1.0)
+    )
+    assert float(rel_T) <= 0.05
 
 
 def test_variable_gamma_changes_mdot():

tests/test_validity_and_units.py

@@ -11,11 +11,12 @@
 from venting.solver import solve_case
 
 
-def test_profile_table_rejects_mmhg_looking_values_when_unit_pa(tmp_path: Path):
+def test_profile_table_warns_mmhg_looking_values_when_unit_pa(tmp_path: Path):
     p = tmp_path / "profile.csv"
     p.write_text("0,760\n1,740\n", encoding="utf-8")
-    with pytest.raises(ValueError, match="looks like mmHg"):
-        make_profile_from_table("tab", p, pressure_unit="Pa")
+    with pytest.warns(UserWarning, match="looks like mmHg"):
+        prof = make_profile_from_table("tab", p, pressure_unit="Pa")
+    assert prof.P(0.0) > 0
 
 
 def test_profile_table_accepts_mmhg_when_explicit(tmp_path: Path):