diff --git a/src/venting/cli.py b/src/venting/cli.py index dd192bd..d383e8d 100644 --- a/src/venting/cli.py +++ b/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, diff --git a/src/venting/compare.py b/src/venting/compare.py index 69e4c0e..c6dd48c 100644 --- a/src/venting/compare.py +++ b/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), } diff --git a/src/venting/flow.py b/src/venting/flow.py index 8faeeae..33b294d 100644 --- a/src/venting/flow.py +++ b/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( diff --git a/src/venting/io.py b/src/venting/io.py index 56edb72..7e39e51 100644 --- a/src/venting/io.py +++ b/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: diff --git a/src/venting/profiles.py b/src/venting/profiles.py index 1be20f7..7ec0dd1 100644 --- a/src/venting/profiles.py +++ b/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: diff --git a/src/venting/solver.py b/src/venting/solver.py index d8dfe40..2e479fe 100644 --- a/src/venting/solver.py +++ b/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) diff --git a/src/venting/thermo.py b/src/venting/thermo.py index 440a3d7..e24d800 100644 --- a/src/venting/thermo.py +++ b/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: diff --git a/src/venting/validity.py b/src/venting/validity.py index 212b3fc..2538637 100644 --- a/src/venting/validity.py +++ b/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 = ( diff --git a/tests/test_solver_stream.py b/tests/test_solver_stream.py index 00df15a..b793fcb 100644 --- a/tests/test_solver_stream.py +++ b/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(): diff --git a/tests/test_v9_features.py b/tests/test_v9_features.py index f6c9f19..42a1b46 100644 --- a/tests/test_v9_features.py +++ b/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(): diff --git a/tests/test_validity_and_units.py b/tests/test_validity_and_units.py index 4ba8508..76d2bc6 100644 --- a/tests/test_validity_and_units.py +++ b/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):