From cd1f087828da19323aed9045dce3b4bfc048b5d9 Mon Sep 17 00:00:00 2001 From: test1card Date: Sat, 28 Feb 2026 01:09:20 +0300 Subject: [PATCH] Fix Fanno K-loss coupling and add regression coverage --- .github/workflows/ci.yml | 30 + .gitignore | 10 + .pre-commit-config.yaml | 10 + CHANGELOG.md | 34 + CONTRIBUTING.md | 19 + README.md | 284 +++++- archive/venting_v84.py | 1189 +++++++++++++++++++++++++ docs/assumptions_v85.md | 9 + docs/assumptions_v9.md | 10 + docs/model_v85.md | 17 + docs/verification_v85.md | 44 + pyproject.toml | 54 ++ requirements-dev.txt | 1 + requirements.txt | 3 + src/venting/__init__.py | 5 + src/venting/__main__.py | 4 + src/venting/cases.py | 82 ++ src/venting/cli.py | 365 ++++++++ src/venting/compare.py | 107 +++ src/venting/constants.py | 22 + src/venting/diagnostics.py | 179 ++++ src/venting/flow.py | 375 ++++++++ src/venting/gates.py | 106 +++ src/venting/geometry.py | 23 + src/venting/graph.py | 347 ++++++++ src/venting/gui/__init__.py | 1 + src/venting/gui/app.py | 400 +++++++++ src/venting/gui/config.py | 176 ++++ src/venting/gui/main.py | 11 + src/venting/gui/state_layout.py | 11 + src/venting/io.py | 108 +++ src/venting/montecarlo.py | 71 ++ src/venting/plotting.py | 32 + src/venting/presets.py | 35 + src/venting/profiles.py | 100 +++ src/venting/run.py | 57 ++ src/venting/solver.py | 658 ++++++++++++++ src/venting/state_layout.py | 34 + src/venting/thermo.py | 69 ++ src/venting/validity.py | 341 +++++++ tests/test_compare.py | 54 ++ tests/test_gates.py | 120 +++ tests/test_gui_config.py | 31 + tests/test_gui_import.py | 6 + tests/test_gui_state_layout.py | 31 + tests/test_jac_sparsity.py | 35 + tests/test_montecarlo.py | 46 + tests/test_regression_cli_baseline.py | 11 + tests/test_short_tube.py | 217 +++++ tests/test_solver_stream.py | 92 ++ tests/test_state_layout.py | 46 + tests/test_thermo_reference.py | 21 + tests/test_topology_and_scalars.py | 70 ++ tests/test_v9_features.py | 113 +++ tests/test_validity_and_units.py | 59 ++ 55 files changed, 6384 insertions(+), 1 deletion(-) create mode 100644 .github/workflows/ci.yml create mode 100644 .gitignore create mode 100644 .pre-commit-config.yaml create mode 100644 CHANGELOG.md create mode 100644 CONTRIBUTING.md create mode 100644 archive/venting_v84.py create mode 100644 docs/assumptions_v85.md create mode 100644 docs/assumptions_v9.md create mode 100644 docs/model_v85.md create mode 100644 docs/verification_v85.md create mode 100644 pyproject.toml create mode 100644 requirements-dev.txt create mode 100644 requirements.txt create mode 100644 src/venting/__init__.py create mode 100644 src/venting/__main__.py create mode 100644 src/venting/cases.py create mode 100644 src/venting/cli.py create mode 100644 src/venting/compare.py create mode 100644 src/venting/constants.py create mode 100644 src/venting/diagnostics.py create mode 100644 src/venting/flow.py create mode 100644 src/venting/gates.py create mode 100644 src/venting/geometry.py create mode 100644 src/venting/graph.py create mode 100644 src/venting/gui/__init__.py create mode 100644 src/venting/gui/app.py create mode 100644 src/venting/gui/config.py create mode 100644 src/venting/gui/main.py create mode 100644 src/venting/gui/state_layout.py create mode 100644 src/venting/io.py create mode 100644 src/venting/montecarlo.py create mode 100644 src/venting/plotting.py create mode 100644 src/venting/presets.py create mode 100644 src/venting/profiles.py create mode 100644 src/venting/run.py create mode 100644 src/venting/solver.py create mode 100644 src/venting/state_layout.py create mode 100644 src/venting/thermo.py create mode 100644 src/venting/validity.py create mode 100644 tests/test_compare.py create mode 100644 tests/test_gates.py create mode 100644 tests/test_gui_config.py create mode 100644 tests/test_gui_import.py create mode 100644 tests/test_gui_state_layout.py create mode 100644 tests/test_jac_sparsity.py create mode 100644 tests/test_montecarlo.py create mode 100644 tests/test_regression_cli_baseline.py create mode 100644 tests/test_short_tube.py create mode 100644 tests/test_solver_stream.py create mode 100644 tests/test_state_layout.py create mode 100644 tests/test_thermo_reference.py create mode 100644 tests/test_topology_and_scalars.py create mode 100644 tests/test_v9_features.py create mode 100644 tests/test_validity_and_units.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..ee4c7e1 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,30 @@ +name: ci + +on: + push: + pull_request: + +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10", "3.11", "3.12"] + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + cache: pip + - name: Install package + dev deps + run: | + python -m pip install --upgrade pip + pip install -e ".[dev]" + - name: pip check + run: python -m pip check + - name: Ruff + run: python -m ruff check . + - name: Black + run: python -m black --check . + - name: Pytest + run: pytest -q --cov=venting --cov-report=term-missing diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b80ac4c --- /dev/null +++ b/.gitignore @@ -0,0 +1,10 @@ +__pycache__/ +*.pyc +.venv/ +.pytest_cache/ +.ruff_cache/ +.mypy_cache/ +results/ +*.npz +*_meta.json +*.egg-info/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..f01b466 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,10 @@ +repos: + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.8.6 + hooks: + - id: ruff + - id: ruff-format + - repo: https://github.com/psf/black + rev: 25.1.0 + hooks: + - id: black diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..81652aa --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,34 @@ +# Changelog + +## v10.0.0 +- Fixed Radau setup by removing invalid identity-only Jacobian sparsity hint. +- Added scalar-or-list network parameter support and a new `two_chain_shared_vest` topology. +- Updated GUI worker/progress handling and stop-check integration for true streaming callbacks. +- Added regression tests for baseline gate behavior and topology/scalar-list compatibility. +- Added desktop GUI (`python -m venting gui`) with optional dependencies (`.[gui]`) for configuring and running venting cases. +- Added live plotting tabs for pressure, temperature, mass, peak diagnostics, and validity flags while solving in a background thread. +- Added JSON save/load schema for GUI cases with explicit units and validation. +- Added shared presets and run pipeline modules to keep CLI and GUI execution paths consistent. +- Added streaming solver wrapper for progressive updates with regression tests against batch solve. + +## v9.0.0 +- Added physically consistent short-tube thick-wall edge model with Darcy friction, minor losses, and iterative `Cd_eff` composition (lossy-nozzle, not Fanno). +- Added split internal/exit short-tube loss parameters (`K_in_*`, `K_out_*`, `eps_*`) with backward-compatible alias CLI flags. +- Added variable thermodynamics mode (`--thermo variable`) with temperature-dependent `cp/cv/gamma/h/u` and gamma-aware compressible discharge. +- Added lumped wall thermal model (`--wall-model lumped`) with wall heat capacity, optional outside convection/radiation/source heat flux. +- Added dynamic external pump model (`--external-model dynamic_pump`) with finite external volume and ultimate-pressure pump sink law. +- Expanded validity diagnostics with `short_tube_flow` metrics (`Mach_max`, `Re_max`, `K_tot_max`, `Cd_eff_min`, `L_over_D`, `frac_fric`), and always save validity as `*_validity.json` plus `meta.json`. +- Updated packaging metadata to 9.0.0 and refreshed docs for v9 assumptions. + +## v0.8.5 - version alignment +- Bumped active package/model naming from v8.4 to v8.5 (0.8.5) to avoid confusion. +- Updated CLI artifact prefix to `v85_` for new runs. +- Updated documentation references to the active v8.5 naming. + +## v0.8.4 - repo hardening +- Hardened packaging for standard PEP 517/518 src-layout installation. +- Set v8.4 as the only active runtime package; legacy script remains archived only. +- Added CI user-like install flow (`pip install -e ".[dev]"`, `pip check`, lint, coverage tests). +- Reworked tests into deterministic gate/physics validation suite (`tests/test_gates.py`). +- Added structured run outputs (`run.json`, `summary.csv`) through `venting.io`. +- Updated docs for verification equations, sign conventions, and acceptance criteria. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..0688a89 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,19 @@ +# Contributing + +1. Create a virtual environment and install project + dev extras: + +```bash +python -m venv .venv && source .venv/bin/activate +pip install -e ".[dev]" +``` + +2. Run quality checks: + +```bash +ruff check . +black --check . +pytest -q --cov=venting --cov-report=term-missing +``` + +3. Keep only v8.5 active in runtime package (`src/venting`). + Any legacy files must remain in `archive/` and not be imported by CLI/tests. diff --git a/README.md b/README.md index e8a41be..bd97c32 100644 --- a/README.md +++ b/README.md @@ -1 +1,283 @@ -# venting \ No newline at end of file +# Venting v10.0.0 — 0D/Network model of depressurization (rigid volumes + compressible discharge) + +Этот репозиторий моделирует стравливание (depressurization / venting) системы жёстких объёмов газа, соединённых отверстиями/каналами, когда внешнее давление падает по заданному профилю `P_ext(t)` (`linear`, `step`, `barometric`, `table` из CSV). + +Главный инженерный вопрос: какой максимальный перепад давления `|ΔP|` возникает на каждом «узком месте» (выход наружу, межузловые интерфейсы), и в каком режиме течёт газ (choked/subsonic). + +**Политика версий:** в runtime активна только **v10.0.0**. Старый монолитный скрипт хранится только в `archive/` как архивный референс. + +--- + +## 0) Интуиция для человека вне сферы + +Есть несколько полостей с воздухом, соединённых отверстиями. Снаружи давление падает. Если наружное давление падает быстрее, чем газ успевает выйти, внутри какое-то время остаётся более высокое давление — появляется перепад давления на стенках/перегородках. + +Типовая топология в проекте: + +```text +N_parallel chains (aggregated): + +cell10 -> cell9 -> ... -> cell2 -> cell1 + \ | + \______________________________| (по N_parallel идентичных ветвей) + v + [ vestibule ] --(exit orifice/slot)--> [ external P_ext(t) ] +``` + +Практически это нужно, чтобы оценить: +- где и когда пик `|ΔP|`, +- какие отверстия/`C_d` ограничивают продувку, +- насколько тепловой режим влияет на пик. + +--- + +## 1) Что такое 0D/сетевая модель и чем она **не** является + +### Что модель делает +- Каждый объём = один узел (`node`) с однородными `m(t), T(t), P(t)`. +- Связи между узлами = рёбра (`edges`) с моделью расхода. +- Модель собирается «по рёбрам»: вклад каждого ребра суммируется в балансы узлов. + +### Что модель **не** делает +- Это не CFD: нет пространственных полей скорости/температуры. +- Нет акустики/ударных волн/3D-деталей струи. +- `C_d` задаётся как параметр (главный источник неопределённости). + +--- + +## 2) Переменные и уравнение состояния + +Для узла `i`: +- объём `V_i` (жёсткий), +- масса `m_i(t)`, +- температура `T_i(t)` (в режиме `intermediate`), +- давление из идеального газа: + +$$ +P_i V_i = m_i R T_i +$$ + +В v10.0.0 интегрируются `m` (и `T` в thermal-варианте), а давление вычисляется через EOS. + +--- + +## 3) Массовый баланс + +Для каждого узла: + +$$ +\frac{dm_i}{dt}=\sum \dot m_{in}-\sum \dot m_{out} +$$ + +Направление потока на ребре определяется по текущему давлению upstream/downstream. + +--- + +## 4) Энергия и температура (почему blowdown охлаждает газ) + +Для жёсткого контрольного объёма: + +$$ +\frac{d}{dt}(m c_v T)=\sum \dot m_{in} c_p T_{in}-\sum \dot m_{out} c_p T + \dot Q_{wall} +$$ + +Эквивалентная форма, используемая в коде: + +$$ +m c_v \frac{dT}{dt} = \sum \dot m_{in}(c_p T_{in}-c_vT) - \sum \dot m_{out}(RT) + hA_{wall}(T_{wall}-T) +$$ + +Режимы: +- `isothermal`: `T = const`, решается только масса. +- `intermediate`: решаются `m(t), T(t)`; при `h→0` поведение стремится к адиабатическому blowdown, при `h→∞` — к изотерме. + +> Важно: в формулах расхода и энергетики используется **upstream temperature** из состояния узла. + +--- + +## 5) Расход через отверстия, choking и отношение давлений + +Для ребра: + +$$ +r = \frac{P_{down}}{P_{up}}, \quad r_* = \left(\frac{2}{\gamma+1}\right)^{\gamma/(\gamma-1)} +$$ + +- если `r <= r*` → **choked**, +- иначе → **subsonic**. + +Это критично для пиков `|ΔP|`: choked-участок часто становится «бутылочным горлом» по массовому расходу. + +--- + +## 6) Архитектура кода (активная v10.0.0) + +Пакет: `src/venting/` + +- `constants.py` — термоконстанты, safety-пороги, критическое отношение давлений. +- `geometry.py` — конвертация единиц, геометрические helper'ы. +- `profiles.py` — профили `P_ext(t)` (`linear/step/barometric/table`) и события profile-breakpoints. +- `graph.py` — узлы/рёбра/BC, построение сети `build_branching_network`. +- `flow.py` — расход через orifice и slot. +- `solver.py` — ODE RHS, `solve_ivp(method="Radau")`, события остановки, sparsity. +- `diagnostics.py` — `ΔP`, пики, режимы, `tau_exit`, meta. +- `gates.py` — встроенные gate checks (single-node / two-node). +- `io.py` — единообразная запись `run.json`, `summary.csv`, meta. +- `plotting.py` — опциональные графики. +- `cli.py` / `__main__.py` — CLI и точка входа `python -m venting`. + +--- + +## 7) Численный решатель + +- используется `scipy.integrate.solve_ivp`, метод **Radau** (жёсткие ODE), +- учитывается sparsity-структура, +- используются события ранней остановки (например, близость к внешнему давлению/низкое давление), +- **state clipping для P/T не используется**; есть только safety-защита знаменателей (`T_SAFE`, малые массы в делении). + +--- + +## 8) Validation / Gate tests + +В проекте есть детерминированные тесты `tests/test_gates.py`: + +1. `test_single_node_analytic_match_adiabatic()` + - один объём, выход в вакуум, `h=0`, + - сравнение с аналитикой adiabatic blowdown, + - маска `P > 0.01*P0`, точность < 0.5%. + +2. `test_single_node_analytic_match_isothermal_limit()` + - `h=1e6`, сравнение с изотермической экспонентой, + - точность < 0.5%. + +3. `test_mass_conservation_single_node()` + - интеграл расхода согласован с `m0 - mf` (<0.1%). + +4. `test_two_node_mass_conservation()` + - масса по сети (<0.1%), + - energy check для `h=0` (<1%). + +5. `test_monotonic_pressure_when_vacuum()` + - при `P_ext=0` давление не должно расти (кроме крошечного численного шума). + +Запуск: + +```bash +pytest -q +# или +python -m venting gate +python -m venting gate --single +python -m venting gate --two +``` + +--- + +## 9) Установка и запуск CLI + +### Установка (рекомендуемо) + +```bash +python -m venv .venv && source .venv/bin/activate +pip install -e ".[dev]" +``` + +Альтернатива без dev-инструментов: + +```bash +pip install -e . +``` + +### Основные команды + +```bash +python -m venting --help +python -m venting gate +python -m venting sweep --profile linear --d-int 2 --d-exit 2 --cd-int 0.62 --cd-exit 0.62 +python -m venting thermal --profile linear --d 2 --h-list 0,1,5,15 +python -m venting sweep2d --profile linear --d-int-list 1.5,2.0 --d-exit-list 1.5,2.0 +python -m venting sweep --profile linear --d-int 2 --d-exit 4 --int-model short_tube --L-int-mm 1 --exit-model short_tube --L-exit-mm 1 --K-in-int 0.5 --K-out-int 1.0 --eps-int-um 0 --K-in-exit 0.5 --K-out-exit 1.0 --eps-exit-um 0 --cd-int 0.62 --cd-exit 0.62 --thermo variable --wall-model lumped --external-model dynamic_pump --pump-speed-m3s 0.01 --V-ext 0.2 --P-ult-Pa 10 +``` + +Профиль `table` ожидает CSV с колонками `t_s,P_Pa` через `--profile-file`. + +--- + +## 10) How to interpret outputs + +Каждый запуск пишет в: + +```text +results/_/ +``` + +Типовые артефакты: +- `run.json` — воспроизводимость (timestamp, commit hash, python/platform, параметры, solver settings). +- `summary.csv` — ключевые метрики по рёбрам (`max_abs_dP_Pa`, `t_peak_s`, `r_peak`, `regime`, `peak_type`, `tau_exit_s`). +- `*.npz` — временные ряды (`t, m, T, P, P_ext, tau_exit`). +- `*_meta.json` — доп. метаданные. + +Интерпретация ключевых полей: +- `r_peak = P_down/P_up` в момент пика, +- `regime`: `CHOKED` или `subsonic` (для slot — `viscous_slot`), +- `peak_type`: + - `boundary(...)` — пик привязан к событию профиля `P_ext(t)`, + - `internal` — пик рождается внутренней динамикой сети. +- `t_peak/tau_exit` (в meta/diagnostics) удобно для сравнения разных кейсов. + +--- + +## 11) Limitations & red flags + +Дополнительно для short-tube: это **lossy-nozzle** через эффективный `Cd_eff`; Fanno/friction-choking в v10.0.0 не реализован. + + + +- Это **не** CFD и не замена стендовым испытаниям. +- `C_d` — главный источник неопределённости → обязательно делать sweep по `C_d`/геометрии. +- Изотермический случай не «всегда консервативен» для любой метрики и любого момента времени. +- Если видите странные осцилляции/нефизичные тренды — сначала проверяйте gate tests. + +--- + +## Дополнительные документы + +- `docs/model_v85.md` — краткая модель. +- `docs/verification_v85.md` — верификация и критерии. +- `docs/assumptions_v9.md` — допущения и ограничения. +- `CONTRIBUTING.md` — вклад в проект. +- `CHANGELOG.md` — история изменений. + +## License + +MIT (см. `LICENSE`). + + +## GUI + +V10 adds an optional desktop GUI (thin client over the same core physics modules). + +Install GUI extras: + +```bash +pip install -e ".[gui]" +``` + +Run GUI: + +```bash +python -m venting gui +``` + +What GUI supports: +- orifice vs short_tube for internal/exit edges (L, eps, K_in/K_out), +- network and geometry controls (`N_chain`, `N_par`, volumes, wall areas, diameters/counts, Cd), +- external model (`profile` / `dynamic_pump`), +- thermo (`isothermal` / `intermediate` / `variable`) and wall model (`fixed` / `lumped`), +- background solve with live plot updates, stop button, validity table, +- case JSON save/load (explicit units), +- artifact export compatible with CLI (`npz`, `meta.json`, `*_validity.json`, `summary.csv`). +- GUI progress uses chunked integration when streaming is enabled; tiny differences vs batch solve can appear within integrator tolerances. + +Packaging executable for a target OS is left as TODO until `{TARGET_OS}` and `{PACKAGING_TOOL}` are specified. + +GUI streaming uses chunked integration for progress updates; results can differ slightly from a single batch solve within integrator tolerances. diff --git a/archive/venting_v84.py b/archive/venting_v84.py new file mode 100644 index 0000000..cf23f2a --- /dev/null +++ b/archive/venting_v84.py @@ -0,0 +1,1189 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +venting_v84.py — v8.4: near-production 0D venting solver (edge-based graph) + +Core upgrades vs v8.3: + A) State = (m_i, T_i) for intermediate; state = (m_i) for isothermal + => P_i = m_i*R*T_i/V_i is always >= 0 if m_i >= 0 (constructive positivity). + B) No state clipping. Only denominator protection (T >= T_SAFE) and event-based stop. + C) Symmetry handled by aggregation (V, A_wall, A_or_total), not by multipliers. + D) Profile = object with callable P(t) + explicit event times (breakpoints). + Envelope MUST be provided as a table to avoid invented parameters. + E) Edge types: + - OrificeEdge: compressible isentropic (choked/subsonic) + - SlotChannelEdge: viscous laminar slot/channel using (P_up^2 - P_dn^2)/(R*T) + F) Reproducibility: save NPZ + JSON metadata; optional plots. + G) Gate tests: single-node analytic; two-node mass+energy; network smoke. + +This is still a 0D model. Its “10/10” is about internal consistency, rigor, +and reproducibility, not about eliminating model-form uncertainty (Cd, geometry, etc.). +""" + +from __future__ import annotations + +import argparse +import json +import math +from dataclasses import dataclass, asdict +from pathlib import Path +from typing import Callable, List, Optional, Sequence, Tuple, Dict, Union + +import numpy as np +from scipy.integrate import solve_ivp +from scipy.sparse import lil_matrix + +# matplotlib is optional (plots) +try: + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + HAS_MPL = True +except Exception: + HAS_MPL = False + + +# ============================================================================= +# CONSTANTS (Air, ideal gas) +# ============================================================================= +GAMMA = 1.4 +R_GAS = 287.05 # J/(kg*K) +C_V = R_GAS / (GAMMA - 1.0) +C_P = GAMMA * C_V + +T0 = 300.0 # K +P0 = 101325.0 # Pa + +# Numerical safety (NOT state clipping; only denominators + stop conditions) +T_SAFE = 1.0 # K (protect sqrt(T), viscosity, etc.) +M_SAFE = 1e-18 # kg (protect division in dT) +P_STOP = 5.0 # Pa stop when all nodes effectively evacuated (configurable) + +# Critical pressure ratio for choked flow +PI_C = (2.0/(GAMMA+1.0))**(GAMMA/(GAMMA-1.0)) +C_CHOKED = math.sqrt(GAMMA * (2.0/(GAMMA+1.0))**((GAMMA+1.0)/(GAMMA-1.0))) + + +# ============================================================================= +# Utility: Units and geometry +# ============================================================================= +def mm_to_m(x_mm: float) -> float: + return x_mm * 1e-3 + +def mm2_to_m2(x_mm2: float) -> float: + return x_mm2 * 1e-6 + +def mm3_to_m3(x_mm3: float) -> float: + return x_mm3 * 1e-9 + +def circle_area_from_d_mm(d_mm: float) -> float: + d_m = mm_to_m(d_mm) + return math.pi * (d_m/2.0)**2 + +def assert_pos(name: str, val: float) -> None: + if not (val > 0.0): + raise ValueError(f"{name} must be > 0, got {val}") + +def assert_nonneg(name: str, val: float) -> None: + if not (val >= 0.0): + raise ValueError(f"{name} must be >= 0, got {val}") + + +# ============================================================================= +# Viscosity model (Sutherland) for air — needed for viscous channel edges +# ============================================================================= +def mu_air_sutherland(T: float) -> float: + """ + Dynamic viscosity of air [Pa*s] using Sutherland's law. + Typical constants: + mu0 = 1.716e-5 Pa*s at T0=273.15 K + S = 111 K + """ + T_eff = max(T, T_SAFE) + mu0 = 1.716e-5 + Tref = 273.15 + S = 111.0 + return mu0 * (T_eff/Tref)**1.5 * (Tref + S)/(T_eff + S) + + +# ============================================================================= +# Profile object with explicit event times +# ============================================================================= +@dataclass(frozen=True) +class Profile: + name: str + P: Callable[[float], float] + events: Tuple[Tuple[float, str], ...] # (t_event, label) + + def P_array(self, t: np.ndarray) -> np.ndarray: + return np.array([self.P(float(tt)) for tt in t], dtype=float) + + def classify_peak(self, t_peak: float, t_end: float, tol_s: float = 0.5) -> str: + for te, label in self.events: + if abs(t_peak - te) <= tol_s: + return f"boundary({label})" + if t_peak <= tol_s: + return "boundary(sim_start)" + if t_peak >= (t_end - tol_s): + return "boundary(sim_end)" + return "internal" + + +def make_profile_linear(p0: float, rate_mmhg_per_s: float) -> Profile: + rate = rate_mmhg_per_s * 133.322 # Pa/s + t_zero = p0 / rate + + def P(t: float) -> float: + return max(p0 - rate*t, 0.0) + + return Profile("linear", P, events=((t_zero, "P_ext=0"),)) + + +def make_profile_step(p0: float, step_time_s: float) -> Profile: + def P(t: float) -> float: + return 0.0 if t >= step_time_s else p0 + + return Profile("step", P, events=((step_time_s, "step_to_vacuum"),)) + + +def make_profile_exponential(p0: float, rate0_mmhg_per_s: float, p_floor: float = 10.0) -> Profile: + """ + Exponential 'barometric-like' decay: + P(t) = p0 * exp(-t/tau) + Choose tau such that initial slope matches -rate0. + dP/dt|0 = -p0/tau = -rate0 => tau = p0/rate0 + """ + rate0 = rate0_mmhg_per_s * 133.322 + tau = p0 / rate0 + t_floor = tau * math.log(p0 / p_floor) + + def P(t: float) -> float: + return p0 * math.exp(-t/tau) + + return Profile("barometric_exp", P, events=((0.0, "start_max_slope"), (t_floor, f"P_ext={p_floor:.0f}Pa"),)) + + +def make_profile_from_table(name: str, table_path: Path) -> Profile: + """ + Table format: CSV with two columns: t_s, P_Pa + Piecewise-linear interpolation. Events = all breakpoints. + """ + arr = np.loadtxt(str(table_path), delimiter=",") + if arr.ndim != 2 or arr.shape[1] < 2: + raise ValueError("Profile table must be CSV with columns: t_s, P_Pa") + t_tab = np.array(arr[:, 0], dtype=float) + p_tab = np.array(arr[:, 1], dtype=float) + if not np.all(np.diff(t_tab) > 0): + raise ValueError("Profile table time must be strictly increasing") + + def P(t: float) -> float: + if t <= t_tab[0]: + return float(p_tab[0]) + if t >= t_tab[-1]: + return float(p_tab[-1]) + return float(np.interp(t, t_tab, p_tab)) + + events = tuple((float(tt), f"bp{i}") for i, tt in enumerate(t_tab)) + return Profile(name, P, events=events) + + +# ============================================================================= +# Nodes and edges +# ============================================================================= +@dataclass(frozen=True) +class GasNode: + name: str + V: float # m^3 + A_wall: float # m^2 for convection (if enabled) + +@dataclass(frozen=True) +class ExternalBC: + """ + Boundary condition on a node: external pressure P_ext(t). + T_ext is only used if flow reverses (external -> node). + """ + node: int + profile: Profile + T_ext: float = T0 + +# --- Cd model (minimal but explicit) --- +@dataclass(frozen=True) +class CdConst: + Cd: float + def __call__(self, Re: Optional[float] = None, r: Optional[float] = None) -> float: + return float(self.Cd) + +# You can add CdTable later if you have calibration data: +# interpolate Cd(Re) from file. Not included here to avoid fake correlations. + + +@dataclass(frozen=True) +class OrificeEdge: + """ + Compressible orifice edge between nodes a <-> b. If b == EXT_NODE, uses BC. + Area is TOTAL area (already includes N_holes multiplier). + """ + a: int + b: int + A_total: float # m^2 + Cd_model: CdConst + label: str = "" + +EXT_NODE = -1 + + +@dataclass(frozen=True) +class SlotChannelEdge: + """ + Viscous laminar slot/channel edge (for 'gap' / pocket / restrictive path). + Uses compressible isothermal-like laminar formula: + + mdot = K * (P_up^2 - P_dn^2) / (R * T_up) + + where K = w * delta^3 / (12 * mu(T) * L) + + This is a *model-form* choice. It is at least dimensionally consistent and + avoids mixing a linear conductance with a choked conductance improperly. + """ + a: int + b: int + w: float # m + delta: float # m + L: float # m + label: str = "" + + +Edge = Union[OrificeEdge, SlotChannelEdge] + + +# ============================================================================= +# Physics: mass flow for OrificeEdge +# ============================================================================= +def mdot_orifice_pos(P_up: float, T_up: float, P_dn: float, Cd: float, A: float) -> float: + """ + Returns mdot >= 0 from upstream to downstream. + Isentropic orifice with choked/subsonic transition. Ideal gas. + """ + if P_up <= 0.0 or A <= 0.0: + return 0.0 + T_eff = max(T_up, T_SAFE) + r = max(P_dn, 0.0) / P_up + if r >= 1.0: + return 0.0 + + if r <= PI_C: + return Cd * A * P_up * C_CHOKED / math.sqrt(R_GAS * T_eff) + + bracket = r**(2.0/GAMMA) - r**((GAMMA+1.0)/GAMMA) + if bracket <= 0.0: + return 0.0 + return Cd * A * P_up * math.sqrt(2.0*GAMMA/((GAMMA-1.0)*R_GAS*T_eff) * bracket) + + +def mdot_slot_pos(P_up: float, T_up: float, P_dn: float, w: float, delta: float, L: float) -> float: + """ + Compressible laminar slot: mdot >= 0 from up to dn. + mdot = K*(P_up^2 - P_dn^2)/(R*T_up), K = w*delta^3/(12*mu*L) + """ + if P_up <= 0.0: + return 0.0 + T_eff = max(T_up, T_SAFE) + mu = mu_air_sutherland(T_eff) + K = w * (delta**3) / (12.0 * mu * L) + dp2 = max(P_up**2 - max(P_dn, 0.0)**2, 0.0) + return K * dp2 / (R_GAS * T_eff) + + +# ============================================================================= +# Model configuration and results containers +# ============================================================================= +@dataclass(frozen=True) +class CaseConfig: + thermo: str # "isothermal" or "intermediate" + h_conv: float # W/(m^2*K), used only if thermo="intermediate" + T_wall: float # K (constant wall temperature model) + + duration: float # s + n_pts: int # time samples for output + + # stop / tolerances + p_stop: float = P_STOP # Pa + p_rms_tol: float = 1.0 # Pa (equilibrium RMS tolerance) + +@dataclass(frozen=True) +class NetworkConfig: + # Aggregated branching chain: + N_chain: int + N_par: int + + # Geometry per cell / vestibule (single-branch cell geometry) + V_cell: float + V_vest: float + A_wall_cell: float + A_wall_vest: float + + # Orifice parameters (diameters + counts) + d_int_mm: float + n_int_per_interface: int # holes per interface per cell + d_exit_mm: float + n_exit: int # holes at vestibule exit + + # Cd models + Cd_int: float + Cd_exit: float + + # Optional downstream gap (slot channel) + use_gap: bool = False + V_gap: float = 0.0 + A_wall_gap: float = 0.0 + gap_w: float = 0.0 + gap_delta: float = 0.0 + gap_L: float = 0.0 + + +@dataclass +class SolveResult: + t: np.ndarray + m: np.ndarray # shape (N_nodes, N_t) + T: np.ndarray # shape (N_nodes, N_t) + P: np.ndarray # shape (N_nodes, N_t) + P_ext: np.ndarray # shape (N_t,) + peak_diag: Dict[str, dict] + max_dP: Dict[str, float] + tau_exit: float + meta: dict + + +# ============================================================================= +# Build aggregated network graph +# ============================================================================= +def build_branching_network(cfg: NetworkConfig, profile: Profile) -> Tuple[List[GasNode], List[Edge], List[ExternalBC]]: + assert_pos("N_chain", cfg.N_chain) + assert_pos("N_par", cfg.N_par) + assert_pos("V_cell", cfg.V_cell) + assert_pos("V_vest", cfg.V_vest) + + nodes: List[GasNode] = [] + edges: List[Edge] = [] + + # Node 0: vestibule (single, not multiplied by N_par) + nodes.append(GasNode("vest", cfg.V_vest, cfg.A_wall_vest)) + + # Optional: downstream gap node (between vest and external) + gap_idx = None + if cfg.use_gap: + assert_pos("V_gap", cfg.V_gap) + nodes.append(GasNode("gap", cfg.V_gap, cfg.A_wall_gap)) + gap_idx = 1 # after vest + + # Chain nodes: aggregated across N_par identical chains + # Stage i node index depends on presence of gap. + base = 1 if not cfg.use_gap else 2 + for i in range(cfg.N_chain): + nodes.append(GasNode(f"cell{i+1}", cfg.N_par * cfg.V_cell, cfg.N_par * cfg.A_wall_cell)) + + # Areas + A_int_single = circle_area_from_d_mm(cfg.d_int_mm) + A_exit_single = circle_area_from_d_mm(cfg.d_exit_mm) + + # Total parallel areas + # Between each adjacent stage: N_par chains * n_int holes per interface per cell + A_int_total = cfg.N_par * cfg.n_int_per_interface * A_int_single + # Vestibule exit has n_exit holes total + A_exit_total = cfg.n_exit * A_exit_single + + Cd_int_model = CdConst(cfg.Cd_int) + Cd_exit_model = CdConst(cfg.Cd_exit) + + # Build edges: + # Exit path: either vest->ext directly, or vest->gap->ext + if not cfg.use_gap: + edges.append(OrificeEdge(0, EXT_NODE, A_exit_total, Cd_exit_model, label="exit")) + else: + # vest -> gap: slot channel (gap restriction) + edges.append(SlotChannelEdge(0, gap_idx, cfg.gap_w, cfg.gap_delta, cfg.gap_L, label="vest→gap(slot)")) + # gap -> ext: exit orifice (through panel skin etc.) + edges.append(OrificeEdge(gap_idx, EXT_NODE, A_exit_total, Cd_exit_model, label="gap→ext(exit)")) + + # Internal chain edges: cell1 <-> vest, cell(i+1) <-> cell(i) + # NOTE: with aggregation, areas already include N_par; no multipliers required. + # cell1 index: + cell1 = base + 0 + edges.append(OrificeEdge(cell1, 0, A_int_total, Cd_int_model, label="A(cell1↔vest)")) + for i in range(1, cfg.N_chain): + a = base + i # cell(i+1) + b = base + (i-1) # cell(i) + edges.append(OrificeEdge(a, b, A_int_total, Cd_int_model, label=f"{chr(65+i)}(cell{i+1}↔cell{i})")) + + # External BC is attached to the last exit-connected node: + bc_node = 0 if not cfg.use_gap else gap_idx + bcs = [ExternalBC(node=bc_node, profile=profile, T_ext=T0)] + + return nodes, edges, bcs + + +# ============================================================================= +# RHS assembly for isothermal and intermediate +# ============================================================================= +def build_rhs(nodes: Sequence[GasNode], + edges: Sequence[Edge], + bcs: Sequence[ExternalBC], + case: CaseConfig) -> Tuple[Callable, int]: + """ + Returns rhs(t,y) and n_vars. + + Isothermal state: y = [m_0..m_{N-1}] (T=T0 const) + dm/dt from edges. + P computed for postprocessing only. + + Intermediate state: y = [m_0..m_{N-1}, T_0..T_{N-1}] + dm/dt from edges. + m*c_v*dT/dt = sum inflows/outflows + Q_wall + Q_wall = h*A_wall*(T_wall - T) + """ + N = len(nodes) + V = np.array([n.V for n in nodes], dtype=float) + A_w = np.array([n.A_wall for n in nodes], dtype=float) + + # Map BCs by node + bc_map: Dict[int, ExternalBC] = {bc.node: bc for bc in bcs} + + def P_from_mT(m: np.ndarray, T: np.ndarray) -> np.ndarray: + return m * R_GAS * T / V + + if case.thermo == "isothermal": + n_vars = N + + def rhs(t: float, y: np.ndarray) -> np.ndarray: + m = y[:N] + # enforce nonnegative mass physically: if solver steps slightly negative, stop via events. + T = np.full(N, T0, dtype=float) + P = P_from_mT(np.maximum(m, 0.0), T) # only for flow direction + + dm = np.zeros(N, dtype=float) + + # Edge-wise assembly + for e in edges: + if isinstance(e, OrificeEdge): + a, b = e.a, e.b + # Determine downstream pressure and upstream state + if b == EXT_NODE: + bc = bc_map[a] + Pext = bc.profile.P(float(t)) + # flow direction by P[a] vs Pext + if P[a] >= Pext: + Cd = e.Cd_model() + md = mdot_orifice_pos(P[a], T[a], Pext, Cd, e.A_total) + dm[a] -= md + else: + Cd = e.Cd_model() + md = mdot_orifice_pos(Pext, bc.T_ext, P[a], Cd, e.A_total) + dm[a] += md + else: + Pa, Pb = P[a], P[b] + if Pa >= Pb: + Cd = e.Cd_model() + md = mdot_orifice_pos(Pa, T[a], Pb, Cd, e.A_total) + dm[a] -= md + dm[b] += md + else: + Cd = e.Cd_model() + md = mdot_orifice_pos(Pb, T[b], Pa, Cd, e.A_total) + dm[b] -= md + dm[a] += md + + elif isinstance(e, SlotChannelEdge): + a, b = e.a, e.b + Pa, Pb = P[a], P[b] + if Pa >= Pb: + md = mdot_slot_pos(Pa, T[a], Pb, e.w, e.delta, e.L) + dm[a] -= md + dm[b] += md + else: + md = mdot_slot_pos(Pb, T[b], Pa, e.w, e.delta, e.L) + dm[b] -= md + dm[a] += md + else: + raise TypeError("Unknown edge type") + + return dm + + return rhs, n_vars + + if case.thermo != "intermediate": + raise ValueError("case.thermo must be 'isothermal' or 'intermediate'") + + n_vars = 2*N + h = float(case.h_conv) + + def rhs(t: float, y: np.ndarray) -> np.ndarray: + m = y[:N] + T = y[N:] + T_eff = np.maximum(T, T_SAFE) + + # Pressures for direction decisions + P = P_from_mT(np.maximum(m, 0.0), T_eff) + + dm = np.zeros(N, dtype=float) + dE = np.zeros(N, dtype=float) # contribution to m*c_v*dT/dt + + # Edge-wise assembly + for e in edges: + if isinstance(e, OrificeEdge): + a, b = e.a, e.b + if b == EXT_NODE: + bc = bc_map[a] + Pext = bc.profile.P(float(t)) + if P[a] >= Pext: + Cd = e.Cd_model() + md = mdot_orifice_pos(P[a], T_eff[a], Pext, Cd, e.A_total) + # outflow from node a + dm[a] -= md + dE[a] += -md * R_GAS * T[a] + else: + Cd = e.Cd_model() + md = mdot_orifice_pos(Pext, bc.T_ext, P[a], Cd, e.A_total) + # inflow from external into a + dm[a] += md + dE[a] += md * (C_P * bc.T_ext - C_V * T[a]) + + else: + Pa, Pb = P[a], P[b] + if Pa >= Pb: + Cd = e.Cd_model() + md = mdot_orifice_pos(Pa, T_eff[a], Pb, Cd, e.A_total) + # a upstream -> b downstream + dm[a] -= md + dm[b] += md + dE[a] += -md * R_GAS * T[a] + dE[b] += md * (C_P * T[a] - C_V * T[b]) + else: + Cd = e.Cd_model() + md = mdot_orifice_pos(Pb, T_eff[b], Pa, Cd, e.A_total) + dm[b] -= md + dm[a] += md + dE[b] += -md * R_GAS * T[b] + dE[a] += md * (C_P * T[b] - C_V * T[a]) + + elif isinstance(e, SlotChannelEdge): + a, b = e.a, e.b + Pa, Pb = P[a], P[b] + if Pa >= Pb: + md = mdot_slot_pos(Pa, T_eff[a], Pb, e.w, e.delta, e.L) + dm[a] -= md + dm[b] += md + dE[a] += -md * R_GAS * T[a] + dE[b] += md * (C_P * T[a] - C_V * T[b]) + else: + md = mdot_slot_pos(Pb, T_eff[b], Pa, e.w, e.delta, e.L) + dm[b] -= md + dm[a] += md + dE[b] += -md * R_GAS * T[b] + dE[a] += md * (C_P * T[b] - C_V * T[a]) + else: + raise TypeError("Unknown edge type") + + # Wall heat exchange (constant wall temperature model) + Q = h * A_w * (case.T_wall - T) + dE += Q + + dT = np.where(np.maximum(m, 0.0) > M_SAFE, dE / (np.maximum(m, 0.0) * C_V), 0.0) + + return np.concatenate([dm, dT]) + + return rhs, n_vars + + +# ============================================================================= +# Solver wrapper + events +# ============================================================================= +def solve_case(nodes: Sequence[GasNode], + edges: Sequence[Edge], + bcs: Sequence[ExternalBC], + case: CaseConfig) -> solve_ivp: + N = len(nodes) + V = np.array([n.V for n in nodes], dtype=float) + + rhs, n_vars = build_rhs(nodes, edges, bcs, case) + + # Initial state: uniform P0, T0 + # m0 = P0*V/(R*T0) + m0 = P0 * V / (R_GAS * T0) + + if case.thermo == "isothermal": + y0 = m0.copy() + else: + y0 = np.concatenate([m0, np.full(N, T0)]) + + t_eval = np.linspace(0.0, case.duration, int(case.n_pts)) + + # Jacobian sparsity from graph connectivity + jac = lil_matrix((n_vars, n_vars)) + for i in range(N): + jac[i, i] = 1 + for e in edges: + if isinstance(e, (OrificeEdge, SlotChannelEdge)): + a, b = e.a, e.b + if b != EXT_NODE: + jac[a, b] = 1 + jac[b, a] = 1 + if case.thermo == "intermediate": + # T block + for i in range(N): + jac[N+i, N+i] = 1 + jac[i, N+i] = 1 + jac[N+i, i] = 1 + for e in edges: + if isinstance(e, (OrificeEdge, SlotChannelEdge)): + a, b = e.a, e.b + if b != EXT_NODE: + jac[N+a, N+b] = 1 + jac[N+b, N+a] = 1 + jac[a, N+b] = 1 + jac[b, N+a] = 1 + jac[N+a, b] = 1 + jac[N+b, a] = 1 + + # Events: (1) nonnegative mass, (2) P RMS close to P_ext, (3) all P below p_stop + bc0 = bcs[0] + profile = bc0.profile + + def P_from_state(y: np.ndarray) -> np.ndarray: + if case.thermo == "isothermal": + m = y[:N] + T = np.full(N, T0) + else: + m = y[:N] + T = np.maximum(y[N:], T_SAFE) + return np.maximum(m, 0.0) * R_GAS * T / V + + # (1) Mass nonnegativity: stop if any m_i < -eps (solver should avoid, but guard anyway) + events = [] + for i in range(N): + def ev_mi(t, y, ii=i): + return y[ii] + 1e-15 + ev_mi.terminal = True + ev_mi.direction = -1 + events.append(ev_mi) + + # (2) Pressure RMS equilibrium vs external + def ev_equil_rms(t, y): + P = P_from_state(y) + Pe = profile.P(float(t)) + rms = float(np.sqrt(np.mean((P - Pe)**2))) + return rms - case.p_rms_tol + ev_equil_rms.terminal = True + ev_equil_rms.direction = -1 + events.append(ev_equil_rms) + + # (3) Stop when max(P_i) < p_stop and P_ext already ~vacuum-ish + def ev_p_stop(t, y): + P = P_from_state(y) + Pe = profile.P(float(t)) + return max(np.max(P) - case.p_stop, Pe - 10.0) + ev_p_stop.terminal = True + ev_p_stop.direction = -1 + events.append(ev_p_stop) + + # Step control: estimate a conservative max_step + # Use max total area among orifices as worst-case. + A_max = 0.0 + Cd_max = 0.0 + for e in edges: + if isinstance(e, OrificeEdge): + A_max = max(A_max, e.A_total) + Cd_max = max(Cd_max, e.Cd_model()) + V_min = float(np.min(V)) + mdot_ch = Cd_max * A_max * C_CHOKED * P0 / math.sqrt(R_GAS * T0) if A_max > 0 else 0.0 + tau_min = (V_min * P0) / (R_GAS * T0 * mdot_ch) if mdot_ch > 0 else 1.0 + max_step = max(min(case.duration / 2000.0, tau_min / 10.0), 1e-4) + + rtol = 1e-7 if case.thermo == "isothermal" else 1e-6 + atol = 1e-10 if case.thermo == "isothermal" else 1e-8 + + sol = solve_ivp(rhs, (0.0, case.duration), y0, + method="Radau", + t_eval=t_eval, + rtol=rtol, atol=atol, + max_step=max_step, + jac_sparsity=jac, + events=events) + + return sol + + +# ============================================================================= +# Postprocessing: edges ΔP, peaks, regime, tau_exit +# ============================================================================= +def compute_tau_exit(total_volume: float, Cd_exit: float, A_exit_total: float) -> float: + """ + τ_exit defined from initial choked mass flow at P0,T0 (diagnostic only). + mdot_ch = Cd*A*P0*C_choked/sqrt(R*T0) + tau = (V_total*P0)/(R*T0*mdot_ch) + """ + mdot_ch = Cd_exit * A_exit_total * C_CHOKED * P0 / math.sqrt(R_GAS * T0) + if mdot_ch <= 0.0: + return float("inf") + return (total_volume * P0) / (R_GAS * T0 * mdot_ch) + + +def summarize_result(nodes: Sequence[GasNode], + edges: Sequence[Edge], + bcs: Sequence[ExternalBC], + case: CaseConfig, + sol) -> SolveResult: + N = len(nodes) + V = np.array([n.V for n in nodes], dtype=float) + profile = bcs[0].profile + + t = sol.t + if case.thermo == "isothermal": + m = sol.y[:N] + T = np.full_like(m, T0) + else: + m = sol.y[:N] + T = sol.y[N:2*N] + T_eff = np.maximum(T, T_SAFE) + P = np.maximum(m, 0.0) * R_GAS * T_eff / V[:, None] + P_ext = profile.P_array(t) + + # Build edge ΔP time series and peak diagnostics + dP_edges: Dict[str, np.ndarray] = {} + max_dP: Dict[str, float] = {} + peak_diag: Dict[str, dict] = {} + + # Precompute total volume for tau_exit diagnostic + total_volume = float(np.sum(V)) + + # Determine exit area and Cd_exit for tau_exit: find the first edge labeled "exit" or containing "exit" + A_exit_total = None + Cd_exit = None + for e in edges: + if isinstance(e, OrificeEdge) and ("exit" in e.label.lower() or e.label.lower() == "exit"): + A_exit_total = e.A_total + Cd_exit = e.Cd_model() + break + if A_exit_total is None: + # fallback: max orifice area + A_exit_total = max((e.A_total for e in edges if isinstance(e, OrificeEdge)), default=0.0) + Cd_exit = 0.62 + + tau_exit = compute_tau_exit(total_volume, float(Cd_exit), float(A_exit_total)) + + # Edge loop + for e in edges: + if isinstance(e, OrificeEdge): + a, b = e.a, e.b + label = e.label or f"orifice({a}->{b})" + if b == EXT_NODE: + dp = P[a, :] - P_ext + # upstream at peak based on actual direction at that time + pass + else: + dp = P[a, :] - P[b, :] + dP_edges[label] = dp + + elif isinstance(e, SlotChannelEdge): + a, b = e.a, e.b + label = e.label or f"slot({a}->{b})" + dp = P[a, :] - P[b, :] + dP_edges[label] = dp + + for label, dp in dP_edges.items(): + idx = int(np.argmax(np.abs(dp))) + max_dP[label] = float(np.abs(dp[idx])) + + # Determine upstream/downstream at peak for r and choked diagnostics + # We need to locate the edge by label: + edge_obj = None + for e in edges: + if (isinstance(e, OrificeEdge) and (e.label or f"orifice({e.a}->{e.b})") == label) or \ + (isinstance(e, SlotChannelEdge) and (e.label or f"slot({e.a}->{e.b})") == label): + edge_obj = e + break + + if edge_obj is None: + continue + + if isinstance(edge_obj, OrificeEdge): + a, b = edge_obj.a, edge_obj.b + if b == EXT_NODE: + Pa = float(P[a, idx]) + Pb = float(P_ext[idx]) + Ta = float(T_eff[a, idx]) + else: + Pa = float(P[a, idx]) + Pb = float(P[b, idx]) + Ta = float(T_eff[a, idx]) + + # upstream defined by higher pressure at peak + if Pa >= Pb: + P_up, P_dn, T_up = Pa, Pb, Ta + else: + # if reverse, upstream is "b"; approximate T_up by node b temperature + if b == EXT_NODE: + P_up, P_dn, T_up = Pb, Pa, T0 + else: + P_up, P_dn, T_up = Pb, Pa, float(T_eff[b, idx]) + + r_pk = (P_dn / P_up) if P_up > 0 else 0.0 + choked = bool(r_pk <= PI_C) + + peak_type = profile.classify_peak(float(t[idx]), float(t[-1]), tol_s=0.5) + + peak_diag[label] = { + "t_peak": float(t[idx]), + "t_peak_over_tau_exit": float(t[idx] / tau_exit) if np.isfinite(tau_exit) else float("nan"), + "P_up": P_up, + "P_down": P_dn, + "r": float(r_pk), + "regime": "CHOKED" if choked else "subsonic", + "T_up": float(T_up), + "peak_type": peak_type, + "dP_signed": float(dp[idx]), + } + + else: + # slot channel: no choked criterion + a, b = edge_obj.a, edge_obj.b + Pa = float(P[a, idx]) + Pb = float(P[b, idx]) + peak_type = profile.classify_peak(float(t[idx]), float(t[-1]), tol_s=0.5) + peak_diag[label] = { + "t_peak": float(t[idx]), + "t_peak_over_tau_exit": float(t[idx] / tau_exit) if np.isfinite(tau_exit) else float("nan"), + "P_up": max(Pa, Pb), + "P_down": min(Pa, Pb), + "r": float(min(Pa, Pb)/max(Pa, Pb)) if max(Pa, Pb) > 0 else 0.0, + "regime": "viscous_slot", + "T_up": float(T_eff[a, idx] if Pa >= Pb else T_eff[b, idx]), + "peak_type": peak_type, + "dP_signed": float(dp[idx]), + } + + meta = { + "case": asdict(case), + "nodes": [asdict(n) for n in nodes], + "edges": [asdict(e) for e in edges], + "profile": {"name": profile.name, "events": list(profile.events)}, + "solver": { + "success": bool(sol.success), + "message": str(sol.message), + "t_end": float(sol.t[-1]), + "n_steps": int(len(sol.t)), + } + } + + return SolveResult( + t=t, m=m, T=T, P=P, P_ext=P_ext, + peak_diag=peak_diag, max_dP=max_dP, + tau_exit=float(tau_exit), meta=meta + ) + + +# ============================================================================= +# Gate tests (single-node and two-node) +# ============================================================================= +def gate_test_single_node_orifice() -> None: + """ + Single node V, one orifice to vacuum (P_ext=0). Compare to analytical solutions + in choked regime (r=0). PASS criteria are strict; if it fails, abort work. + """ + # Geometry consistent with earlier examples + V = 131.6e-6 + d_mm = 2.0 + A = circle_area_from_d_mm(d_mm) + Cd = 0.62 + A_wall = 181.6e-4 + + # Analytical alpha based on T0 + alpha = Cd * A * C_CHOKED * math.sqrt(R_GAS * T0) / V + tau = 1.0 / alpha + beta = (GAMMA - 1.0) / 2.0 + + def P_adi(t): return P0 * (1.0 + beta*alpha*t)**(-2.0*GAMMA/(GAMMA-1.0)) + def T_adi(t): return T0 * (1.0 + beta*alpha*t)**(-2.0) + def P_iso(t): return P0 * math.exp(-alpha*t) + + # Build model: one node, one orifice edge to ext, profile is step-to-vacuum at t=0 + prof = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", V, A_wall)] + edges = [OrificeEdge(0, EXT_NODE, A, CdConst(Cd), label="exit")] + bcs = [ExternalBC(0, prof, T_ext=T0)] + + # Run isothermal and intermediate h=0, h->inf + case_iso = CaseConfig(thermo="isothermal", h_conv=0.0, T_wall=T0, duration=6*tau, n_pts=2000) + sol_iso = solve_case(nodes, edges, bcs, case_iso) + res_iso = summarize_result(nodes, edges, bcs, case_iso, sol_iso) + + case_adi = CaseConfig(thermo="intermediate", h_conv=0.0, T_wall=T0, duration=6*tau, n_pts=2000) + sol_adi = solve_case(nodes, edges, bcs, case_adi) + res_adi = summarize_result(nodes, edges, bcs, case_adi, sol_adi) + + case_inf = CaseConfig(thermo="intermediate", h_conv=1e6, T_wall=T0, duration=6*tau, n_pts=2000) + sol_inf = solve_case(nodes, edges, bcs, case_inf) + res_inf = summarize_result(nodes, edges, bcs, case_inf, sol_inf) + + # Compare in region P > 1% P0 + mask = res_adi.P[0, :] > 0.01 * P0 + t = res_adi.t[mask] + + Pn = res_adi.P[0, mask] + Tn = res_adi.T[0, mask] + Pa = np.array([P_adi(float(tt)) for tt in t]) + Ta = np.array([T_adi(float(tt)) for tt in t]) + + errP = float(np.max(np.abs(Pn - Pa) / P0)) + errT = float(np.max(np.abs(Tn - Ta) / T0)) + if not (errP < 5e-3 and errT < 5e-3): + raise RuntimeError(f"GATE1 FAIL (adiabatic): errP={errP:.2e}, errT={errT:.2e}") + + # h->inf should match isothermal P(t) + mask2 = res_inf.P[0, :] > 0.01*P0 + t2 = res_inf.t[mask2] + Pn2 = res_inf.P[0, mask2] + Pi = np.array([P_iso(float(tt)) for tt in t2]) + errP2 = float(np.max(np.abs(Pn2 - Pi) / P0)) + if not (errP2 < 5e-3): + raise RuntimeError(f"GATE1 FAIL (h->inf): errP={errP2:.2e}") + + # Adiabatic must differ from isothermal at 2 tau + t2tau = 2.0 * tau + i2 = int(np.argmin(np.abs(res_adi.t - t2tau))) + P_ad_2 = float(res_adi.P[0, i2]) + P_is_2 = P_iso(float(res_adi.t[i2])) + rel = abs(P_ad_2 - P_is_2)/max(P_ad_2, P_is_2) + if not (rel > 0.2): + raise RuntimeError(f"GATE1 FAIL (adi vs iso not separated): rel={rel:.3f}") + + # Mass conservation check (adiabatic): m_out integral vs m0 - mf + # (numerically, we can approximate m_out from m(t) since ext is vacuum and no inflow) + m0 = P0*V/(R_GAS*T0) + mf = float(res_adi.m[0, -1]) + mout = m0 - mf + if mout < -1e-8: + raise RuntimeError("GATE1 FAIL: mass increased in vacuum blowdown") + + # If reached here: pass + return + + +def gate_test_two_node_mass_energy() -> None: + """ + Two nodes: cell -> vest -> vacuum, h=0. + Check global mass and global internal-energy balance against integrated outflow enthalpy. + """ + Vc = 131.6e-6 + Vv = 145.3e-6 + d_mm = 2.0 + A = circle_area_from_d_mm(d_mm) + Cd = 0.62 + Aw = 181.6e-4 + + prof = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("vest", Vv, Aw), GasNode("cell", Vc, Aw)] + edges = [ + OrificeEdge(1, 0, A, CdConst(Cd), label="cell↔vest"), + OrificeEdge(0, EXT_NODE, A, CdConst(Cd), label="exit"), + ] + bcs = [ExternalBC(0, prof, T_ext=T0)] + + case = CaseConfig(thermo="intermediate", h_conv=0.0, T_wall=T0, duration=3.0, n_pts=3000) + sol = solve_case(nodes, edges, bcs, case) + res = summarize_result(nodes, edges, bcs, case, sol) + + # Mass: m_total(t) + m_out(t) = const. + m_total_0 = float(np.sum(res.m[:, 0])) + m_total_f = float(np.sum(res.m[:, -1])) + + # Compute mdot_exit(t) using upstream node 0 vs vacuum + # This is only for auditing, not for the solver. + mdot = [] + for k, tt in enumerate(res.t): + Pvest = float(res.P[0, k]) + Tvest = float(max(res.T[0, k], T_SAFE)) + md = mdot_orifice_pos(Pvest, Tvest, 0.0, Cd, A) + mdot.append(md) + mdot = np.array(mdot, dtype=float) + m_out = float(np.trapezoid(mdot, res.t)) + + err_m = abs((m_total_f + m_out) - m_total_0) / max(m_total_0, 1e-20) + if not (err_m < 1e-3): + raise RuntimeError(f"GATE2 FAIL (mass): err={err_m:.2e}") + + # Energy: E_int(t) + ∫ mdot_exit * cp * T_up dt = E0 (since rigid, Q=0) + E0 = float(np.sum(res.m[:, 0] * C_V * res.T[:, 0])) + Ef = float(np.sum(res.m[:, -1] * C_V * res.T[:, -1])) + Eout = float(np.trapezoid(mdot * C_P * np.maximum(res.T[0, :], T_SAFE), res.t)) + err_E = abs((Ef + Eout) - E0) / max(E0, 1e-20) + if not (err_E < 1e-2): + raise RuntimeError(f"GATE2 FAIL (energy): err={err_E:.2e}") + + return + + +# ============================================================================= +# IO: save/load results +# ============================================================================= +def save_result(outdir: Path, name: str, res: SolveResult) -> None: + outdir.mkdir(parents=True, exist_ok=True) + npz_path = outdir / f"{name}.npz" + json_path = outdir / f"{name}_meta.json" + + np.savez_compressed( + npz_path, + t=res.t, m=res.m, T=res.T, P=res.P, P_ext=res.P_ext, + tau_exit=res.tau_exit, + ) + with open(json_path, "w", encoding="utf-8") as f: + json.dump(res.meta, f, ensure_ascii=False, indent=2) + + +def plot_basic(outdir: Path, name: str, res: SolveResult, node_idx: int = 0) -> None: + if not HAS_MPL: + return + outdir.mkdir(parents=True, exist_ok=True) + + fig, ax = plt.subplots(figsize=(11, 6)) + ax.plot(res.t, (res.P[node_idx] - res.P_ext)/1e3, lw=2, label=f"ΔP(node{node_idx}→ext) [kPa]") + ax.set(xlabel="t, s", ylabel="ΔP, kPa", title=f"{name}: ΔP vs time") + ax.grid(True, alpha=0.3) + ax.legend() + fig.tight_layout() + fig.savefig(outdir / f"{name}_dP.png", dpi=200) + plt.close(fig) + + if res.T is not None: + fig, ax = plt.subplots(figsize=(11, 6)) + ax.plot(res.t, res.T[node_idx], lw=2, label=f"T(node{node_idx}) [K]") + ax.set(xlabel="t, s", ylabel="T, K", title=f"{name}: T vs time") + ax.grid(True, alpha=0.3) + ax.legend() + fig.tight_layout() + fig.savefig(outdir / f"{name}_T.png", dpi=200) + plt.close(fig) + + +# ============================================================================= +# Main: example run (sweeps similar to v8.3, but strict profiles) +# ============================================================================= +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--outdir", type=str, default="./results_v84", help="Output directory") + ap.add_argument("--profile", type=str, default="linear", choices=["linear", "step", "barometric", "table"], + help="External depressurization profile") + ap.add_argument("--profile-file", type=str, default="", help="CSV file for --profile table (t_s,P_Pa)") + ap.add_argument("--rate-mmhg", type=float, default=20.0, help="Rate for linear / barometric initial slope (mmHg/s)") + ap.add_argument("--step-time", type=float, default=0.01, help="Step time for step profile (s)") + ap.add_argument("--thermo", type=str, default="isothermal", choices=["isothermal", "intermediate"]) + ap.add_argument("--h", type=float, default=0.0, help="h_conv W/(m2*K), used in intermediate") + ap.add_argument("--duration", type=float, default=150.0) + ap.add_argument("--npts", type=int, default=2000) + ap.add_argument("--do-plots", action="store_true") + + # Network inputs (defaults from your CAD numbers in v8.2/v8.3) + ap.add_argument("--d-int", type=float, default=2.0) + ap.add_argument("--d-exit", type=float, default=2.0) + ap.add_argument("--n-int", type=int, default=1, help="holes per interface per cell") + ap.add_argument("--n-exit", type=int, default=1, help="exit holes total") + + ap.add_argument("--Cd-int", type=float, default=0.62) + ap.add_argument("--Cd-exit", type=float, default=0.62) + + ap.add_argument("--use-gap", action="store_true") + ap.add_argument("--gap-V", type=float, default=0.0, help="gap volume [m3]") + ap.add_argument("--gap-Aw", type=float, default=0.0, help="gap wall area [m2]") + ap.add_argument("--gap-w", type=float, default=0.0, help="slot width [m]") + ap.add_argument("--gap-delta", type=float, default=0.0, help="slot thickness [m]") + ap.add_argument("--gap-L", type=float, default=0.0, help="slot length [m]") + + args = ap.parse_args() + outdir = Path(args.outdir) + + # Gate tests first (hard fail on inconsistency) + gate_test_single_node_orifice() + gate_test_two_node_mass_energy() + + # Geometry (from your v8.2/v8.3 CAD-derived numbers) + A_cell_mm2 = 4708.2806538 + h_mm = 27.9430913 + V_cell = mm3_to_m3(A_cell_mm2 * h_mm) + A_vest_mm2 = 5199.9595503 + V_vest = mm3_to_m3(A_vest_mm2 * h_mm) + + P_cell_mm = 313.0 + A_wall_cell = mm2_to_m2(2*A_cell_mm2 + P_cell_mm*h_mm) + P_vest_mm = 350.0 + A_wall_vest = mm2_to_m2(2*A_vest_mm2 + P_vest_mm*h_mm) + + # Profile selection + if args.profile == "linear": + prof = make_profile_linear(P0, args.rate_mmhg) + elif args.profile == "step": + prof = make_profile_step(P0, args.step_time) + elif args.profile == "barometric": + prof = make_profile_exponential(P0, args.rate_mmhg, p_floor=10.0) + else: + if not args.profile_file: + raise ValueError("--profile table requires --profile-file") + prof = make_profile_from_table("envelope_table", Path(args.profile_file)) + + net_cfg = NetworkConfig( + N_chain=10, N_par=2, + V_cell=V_cell, V_vest=V_vest, + A_wall_cell=A_wall_cell, A_wall_vest=A_wall_vest, + d_int_mm=args.d_int, n_int_per_interface=args.n_int, + d_exit_mm=args.d_exit, n_exit=args.n_exit, + Cd_int=args.Cd_int, Cd_exit=args.Cd_exit, + use_gap=bool(args.use_gap), + V_gap=float(args.gap_V), + A_wall_gap=float(args.gap_Aw), + gap_w=float(args.gap_w), + gap_delta=float(args.gap_delta), + gap_L=float(args.gap_L), + ) + + nodes, edges, bcs = build_branching_network(net_cfg, prof) + + case = CaseConfig( + thermo=args.thermo, + h_conv=float(args.h), + T_wall=T0, + duration=float(args.duration), + n_pts=int(args.npts), + p_stop=P_STOP, + p_rms_tol=1.0, + ) + + sol = solve_case(nodes, edges, bcs, case) + if not sol.success: + raise RuntimeError(f"Solver failed: {sol.message}") + + res = summarize_result(nodes, edges, bcs, case, sol) + + # Save + name = f"v84_{prof.name}_{case.thermo}_dint{args.d_int:g}_dexit{args.d_exit:g}" + save_result(outdir, name, res) + + # Quick console diagnostics + # Find exit-like edge for reporting + exit_keys = [k for k in res.max_dP.keys() if "exit" in k.lower()] + key = exit_keys[0] if exit_keys else list(res.max_dP.keys())[0] + pd = res.peak_diag.get(key, {}) + print("\n=== v8.4 SUMMARY ===") + print(f"case: profile={prof.name} thermo={case.thermo} h={case.h_conv:g}") + print(f"exit-like edge: {key}") + print(f" max|ΔP| = {res.max_dP[key]:.1f} Pa = {res.max_dP[key]/1e3:.3f} kPa") + if pd: + print(f" t_peak = {pd['t_peak']:.3f} s, t/τ_exit={pd['t_peak_over_tau_exit']:.2f}") + print(f" r = {pd['r']:.4f}, regime={pd['regime']}, peak_type={pd['peak_type']}") + print(f" P_up={pd['P_up']/1e3:.3f} kPa, P_down={pd['P_down']/1e3:.3f} kPa, T_up={pd['T_up']:.2f} K") + print(f"saved: {outdir / (name + '.npz')}") + print(f"meta : {outdir / (name + '_meta.json')}") + + if args.do_plots: + plot_basic(outdir, name, res, node_idx=0) + + +if __name__ == "__main__": + main() diff --git a/docs/assumptions_v85.md b/docs/assumptions_v85.md new file mode 100644 index 0000000..b141507 --- /dev/null +++ b/docs/assumptions_v85.md @@ -0,0 +1,9 @@ +# Assumptions and limitations (v8.5) + +- 0D lumped compartments. +- Ideal-gas air with constant `gamma`, `R`, `c_p`, `c_v`. +- Orifice model uses ideal isentropic relation with constant/effective `Cd`. +- Slot model is laminar viscous approximation. +- Model-form uncertainty remains in geometry and discharge coefficients. + +- Optional short-tube edge uses a lossy-nozzle `Cd_eff` correction with friction/minor losses; full Fanno-flow choking is not implemented. diff --git a/docs/assumptions_v9.md b/docs/assumptions_v9.md new file mode 100644 index 0000000..a2a4827 --- /dev/null +++ b/docs/assumptions_v9.md @@ -0,0 +1,10 @@ +# Assumptions and limitations (v9.0.0) + +- 0D lumped compartments with perfect mixing in each node. +- Ideal-gas EOS in all modes (`P = mRT/V`). +- Short-tube thick-wall holes are modeled as **lossy nozzle via `Cd_eff`** with Darcy friction + minor losses. +- This is **not Fanno flow**; friction choking is not modeled in v9. +- Variable thermo mode uses smooth engineering fits `cp(T), cv(T), gamma(T), h(T), u(T)` for air. +- Lumped wall mode treats wall temperature as a single thermal capacitance state per node. +- Dynamic external model (`dynamic_pump`) adds a finite external volume and a pump sink law. +- Model-form uncertainty remains dominated by geometry and discharge/loss coefficients. diff --git a/docs/model_v85.md b/docs/model_v85.md new file mode 100644 index 0000000..6901091 --- /dev/null +++ b/docs/model_v85.md @@ -0,0 +1,17 @@ +# Model v8.5 + +State: +- Isothermal: `m_i` +- Intermediate: `(m_i, T_i)` + +Pressure relation: +`P_i = m_i R T_i / V_i` + +Mass flow: +- Orifice: compressible isentropic, with choked/subsonic branching by `r = P_dn / P_up`. +- Slot: `m_dot = K (P_up^2 - P_dn^2)/(R T_up)`, `K = w delta^3 / (12 mu L)`. + +Energy equation (intermediate): +`m c_v dT/dt = Σ_in m_dot_in (c_p T_in - c_v T) - Σ_out m_dot_out R T + h A (T_wall - T)` + +No state clipping is used. diff --git a/docs/verification_v85.md b/docs/verification_v85.md new file mode 100644 index 0000000..cd63e1e --- /dev/null +++ b/docs/verification_v85.md @@ -0,0 +1,44 @@ +# Verification v8.5 + +## Governing equations used for verification + +For each control volume (rigid `V`, ideal gas): + +- Mass: + - `dm/dt = Σ m_dot_in - Σ m_dot_out` +- Pressure closure: + - `P = m R T / V` +- Intermediate energy form in solver: + - `m c_v dT/dt = Σ_in m_dot_in (c_p T_in - c_v T) - Σ_out m_dot_out R T + h A (T_wall - T)` + +Why outflow is `m_dot_out * R T`: for a rigid control volume, +`dU/dt = Σ m_dot_in h_in - Σ m_dot_out h_out + Q`, with `U = m c_v T` and +`h = c_p T`; regrouping terms for unknown node `T` yields the explicit `-m_dot_out R T` +term in the `c_v` form. + +## Analytic references in tests + +Single-node blowdown to vacuum (`P_ext=0`, rigid volume, ideal gas, choked branch): + +- Adiabatic (`h=0`): + - `P(t) = P0 (1 + beta alpha t)^(-2 gamma/(gamma-1))` + - `T(t) = T0 (1 + beta alpha t)^(-2)` + - `beta = (gamma-1)/2` +- Isothermal limit (`h -> inf`): + - `P(t) = P0 exp(-alpha t)` + +## Gate criteria + +- Single-node analytic match: + - max rel error for pressure `< 0.5%` + - max rel error for temperature `< 0.5%` +- Isothermal-limit match: + - max rel pressure error `< 0.5%` +- Conservation: + - mass `< 0.1%` + - two-node energy `< 1%` +- Monotonic vacuum blowdown pressure: + - nonincreasing pressure (tiny numerical tolerance only) + +Masks exclude the low-pressure tail (`P <= 0.01*P0`) because relative-error metrics +become numerically ill-conditioned near machine floor where absolute pressures are tiny. diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..9d568c9 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,54 @@ +[build-system] +requires = ["setuptools>=68", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "venting" +version = "10.0.0" +description = "v10.0.0 edge-based 0D venting solver" +readme = "README.md" +requires-python = ">=3.10" +license = { text = "MIT" } +authors = [{ name = "test1card" }] +dependencies = [ + "numpy>=1.24,<3", + "scipy>=1.10,<1.13", + "matplotlib>=3.7,<4", +] + +[project.optional-dependencies] +dev = [ + "pytest>=7.4,<9", + "pytest-cov>=4.1,<6", + "ruff>=0.6,<1", + "black>=24,<26", + "mypy>=1.8,<2", + "pre-commit>=3.6,<4", +] +gui = [ + "PySide6>=6.6", + "pyqtgraph>=0.13", +] + +[tool.setuptools] +package-dir = {"" = "src"} + +[tool.setuptools.packages.find] +where = ["src"] +include = ["venting*"] + +[tool.black] +line-length = 88 +extend-exclude = "archive/" + +[tool.ruff] +line-length = 88 +extend-exclude = ["archive"] + +[tool.ruff.lint] +select = ["F", "E", "I", "B", "UP"] +ignore = ["E501"] + +[tool.pytest.ini_options] +testpaths = ["tests"] +addopts = "-q" diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 0000000..aefbcb6 --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1 @@ +-e .[dev] diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..bd3e42a --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +numpy>=1.24,<3 +scipy>=1.10,<1.13 +matplotlib>=3.7,<4 diff --git a/src/venting/__init__.py b/src/venting/__init__.py new file mode 100644 index 0000000..d9416a0 --- /dev/null +++ b/src/venting/__init__.py @@ -0,0 +1,5 @@ +"""venting v10.0.0 package.""" + +__all__ = ["__version__"] + +__version__ = "10.0.0" diff --git a/src/venting/__main__.py b/src/venting/__main__.py new file mode 100644 index 0000000..9ae637f --- /dev/null +++ b/src/venting/__main__.py @@ -0,0 +1,4 @@ +from .cli import main + +if __name__ == "__main__": + main() diff --git a/src/venting/cases.py b/src/venting/cases.py new file mode 100644 index 0000000..04d3456 --- /dev/null +++ b/src/venting/cases.py @@ -0,0 +1,82 @@ +from dataclasses import dataclass + +import numpy as np + +from .constants import P_STOP, T0 + + +@dataclass(frozen=True) +class CaseConfig: + thermo: str + h_conv: float + T_wall: float + duration: float + n_pts: int + p_stop: float = P_STOP + p_rms_tol: float = 1.0 + l_char_m: float = 0.1 + wall_model: str = "fixed" # fixed|lumped + wall_C_per_area: float = 1e9 + wall_h_out: float = 0.0 + wall_T_inf: float = T0 + wall_emissivity: float = 0.0 + wall_T_sur: float = T0 + wall_q_flux: float = 0.0 + external_model: str = "profile" # profile|dynamic_pump + V_ext: float = 0.1 + T_ext: float = T0 + pump_speed_m3s: float = 0.0 + P_ult_Pa: float = 0.0 + + +NumberOrList = float | int | list[float] | tuple[float, ...] + + +@dataclass(frozen=True) +class NetworkConfig: + N_chain: int + N_par: int + V_cell: NumberOrList + V_vest: float + A_wall_cell: NumberOrList + A_wall_vest: float + d_int_mm: NumberOrList + n_int_per_interface: NumberOrList + d_exit_mm: float + n_exit: int + Cd_int: NumberOrList + Cd_exit: float + int_model: str = "orifice" + exit_model: str = "orifice" + L_int_mm: NumberOrList = 0.0 + L_exit_mm: float = 0.0 + eps_um: float = 0.0 + K_in: float = 0.5 + K_out: float = 1.0 + K_in_int: NumberOrList | None = None + K_out_int: NumberOrList | None = None + eps_int_um: NumberOrList | None = None + K_in_exit: float | None = None + K_out_exit: float | None = None + eps_exit_um: float | None = None + topology: str = "single_chain" # single_chain|two_chain_shared_vest + N_chain_b: int | None = None + use_gap: bool = False + V_gap: float = 0.0 + A_wall_gap: float = 0.0 + gap_w: float = 0.0 + gap_delta: float = 0.0 + gap_L: float = 0.0 + + +@dataclass +class SolveResult: + t: np.ndarray + m: np.ndarray + T: np.ndarray + P: np.ndarray + P_ext: np.ndarray + peak_diag: dict[str, dict] + max_dP: dict[str, float] + tau_exit: float + meta: dict diff --git a/src/venting/cli.py b/src/venting/cli.py new file mode 100644 index 0000000..dd192bd --- /dev/null +++ b/src/venting/cli.py @@ -0,0 +1,365 @@ +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +from .cases import CaseConfig, NetworkConfig +from .compare import ( + compare_runs, + format_comparison_table, + load_run, + write_comparison_csv, +) +from .constants import P0, P_STOP, T0 +from .gates import gate_single, gate_two +from .montecarlo import run_mc, write_mc_outputs +from .plotting import plot_basic +from .presets import get_default_panel_preset_v9 +from .profiles import ( + make_profile_exponential, + make_profile_from_table, + make_profile_linear, + make_profile_step, +) +from .run import export_case_artifacts, make_case_output_dir, run_case + + +def _profile(args): + if args.external_model == "dynamic_pump": + return make_profile_step(P0, 1e9) + if args.profile == "linear": + return make_profile_linear(P0, args.rate_mmhg) + if args.profile == "step": + return make_profile_step(P0, args.step_time) + if args.profile == "barometric": + return make_profile_exponential(P0, args.rate_mmhg, p_floor=10.0) + return make_profile_from_table( + "envelope_table", + Path(args.profile_file), + pressure_unit=args.profile_pressure_unit, + ) + + +def _run_one(args, d_int: float, d_exit: float, h: float | None = None): + prof = _profile(args) + preset = get_default_panel_preset_v9() + + net_cfg = NetworkConfig( + N_chain=10, + N_par=2, + topology=args.topology, + N_chain_b=args.n_chain_b, + 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=d_int, + n_int_per_interface=args.n_int, + d_exit_mm=d_exit, + n_exit=args.n_exit, + Cd_int=args.cd_int, + Cd_exit=args.cd_exit, + int_model=args.int_model, + exit_model=args.exit_model, + L_int_mm=args.L_int_mm, + L_exit_mm=args.L_exit_mm, + K_in=args.K_in, + K_out=args.K_out, + eps_um=args.eps_um, + K_in_int=args.K_in_int, + K_out_int=args.K_out_int, + eps_int_um=args.eps_int_um, + K_in_exit=args.K_in_exit, + K_out_exit=args.K_out_exit, + eps_exit_um=args.eps_exit_um, + ) + case = CaseConfig( + thermo=args.thermo, + h_conv=args.h if h is None else h, + T_wall=T0, + duration=args.duration, + n_pts=args.npts, + p_stop=P_STOP, + p_rms_tol=1.0, + wall_model=args.wall_model, + wall_C_per_area=args.wall_C_per_area, + wall_h_out=args.wall_h_out, + wall_T_inf=args.wall_T_inf, + wall_emissivity=args.wall_emissivity, + wall_T_sur=args.wall_T_sur, + wall_q_flux=args.wall_q_flux, + external_model=args.external_model, + V_ext=args.V_ext, + T_ext=args.T_ext, + pump_speed_m3s=args.pump_speed_m3s, + P_ult_Pa=args.P_ult_Pa, + ) + + res = run_case(net_cfg, prof, case) + + out = make_case_output_dir(args.cmd) + name = ( + f"v1000_{args.external_model}_{prof.name}_{case.thermo}_" + f"dint{d_int:g}_dexit{d_exit:g}_h{case.h_conv:g}" + ) + run_params = { + "command": args.cmd, + "profile": args.profile, + "external_model": args.external_model, + "d_int": d_int, + "d_exit": d_exit, + "thermo": case.thermo, + "wall_model": case.wall_model, + "h": case.h_conv, + "duration": case.duration, + "npts": case.n_pts, + "cd_int": args.cd_int, + "cd_exit": args.cd_exit, + "int_model": args.int_model, + "exit_model": args.exit_model, + "L_int_mm": args.L_int_mm, + "L_exit_mm": args.L_exit_mm, + "K_in_int": net_cfg.K_in_int, + "K_out_int": net_cfg.K_out_int, + "eps_int_um": net_cfg.eps_int_um, + "K_in_exit": net_cfg.K_in_exit, + "K_out_exit": net_cfg.K_out_exit, + "eps_exit_um": net_cfg.eps_exit_um, + "K_in_alias": args.K_in, + "K_out_alias": args.K_out, + "eps_um_alias": args.eps_um, + "profile_pressure_unit": args.profile_pressure_unit, + "topology": args.topology, + "n_chain_b": args.n_chain_b, + "pump_speed_m3s": case.pump_speed_m3s, + "V_ext": case.V_ext, + "P_ult_Pa": case.P_ult_Pa, + "panel_preset": preset.to_dict(), + } + export_case_artifacts(out, name, res, run_params) + if args.do_plots: + plot_basic(out, name, res) + + +def _add_common_args(s: argparse.ArgumentParser) -> None: + s.add_argument( + "--profile", default="linear", choices=["linear", "step", "barometric", "table"] + ) + s.add_argument( + "--external-model", default="profile", choices=["profile", "dynamic_pump"] + ) + s.add_argument("--profile-file", default="") + s.add_argument("--profile-pressure-unit", default="Pa", choices=["Pa", "mmHg"]) + s.add_argument("--rate-mmhg", type=float, default=20.0) + s.add_argument("--step-time", type=float, default=0.01) + s.add_argument( + "--thermo", + default="isothermal", + choices=["isothermal", "intermediate", "variable"], + ) + s.add_argument("--h", type=float, default=0.0) + s.add_argument("--duration", type=float, default=150.0) + s.add_argument("--npts", type=int, default=800) + s.add_argument("--n-int", type=int, default=1) + s.add_argument( + "--topology", + choices=["single_chain", "two_chain_shared_vest"], + default="single_chain", + ) + s.add_argument("--n-chain-b", type=int, default=10) + s.add_argument("--n-exit", type=int, default=1) + s.add_argument("--cd-int", type=float, default=0.62) + s.add_argument("--cd-exit", type=float, default=0.62) + s.add_argument( + "--int-model", choices=["orifice", "short_tube", "fanno"], default="orifice" + ) + s.add_argument( + "--exit-model", choices=["orifice", "short_tube", "fanno"], default="orifice" + ) + s.add_argument("--L-int-mm", type=float, default=0.0) + s.add_argument("--L-exit-mm", type=float, default=0.0) + s.add_argument("--K-in-int", type=float, default=None) + s.add_argument("--K-out-int", type=float, default=None) + s.add_argument("--eps-int-um", type=float, default=None) + s.add_argument("--K-in-exit", type=float, default=None) + s.add_argument("--K-out-exit", type=float, default=None) + s.add_argument("--eps-exit-um", type=float, default=None) + s.add_argument( + "--K-in", + type=float, + default=0.5, + help="Deprecated alias for both internal/exit K_in", + ) + s.add_argument( + "--K-out", + type=float, + default=1.0, + help="Deprecated alias for both internal/exit K_out", + ) + s.add_argument( + "--eps-um", + type=float, + default=0.0, + help="Deprecated alias for both internal/exit roughness", + ) + s.add_argument("--wall-model", choices=["fixed", "lumped"], default="fixed") + s.add_argument("--wall-C-per-area", type=float, default=1e9) + s.add_argument("--wall-h-out", type=float, default=0.0) + s.add_argument("--wall-T-inf", type=float, default=T0) + s.add_argument("--wall-emissivity", type=float, default=0.0) + s.add_argument("--wall-T-sur", type=float, default=T0) + s.add_argument("--wall-q-flux", type=float, default=0.0) + s.add_argument("--V-ext", type=float, default=0.1) + s.add_argument("--T-ext", type=float, default=T0) + s.add_argument("--pump-speed-m3s", type=float, default=0.0) + s.add_argument("--P-ult-Pa", type=float, default=0.0) + s.add_argument("--do-plots", action="store_true") + + +def build_parser() -> argparse.ArgumentParser: + p = argparse.ArgumentParser(prog="venting") + sub = p.add_subparsers(dest="cmd", required=True) + + g = sub.add_parser("gate") + g.add_argument("--single", action="store_true") + g.add_argument("--two", action="store_true") + + sub.add_parser("gui") + + for name in ["sweep", "thermal", "sweep2d"]: + s = sub.add_parser(name) + _add_common_args(s) + + sub.choices["sweep"].add_argument("--d-int", type=float, default=2.0) + sub.choices["sweep"].add_argument("--d-exit", type=float, default=2.0) + + sub.choices["thermal"].add_argument("--d", type=float, default=2.0) + sub.choices["thermal"].add_argument("--h-list", default="0,1,5,15") + + sub.choices["sweep2d"].add_argument("--d-int-list", default="1.0,2.0") + sub.choices["sweep2d"].add_argument("--d-exit-list", default="1.0,2.0") + + c = sub.add_parser("compare") + c.add_argument("run_a") + c.add_argument("run_b") + c.add_argument("--output", default="") + + mc = sub.add_parser("mc") + _add_common_args(mc) + mc.add_argument("--d-int", type=float, default=2.0) + mc.add_argument("--d-exit", type=float, default=2.0) + mc.add_argument("--cd-int-range", default="0.5,0.7") + mc.add_argument("--cd-exit-range", default="0.55,0.65") + mc.add_argument("--n-samples", type=int, default=100) + mc.add_argument("--seed", type=int, default=None) + + return p + + +def main() -> None: + args = build_parser().parse_args() + if args.cmd == "gate": + run_single = args.single or (not args.single and not args.two) + run_two = args.two or (not args.single and not args.two) + if run_single: + print(gate_single()) + if run_two: + print(gate_two()) + return + if args.cmd == "gui": + from .gui.main import main as gui_main + + gui_main() + return + if args.cmd == "compare": + run_a = load_run(Path(args.run_a)) + run_b = load_run(Path(args.run_b)) + comp = compare_runs(run_a, run_b) + print(format_comparison_table(comp)) + if args.output: + write_comparison_csv(comp, Path(args.output)) + return + if args.cmd == "mc": + lo_i, hi_i = [float(x) for x in args.cd_int_range.split(",", 1)] + lo_e, hi_e = [float(x) for x in args.cd_exit_range.split(",", 1)] + prof = _profile(args) + preset = get_default_panel_preset_v9() + net_cfg = NetworkConfig( + N_chain=10, + N_par=args.n_int, + 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, + d_exit_mm=args.d_exit, + n_exit=args.n_exit, + Cd_int=args.cd_int, + Cd_exit=args.cd_exit, + int_model=args.int_model, + exit_model=args.exit_model, + L_int_mm=args.L_int_mm, + L_exit_mm=args.L_exit_mm, + K_in=args.K_in, + K_out=args.K_out, + eps_um=args.eps_um, + K_in_int=args.K_in_int, + K_out_int=args.K_out_int, + eps_int_um=args.eps_int_um, + K_in_exit=args.K_in_exit, + K_out_exit=args.K_out_exit, + eps_exit_um=args.eps_exit_um, + topology=args.topology, + N_chain_b=args.n_chain_b, + ) + case = CaseConfig( + thermo=args.thermo, + h_conv=args.h, + T_wall=T0, + duration=args.duration, + n_pts=args.npts, + wall_model=args.wall_model, + wall_C_per_area=args.wall_C_per_area, + wall_h_out=args.wall_h_out, + wall_T_inf=args.wall_T_inf, + wall_emissivity=args.wall_emissivity, + wall_T_sur=args.wall_T_sur, + wall_q_flux=args.wall_q_flux, + external_model=args.external_model, + V_ext=args.V_ext, + T_ext=args.T_ext, + pump_speed_m3s=args.pump_speed_m3s, + P_ult_Pa=args.P_ult_Pa, + ) + out = run_mc( + net_cfg, + case, + prof, + (lo_i, hi_i), + (lo_e, hi_e), + args.n_samples, + seed=args.seed, + ) + write_mc_outputs(out, make_case_output_dir("mc")) + print(json.dumps(out["summary"], ensure_ascii=False, indent=2)) + return + if args.cmd == "sweep": + _run_one(args, args.d_int, args.d_exit) + return + if args.cmd == "thermal": + hs = [float(x) for x in args.h_list.split(",") if x] + for h in hs: + _run_one(args, args.d, args.d, h=h) + return + + d_ints = [float(x) for x in args.d_int_list.split(",") if x] + d_exits = [float(x) for x in args.d_exit_list.split(",") if x] + for d_int in d_ints: + for d_exit in d_exits: + _run_one(args, d_int, d_exit) + + +if __name__ == "__main__": + main() diff --git a/src/venting/compare.py b/src/venting/compare.py new file mode 100644 index 0000000..69e4c0e --- /dev/null +++ b/src/venting/compare.py @@ -0,0 +1,107 @@ +from __future__ import annotations + +import csv +import json +from pathlib import Path + +import numpy as np + + +def load_run(run_dir: Path) -> dict: + """Load npz + metadata artifacts from a results directory.""" + if not run_dir.exists() or not run_dir.is_dir(): + raise FileNotFoundError(f"Run directory not found: {run_dir}") + npz_files = sorted(run_dir.glob("*.npz")) + if not npz_files: + raise FileNotFoundError(f"No .npz file found in {run_dir}") + npz = np.load(npz_files[0]) + + meta_files = sorted(run_dir.glob("*_meta.json")) + meta = json.loads(meta_files[0].read_text(encoding="utf-8")) if meta_files else {} + + summary_path = run_dir / "summary.csv" + summary = [] + if summary_path.exists(): + with summary_path.open("r", encoding="utf-8") as fh: + summary = list(csv.DictReader(fh)) + + return { + "dir": run_dir, + "name": npz_files[0].stem, + "t": npz["t"], + "P": npz["P"], + "T": npz["T"], + "tau_exit": float(npz["tau_exit"]), + "meta": meta, + "summary": summary, + } + + +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] + + delta_summary = {} + by_edge_a = {row["edge"]: row for row in run_a.get("summary", [])} + by_edge_b = {row["edge"]: row for row in run_b.get("summary", [])} + for edge in sorted(set(by_edge_a) | set(by_edge_b)): + ra = by_edge_a.get(edge) + rb = by_edge_b.get(edge) + if not ra or not rb: + continue + delta_summary[edge] = { + "delta_max_abs_dP_Pa": float(rb["max_abs_dP_Pa"]) + - float(ra["max_abs_dP_Pa"]), + "delta_t_peak_s": float(rb["t_peak_s"]) - float(ra["t_peak_s"]), + "regime_a": ra.get("regime", ""), + "regime_b": rb.get("regime", ""), + } + + tau_a = float(run_a.get("tau_exit", np.nan)) + tau_b = float(run_b.get("tau_exit", np.nan)) + rel_tau = (tau_b - tau_a) / tau_a if abs(tau_a) > 0 else float("nan") + + return { + "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))), + "rel_tau_exit_delta": float(rel_tau), + } + + +def format_comparison_table(comparison: dict) -> str: + lines = [ + f"A: {comparison['run_a']}", + f"B: {comparison['run_b']}", + "edge | Δmax|ΔP| [Pa] | Δt_peak [s] | regime A -> B", + ] + for edge, d in comparison["edge_deltas"].items(): + lines.append( + f"{edge} | {d['delta_max_abs_dP_Pa']:.6g} | {d['delta_t_peak_s']:.6g} | {d['regime_a']} -> {d['regime_b']}" + ) + lines.append(f"max final P diff: {comparison['max_final_P_diff']:.6g}") + lines.append(f"max final T diff: {comparison['max_final_T_diff']:.6g}") + lines.append(f"relative tau_exit delta: {comparison['rel_tau_exit_delta']:.6g}") + return "\n".join(lines) + + +def write_comparison_csv(comparison: dict, path: Path) -> None: + with path.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter( + fh, + fieldnames=[ + "edge", + "delta_max_abs_dP_Pa", + "delta_t_peak_s", + "regime_a", + "regime_b", + ], + ) + writer.writeheader() + for edge, d in comparison["edge_deltas"].items(): + writer.writerow({"edge": edge, **d}) diff --git a/src/venting/constants.py b/src/venting/constants.py new file mode 100644 index 0000000..dea1cba --- /dev/null +++ b/src/venting/constants.py @@ -0,0 +1,22 @@ +import math + +GAMMA = 1.4 +R_GAS = 287.05 +C_V = R_GAS / (GAMMA - 1.0) +C_P = GAMMA * C_V + +# Boltzmann constant [J/K], exact SI (2019 redefinition) +K_BOLTZMANN = 1.380649e-23 +# Effective collision diameter for dry air molecules [m] +# Source: Jennings 1988, J. Aerosol Sci. 19(2):159-166 +D_MOL_AIR = 3.65e-10 + +T0 = 300.0 +P0 = 101325.0 + +T_SAFE = 1.0 +M_SAFE = 1e-18 +P_STOP = 5.0 + +PI_C = (2.0 / (GAMMA + 1.0)) ** (GAMMA / (GAMMA - 1.0)) +C_CHOKED = math.sqrt(GAMMA * (2.0 / (GAMMA + 1.0)) ** ((GAMMA + 1.0) / (GAMMA - 1.0))) diff --git a/src/venting/diagnostics.py b/src/venting/diagnostics.py new file mode 100644 index 0000000..e35a223 --- /dev/null +++ b/src/venting/diagnostics.py @@ -0,0 +1,179 @@ +import math +from dataclasses import asdict + +import numpy as np + +from .cases import CaseConfig, SolveResult +from .constants import C_CHOKED, P0, PI_C, R_GAS, T0, T_SAFE +from .graph import EXT_NODE, GasNode, OrificeEdge, ShortTubeEdge, SlotChannelEdge +from .validity import evaluate_validity_flags + + +def compute_tau_exit(total_volume: float, Cd_exit: float, A_exit_total: float) -> float: + mdot_ch = Cd_exit * A_exit_total * C_CHOKED * P0 / math.sqrt(R_GAS * T0) + if mdot_ch <= 0.0: + return float("inf") + return (total_volume * P0) / (R_GAS * T0 * mdot_ch) + + +def summarize_result( + nodes: list[GasNode], edges: list, bcs: list, case: CaseConfig, sol +) -> SolveResult: + nodes_use = getattr(sol, "nodes_local", nodes) + edges_use = getattr(sol, "edges_local", edges) + bcs_use = getattr(sol, "bcs_local", bcs) + N = len(nodes_use) + V = np.array([n.V for n in nodes_use], dtype=float) + t = sol.t + + if case.thermo == "isothermal": + m = sol.y[:N] + T = np.full_like(m, T0) + if getattr(sol, "ext_idx", None) is not None: + T[sol.ext_idx, :] = case.T_ext + else: + m = sol.y[:N] + T = sol.y[N : 2 * N] + if getattr(sol, "ext_idx", None) is not None: + T[sol.ext_idx, :] = case.T_ext + T_eff = np.maximum(T, T_SAFE) + P = np.maximum(m, 0.0) * R_GAS * T_eff / V[:, None] + + ext_idx = getattr(sol, "ext_idx", None) + if case.external_model == "dynamic_pump" and ext_idx is not None: + P_ext = P[ext_idx, :] + else: + profile = bcs_use[0].profile + P_ext = profile.P_array(t) + + dP_edges: dict[str, np.ndarray] = {} + max_dP: dict[str, float] = {} + peak_diag: dict[str, dict] = {} + + total_volume = float(np.sum(V if ext_idx is None else np.delete(V, ext_idx))) + A_exit_total = None + Cd_exit = None + for e in edges_use: + if isinstance(e, (OrificeEdge, ShortTubeEdge)) and "exit" in e.label.lower(): + A_exit_total = e.A_total + Cd_exit = e.Cd_model() + break + if A_exit_total is None: + A_exit_total = max( + ( + e.A_total + for e in edges_use + if isinstance(e, (OrificeEdge, ShortTubeEdge)) + ), + default=0.0, + ) + Cd_exit = 0.62 + tau_exit = compute_tau_exit(total_volume, float(Cd_exit), float(A_exit_total)) + + for e in edges_use: + if isinstance(e, (OrificeEdge, ShortTubeEdge)): + label = e.label or f"edge({e.a}->{e.b})" + if e.b == EXT_NODE: + dp = P[e.a, :] - P_ext + elif ext_idx is not None and e.b == ext_idx: + dp = P[e.a, :] - P[e.b, :] + else: + dp = P[e.a, :] - P[e.b, :] + dP_edges[label] = dp + elif isinstance(e, SlotChannelEdge): + label = e.label or f"slot({e.a}->{e.b})" + dP_edges[label] = P[e.a, :] - P[e.b, :] + + for label, dp in dP_edges.items(): + idx = int(np.argmax(np.abs(dp))) + max_dP[label] = float(np.abs(dp[idx])) + edge_obj = next( + e + for e in edges_use + if ((e.label or f"edge({e.a}->{e.b})") == label) + or ((e.label or f"slot({e.a}->{e.b})") == label) + ) + + a, b = edge_obj.a, edge_obj.b + Pa = float(P[a, idx]) + if b == EXT_NODE: + Pb = float(P_ext[idx]) + else: + Pb = float(P[b, idx]) + + if Pa >= Pb: + p_up, p_dn = Pa, Pb + t_up = float(T_eff[a, idx]) + else: + p_up, p_dn = Pb, Pa + if b == EXT_NODE: + t_up = case.T_ext + else: + t_up = float(T_eff[b, idx]) + + regime = "viscous_slot" + if isinstance(edge_obj, (OrificeEdge, ShortTubeEdge)): + r_pk = (p_dn / p_up) if p_up > 0 else 0.0 + regime = "CHOKED" if r_pk <= PI_C else "subsonic" + else: + r_pk = float(min(Pa, Pb) / max(Pa, Pb)) if max(Pa, Pb) > 0 else 0.0 + + if case.external_model == "dynamic_pump": + peak_type = "internal" + else: + profile = bcs_use[0].profile + peak_type = profile.classify_peak(float(t[idx]), float(t[-1])) + + peak_diag[label] = { + "t_peak": float(t[idx]), + "t_peak_over_tau_exit": ( + float(t[idx] / tau_exit) if np.isfinite(tau_exit) else float("nan") + ), + "P_up": p_up, + "P_down": p_dn, + "r": float(r_pk), + "regime": regime, + "T_up": t_up, + "peak_type": peak_type, + "dP_signed": float(dp[idx]), + } + + validity_flags = evaluate_validity_flags( + nodes=nodes_use, + edges=edges_use, + P=P, + T=T_eff, + m=m, + P_ext=P_ext, + t=t, + l_char_m=case.l_char_m, + ) + + meta = { + "case": asdict(case), + "nodes": [asdict(n) for n in nodes_use], + "edges": [asdict(e) for e in edges_use], + "profile": ( + {"name": bcs_use[0].profile.name, "events": list(bcs_use[0].profile.events)} + if bcs_use + else {"name": "dynamic_pump", "events": []} + ), + "solver": { + "success": bool(sol.success), + "message": str(sol.message), + "t_end": float(sol.t[-1]), + "n_steps": int(len(sol.t)), + }, + "validity_flags": validity_flags, + } + return SolveResult( + t=t, + m=m, + T=T, + P=P, + P_ext=P_ext, + peak_diag=peak_diag, + max_dP=max_dP, + tau_exit=tau_exit, + meta=meta, + ) diff --git a/src/venting/flow.py b/src/venting/flow.py new file mode 100644 index 0000000..9e47b47 --- /dev/null +++ b/src/venting/flow.py @@ -0,0 +1,375 @@ +import math + +from .constants import GAMMA, R_GAS, T_SAFE + + +def mu_air_sutherland(T: float) -> float: + t_eff = max(T, T_SAFE) + mu0 = 1.716e-5 + t_ref = 273.15 + s = 111.0 + return mu0 * (t_eff / t_ref) ** 1.5 * (t_ref + s) / (t_eff + s) + + +def mdot_orifice_pos_props( + P_up: float, + T_up: float, + P_dn: float, + Cd: float, + A: float, + gamma: float = GAMMA, + r_gas: float = R_GAS, +) -> float: + if P_up <= 0.0 or A <= 0.0: + return 0.0 + t_eff = max(T_up, T_SAFE) + r = max(P_dn, 0.0) / P_up + if r >= 1.0: + return 0.0 + + pi_c = (2.0 / (gamma + 1.0)) ** (gamma / (gamma - 1.0)) + c_choked = math.sqrt( + gamma * (2.0 / (gamma + 1.0)) ** ((gamma + 1.0) / (gamma - 1.0)) + ) + + if r <= pi_c: + return Cd * A * P_up * c_choked / math.sqrt(r_gas * t_eff) + + bracket = r ** (2.0 / gamma) - r ** ((gamma + 1.0) / gamma) + if bracket <= 0.0: + return 0.0 + return ( + Cd + * A + * P_up + * math.sqrt(2.0 * gamma / ((gamma - 1.0) * r_gas * t_eff) * bracket) + ) + + +def mdot_orifice_pos( + P_up: float, T_up: float, P_dn: float, Cd: float, A: float +) -> float: + return mdot_orifice_pos_props(P_up, T_up, P_dn, Cd, A, gamma=GAMMA, r_gas=R_GAS) + + +def mdot_slot_pos( + P_up: float, T_up: float, P_dn: float, w: float, delta: float, L: float +) -> float: + if P_up <= 0.0: + return 0.0 + t_eff = max(T_up, T_SAFE) + mu = mu_air_sutherland(t_eff) + K = w * (delta**3) / (12.0 * mu * L) + dp2 = max(P_up**2 - max(P_dn, 0.0) ** 2, 0.0) + return K * dp2 / (R_GAS * t_eff) + + +def friction_factor(Re: float, eps_over_D: float) -> float: + """Darcy friction factor (not Fanning).""" + if Re <= 0.0: + return 0.0 + if Re < 2300.0: + return 64.0 / Re + term = eps_over_D / 3.7 + 5.74 / (Re**0.9) + return 0.25 / (math.log10(max(term, 1e-20)) ** 2) + + +def mdot_short_tube_pos( + P_up: float, + T_up: float, + P_dn: float, + Cd0: float, + A_total: float, + D: float, + L: float, + eps: float, + K_in: float, + K_out: float, + gamma: float = GAMMA, + r_gas: float = R_GAS, +) -> float: + """Lossy nozzle via effective Cd for short-tube (thick wall). + + This model applies additional minor/friction losses via Cd_eff and does not + model Fanno friction choking. + """ + if P_up <= 0.0 or A_total <= 0.0 or D <= 0.0: + return 0.0 + + t_eff = max(T_up, T_SAFE) + if L <= 0.0 and K_in <= 0.0 and K_out <= 0.0 and eps <= 0.0: + return mdot_orifice_pos_props( + P_up, t_eff, P_dn, Cd0, A_total, gamma=gamma, r_gas=r_gas + ) + + Cd_eff = float(Cd0) + mdot = 0.0 + for _ in range(5): + mdot_new = mdot_orifice_pos_props( + P_up, t_eff, P_dn, Cd_eff, A_total, gamma=gamma, r_gas=r_gas + ) + rho = P_up / (r_gas * t_eff) + u = mdot_new / max(rho * A_total, 1e-18) + mu = mu_air_sutherland(t_eff) + Re = rho * u * D / max(mu, 1e-18) + f_D = friction_factor(Re, eps / max(D, 1e-12)) + K_fric = f_D * (L / max(D, 1e-12)) + K_tot = K_in + K_out + K_fric + Cd_eff_new = 1.0 / math.sqrt(max(Cd0 ** (-2.0) + K_tot, 1e-18)) + + if mdot > 0 and abs(mdot_new - mdot) / max(mdot_new, 1e-18) < 0.01: + mdot = mdot_new + break + mdot = 0.5 * mdot + 0.5 * mdot_new + Cd_eff = 0.5 * Cd_eff + 0.5 * Cd_eff_new + + return max(mdot, 0.0) + + +def _fanno_length_to_sonic(M: float, gamma: float) -> float: + m2 = max(M * M, 1e-12) + term1 = (1.0 - m2) / (gamma * m2) + term2 = ((gamma + 1.0) / (2.0 * gamma)) * math.log( + ((gamma + 1.0) * m2) / (2.0 + (gamma - 1.0) * m2) + ) + return term1 + term2 + + +def _p_over_pstar(M: float, gamma: float) -> float: + return (1.0 / max(M, 1e-12)) * math.sqrt( + (gamma + 1.0) / (2.0 + (gamma - 1.0) * M * M) + ) + + +def _solve_m2_from_m1(m1: float, f4ld: float, gamma: float) -> tuple[float, bool]: + f1 = _fanno_length_to_sonic(m1, gamma) + target = f1 - f4ld + if target <= 0.0: + return 1.0, True + lo = m1 + hi = 1.0 - 1e-7 + for _ in range(60): + mid = 0.5 * (lo + hi) + if _fanno_length_to_sonic(mid, gamma) > target: + lo = mid + else: + hi = mid + return 0.5 * (lo + hi), False + + +def _fanno_state( + P_up: float, + T_up: float, + P_dn: float, + Cd0: float, + A_total: float, + D: float, + L: float, + eps: float, + K_in: float, + K_out: float, + gamma: float, + r_gas: float, +): + if P_up <= 0.0 or A_total <= 0.0 or D <= 0.0: + return {"mdot": 0.0, "choked": False, "mach_exit": 0.0, "f_D": 0.0} + + t_eff = max(T_up, T_SAFE) + if L <= 0.0: + md = mdot_orifice_pos_props(P_up, t_eff, P_dn, Cd0, A_total, gamma, r_gas) + return {"mdot": md, "choked": False, "mach_exit": 0.0, "f_D": 0.0} + + if max(P_dn, 0.0) / P_up >= 0.6: + md = mdot_short_tube_pos( + P_up, + t_eff, + P_dn, + Cd0, + A_total, + D, + L, + eps, + K_in, + K_out, + gamma, + r_gas, + ) + return {"mdot": md, "choked": False, "mach_exit": 0.0, "f_D": 0.0} + + a_eff = Cd0 * A_total + re_guess = 1e5 + m1_solution = 0.1 + m2_solution = 0.1 + choked_solution = False + f_D = 0.03 + + for _ in range(5): + f_D = max(friction_factor(re_guess, eps / max(D, 1e-12)), 1e-4) + f4ld = 4.0 * f_D * L / max(D, 1e-12) + K_in + K_out + + if f4ld < 1e-3: + md = mdot_short_tube_pos( + P_up, + t_eff, + P_dn, + Cd0, + A_total, + D, + L, + eps, + K_in, + K_out, + gamma, + r_gas, + ) + return {"mdot": md, "choked": False, "mach_exit": 0.0, "f_D": f_D} + + lo = 1e-5 + hi = 1.0 - 1e-7 + for _i in range(80): + mid = 0.5 * (lo + hi) + if _fanno_length_to_sonic(mid, gamma) > f4ld: + lo = mid + else: + hi = mid + m1_cap = 0.5 * (lo + hi) + + def outlet_pressure( + m1: float, f4ld_local: float = f4ld + ) -> tuple[float, float, bool]: + m2, choked = _solve_m2_from_m1(m1, f4ld_local, gamma) + p1 = P_up / (1.0 + 0.5 * (gamma - 1.0) * m1 * m1) ** (gamma / (gamma - 1.0)) + t1 = t_eff / (1.0 + 0.5 * (gamma - 1.0) * m1 * m1) + p2 = p1 * (_p_over_pstar(m2, gamma) / _p_over_pstar(m1, gamma)) + p_out = p2 / (1.0 + K_out) + if K_in > 0.0: + p_out /= 1.0 + K_in + return p_out, t1, choked + + p_cap, _t1_cap, _ = outlet_pressure(m1_cap) + if max(P_dn, 0.0) <= p_cap: + m1_solution = m1_cap + m2_solution = 1.0 + choked_solution = True + else: + lo = 1e-5 + hi = m1_cap + for _j in range(60): + mid = 0.5 * (lo + hi) + p_mid, _, _ = outlet_pressure(mid) + if p_mid > P_dn: + lo = mid + else: + hi = mid + m1_solution = 0.5 * (lo + hi) + m2_solution, choked_solution = _solve_m2_from_m1(m1_solution, f4ld, gamma) + + t1 = t_eff / (1.0 + 0.5 * (gamma - 1.0) * m1_solution * m1_solution) + mdot = P_up * m1_solution * a_eff * math.sqrt(gamma / (r_gas * max(t1, T_SAFE))) + rho1 = P_up / (r_gas * max(t_eff, T_SAFE)) + u1 = mdot / max(rho1 * A_total, 1e-18) + re_new = rho1 * u1 * D / max(mu_air_sutherland(t_eff), 1e-18) + if abs(re_new - re_guess) / max(re_guess, 1.0) < 0.05: + re_guess = re_new + break + re_guess = 0.5 * re_guess + 0.5 * re_new + + mdot = ( + P_up + * m1_solution + * a_eff + * math.sqrt( + gamma + / ( + r_gas + * max( + t_eff / (1.0 + 0.5 * (gamma - 1.0) * m1_solution * m1_solution), + T_SAFE, + ) + ) + ) + ) + return { + "mdot": max(mdot, 0.0), + "choked": bool(choked_solution), + "mach_exit": float(m2_solution), + "f_D": float(f_D), + } + + +def mdot_fanno_tube( + P_up: float, + T_up: float, + P_dn: float, + Cd0: float, + A_total: float, + D: float, + L: float, + eps: float, + K_in: float, + K_out: float, + gamma: float = GAMMA, + r_gas: float = R_GAS, +) -> float: + """Mass flow in a short tube including Fanno friction choking. + + Reference: Shapiro, Dynamics and Thermodynamics of Compressible Fluid Flow, + 1953, Chapter 6. + """ + if L <= 0.0: + return mdot_orifice_pos_props(P_up, T_up, P_dn, Cd0, A_total, gamma, r_gas) + + md_f = _fanno_state( + P_up, + T_up, + P_dn, + Cd0, + A_total, + D, + L, + eps, + K_in, + K_out, + gamma, + r_gas, + )["mdot"] + md_l = mdot_short_tube_pos( + P_up, T_up, P_dn, Cd0, A_total, D, L, eps, K_in, K_out, gamma, r_gas + ) + if md_f < 1e-3 * max(md_l, 1e-18): + md_f = md_l + md = min(md_f, md_l) + if P_up > 0 and (P_dn / P_up) <= 0.5: + md = min(md_l, max(md, 0.99 * md_l)) + return md + + +def fanno_choked_state( + P_up: float, + T_up: float, + P_dn: float, + Cd0: float, + A_total: float, + D: float, + L: float, + eps: float, + K_in: float, + K_out: float, + gamma: float = GAMMA, + r_gas: float = R_GAS, +) -> tuple[bool, float]: + state = _fanno_state( + P_up, + T_up, + P_dn, + Cd0, + A_total, + D, + L, + eps, + K_in, + K_out, + gamma, + r_gas, + ) + return bool(state["choked"]), float(state["mach_exit"]) diff --git a/src/venting/gates.py b/src/venting/gates.py new file mode 100644 index 0000000..bdff9ff --- /dev/null +++ b/src/venting/gates.py @@ -0,0 +1,106 @@ +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np + +from .cases import CaseConfig +from .constants import C_CHOKED, C_P, C_V, GAMMA, P0, R_GAS, T0, T_SAFE +from .diagnostics import summarize_result +from .flow import mdot_orifice_pos +from .geometry import circle_area_from_d_mm +from .graph import EXT_NODE, CdConst, ExternalBC, GasNode, OrificeEdge +from .profiles import Profile +from .solver import solve_case + + +@dataclass(frozen=True) +class GateMetrics: + errP: float + errT: float + errMass: float + errEnergy: float | None = None + + +def gate_single() -> GateMetrics: + V = 131.6e-6 + A = circle_area_from_d_mm(2.0) + Cd = 0.62 + alpha = Cd * A * C_CHOKED * math.sqrt(R_GAS * T0) / V + beta = (GAMMA - 1.0) / 2.0 + + def p_adi(t: float) -> float: + return P0 * (1.0 + beta * alpha * t) ** (-2.0 * GAMMA / (GAMMA - 1.0)) + + def t_adi(t: float) -> float: + return T0 * (1.0 + beta * alpha * t) ** (-2.0) + + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", V, 181.6e-4)] + edges = [OrificeEdge(0, EXT_NODE, A, CdConst(Cd), label="exit")] + bcs = [ExternalBC(0, profile, T_ext=T0)] + + case_adi = CaseConfig("intermediate", 0.0, T0, 1.5, 700) + ad = summarize_result( + nodes, edges, bcs, case_adi, solve_case(nodes, edges, bcs, case_adi) + ) + + mask = ad.P[0, :] > 0.01 * P0 + t = ad.t[mask] + errP = float( + np.max(np.abs(ad.P[0, mask] - np.array([p_adi(float(tt)) for tt in t])) / P0) + ) + errT = float( + np.max(np.abs(ad.T[0, mask] - np.array([t_adi(float(tt)) for tt in t])) / T0) + ) + + m0 = P0 * V / (R_GAS * T0) + mf = float(ad.m[0, -1]) + mdot = np.array( + [ + mdot_orifice_pos( + float(ad.P[0, k]), float(max(ad.T[0, k], T_SAFE)), 0.0, Cd, A + ) + for k in range(len(ad.t)) + ] + ) + m_out = float(np.trapz(mdot, ad.t)) + errMass = abs((mf + m_out) - m0) / m0 + + return GateMetrics(errP=errP, errT=errT, errMass=errMass) + + +def gate_two() -> GateMetrics: + Vc, Vv = 131.6e-6, 145.3e-6 + A = circle_area_from_d_mm(2.0) + Cd = 0.62 + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("vest", Vv, 181.6e-4), GasNode("cell", Vc, 181.6e-4)] + edges = [ + OrificeEdge(1, 0, A, CdConst(Cd), label="cell↔vest"), + OrificeEdge(0, EXT_NODE, A, CdConst(Cd), label="exit"), + ] + bcs = [ExternalBC(0, profile, T_ext=T0)] + case = CaseConfig("intermediate", 0.0, T0, 2.0, 900) + res = summarize_result(nodes, edges, bcs, case, solve_case(nodes, edges, bcs, case)) + + m_total_0 = float(np.sum(res.m[:, 0])) + m_total_f = float(np.sum(res.m[:, -1])) + mdot = np.array( + [ + mdot_orifice_pos( + float(res.P[0, k]), float(max(res.T[0, k], T_SAFE)), 0.0, Cd, A + ) + for k in range(len(res.t)) + ], + dtype=float, + ) + m_out = float(np.trapz(mdot, res.t)) + errMass = abs((m_total_f + m_out) - m_total_0) / m_total_0 + + E0 = float(np.sum(res.m[:, 0] * C_V * res.T[:, 0])) + Ef = float(np.sum(res.m[:, -1] * C_V * res.T[:, -1])) + Eout = float(np.trapz(mdot * C_P * np.maximum(res.T[0, :], T_SAFE), res.t)) + errE = abs((Ef + Eout) - E0) / E0 + return GateMetrics(errP=0.0, errT=0.0, errMass=errMass, errEnergy=errE) diff --git a/src/venting/geometry.py b/src/venting/geometry.py new file mode 100644 index 0000000..8321327 --- /dev/null +++ b/src/venting/geometry.py @@ -0,0 +1,23 @@ +import math + + +def mm_to_m(x_mm: float) -> float: + return x_mm * 1e-3 + + +def mm2_to_m2(x_mm2: float) -> float: + return x_mm2 * 1e-6 + + +def mm3_to_m3(x_mm3: float) -> float: + return x_mm3 * 1e-9 + + +def circle_area_from_d_mm(d_mm: float) -> float: + d_m = mm_to_m(d_mm) + return math.pi * (d_m / 2.0) ** 2 + + +def assert_pos(name: str, val: float) -> None: + if not (val > 0.0): + raise ValueError(f"{name} must be > 0, got {val}") diff --git a/src/venting/graph.py b/src/venting/graph.py new file mode 100644 index 0000000..15f95c7 --- /dev/null +++ b/src/venting/graph.py @@ -0,0 +1,347 @@ +from dataclasses import dataclass + +from .cases import NetworkConfig +from .constants import T0 +from .geometry import assert_pos, circle_area_from_d_mm +from .profiles import Profile + + +@dataclass(frozen=True) +class GasNode: + name: str + V: float + A_wall: float + + +@dataclass(frozen=True) +class ExternalBC: + node: int + profile: Profile + T_ext: float = T0 + + +@dataclass(frozen=True) +class CdConst: + Cd: float + + def __call__(self, Re: float | None = None, r: float | None = None) -> float: + return float(self.Cd) + + +@dataclass(frozen=True) +class OrificeEdge: + a: int + b: int + A_total: float + Cd_model: CdConst + label: str = "" + + +@dataclass(frozen=True) +class ShortTubeEdge: + a: int + b: int + A_total: float + D: float + L: float + eps: float + K_in: float + K_out: float + Cd_model: CdConst + label: str = "" + fanno: bool = False + + +@dataclass(frozen=True) +class SlotChannelEdge: + a: int + b: int + w: float + delta: float + L: float + label: str = "" + + +EXT_NODE = -1 +Edge = OrificeEdge | SlotChannelEdge | ShortTubeEdge + + +def _expand_scalar_or_list(value, expected_len: int, name: str) -> list[float]: + if isinstance(value, (int, float)): + return [float(value)] * expected_len + seq = list(value) + if len(seq) != expected_len: + raise ValueError( + f"{name} length mismatch: expected {expected_len}, got {len(seq)}" + ) + return [float(v) for v in seq] + + +def _resolve_int_losses(cfg: NetworkConfig, n_ifaces: int): + k_in_src = cfg.K_in_int if cfg.K_in_int is not None else cfg.K_in + k_out_src = cfg.K_out_int if cfg.K_out_int is not None else cfg.K_out + eps_src = cfg.eps_int_um if cfg.eps_int_um is not None else cfg.eps_um + k_in = _expand_scalar_or_list(k_in_src, n_ifaces, "K_in_int") + k_out = _expand_scalar_or_list(k_out_src, n_ifaces, "K_out_int") + eps = [v * 1e-6 for v in _expand_scalar_or_list(eps_src, n_ifaces, "eps_int_um")] + return k_in, k_out, eps + + +def _make_exit_edge( + cfg: NetworkConfig, a: int, A_exit_total: float, cd_exit: CdConst, label: str +) -> OrificeEdge | ShortTubeEdge: + if cfg.exit_model == "orifice" or cfg.L_exit_mm <= 0.0: + return OrificeEdge(a, EXT_NODE, A_exit_total, cd_exit, label=label) + fanno = cfg.exit_model == "fanno" + d = cfg.d_exit_mm * 1e-3 + length = cfg.L_exit_mm * 1e-3 + k_in = float(cfg.K_in_exit if cfg.K_in_exit is not None else cfg.K_in) + k_out = float(cfg.K_out_exit if cfg.K_out_exit is not None else cfg.K_out) + eps_um = float(cfg.eps_exit_um if cfg.eps_exit_um is not None else cfg.eps_um) + return ShortTubeEdge( + a=a, + b=EXT_NODE, + A_total=A_exit_total, + D=d, + L=length, + eps=eps_um * 1e-6, + K_in=k_in, + K_out=k_out, + Cd_model=cd_exit, + label=label, + fanno=fanno, + ) + + +def _make_int_edge( + cfg: NetworkConfig, + a: int, + b: int, + d_int_mm: float, + n_int_per_interface: float, + cd_int: float, + l_int_mm: float, + eps_int_m: float, + k_in: float, + k_out: float, + label: str, +) -> OrificeEdge | ShortTubeEdge: + A_int_total = cfg.N_par * n_int_per_interface * circle_area_from_d_mm(d_int_mm) + cd_model = CdConst(cd_int) + if cfg.int_model == "orifice" or l_int_mm <= 0.0: + return OrificeEdge(a, b, A_int_total, cd_model, label=label) + fanno = cfg.int_model == "fanno" + return ShortTubeEdge( + a=a, + b=b, + A_total=A_int_total, + D=d_int_mm * 1e-3, + L=l_int_mm * 1e-3, + eps=eps_int_m, + K_in=k_in, + K_out=k_out, + Cd_model=cd_model, + label=label, + fanno=fanno, + ) + + +def _build_single_chain( + cfg: NetworkConfig, profile: Profile +) -> tuple[list[GasNode], list[Edge], list[ExternalBC]]: + assert_pos("N_chain", cfg.N_chain) + assert_pos("N_par", cfg.N_par) + assert_pos("V_vest", cfg.V_vest) + + n_ifaces = cfg.N_chain + v_cells = _expand_scalar_or_list(cfg.V_cell, cfg.N_chain, "V_cell") + a_wall_cells = _expand_scalar_or_list(cfg.A_wall_cell, cfg.N_chain, "A_wall_cell") + d_ints = _expand_scalar_or_list(cfg.d_int_mm, n_ifaces, "d_int_mm") + n_ints = _expand_scalar_or_list( + cfg.n_int_per_interface, n_ifaces, "n_int_per_interface" + ) + cd_ints = _expand_scalar_or_list(cfg.Cd_int, n_ifaces, "Cd_int") + l_ints = _expand_scalar_or_list(cfg.L_int_mm, n_ifaces, "L_int_mm") + k_in_int, k_out_int, eps_int = _resolve_int_losses(cfg, n_ifaces) + + nodes: list[GasNode] = [GasNode("vest", cfg.V_vest, cfg.A_wall_vest)] + edges: list[Edge] = [] + + gap_idx = None + if cfg.use_gap: + assert_pos("V_gap", cfg.V_gap) + nodes.append(GasNode("gap", cfg.V_gap, cfg.A_wall_gap)) + gap_idx = 1 + + base = 1 if not cfg.use_gap else 2 + for i in range(cfg.N_chain): + nodes.append( + GasNode( + f"cell{i + 1}", + cfg.N_par * v_cells[i], + cfg.N_par * a_wall_cells[i], + ) + ) + + A_exit_total = cfg.n_exit * circle_area_from_d_mm(cfg.d_exit_mm) + cd_exit = CdConst(cfg.Cd_exit) + if not cfg.use_gap: + edges.append(_make_exit_edge(cfg, 0, A_exit_total, cd_exit, label="exit")) + else: + edges.append( + SlotChannelEdge( + 0, gap_idx, cfg.gap_w, cfg.gap_delta, cfg.gap_L, label="vest→gap(slot)" + ) + ) + edges.append( + _make_exit_edge(cfg, gap_idx, A_exit_total, cd_exit, label="gap→ext(exit)") + ) + + cell1 = base + edges.append( + _make_int_edge( + cfg, + cell1, + 0, + d_ints[0], + n_ints[0], + cd_ints[0], + l_ints[0], + eps_int[0], + k_in_int[0], + k_out_int[0], + label="A(cell1↔vest)", + ) + ) + for i in range(1, cfg.N_chain): + a = base + i + b = base + (i - 1) + edges.append( + _make_int_edge( + cfg, + a, + b, + d_ints[i], + n_ints[i], + cd_ints[i], + l_ints[i], + eps_int[i], + k_in_int[i], + k_out_int[i], + label=f"{chr(65 + i)}(cell{i + 1}↔cell{i})", + ) + ) + + bc_node = 0 if not cfg.use_gap else gap_idx + bcs = [ExternalBC(node=bc_node, profile=profile, T_ext=T0)] + return nodes, edges, bcs + + +def _build_two_chain_shared_vest( + cfg: NetworkConfig, profile: Profile +) -> tuple[list[GasNode], list[Edge], list[ExternalBC]]: + n_chain_a = cfg.N_chain + n_chain_b = cfg.N_chain_b if cfg.N_chain_b is not None else cfg.N_chain + if n_chain_a <= 0 or n_chain_b <= 0: + raise ValueError("N_chain and N_chain_b must be > 0") + + n_ifaces = max(n_chain_a, n_chain_b) + v_cells = _expand_scalar_or_list(cfg.V_cell, n_ifaces, "V_cell") + a_wall_cells = _expand_scalar_or_list(cfg.A_wall_cell, n_ifaces, "A_wall_cell") + d_ints = _expand_scalar_or_list(cfg.d_int_mm, n_ifaces, "d_int_mm") + n_ints = _expand_scalar_or_list( + cfg.n_int_per_interface, n_ifaces, "n_int_per_interface" + ) + cd_ints = _expand_scalar_or_list(cfg.Cd_int, n_ifaces, "Cd_int") + l_ints = _expand_scalar_or_list(cfg.L_int_mm, n_ifaces, "L_int_mm") + k_in_int, k_out_int, eps_int = _resolve_int_losses(cfg, n_ifaces) + + nodes: list[GasNode] = [GasNode("vest", cfg.V_vest, cfg.A_wall_vest)] + edges: list[Edge] = [] + + A_exit_total = cfg.n_exit * circle_area_from_d_mm(cfg.d_exit_mm) + edges.append( + _make_exit_edge(cfg, 0, A_exit_total, CdConst(cfg.Cd_exit), label="exit") + ) + + base_a = len(nodes) + for i in range(n_chain_a): + nodes.append( + GasNode( + f"A_cell{i + 1}", + cfg.N_par * v_cells[i], + cfg.N_par * a_wall_cells[i], + ) + ) + base_b = len(nodes) + for i in range(n_chain_b): + nodes.append( + GasNode( + f"B_cell{i + 1}", + cfg.N_par * v_cells[i], + cfg.N_par * a_wall_cells[i], + ) + ) + + # chain A + for i in range(n_chain_a): + a = base_a + i + b = 0 if i == 0 else (base_a + i - 1) + edges.append( + _make_int_edge( + cfg, + a, + b, + d_ints[i], + n_ints[i], + cd_ints[i], + l_ints[i], + eps_int[i], + k_in_int[i], + k_out_int[i], + label=f"A{i + 1}", + ) + ) + + # chain B + for i in range(n_chain_b): + a = base_b + i + b = 0 if i == 0 else (base_b + i - 1) + edges.append( + _make_int_edge( + cfg, + a, + b, + d_ints[i], + n_ints[i], + cd_ints[i], + l_ints[i], + eps_int[i], + k_in_int[i], + k_out_int[i], + label=f"B{i + 1}", + ) + ) + + return nodes, edges, [ExternalBC(node=0, profile=profile, T_ext=T0)] + + +def build_branching_network( + cfg: NetworkConfig, profile: Profile +) -> tuple[list[GasNode], list[Edge], list[ExternalBC]]: + if cfg.int_model not in {"orifice", "short_tube", "fanno"}: + raise ValueError("int_model must be 'orifice', 'short_tube', or 'fanno'") + if cfg.exit_model not in {"orifice", "short_tube", "fanno"}: + raise ValueError("exit_model must be 'orifice', 'short_tube', or 'fanno'") + + if cfg.topology == "single_chain": + return _build_single_chain(cfg, profile) + if cfg.topology == "two_chain_shared_vest": + return _build_two_chain_shared_vest(cfg, profile) + raise ValueError("topology must be 'single_chain' or 'two_chain_shared_vest'") + + +def build_network( + cfg: NetworkConfig, profile: Profile +) -> tuple[list[GasNode], list[Edge], list[ExternalBC]]: + return build_branching_network(cfg, profile) diff --git a/src/venting/gui/__init__.py b/src/venting/gui/__init__.py new file mode 100644 index 0000000..a15240e --- /dev/null +++ b/src/venting/gui/__init__.py @@ -0,0 +1 @@ +"""GUI package for venting v10.""" diff --git a/src/venting/gui/app.py b/src/venting/gui/app.py new file mode 100644 index 0000000..7e9822c --- /dev/null +++ b/src/venting/gui/app.py @@ -0,0 +1,400 @@ +from __future__ import annotations + +from dataclasses import asdict +from pathlib import Path + +from venting.diagnostics import summarize_result +from venting.graph import build_branching_network +from venting.presets import get_default_panel_preset_v9 +from venting.profiles import ( + make_profile_exponential, + make_profile_from_table, + make_profile_linear, + make_profile_step, +) +from venting.run import export_case_artifacts, make_case_output_dir +from venting.solver import solve_case_stream + +from .config import GuiCaseConfig +from .state_layout import infer_layout_from_modes + + +def load_gui_deps(): + import pyqtgraph as pg + from PySide6 import QtCore, QtWidgets + + return QtCore, QtWidgets, pg + + +def _profile_from_cfg(cfg: GuiCaseConfig): + if cfg.external_model == "dynamic_pump": + return make_profile_step(101325.0, 1e9) + if cfg.profile_kind == "linear": + return make_profile_linear(101325.0, cfg.rate_mmhg_per_s) + if cfg.profile_kind == "step": + return make_profile_step(101325.0, cfg.step_time_s) + if cfg.profile_kind == "barometric": + return make_profile_exponential(101325.0, cfg.rate_mmhg_per_s, p_floor=10.0) + return make_profile_from_table( + "gui_table", + Path(cfg.profile_file), + pressure_unit=cfg.profile_pressure_unit, + ) + + +def create_main_window(): + QtCore, QtWidgets, pg = load_gui_deps() + + class SolverWorker(QtCore.QThread): + progress = QtCore.Signal(object) + finished_result = QtCore.Signal(object) + failed = QtCore.Signal(str) + + def __init__(self, cfg: GuiCaseConfig): + super().__init__() + self.cfg = cfg + self.stop_requested = False + + def request_stop(self): + self.stop_requested = True + + def run(self): + try: + profile = _profile_from_cfg(self.cfg) + net = self.cfg.to_network_config() + case = self.cfg.to_case_config() + nodes, edges, bcs = build_branching_network(net, profile) + + def on_chunk(payload): + if self.stop_requested: + return + self.progress.emit(payload) + + sol = solve_case_stream( + nodes, + edges, + bcs, + case, + callback=on_chunk, + dt_chunk_s=1.0, + should_stop=lambda: self.stop_requested, + ) + if self.stop_requested or not sol.success: + return + res = summarize_result(nodes, edges, bcs, case, sol) + self.finished_result.emit( + {"res": res, "cfg": self.cfg, "profile": profile} + ) + except Exception as exc: # pragma: no cover - UI path + self.failed.emit(str(exc)) + + class MainWindow(QtWidgets.QMainWindow): + def __init__(self) -> None: + super().__init__() + self.setWindowTitle("Venting v10 GUI") + self.resize(1400, 840) + self.worker = None + self.latest_res = None + + preset = get_default_panel_preset_v9() + self.cfg = GuiCaseConfig( + V_cell_m3=preset.V_cell, + V_vest_m3=preset.V_vest, + A_wall_cell_m2=preset.A_wall_cell, + A_wall_vest_m2=preset.A_wall_vest, + ) + + root = QtWidgets.QWidget() + root_layout = QtWidgets.QVBoxLayout(root) + self.setCentralWidget(root) + + toolbar = QtWidgets.QHBoxLayout() + self.btn_run = QtWidgets.QPushButton("Run") + self.btn_stop = QtWidgets.QPushButton("Stop") + self.btn_open = QtWidgets.QPushButton("Open config") + self.btn_save = QtWidgets.QPushButton("Save config") + self.status = QtWidgets.QLabel("Ready") + for w in [ + self.btn_run, + self.btn_stop, + self.btn_open, + self.btn_save, + self.status, + ]: + toolbar.addWidget(w) + toolbar.addStretch(1) + root_layout.addLayout(toolbar) + + split = QtWidgets.QSplitter() + root_layout.addWidget(split, 1) + + # left config panel + left_scroll = QtWidgets.QScrollArea() + left_scroll.setWidgetResizable(True) + left_body = QtWidgets.QWidget() + form = QtWidgets.QFormLayout(left_body) + left_scroll.setWidget(left_body) + split.addWidget(left_scroll) + + self.fields = {} + + def add_line(key, value): + w = QtWidgets.QLineEdit(str(value)) + self.fields[key] = w + form.addRow(key, w) + + def add_combo(key, items, value): + w = QtWidgets.QComboBox() + w.addItems(items) + w.setCurrentText(value) + self.fields[key] = w + form.addRow(key, w) + + # required controls + add_combo( + "int_model", ["orifice", "short_tube", "fanno"], self.cfg.int_model + ) + add_combo( + "exit_model", ["orifice", "short_tube", "fanno"], self.cfg.exit_model + ) + add_line("L_int_mm", self.cfg.L_int_mm) + add_line("L_exit_mm", self.cfg.L_exit_mm) + add_line("eps_int_um", self.cfg.eps_int_um) + add_line("eps_exit_um", self.cfg.eps_exit_um) + add_line("K_in_int", self.cfg.K_in_int) + add_line("K_out_int", self.cfg.K_out_int) + add_line("K_in_exit", self.cfg.K_in_exit) + add_line("K_out_exit", self.cfg.K_out_exit) + add_combo( + "topology", ["single_chain", "two_chain_shared_vest"], self.cfg.topology + ) + add_line("N_chain", self.cfg.N_chain) + add_line("N_chain_b", self.cfg.N_chain_b) + add_line("N_par", self.cfg.N_par) + add_line("V_cell_m3", self.cfg.V_cell_m3) + add_line("V_vest_m3", self.cfg.V_vest_m3) + add_line("A_wall_cell_m2", self.cfg.A_wall_cell_m2) + add_line("A_wall_vest_m2", self.cfg.A_wall_vest_m2) + add_line("d_int_mm", self.cfg.d_int_mm) + add_line("n_int_per_interface", self.cfg.n_int_per_interface) + add_line("d_exit_mm", self.cfg.d_exit_mm) + add_line("n_exit", self.cfg.n_exit) + add_line("Cd_int", self.cfg.Cd_int) + add_line("Cd_exit", self.cfg.Cd_exit) + add_combo( + "external_model", ["profile", "dynamic_pump"], self.cfg.external_model + ) + add_combo( + "profile_kind", + ["linear", "step", "barometric", "table"], + self.cfg.profile_kind, + ) + add_line("profile_file", self.cfg.profile_file) + add_combo( + "profile_pressure_unit", ["Pa", "mmHg"], self.cfg.profile_pressure_unit + ) + add_line("rate_mmhg_per_s", self.cfg.rate_mmhg_per_s) + add_line("step_time_s", self.cfg.step_time_s) + add_combo( + "thermo", ["isothermal", "intermediate", "variable"], self.cfg.thermo + ) + add_combo("wall_model", ["fixed", "lumped"], self.cfg.wall_model) + add_line("h_conv_W_m2K", self.cfg.h_conv_W_m2K) + add_line("wall_C_per_area_J_m2K", self.cfg.wall_C_per_area_J_m2K) + add_line("wall_h_out_W_m2K", self.cfg.wall_h_out_W_m2K) + add_line("wall_T_inf_K", self.cfg.wall_T_inf_K) + add_line("wall_emissivity", self.cfg.wall_emissivity) + add_line("wall_T_sur_K", self.cfg.wall_T_sur_K) + add_line("wall_q_flux_W_m2", self.cfg.wall_q_flux_W_m2) + add_line("duration_s", self.cfg.duration_s) + add_line("n_pts", self.cfg.n_pts) + add_line("V_ext_m3", self.cfg.V_ext_m3) + add_line("pump_speed_m3s", self.cfg.pump_speed_m3s) + add_line("P_ult_Pa", self.cfg.P_ult_Pa) + add_line("T_ext_K", self.cfg.T_ext_K) + add_line("output_case_name", self.cfg.output_case_name) + + # right tabs/plots + tabs = QtWidgets.QTabWidget() + split.addWidget(tabs) + + self.plot_p = pg.PlotWidget(title="Pressure") + self.plot_t = pg.PlotWidget(title="Temperature") + self.plot_m = pg.PlotWidget(title="Mass") + self.table_peaks = QtWidgets.QTableWidget(0, 4) + self.table_peaks.setHorizontalHeaderLabels( + ["edge", "|ΔP|max", "t_peak", "regime"] + ) + self.table_valid = QtWidgets.QTableWidget(0, 3) + self.table_valid.setHorizontalHeaderLabels(["flag", "status", "message"]) + + tabs.addTab(self.plot_p, "P(t)") + tabs.addTab(self.plot_t, "T(t)") + tabs.addTab(self.plot_m, "m(t)") + tabs.addTab(self.table_peaks, "ΔP/peaks") + tabs.addTab(self.table_valid, "Validity") + + self.btn_run.clicked.connect(self.on_run) + self.btn_stop.clicked.connect(self.on_stop) + self.btn_open.clicked.connect(self.on_open) + self.btn_save.clicked.connect(self.on_save) + + def _read_cfg(self) -> GuiCaseConfig: + kwargs = {} + for key, widget in self.fields.items(): + if hasattr(widget, "currentText"): + value = widget.currentText() + else: + value = widget.text() + kwargs[key] = value + + int_fields = { + "N_chain", + "N_chain_b", + "N_par", + "n_int_per_interface", + "n_exit", + "n_pts", + } + float_fields = ( + set(kwargs) + - int_fields + - { + "int_model", + "exit_model", + "external_model", + "profile_kind", + "profile_pressure_unit", + "thermo", + "wall_model", + "profile_file", + "output_case_name", + } + ) + for k in int_fields: + kwargs[k] = int(kwargs[k]) + for k in float_fields: + kwargs[k] = float(kwargs[k]) + cfg = GuiCaseConfig(**kwargs) + cfg.validate() + return cfg + + def on_run(self): + try: + cfg = self._read_cfg() + except Exception as exc: + self.status.setText(f"Invalid config: {exc}") + return + self.status.setText("Running...") + self.plot_p.clear() + self.plot_t.clear() + self.plot_m.clear() + self.worker = SolverWorker(cfg) + self.worker.progress.connect(self.on_progress) + self.worker.finished_result.connect(self.on_finished) + self.worker.failed.connect( + lambda msg: self.status.setText(f"Failed: {msg}") + ) + self.worker.start() + + def on_stop(self): + if self.worker is not None: + self.worker.request_stop() + self.status.setText("Stopping requested...") + + def on_progress(self, payload): + t = payload["t"] + y = payload["y"] + n_nodes = int(payload.get("node_count", payload.get("n_nodes", 1))) + layout = infer_layout_from_modes( + payload.get("thermo", "isothermal"), + payload.get("wall_model", "fixed"), + n_nodes, + ) + + m = y[layout.m_slice, :] + self.plot_m.clear() + self.plot_m.plot(t, m[0], pen="y") + + self.plot_t.clear() + if layout.t_slice is not None: + t_arr = y[layout.t_slice, :] + self.plot_t.plot(t, t_arr[0], pen="c") + + self.status.setText(f"Running... {100 * payload['progress']:.0f}%") + + def on_finished(self, payload): + self.latest_res = payload["res"] + cfg = payload["cfg"] + res = payload["res"] + self.plot_p.clear() + self.plot_p.plot(res.t, res.P[0], pen="g") + self.plot_t.clear() + self.plot_t.plot(res.t, res.T[0], pen="c") + self.plot_m.clear() + self.plot_m.plot(res.t, res.m[0], pen="y") + + self._fill_tables(res) + + out = make_case_output_dir(cfg.output_case_name) + run_params = asdict(cfg) + stem = f"v1000_gui_{cfg.external_model}_{cfg.profile_kind}_{cfg.thermo}" + export_case_artifacts(out, stem, res, run_params) + self.status.setText(f"Done. Saved to {out}") + + def _fill_tables(self, res): + peaks = res.peak_diag + self.table_peaks.setRowCount(len(peaks)) + for i, (edge, peak) in enumerate(peaks.items()): + self.table_peaks.setItem(i, 0, QtWidgets.QTableWidgetItem(edge)) + self.table_peaks.setItem( + i, 1, QtWidgets.QTableWidgetItem(f"{res.max_dP.get(edge, 0.0):.2f}") + ) + self.table_peaks.setItem( + i, 2, QtWidgets.QTableWidgetItem(f"{peak.get('t_peak', 0.0):.4g}") + ) + self.table_peaks.setItem( + i, 3, QtWidgets.QTableWidgetItem(str(peak.get("regime", ""))) + ) + + flags = res.meta.get("validity_flags", {}) + self.table_valid.setRowCount(len(flags)) + for i, (name, flag) in enumerate(flags.items()): + self.table_valid.setItem(i, 0, QtWidgets.QTableWidgetItem(name)) + self.table_valid.setItem( + i, 1, QtWidgets.QTableWidgetItem(str(flag.get("status", ""))) + ) + self.table_valid.setItem( + i, 2, QtWidgets.QTableWidgetItem(str(flag.get("message", ""))) + ) + + def on_open(self): + path, _ = QtWidgets.QFileDialog.getOpenFileName( + self, "Open config", "", "JSON (*.json)" + ) + if not path: + return + cfg = GuiCaseConfig.load_json(path) + for key, widget in self.fields.items(): + val = getattr(cfg, key) + if hasattr(widget, "setCurrentText"): + widget.setCurrentText(str(val)) + else: + widget.setText(str(val)) + self.status.setText(f"Loaded {path}") + + def on_save(self): + try: + cfg = self._read_cfg() + except Exception as exc: + self.status.setText(f"Cannot save invalid config: {exc}") + return + path, _ = QtWidgets.QFileDialog.getSaveFileName( + self, "Save config", "case.json", "JSON (*.json)" + ) + if not path: + return + cfg.save_json(path) + self.status.setText(f"Saved {path}") + + return MainWindow() diff --git a/src/venting/gui/config.py b/src/venting/gui/config.py new file mode 100644 index 0000000..bb87593 --- /dev/null +++ b/src/venting/gui/config.py @@ -0,0 +1,176 @@ +from __future__ import annotations + +import json +from dataclasses import asdict, dataclass +from pathlib import Path + +from venting.cases import CaseConfig, NetworkConfig + +ALLOWED_THERMO = {"isothermal", "intermediate", "variable"} +ALLOWED_WALL = {"fixed", "lumped"} +ALLOWED_EXTERNAL = {"profile", "dynamic_pump"} +ALLOWED_EDGE_MODEL = {"orifice", "short_tube", "fanno"} +ALLOWED_PROFILE = {"linear", "step", "barometric", "table"} + + +@dataclass +class GuiCaseConfig: + # geometry/network (SI unless suffix says otherwise) + topology: str = "single_chain" + N_chain: int = 10 + N_chain_b: int = 10 + N_par: int = 2 + V_cell_m3: float = 1.0e-4 + V_vest_m3: float = 1.0e-4 + A_wall_cell_m2: float = 0.01 + A_wall_vest_m2: float = 0.01 + d_int_mm: float = 2.0 + n_int_per_interface: int = 1 + d_exit_mm: float = 2.0 + n_exit: int = 1 + Cd_int: float = 0.62 + Cd_exit: float = 0.62 + int_model: str = "orifice" + exit_model: str = "orifice" + L_int_mm: float = 0.0 + L_exit_mm: float = 0.0 + K_in_int: float = 0.5 + K_out_int: float = 1.0 + eps_int_um: float = 0.0 + K_in_exit: float = 0.5 + K_out_exit: float = 1.0 + eps_exit_um: float = 0.0 + + # boundary/profile + external_model: str = "profile" + profile_kind: str = "linear" + profile_file: str = "" + profile_pressure_unit: str = "Pa" + rate_mmhg_per_s: float = 20.0 + step_time_s: float = 0.01 + + # thermal/case + thermo: str = "isothermal" + duration_s: float = 150.0 + n_pts: int = 800 + h_conv_W_m2K: float = 0.0 + T_wall_K: float = 300.0 + wall_model: str = "fixed" + wall_C_per_area_J_m2K: float = 1e9 + wall_h_out_W_m2K: float = 0.0 + wall_T_inf_K: float = 300.0 + wall_emissivity: float = 0.0 + wall_T_sur_K: float = 300.0 + wall_q_flux_W_m2: float = 0.0 + + # dynamic pump + V_ext_m3: float = 0.1 + T_ext_K: float = 300.0 + pump_speed_m3s: float = 0.0 + P_ult_Pa: float = 0.0 + + # output + output_case_name: str = "gui" + + def validate(self) -> None: + if self.thermo not in ALLOWED_THERMO: + raise ValueError("thermo is invalid") + if self.wall_model not in ALLOWED_WALL: + raise ValueError("wall_model is invalid") + if self.external_model not in ALLOWED_EXTERNAL: + raise ValueError("external_model is invalid") + if self.topology not in {"single_chain", "two_chain_shared_vest"}: + raise ValueError("topology is invalid") + if ( + self.int_model not in ALLOWED_EDGE_MODEL + or self.exit_model not in ALLOWED_EDGE_MODEL + ): + raise ValueError("edge model must be orifice|short_tube|fanno") + if self.profile_kind not in ALLOWED_PROFILE: + raise ValueError("profile_kind is invalid") + for name in [ + "N_chain", + "N_par", + "N_chain_b", + "n_int_per_interface", + "n_exit", + "n_pts", + ]: + if getattr(self, name) <= 0: + raise ValueError(f"{name} must be > 0") + for name in [ + "V_cell_m3", + "V_vest_m3", + "A_wall_cell_m2", + "A_wall_vest_m2", + "d_int_mm", + "d_exit_mm", + "duration_s", + ]: + if getattr(self, name) <= 0.0: + raise ValueError(f"{name} must be > 0") + + def to_network_config(self) -> NetworkConfig: + self.validate() + return NetworkConfig( + N_chain=self.N_chain, + N_chain_b=self.N_chain_b, + topology=self.topology, + N_par=self.N_par, + V_cell=self.V_cell_m3, + V_vest=self.V_vest_m3, + A_wall_cell=self.A_wall_cell_m2, + A_wall_vest=self.A_wall_vest_m2, + d_int_mm=self.d_int_mm, + n_int_per_interface=self.n_int_per_interface, + d_exit_mm=self.d_exit_mm, + n_exit=self.n_exit, + Cd_int=self.Cd_int, + Cd_exit=self.Cd_exit, + int_model=self.int_model, + exit_model=self.exit_model, + L_int_mm=self.L_int_mm, + L_exit_mm=self.L_exit_mm, + K_in_int=self.K_in_int, + K_out_int=self.K_out_int, + eps_int_um=self.eps_int_um, + K_in_exit=self.K_in_exit, + K_out_exit=self.K_out_exit, + eps_exit_um=self.eps_exit_um, + ) + + def to_case_config(self) -> CaseConfig: + self.validate() + return CaseConfig( + thermo=self.thermo, + h_conv=self.h_conv_W_m2K, + T_wall=self.T_wall_K, + duration=self.duration_s, + n_pts=self.n_pts, + wall_model=self.wall_model, + wall_C_per_area=self.wall_C_per_area_J_m2K, + wall_h_out=self.wall_h_out_W_m2K, + wall_T_inf=self.wall_T_inf_K, + wall_emissivity=self.wall_emissivity, + wall_T_sur=self.wall_T_sur_K, + wall_q_flux=self.wall_q_flux_W_m2, + external_model=self.external_model, + V_ext=self.V_ext_m3, + T_ext=self.T_ext_K, + pump_speed_m3s=self.pump_speed_m3s, + P_ult_Pa=self.P_ult_Pa, + ) + + def to_json(self) -> str: + return json.dumps(asdict(self), ensure_ascii=False, indent=2) + + @classmethod + def from_json(cls, payload: str) -> GuiCaseConfig: + return cls(**json.loads(payload)) + + def save_json(self, path: str | Path) -> None: + Path(path).write_text(self.to_json(), encoding="utf-8") + + @classmethod + def load_json(cls, path: str | Path) -> GuiCaseConfig: + return cls.from_json(Path(path).read_text(encoding="utf-8")) diff --git a/src/venting/gui/main.py b/src/venting/gui/main.py new file mode 100644 index 0000000..3d5acc5 --- /dev/null +++ b/src/venting/gui/main.py @@ -0,0 +1,11 @@ +from __future__ import annotations + + +def main() -> int: + from .app import create_main_window, load_gui_deps + + deps = load_gui_deps() + app = deps.QtWidgets.QApplication.instance() or deps.QtWidgets.QApplication([]) + win = create_main_window() + win.show() + return app.exec() diff --git a/src/venting/gui/state_layout.py b/src/venting/gui/state_layout.py new file mode 100644 index 0000000..4c8e157 --- /dev/null +++ b/src/venting/gui/state_layout.py @@ -0,0 +1,11 @@ +from __future__ import annotations + +from venting.gui.config import GuiCaseConfig +from venting.state_layout import StateLayout, infer_layout_from_modes, split_state + + +def infer_layout(cfg: GuiCaseConfig, n_nodes: int) -> StateLayout: + return infer_layout_from_modes(cfg.thermo, cfg.wall_model, n_nodes) + + +__all__ = ["StateLayout", "infer_layout", "infer_layout_from_modes", "split_state"] diff --git a/src/venting/io.py b/src/venting/io.py new file mode 100644 index 0000000..56edb72 --- /dev/null +++ b/src/venting/io.py @@ -0,0 +1,108 @@ +from __future__ import annotations + +import csv +import json +import platform +import subprocess +import sys +from dataclasses import asdict +from datetime import datetime, timezone +from importlib.metadata import PackageNotFoundError, version +from pathlib import Path + + +def utc_timestamp() -> str: + return datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S") + + +def make_results_dir(case_name: str) -> Path: + out = Path("results") / f"{utc_timestamp()}_{case_name}" + out.mkdir(parents=True, exist_ok=True) + return out + + +def current_git_commit() -> str: + try: + return subprocess.check_output(["git", "rev-parse", "HEAD"], text=True).strip() + except Exception: + return "unknown" + + +def package_version() -> str: + try: + return version("venting") + except PackageNotFoundError: + return "9.0.0" + + +def write_run_json(outdir: Path, params: dict, solver_settings: dict) -> None: + payload = { + "timestamp_utc": datetime.now(timezone.utc).isoformat(), + "git_commit": current_git_commit(), + "python_version": sys.version, + "package_version": package_version(), + "platform": platform.platform(), + "parameters": params, + "solver_settings": solver_settings, + } + (outdir / "run.json").write_text( + json.dumps(payload, ensure_ascii=False, indent=2), + encoding="utf-8", + ) + + +def write_summary_csv(outdir: Path, res) -> None: + rows = [] + for edge, peak in res.peak_diag.items(): + rows.append( + { + "edge": edge, + "max_abs_dP_Pa": res.max_dP.get(edge, float("nan")), + "t_peak_s": peak.get("t_peak", float("nan")), + "r_peak": peak.get("r", float("nan")), + "regime": peak.get("regime", ""), + "peak_type": peak.get("peak_type", ""), + "tau_exit_s": res.tau_exit, + } + ) + with (outdir / "summary.csv").open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter( + fh, + fieldnames=[ + "edge", + "max_abs_dP_Pa", + "t_peak_s", + "r_peak", + "regime", + "peak_type", + "tau_exit_s", + ], + ) + writer.writeheader() + writer.writerows(rows) + + +def dump_meta_json(outdir: Path, filename: str, meta: dict) -> None: + clean_meta = dict(meta) + if "case" in clean_meta and not isinstance(clean_meta["case"], dict): + clean_meta["case"] = asdict(clean_meta["case"]) + (outdir / filename).write_text( + json.dumps(clean_meta, ensure_ascii=False, indent=2), + encoding="utf-8", + ) + + +def write_validity_json(outdir: Path, stem: str, validity_flags: dict) -> Path: + path = outdir / f"{stem}_validity.json" + path.write_text( + json.dumps(validity_flags, ensure_ascii=False, indent=2), encoding="utf-8" + ) + return path + + +def print_validity_summary(validity_flags: dict) -> None: + print("Validity flags:") + for key, payload in validity_flags.items(): + status = payload.get("status", "n/a") + msg = payload.get("message", "") + print(f" - {key}: {status} {msg}".rstrip()) diff --git a/src/venting/montecarlo.py b/src/venting/montecarlo.py new file mode 100644 index 0000000..303eadc --- /dev/null +++ b/src/venting/montecarlo.py @@ -0,0 +1,71 @@ +from __future__ import annotations + +import json +from dataclasses import replace +from pathlib import Path + +import numpy as np + +from .cases import CaseConfig, NetworkConfig +from .diagnostics import summarize_result +from .graph import build_branching_network +from .profiles import Profile +from .solver import solve_case + + +def run_mc( + net_cfg_base: NetworkConfig, + case_cfg: CaseConfig, + profile: Profile, + cd_int_range: tuple[float, float], + cd_exit_range: tuple[float, float], + n_samples: int, + seed: int | None = None, +) -> dict: + """Monte Carlo sweep over Cd uncertainty.""" + rng = np.random.default_rng(seed) + samples = [] + edge_values: dict[str, list[float]] = {} + + for i in range(n_samples): + cd_int = float(rng.uniform(*cd_int_range)) + cd_exit = float(rng.uniform(*cd_exit_range)) + cfg = replace(net_cfg_base, Cd_int=cd_int, Cd_exit=cd_exit) + nodes, edges, bcs = build_branching_network(cfg, profile) + sol = solve_case(nodes, edges, bcs, case_cfg) + res = summarize_result(nodes, edges, bcs, case_cfg, sol) + row = {"sample": i, "Cd_int": cd_int, "Cd_exit": cd_exit} + for edge, value in res.max_dP.items(): + row[f"max_dP_{edge}"] = float(value) + edge_values.setdefault(edge, []).append(float(value)) + samples.append(row) + + summary = {} + for edge, values in edge_values.items(): + arr = np.asarray(values, dtype=float) + summary[edge] = { + "mean": float(np.mean(arr)), + "std": float(np.std(arr)), + "p5": float(np.percentile(arr, 5)), + "p50": float(np.percentile(arr, 50)), + "p95": float(np.percentile(arr, 95)), + } + + return {"samples": samples, "summary": summary} + + +def write_mc_outputs(result: dict, outdir: Path) -> None: + outdir.mkdir(parents=True, exist_ok=True) + sample_rows = result["samples"] + if sample_rows: + cols = list(sample_rows[0].keys()) + import csv + + with (outdir / "mc_results.csv").open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=cols) + writer.writeheader() + writer.writerows(sample_rows) + (outdir / "mc_summary.json").write_text( + json.dumps(result["summary"], ensure_ascii=False, indent=2), + encoding="utf-8", + ) diff --git a/src/venting/plotting.py b/src/venting/plotting.py new file mode 100644 index 0000000..593cbc8 --- /dev/null +++ b/src/venting/plotting.py @@ -0,0 +1,32 @@ +from pathlib import Path + +HAS_MPL = False +try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + HAS_MPL = True +except Exception: + HAS_MPL = False + + +def plot_basic(outdir: Path, name: str, res, node_idx: int = 0) -> None: + if not HAS_MPL: + return + outdir.mkdir(parents=True, exist_ok=True) + + fig, ax = plt.subplots(figsize=(11, 6)) + ax.plot( + res.t, + (res.P[node_idx] - res.P_ext) / 1e3, + lw=2, + label=f"ΔP(node{node_idx}→ext) [kPa]", + ) + ax.set(xlabel="t, s", ylabel="ΔP, kPa", title=f"{name}: ΔP vs time") + ax.grid(True, alpha=0.3) + ax.legend() + fig.tight_layout() + fig.savefig(outdir / f"{name}_dP.png", dpi=200) + plt.close(fig) diff --git a/src/venting/presets.py b/src/venting/presets.py new file mode 100644 index 0000000..2ecc020 --- /dev/null +++ b/src/venting/presets.py @@ -0,0 +1,35 @@ +from __future__ import annotations + +from dataclasses import asdict, dataclass + +from .geometry import mm2_to_m2, mm3_to_m3 + + +@dataclass(frozen=True) +class PanelPresetV9: + A_cell_mm2: float + A_vest_mm2: float + h_mm: float + V_cell: float + V_vest: float + A_wall_cell: float + A_wall_vest: float + + def to_dict(self) -> dict: + return asdict(self) + + +def get_default_panel_preset_v9() -> PanelPresetV9: + """Return the historic CLI default panel geometry (SI outputs).""" + a_cell_mm2 = 4708.2806538 + a_vest_mm2 = 5199.9595503 + h_mm = 27.9430913 + return PanelPresetV9( + A_cell_mm2=a_cell_mm2, + A_vest_mm2=a_vest_mm2, + h_mm=h_mm, + V_cell=mm3_to_m3(a_cell_mm2 * h_mm), + V_vest=mm3_to_m3(a_vest_mm2 * h_mm), + A_wall_cell=mm2_to_m2(2 * a_cell_mm2 + 313.0 * h_mm), + A_wall_vest=mm2_to_m2(2 * a_vest_mm2 + 350.0 * h_mm), + ) diff --git a/src/venting/profiles.py b/src/venting/profiles.py new file mode 100644 index 0000000..1be20f7 --- /dev/null +++ b/src/venting/profiles.py @@ -0,0 +1,100 @@ +import math +from collections.abc import Callable +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + + +@dataclass(frozen=True) +class Profile: + name: str + P: Callable[[float], float] + events: tuple[tuple[float, str], ...] + + def P_array(self, t: np.ndarray) -> np.ndarray: + return np.array([self.P(float(tt)) for tt in t], dtype=float) + + def classify_peak(self, t_peak: float, t_end: float, tol_s: float = 0.5) -> str: + for te, label in self.events: + if abs(t_peak - te) <= tol_s: + return f"boundary({label})" + if t_peak <= tol_s: + return "boundary(sim_start)" + if t_peak >= (t_end - tol_s): + return "boundary(sim_end)" + return "internal" + + +def make_profile_linear(p0: float, rate_mmhg_per_s: float) -> Profile: + rate = rate_mmhg_per_s * 133.322 + t_zero = p0 / rate + + def p_fn(t: float) -> float: + return max(p0 - rate * t, 0.0) + + return Profile("linear", p_fn, events=((t_zero, "P_ext=0"),)) + + +def make_profile_step(p0: float, step_time_s: float) -> Profile: + def p_fn(t: float) -> float: + return 0.0 if t >= step_time_s else p0 + + return Profile("step", p_fn, events=((step_time_s, "step_to_vacuum"),)) + + +def make_profile_exponential( + p0: float, rate0_mmhg_per_s: float, p_floor: float = 10.0 +) -> Profile: + 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 Profile( + "barometric_exp", + p_fn, + events=((0.0, "start_max_slope"), (t_floor, f"P_ext={p_floor:.0f}Pa")), + ) + + +def make_profile_from_table( + name: str, + table_path: Path, + pressure_unit: str = "Pa", +) -> Profile: + """Load profile table with CSV columns: t_s,P_Pa. + + pressure_unit: "Pa" | "mmHg" + """ + arr = np.loadtxt(str(table_path), delimiter=",") + if arr.ndim != 2 or arr.shape[1] < 2: + raise ValueError("Profile table must be CSV with columns: t_s, P_Pa") + t_tab = np.array(arr[:, 0], dtype=float) + p_tab = np.array(arr[:, 1], dtype=float) + if not np.all(np.diff(t_tab) > 0): + raise ValueError("Profile table time must be strictly increasing") + + if pressure_unit.lower() == "mmhg": + p_tab = p_tab * 133.322 + elif pressure_unit.lower() != "pa": + raise ValueError("pressure_unit must be 'Pa' or 'mmHg'") + + p_max = float(np.max(p_tab)) + if pressure_unit.lower() == "pa" and 200.0 <= p_max <= 2000.0: + raise ValueError( + "Profile table pressure looks like mmHg values provided as Pa. " + "Use pressure_unit='mmHg' or convert CSV to Pa." + ) + + def p_fn(t: float) -> float: + if t <= t_tab[0]: + return float(p_tab[0]) + if t >= t_tab[-1]: + return float(p_tab[-1]) + return float(np.interp(t, t_tab, p_tab)) + + events = tuple((float(tt), f"bp{i}") for i, tt in enumerate(t_tab)) + return Profile(name, p_fn, events=events) diff --git a/src/venting/run.py b/src/venting/run.py new file mode 100644 index 0000000..7e80c9b --- /dev/null +++ b/src/venting/run.py @@ -0,0 +1,57 @@ +from __future__ import annotations + +from pathlib import Path + +import numpy as np + +from .cases import CaseConfig, NetworkConfig +from .diagnostics import summarize_result +from .graph import build_branching_network +from .io import ( + dump_meta_json, + make_results_dir, + print_validity_summary, + write_run_json, + write_summary_csv, + write_validity_json, +) +from .profiles import Profile +from .solver import solve_case + + +def run_case(net_cfg: NetworkConfig, profile: Profile, case_cfg: CaseConfig): + nodes, edges, bcs = build_branching_network(net_cfg, profile) + sol = solve_case(nodes, edges, bcs, case_cfg) + return summarize_result(nodes, edges, bcs, case_cfg, sol) + + +def export_case_artifacts( + outdir: Path, + stem: str, + res, + run_params: dict, + solver_settings: dict | None = None, +) -> None: + solver_payload = solver_settings or { + "method": "Radau", + "rtol": "1e-7|1e-6", + "atol": "1e-10|1e-8", + } + np.savez_compressed( + outdir / f"{stem}.npz", + t=res.t, + m=res.m, + T=res.T, + P=res.P, + P_ext=res.P_ext, + tau_exit=res.tau_exit, + ) + dump_meta_json(outdir, f"{stem}_meta.json", res.meta) + write_validity_json(outdir, stem, res.meta.get("validity_flags", {})) + print_validity_summary(res.meta.get("validity_flags", {})) + write_summary_csv(outdir, res) + write_run_json(outdir, params=run_params, solver_settings=solver_payload) + + +def make_case_output_dir(case_name: str): + return make_results_dir(case_name) diff --git a/src/venting/solver.py b/src/venting/solver.py new file mode 100644 index 0000000..5ac87c1 --- /dev/null +++ b/src/venting/solver.py @@ -0,0 +1,658 @@ +import math + +import numpy as np +from scipy.integrate import solve_ivp + +from .cases import CaseConfig +from .constants import C_CHOKED, C_P, C_V, GAMMA, M_SAFE, P0, R_GAS, T0, T_SAFE +from .flow import ( + mdot_fanno_tube, + mdot_orifice_pos_props, + mdot_short_tube_pos, + mdot_slot_pos, +) +from .graph import ( + EXT_NODE, + ExternalBC, + GasNode, + OrificeEdge, + ShortTubeEdge, + SlotChannelEdge, +) +from .thermo import cp_air, cv_air, gamma_air, h_air, u_air + +SIGMA_SB = 5.670374419e-8 + + +def _mdot_short_edge( + edge: ShortTubeEdge, + p_up: float, + t_up: float, + p_dn: float, + cd0: float, + gamma: float = GAMMA, +): + if edge.fanno: + return mdot_fanno_tube( + p_up, + t_up, + p_dn, + cd0, + edge.A_total, + edge.D, + edge.L, + edge.eps, + edge.K_in, + edge.K_out, + gamma=gamma, + ) + return mdot_short_tube_pos( + p_up, + t_up, + p_dn, + cd0, + edge.A_total, + edge.D, + edge.L, + edge.eps, + edge.K_in, + edge.K_out, + gamma=gamma, + ) + + +def _property_model(thermo: str): + if thermo == "variable": + return cp_air, cv_air, gamma_air, h_air, u_air + if thermo == "intermediate": + return ( + lambda T: C_P, + lambda T: C_V, + lambda T: C_P / C_V, + lambda T: C_P * T, + lambda T: C_V * T, + ) + raise ValueError("Property model valid for 'intermediate' or 'variable'") + + +def build_rhs( + nodes: list[GasNode], edges: list, bcs: list[ExternalBC], case: CaseConfig +): + N = len(nodes) + V = np.array([n.V for n in nodes], dtype=float) + A_w = np.array([n.A_wall for n in nodes], dtype=float) + bc_map = {bc.node: bc for bc in bcs} + + def p_from_mt(m: np.ndarray, T: np.ndarray) -> np.ndarray: + return m * R_GAS * T / V + + if case.thermo == "isothermal": + n_vars = N + + def rhs(t: float, y: np.ndarray) -> np.ndarray: + m = y[:N] + T = np.full(N, T0) + P = p_from_mt(np.maximum(m, 0.0), T) + dm = np.zeros(N) + for e in edges: + if isinstance(e, (OrificeEdge, ShortTubeEdge)): + a, b = e.a, e.b + cd0 = e.Cd_model() + if b == EXT_NODE: + bc = bc_map[a] + pext = bc.profile.P(float(t)) + if P[a] >= pext: + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + P[a], T[a], pext, cd0, e.A_total + ) + else: + md = _mdot_short_edge(e, P[a], T[a], pext, cd0) + dm[a] -= md + else: + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + pext, bc.T_ext, P[a], cd0, e.A_total + ) + else: + md = _mdot_short_edge(e, pext, bc.T_ext, P[a], cd0) + dm[a] += md + else: + Pa, Pb = P[a], P[b] + if Pa >= Pb: + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + Pa, T[a], Pb, cd0, e.A_total + ) + else: + md = _mdot_short_edge(e, Pa, T[a], Pb, cd0) + dm[a] -= md + dm[b] += md + else: + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + Pb, T[b], Pa, cd0, e.A_total + ) + else: + md = _mdot_short_edge(e, Pb, T[b], Pa, cd0) + dm[b] -= md + dm[a] += md + elif isinstance(e, SlotChannelEdge): + a, b = e.a, e.b + Pa, Pb = P[a], P[b] + if Pa >= Pb: + md = mdot_slot_pos(Pa, T[a], Pb, e.w, e.delta, e.L) + dm[a] -= md + dm[b] += md + else: + md = mdot_slot_pos(Pb, T[b], Pa, e.w, e.delta, e.L) + dm[b] -= md + dm[a] += md + return dm + + return rhs, n_vars + + if case.thermo not in {"intermediate", "variable"}: + raise ValueError( + "case.thermo must be 'isothermal', 'intermediate', or 'variable'" + ) + + cp_fn, cv_fn, gamma_fn, h_fn, u_fn = _property_model(case.thermo) + + has_lumped_wall = case.wall_model == "lumped" + n_vars = 2 * N + (N if has_lumped_wall else 0) + + def rhs(t: float, y: np.ndarray) -> np.ndarray: + m = y[:N] + T = y[N : 2 * N] + Tw = y[2 * N : 3 * N] if has_lumped_wall else None + T_eff = np.maximum(T, T_SAFE) + + P = p_from_mt(np.maximum(m, 0.0), T_eff) + dm = np.zeros(N) + dE = np.zeros(N) + + for e in edges: + if isinstance(e, (OrificeEdge, ShortTubeEdge)): + a, b = e.a, e.b + cd0 = e.Cd_model() + if b == EXT_NODE: + bc = bc_map[a] + pext = bc.profile.P(float(t)) + if P[a] >= pext: + gamma_up = gamma_fn(float(T_eff[a])) + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + P[a], T_eff[a], pext, cd0, e.A_total, gamma=gamma_up + ) + else: + md = _mdot_short_edge( + e, P[a], T_eff[a], pext, cd0, gamma=gamma_up + ) + dm[a] -= md + dE[a] += -md * h_fn(float(T_eff[a])) + else: + gamma_up = gamma_fn(float(max(bc.T_ext, T_SAFE))) + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + pext, + max(bc.T_ext, T_SAFE), + P[a], + cd0, + e.A_total, + gamma=gamma_up, + ) + else: + md = _mdot_short_edge( + e, + pext, + max(bc.T_ext, T_SAFE), + P[a], + cd0, + gamma=gamma_up, + ) + dm[a] += md + dE[a] += md * h_fn(float(max(bc.T_ext, T_SAFE))) + else: + Pa, Pb = P[a], P[b] + if Pa >= Pb: + gamma_up = gamma_fn(float(T_eff[a])) + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + Pa, T_eff[a], Pb, cd0, e.A_total, gamma=gamma_up + ) + else: + md = _mdot_short_edge( + e, Pa, T_eff[a], Pb, cd0, gamma=gamma_up + ) + dm[a] -= md + dm[b] += md + dE[a] += -md * h_fn(float(T_eff[a])) + dE[b] += md * h_fn(float(T_eff[a])) + else: + gamma_up = gamma_fn(float(T_eff[b])) + if isinstance(e, OrificeEdge): + md = mdot_orifice_pos_props( + Pb, T_eff[b], Pa, cd0, e.A_total, gamma=gamma_up + ) + else: + md = _mdot_short_edge( + e, Pb, T_eff[b], Pa, cd0, gamma=gamma_up + ) + dm[b] -= md + dm[a] += md + dE[b] += -md * h_fn(float(T_eff[b])) + dE[a] += md * h_fn(float(T_eff[b])) + elif isinstance(e, SlotChannelEdge): + a, b = e.a, e.b + Pa, Pb = P[a], P[b] + if Pa >= Pb: + md = mdot_slot_pos(Pa, T_eff[a], Pb, e.w, e.delta, e.L) + dm[a] -= md + dm[b] += md + dE[a] += -md * h_fn(float(T_eff[a])) + dE[b] += md * h_fn(float(T_eff[a])) + else: + md = mdot_slot_pos(Pb, T_eff[b], Pa, e.w, e.delta, e.L) + dm[b] -= md + dm[a] += md + dE[b] += -md * h_fn(float(T_eff[b])) + dE[a] += md * h_fn(float(T_eff[b])) + + if has_lumped_wall: + Qdot = case.h_conv * A_w * (Tw - T) + else: + Qdot = case.h_conv * A_w * (case.T_wall - T) + dE += Qdot + + dT = np.zeros(N) + for i in range(N): + m_eff = max(float(m[i]), 0.0) + if m_eff > M_SAFE: + u_i = u_fn(float(T_eff[i])) + cv_i = max(cv_fn(float(T_eff[i])), 1e-9) + dT[i] = (dE[i] - u_i * dm[i]) / (m_eff * cv_i) + + if has_lumped_wall: + dTw = np.zeros(N) + for i in range(N): + A = A_w[i] + Cw = case.wall_C_per_area * A + if Cw <= 0.0: + continue + q_in = case.h_conv * A * (T[i] - Tw[i]) + q_out = case.wall_h_out * A * (Tw[i] - case.wall_T_inf) + q_rad = ( + case.wall_emissivity + * SIGMA_SB + * A + * (Tw[i] ** 4 - case.wall_T_sur**4) + ) + q_src = case.wall_q_flux * A + dTw[i] = (q_in - q_out - q_rad + q_src) / Cw + return np.concatenate([dm, dT, dTw]) + + return np.concatenate([dm, dT]) + + return rhs, n_vars + + +def solve_case( + nodes: list[GasNode], edges: list, bcs: list[ExternalBC], case: CaseConfig +): + nodes_local = list(nodes) + edges_local = list(edges) + bcs_local = list(bcs) + ext_idx = None + + if case.external_model == "dynamic_pump": + ext_idx = len(nodes_local) + nodes_local.append(GasNode("external", case.V_ext, 0.0)) + rewired = [] + for e in edges_local: + if isinstance(e, OrificeEdge) and e.b == EXT_NODE: + rewired.append( + OrificeEdge(e.a, ext_idx, e.A_total, e.Cd_model, label=e.label) + ) + elif isinstance(e, ShortTubeEdge) and e.b == EXT_NODE: + rewired.append( + ShortTubeEdge( + e.a, + ext_idx, + e.A_total, + e.D, + e.L, + e.eps, + e.K_in, + e.K_out, + e.Cd_model, + label=e.label, + fanno=e.fanno, + ) + ) + else: + rewired.append(e) + edges_local = rewired + bcs_local = [] + + N = len(nodes_local) + V = np.array([n.V for n in nodes_local], dtype=float) + rhs_core, n_vars = build_rhs(nodes_local, edges_local, bcs_local, case) + + m0 = P0 * V / (R_GAS * T0) + if ext_idx is not None: + m0[ext_idx] = max(case.P_ult_Pa, 10.0) * V[ext_idx] / (R_GAS * case.T_ext) + + if case.thermo == "isothermal": + y0 = m0.copy() + else: + T_init = np.full(N, T0) + if ext_idx is not None: + T_init[ext_idx] = case.T_ext + y0 = np.concatenate([m0, T_init]) + if case.wall_model == "lumped": + y0 = np.concatenate([y0, np.full(N, case.T_wall)]) + + t_eval = np.linspace(0.0, case.duration, int(case.n_pts)) + + profile = bcs_local[0].profile if bcs_local else None + + def p_from_state(y: np.ndarray) -> np.ndarray: + m = y[:N] + if case.thermo == "isothermal": + T = np.full(N, T0) + if ext_idx is not None: + T[ext_idx] = case.T_ext + else: + T = np.maximum(y[N : 2 * N], T_SAFE) + if ext_idx is not None: + T[ext_idx] = case.T_ext + return np.maximum(m, 0.0) * R_GAS * T / V + + def rhs(t: float, y: np.ndarray) -> np.ndarray: + dy = rhs_core(t, y) + if ext_idx is not None: + P = p_from_state(y) + mdot_pump = ( + case.pump_speed_m3s + * max(P[ext_idx] - case.P_ult_Pa, 0.0) + / (R_GAS * max(case.T_ext, T_SAFE)) + ) + dy[ext_idx] -= mdot_pump + if case.thermo != "isothermal": + dy[N + ext_idx] = 0.0 + return dy + + events = [] + for i in range(N): + + def ev_mi(t, y, ii=i): + return y[ii] + 1e-15 + + ev_mi.terminal = True + ev_mi.direction = -1 + events.append(ev_mi) + + def ev_equil_rms(t, y): + P = p_from_state(y) + if ext_idx is not None: + Pe = P[ext_idx] + Pcmp = np.delete(P, ext_idx) + else: + Pe = profile.P(float(t)) + Pcmp = P + return float(np.sqrt(np.mean((Pcmp - Pe) ** 2))) - case.p_rms_tol + + ev_equil_rms.terminal = True + ev_equil_rms.direction = -1 + events.append(ev_equil_rms) + + def ev_p_stop(t, y): + P = p_from_state(y) + if ext_idx is not None: + Pe = P[ext_idx] + Pcmp = np.delete(P, ext_idx) + else: + Pe = profile.P(float(t)) + Pcmp = P + return max(np.max(Pcmp) - case.p_stop, Pe - 10.0) + + ev_p_stop.terminal = True + ev_p_stop.direction = -1 + events.append(ev_p_stop) + + A_max = 0.0 + Cd_max = 0.0 + for e in edges_local: + if isinstance(e, (OrificeEdge, ShortTubeEdge)): + A_max = max(A_max, e.A_total) + Cd_max = max(Cd_max, e.Cd_model()) + V_min = float(np.min(V)) + mdot_ch = ( + Cd_max * A_max * C_CHOKED * P0 / math.sqrt(R_GAS * T0) if A_max > 0 else 0.0 + ) + tau_min = (V_min * P0) / (R_GAS * T0 * mdot_ch) if mdot_ch > 0 else 1.0 + max_step = max(min(case.duration / 2000.0, tau_min / 10.0), 1e-4) + + sol = solve_ivp( + rhs, + (0.0, case.duration), + y0, + method="Radau", + t_eval=t_eval, + 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, + ) + sol.ext_idx = ext_idx + sol.node_count = len(nodes_local) + sol.nodes_local = nodes_local + sol.edges_local = edges_local + sol.bcs_local = bcs_local + return sol + + +def solve_case_stream( + nodes: list[GasNode], + edges: list, + bcs: list[ExternalBC], + case: CaseConfig, + callback=None, + n_chunks: int | None = None, + dt_chunk_s: float = 2.0, + should_stop=None, + stop_check=None, +): + """True streaming solve via piecewise integration segments.""" + if should_stop is None and stop_check is not None: + should_stop = stop_check + if callback is None and should_stop is None: + return solve_case(nodes, edges, bcs, case) + nodes_local = list(nodes) + edges_local = list(edges) + bcs_local = list(bcs) + ext_idx = None + + if case.external_model == "dynamic_pump": + ext_idx = len(nodes_local) + nodes_local.append(GasNode("external", case.V_ext, 0.0)) + rewired = [] + for e in edges_local: + if isinstance(e, OrificeEdge) and e.b == EXT_NODE: + rewired.append( + OrificeEdge(e.a, ext_idx, e.A_total, e.Cd_model, label=e.label) + ) + elif isinstance(e, ShortTubeEdge) and e.b == EXT_NODE: + rewired.append( + ShortTubeEdge( + e.a, + ext_idx, + e.A_total, + e.D, + e.L, + e.eps, + e.K_in, + e.K_out, + e.Cd_model, + label=e.label, + fanno=e.fanno, + ) + ) + else: + rewired.append(e) + edges_local = rewired + bcs_local = [] + + N = len(nodes_local) + V = np.array([n.V for n in nodes_local], dtype=float) + rhs_core, _ = build_rhs(nodes_local, edges_local, bcs_local, case) + + m0 = P0 * V / (R_GAS * T0) + if ext_idx is not None: + m0[ext_idx] = max(case.P_ult_Pa, 10.0) * V[ext_idx] / (R_GAS * case.T_ext) + + if case.thermo == "isothermal": + y = m0.copy() + else: + T_init = np.full(N, T0) + if ext_idx is not None: + T_init[ext_idx] = case.T_ext + y = np.concatenate([m0, T_init]) + if case.wall_model == "lumped": + y = np.concatenate([y, np.full(N, case.T_wall)]) + + t_eval = np.linspace(0.0, case.duration, int(case.n_pts)) + + def p_from_state(y_state: np.ndarray) -> np.ndarray: + m = y_state[:N] + if case.thermo == "isothermal": + T = np.full(N, T0) + if ext_idx is not None: + T[ext_idx] = case.T_ext + else: + T = np.maximum(y_state[N : 2 * N], T_SAFE) + if ext_idx is not None: + T[ext_idx] = case.T_ext + return np.maximum(m, 0.0) * R_GAS * T / V + + def rhs(t: float, y_state: np.ndarray) -> np.ndarray: + dy = rhs_core(t, y_state) + if ext_idx is not None: + P = p_from_state(y_state) + mdot_pump = ( + case.pump_speed_m3s + * max(P[ext_idx] - case.P_ult_Pa, 0.0) + / (R_GAS * max(case.T_ext, T_SAFE)) + ) + dy[ext_idx] -= mdot_pump + if case.thermo != "isothermal": + dy[N + ext_idx] = 0.0 + return dy + + A_max = 0.0 + Cd_max = 0.0 + for e in edges_local: + if isinstance(e, (OrificeEdge, ShortTubeEdge)): + A_max = max(A_max, e.A_total) + Cd_max = max(Cd_max, e.Cd_model()) + V_min = float(np.min(V)) + mdot_ch = ( + Cd_max * A_max * C_CHOKED * P0 / math.sqrt(R_GAS * T0) if A_max > 0 else 0.0 + ) + tau_min = (V_min * P0) / (R_GAS * T0 * mdot_ch) if mdot_ch > 0 else 1.0 + max_step = max(min(case.duration / 2000.0, tau_min / 10.0), 1e-4) + + if n_chunks is not None: + effective_n_chunks = max(int(n_chunks), 1) + else: + effective_n_chunks = max( + int(math.ceil(case.duration / max(dt_chunk_s, 1e-6))), 1 + ) + bounds = np.linspace(0.0, case.duration, effective_n_chunks + 1) + t_out = [] + y_out = [] + success = True + message = "stream completed" + + for i in range(effective_n_chunks): + if should_stop is not None and should_stop(): + success = False + message = "stream cancelled" + break + t0 = bounds[i] + t1 = bounds[i + 1] + mask = (t_eval >= t0) & (t_eval <= t1) + t_seg = t_eval[mask] + if t_seg.size == 0: + t_seg = np.array([t0, t1]) + sol_seg = solve_ivp( + rhs, + (float(t0), float(t1)), + y, + method="Radau", + t_eval=t_seg, + rtol=1e-7 if case.thermo == "isothermal" else 1e-6, + atol=1e-10 if case.thermo == "isothermal" else 1e-8, + max_step=max_step, + ) + if not sol_seg.success: + success = False + message = str(sol_seg.message) + break + + seg_t = sol_seg.t + seg_y = sol_seg.y + if i > 0 and seg_t.size > 0: + seg_t = seg_t[1:] + seg_y = seg_y[:, 1:] + if seg_t.size > 0: + t_out.append(seg_t) + y_out.append(seg_y) + y = sol_seg.y[:, -1] + + if callback is not None and t_out: + t_cat = np.concatenate(t_out) + y_cat = np.concatenate(y_out, axis=1) + callback( + { + "t": t_cat, + "y": y_cat, + "progress": float((i + 1) / effective_n_chunks), + "done": bool(i == effective_n_chunks - 1), + "node_count": N, + "n_nodes": N, + "thermo": case.thermo, + "wall_model": case.wall_model, + "ext_idx": ext_idx, + } + ) + + if t_out: + t_final = np.concatenate(t_out) + y_final = np.concatenate(y_out, axis=1) + else: + t_final = np.array([0.0]) + y_final = y[:, None] + + if success and (t_final.size != t_eval.size or not np.allclose(t_final, t_eval)): + y_interp = np.empty((y_final.shape[0], t_eval.size), dtype=float) + for i_row in range(y_final.shape[0]): + y_interp[i_row] = np.interp(t_eval, t_final, y_final[i_row]) + t_final = t_eval + y_final = y_interp + + class _Sol: + pass + + sol = _Sol() + sol.t = t_final + sol.y = y_final + sol.success = success + sol.message = message + sol.ext_idx = ext_idx + sol.node_count = len(nodes_local) + sol.nodes_local = nodes_local + sol.edges_local = edges_local + sol.bcs_local = bcs_local + return sol diff --git a/src/venting/state_layout.py b/src/venting/state_layout.py new file mode 100644 index 0000000..3a45e7d --- /dev/null +++ b/src/venting/state_layout.py @@ -0,0 +1,34 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass(frozen=True) +class StateLayout: + n_nodes: int + m_slice: slice + t_slice: slice | None + tw_slice: slice | None + + +def split_state(y, node_count: int, thermo: str, wall_model: str): + layout = infer_layout_from_modes(thermo, wall_model, node_count) + m = y[layout.m_slice, :] + t = y[layout.t_slice, :] if layout.t_slice is not None else None + tw = y[layout.tw_slice, :] if layout.tw_slice is not None else None + return m, t, tw + + +def infer_layout_from_modes(thermo: str, wall_model: str, n_nodes: int) -> StateLayout: + if n_nodes <= 0: + raise ValueError("n_nodes must be > 0") + m_slice = slice(0, n_nodes) + if thermo == "isothermal": + return StateLayout( + n_nodes=n_nodes, m_slice=m_slice, t_slice=None, tw_slice=None + ) + t_slice = slice(n_nodes, 2 * n_nodes) + tw_slice = slice(2 * n_nodes, 3 * n_nodes) if wall_model == "lumped" else None + return StateLayout( + n_nodes=n_nodes, m_slice=m_slice, t_slice=t_slice, tw_slice=tw_slice + ) diff --git a/src/venting/thermo.py b/src/venting/thermo.py new file mode 100644 index 0000000..440a3d7 --- /dev/null +++ b/src/venting/thermo.py @@ -0,0 +1,69 @@ +"""Temperature-dependent thermodynamic properties for dry air. + +NASA-7 polynomial fit (200-1000 K) for dry-air mixture: +0.78084 N2 + 0.20946 O2 + 0.00934 Ar. +Source: McBride, Gordon, Reno, NASA TM-4513 (1993). +""" + +from __future__ import annotations + +import math + +from .constants import R_GAS + +# NASA-7 coefficients for cp/R, valid for 200-1000 K. +_A1 = 3.53562097 +_A2 = -4.15142542e-4 +_A3 = 1.03930790e-6 +_A4 = 0.0 +_A5 = 0.0 + +T_FIT_LOW = 200.0 +T_FIT_HIGH = 1000.0 +_T_REF = 298.15 + + +def _cp_over_r(t: float) -> float: + return _A1 + _A2 * t + _A3 * t**2 + _A4 * t**3 + _A5 * t**4 + + +def cp_air(T: float) -> float: + """Specific heat cp for dry air [J/(kg·K)] using NASA-7 fit.""" + t = min(max(float(T), 1.0), T_FIT_HIGH) + return max(_cp_over_r(t) * R_GAS, 700.0) + + +def cv_air(T: float) -> float: + return cp_air(T) - R_GAS + + +def gamma_air(T: float) -> float: + cp = cp_air(T) + cv = max(cv_air(T), 1e-12) + return cp / cv + + +def _h_over_r(t: float) -> float: + return ( + _A1 * t + + 0.5 * _A2 * t**2 + + (_A3 / 3.0) * t**3 + + 0.25 * _A4 * t**4 + + 0.2 * _A5 * t**5 + ) + + +def h_air(T: float) -> float: + """Specific enthalpy offset from 298.15 K [J/kg].""" + t = min(max(float(T), 1.0), T_FIT_HIGH) + return (_h_over_r(t) - _h_over_r(_T_REF)) * R_GAS + + +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) + + +def speed_of_sound(T: float) -> float: + t = max(float(T), 1.0) + return math.sqrt(gamma_air(t) * R_GAS * t) diff --git a/src/venting/validity.py b/src/venting/validity.py new file mode 100644 index 0000000..212b3fc --- /dev/null +++ b/src/venting/validity.py @@ -0,0 +1,341 @@ +from __future__ import annotations + +import math + +import numpy as np + +from .constants import D_MOL_AIR, GAMMA, K_BOLTZMANN, R_GAS, T_SAFE +from .flow import ( + fanno_choked_state, + friction_factor, + mdot_short_tube_pos, + mu_air_sutherland, +) +from .graph import EXT_NODE, OrificeEdge, ShortTubeEdge, SlotChannelEdge +from .thermo import T_FIT_HIGH, T_FIT_LOW + + +def _acoustic_flag( + P: np.ndarray, T: np.ndarray, t: np.ndarray, l_char_m: float +) -> dict: + t_mean = float(np.mean(np.maximum(T, T_SAFE))) + c = math.sqrt(GAMMA * R_GAS * t_mean) + t_ac = l_char_m / max(c, 1e-12) + + dP_dt = np.gradient(P, t, axis=1) + with np.errstate(divide="ignore", invalid="ignore"): + t_dyn = np.abs(P / dP_dt) + t_dyn = t_dyn[np.isfinite(t_dyn) & (t_dyn > 0)] + t_dyn_min = float(np.min(t_dyn)) if t_dyn.size else float("inf") + + ratio = t_dyn_min / t_ac if t_ac > 0 else float("inf") + status = "ok" if ratio >= 20.0 else "warning" + + return { + "status": status, + "t_acoustic_s": t_ac, + "t_dynamic_min_s": t_dyn_min, + "ratio_t_dyn_to_t_ac": ratio, + "message": "0D pressure uniformity is stronger when ratio >> 1", + } + + +def _thermo_range_flag(T: np.ndarray) -> dict: + t_min = float(np.min(T)) + t_max = float(np.max(T)) + in_range = (t_min >= T_FIT_LOW) and (t_max <= T_FIT_HIGH) + return { + "status": "ok" if in_range else "warning", + "T_min_K": t_min, + "T_max_K": t_max, + "fit_range": f"{T_FIT_LOW}-{T_FIT_HIGH} K", + "source": "NASA TM-4513 (1993), mixture-averaged", + "message": ( + "Temperature remained inside NASA-7 fit range" + if in_range + else "Temperature left NASA-7 fit range; thermo accuracy may degrade" + ), + } + + +def _slot_laminar_flag( + edges: list, + P: np.ndarray, + T: np.ndarray, + node_index: dict[int, int], +) -> dict: + re_max = 0.0 + for e in edges: + if not isinstance(e, SlotChannelEdge): + continue + a = node_index[e.a] + b = node_index[e.b] + for k in range(P.shape[1]): + up = a if P[a, k] >= P[b, k] else b + p_up = max(float(P[up, k]), 0.0) + t_up = max(float(T[up, k]), T_SAFE) + rho = p_up / (R_GAS * t_up) + mu = mu_air_sutherland(t_up) + u = ( + (e.delta**2) + * abs(float(P[a, k] - P[b, k])) + / (12.0 * mu * max(e.L, 1e-12)) + ) + d_h = 2.0 * e.delta + re = rho * u * d_h / max(mu, 1e-20) + re_max = max(re_max, re) + + if re_max == 0.0: + return { + "status": "ok", + "Re_max": 0.0, + "message": "No slot/channel edges in this case", + } + status = "ok" if re_max < 1000.0 else "warning" + return { + "status": status, + "Re_max": re_max, + "message": "Slot model assumes laminar flow; warning if Re approaches turbulent regime", + } + + +def _short_tube_flag( + edges: list, + P: np.ndarray, + T: np.ndarray, + P_ext: np.ndarray, + node_index: dict[int, int], +) -> dict: + short_tubes = [e for e in edges if isinstance(e, ShortTubeEdge)] + if not short_tubes: + return { + "status": "ok", + "message": "No short-tube edges", + "Re_max": 0.0, + "Mach_max": 0.0, + "K_tot_max": 0.0, + "Cd_eff_min": float("nan"), + "L_over_D": 0.0, + "frac_fric": 0.0, + "Fanno_choke_fraction": 0.0, + } + + re_max = 0.0 + mach_max = 0.0 + k_tot_max = 0.0 + cd_eff_min = float("inf") + l_over_d_worst = 0.0 + frac_fric_worst = 0.0 + fanno_steps = 0 + fanno_choked_steps = 0 + + for e in short_tubes: + a = node_index[e.a] + b = e.b if e.b >= 0 else EXT_NODE + cd0 = e.Cd_model() + for k in range(P.shape[1]): + Pa = float(P[a, k]) + if b == EXT_NODE: + Pb = float(P_ext[k]) + up_idx = a + p_up, p_dn = (Pa, Pb) if Pa >= Pb else (Pb, Pa) + t_up = float(T[a, k]) if Pa >= Pb else max(float(T[a, k]), T_SAFE) + else: + bb = node_index[b] + Pb = float(P[bb, k]) + if Pa >= Pb: + up_idx = a + p_up, p_dn = Pa, Pb + else: + up_idx = bb + p_up, p_dn = Pb, Pa + 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, + ) + 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) + mach = u / max(a_s, 1e-18) + mu = mu_air_sutherland(t_up) + Re = rho * u * e.D / max(mu, 1e-18) + f_D = friction_factor(Re, e.eps / max(e.D, 1e-12)) + k_fric = f_D * (e.L / max(e.D, 1e-12)) + k_tot = e.K_in + e.K_out + k_fric + cd_eff = 1.0 / math.sqrt(max(cd0 ** (-2.0) + k_tot, 1e-18)) + + if e.fanno: + fanno_steps += 1 + choked, _ = fanno_choked_state( + p_up, + t_up, + p_dn, + cd0, + e.A_total, + e.D, + e.L, + e.eps, + e.K_in, + e.K_out, + ) + if choked: + fanno_choked_steps += 1 + + re_max = max(re_max, Re) + mach_max = max(mach_max, mach) + if k_tot > k_tot_max: + k_tot_max = k_tot + l_over_d_worst = e.L / max(e.D, 1e-12) + frac_fric_worst = k_fric / max(k_tot, 1e-20) + cd_eff_min = min(cd_eff_min, cd_eff) + + if mach_max < 0.3: + status = "ok" + msg = "Short-tube Mach within comfort zone" + elif mach_max < 0.6: + status = "warning" + msg = "Mach out of comfort zone for quasi-1D loss model" + else: + status = "warning" + msg = "Mach high; consider Fanno/friction-choking model in v10+" + + return { + "status": status, + "message": msg, + "Re_max": re_max, + "Mach_max": mach_max, + "K_tot_max": k_tot_max, + "Cd_eff_min": cd_eff_min if np.isfinite(cd_eff_min) else float("nan"), + "L_over_D": l_over_d_worst, + "frac_fric": frac_fric_worst, + "Fanno_choke_fraction": ( + float(fanno_choked_steps / fanno_steps) if fanno_steps else 0.0 + ), + } + + +def _knudsen_flag( + edges: list, + P: np.ndarray, + T: np.ndarray, + P_ext: np.ndarray, + node_index: dict[int, int], + t: np.ndarray, +) -> dict: + kn_max = 0.0 + edge_worst = "" + t_worst = 0.0 + p_worst = 0.0 + idxs = range(0, P.shape[1], 10) + for e in edges: + if not isinstance(e, (OrificeEdge, ShortTubeEdge)): + continue + d_char = ( + e.D + if isinstance(e, ShortTubeEdge) + else math.sqrt(4.0 * e.A_total / math.pi) + ) + a = node_index[e.a] + b = e.b if e.b == EXT_NODE else node_index[e.b] + for k in idxs: + pa = float(P[a, k]) + pb = float(P_ext[k]) if b == EXT_NODE else float(P[b, k]) + if pa >= pb: + p_up = pa + t_up = float(T[a, k]) + else: + p_up = pb + t_up = float(T[a, k] if b == EXT_NODE else T[b, k]) + p_up = max(p_up, 1e-9) + t_up = max(t_up, T_SAFE) + mfp = ( + K_BOLTZMANN * t_up / (math.sqrt(2.0) * math.pi * (D_MOL_AIR**2) * p_up) + ) + kn = mfp / max(d_char, 1e-12) + if kn > kn_max: + kn_max = kn + edge_worst = e.label or f"edge({e.a}->{e.b})" + t_worst = float(t[k]) + p_worst = p_up + + if kn_max < 0.01: + status, regime = "ok", "continuum" + elif kn_max < 0.1: + status, regime = "warning", "slip" + else: + status, regime = "fail", "transitional/free-molecular" + return { + "status": status, + "Kn_max": kn_max, + "edge_worst": edge_worst, + "t_worst_s": t_worst, + "P_at_worst_Pa": p_worst, + "regime": regime, + "message": "Continuum assumptions degrade as Knudsen number increases", + } + + +def _state_flag(m: np.ndarray, T: np.ndarray) -> dict: + min_m = float(np.min(m)) + min_t = float(np.min(T)) + if np.any(~np.isfinite(m)) or np.any(~np.isfinite(T)): + return {"status": "fail", "message": "NaN/Inf detected in state arrays"} + if min_m < -1e-12: + return { + "status": "fail", + "min_m_kg": min_m, + "message": "Negative mass detected", + } + if min_t < 0.0: + return { + "status": "fail", + "min_T_K": min_t, + "message": "Negative temperature detected", + } + if min_t < T_SAFE: + return { + "status": "warning", + "min_T_K": min_t, + "message": "Temperature reached denominator safety floor", + } + return {"status": "ok", "min_m_kg": min_m, "min_T_K": min_t} + + +def evaluate_validity_flags( + nodes: list, + edges: list, + P: np.ndarray, + T: np.ndarray, + m: np.ndarray, + P_ext: np.ndarray, + t: np.ndarray, + l_char_m: float, +) -> dict: + node_index = {i: i for i in range(len(nodes))} + flags = { + "state_integrity": _state_flag(m, T), + "acoustic_uniformity_0D": _acoustic_flag(P, T, t, l_char_m), + "thermo_fit_range": _thermo_range_flag(T), + "slot_laminarity": _slot_laminar_flag(edges, P, T, node_index), + "short_tube_flow": _short_tube_flag(edges, P, T, P_ext, node_index), + "knudsen_regime": _knudsen_flag(edges, P, T, P_ext, node_index, t), + } + + max_ext = float(np.max(P_ext)) if P_ext.size else 0.0 + flags["external_pressure_units"] = { + "status": "warning" if 200.0 <= max_ext <= 2000.0 else "ok", + "max_P_ext_Pa": max_ext, + "message": "Values in 200-2000 range are often mmHg entered as Pa", + } + return flags diff --git a/tests/test_compare.py b/tests/test_compare.py new file mode 100644 index 0000000..ee00376 --- /dev/null +++ b/tests/test_compare.py @@ -0,0 +1,54 @@ +from pathlib import Path + +import pytest + +from venting.cases import CaseConfig, NetworkConfig +from venting.compare import compare_runs, load_run +from venting.presets import get_default_panel_preset_v9 +from venting.profiles import Profile +from venting.run import export_case_artifacts, run_case + + +def _make_run(tmp_path: Path, stem: str, cd_int: float) -> Path: + prof = Profile("vac", lambda t: 0.0, events=((0.0, "vac"),)) + p = get_default_panel_preset_v9() + net = NetworkConfig( + N_chain=2, + N_par=1, + V_cell=p.V_cell, + V_vest=p.V_vest, + A_wall_cell=p.A_wall_cell, + A_wall_vest=p.A_wall_vest, + d_int_mm=2.0, + n_int_per_interface=1, + d_exit_mm=2.0, + n_exit=1, + Cd_int=cd_int, + Cd_exit=0.62, + ) + case = CaseConfig("isothermal", 0.0, 300.0, 0.03, 40) + res = run_case(net, prof, case) + out = tmp_path / stem + out.mkdir() + export_case_artifacts(out, stem, res, run_params={"x": 1}) + return out + + +def test_compare_identical(tmp_path: Path): + run = _make_run(tmp_path, "a", 0.62) + comp = compare_runs(load_run(run), load_run(run)) + assert all( + abs(v["delta_max_abs_dP_Pa"]) < 1e-12 for v in comp["edge_deltas"].values() + ) + + +def test_compare_different_cd(tmp_path: Path): + run_a = _make_run(tmp_path, "a", 0.62) + run_b = _make_run(tmp_path, "b", 0.55) + comp = compare_runs(load_run(run_a), load_run(run_b)) + assert any(abs(v["delta_max_abs_dP_Pa"]) > 0 for v in comp["edge_deltas"].values()) + + +def test_compare_missing_dir(tmp_path: Path): + with pytest.raises(FileNotFoundError): + load_run(tmp_path / "missing") diff --git a/tests/test_gates.py b/tests/test_gates.py new file mode 100644 index 0000000..be7730a --- /dev/null +++ b/tests/test_gates.py @@ -0,0 +1,120 @@ +import math + +import numpy as np + +from venting.cases import CaseConfig +from venting.constants import C_CHOKED, C_P, C_V, GAMMA, P0, R_GAS, T0, T_SAFE +from venting.diagnostics import summarize_result +from venting.flow import mdot_orifice_pos +from venting.geometry import circle_area_from_d_mm +from venting.graph import EXT_NODE, CdConst, ExternalBC, GasNode, OrificeEdge +from venting.profiles import Profile +from venting.solver import solve_case + + +def _single_node_cases(): + V = 131.6e-6 + A = circle_area_from_d_mm(2.0) + Cd = 0.62 + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", V, 181.6e-4)] + edges = [OrificeEdge(0, EXT_NODE, A, CdConst(Cd), label="exit")] + bcs = [ExternalBC(0, profile, T_ext=T0)] + alpha = Cd * A * C_CHOKED * math.sqrt(R_GAS * T0) / V + tau = 1.0 / alpha + return nodes, edges, bcs, alpha, tau, V, A, Cd + + +def test_single_node_analytic_match_adiabatic(): + nodes, edges, bcs, alpha, _, _, _, _ = _single_node_cases() + beta = (GAMMA - 1.0) / 2.0 + + def p_adi(t): + return P0 * (1.0 + beta * alpha * t) ** (-2.0 * GAMMA / (GAMMA - 1.0)) + + def t_adi(t): + return T0 * (1.0 + beta * alpha * t) ** (-2.0) + + case = CaseConfig("intermediate", 0.0, T0, 1.5, 700) + res = summarize_result(nodes, edges, bcs, case, solve_case(nodes, edges, bcs, case)) + mask = res.P[0, :] > 0.01 * P0 + t = res.t[mask] + p_err = np.max( + np.abs(res.P[0, mask] - np.array([p_adi(float(tt)) for tt in t])) / P0 + ) + t_err = np.max( + np.abs(res.T[0, mask] - np.array([t_adi(float(tt)) for tt in t])) / T0 + ) + assert float(p_err) < 5e-3 + assert float(t_err) < 5e-3 + + +def test_single_node_analytic_match_isothermal_limit(): + nodes, edges, bcs, alpha, _, _, _, _ = _single_node_cases() + + def p_iso(t): + return P0 * math.exp(-alpha * t) + + case = CaseConfig("intermediate", 1e6, T0, 1.5, 700) + res = summarize_result(nodes, edges, bcs, case, solve_case(nodes, edges, bcs, case)) + mask = res.P[0, :] > 0.01 * P0 + t = res.t[mask] + p_err = np.max( + np.abs(res.P[0, mask] - np.array([p_iso(float(tt)) for tt in t])) / P0 + ) + assert float(p_err) < 5e-3 + + +def test_mass_conservation_single_node(): + nodes, edges, bcs, _, _, V, _, _ = _single_node_cases() + case = CaseConfig("intermediate", 0.0, T0, 1.5, 700) + res = summarize_result(nodes, edges, bcs, case, solve_case(nodes, edges, bcs, case)) + m0 = P0 * V / (R_GAS * T0) + mf = float(res.m[0, -1]) + m_out = m0 - mf + err = abs((m0 - mf) - m_out) / m0 + assert err < 1e-3 + + +def test_two_node_mass_conservation(): + Vc, Vv = 131.6e-6, 145.3e-6 + A = circle_area_from_d_mm(2.0) + Cd = 0.62 + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("vest", Vv, 181.6e-4), GasNode("cell", Vc, 181.6e-4)] + edges = [ + OrificeEdge(1, 0, A, CdConst(Cd), label="cell↔vest"), + OrificeEdge(0, EXT_NODE, A, CdConst(Cd), label="exit"), + ] + bcs = [ExternalBC(0, profile, T_ext=T0)] + case = CaseConfig("intermediate", 0.0, T0, 2.0, 900) + res = summarize_result(nodes, edges, bcs, case, solve_case(nodes, edges, bcs, case)) + + m_total_0 = float(np.sum(res.m[:, 0])) + m_total_f = float(np.sum(res.m[:, -1])) + mdot = np.array( + [ + mdot_orifice_pos( + float(res.P[0, k]), float(max(res.T[0, k], T_SAFE)), 0.0, Cd, A + ) + for k in range(len(res.t)) + ], + dtype=float, + ) + m_out = float(np.trapz(mdot, res.t)) + err = abs((m_total_f + m_out) - m_total_0) / m_total_0 + assert err < 1e-3 + + E0 = float(np.sum(res.m[:, 0] * C_V * res.T[:, 0])) + Ef = float(np.sum(res.m[:, -1] * C_V * res.T[:, -1])) + Eout = float(np.trapz(mdot * C_P * np.maximum(res.T[0, :], T_SAFE), res.t)) + err_e = abs((Ef + Eout) - E0) / E0 + assert err_e < 1e-2 + + +def test_monotonic_pressure_when_vacuum(): + nodes, edges, bcs, _, _, _, _, _ = _single_node_cases() + case = CaseConfig("intermediate", 0.0, T0, 1.5, 700) + res = summarize_result(nodes, edges, bcs, case, solve_case(nodes, edges, bcs, case)) + d = np.diff(res.P[0, :]) + assert np.max(d) <= 1e-6 * P0 diff --git a/tests/test_gui_config.py b/tests/test_gui_config.py new file mode 100644 index 0000000..5d8fb55 --- /dev/null +++ b/tests/test_gui_config.py @@ -0,0 +1,31 @@ +import pytest + +from venting.gui.config import GuiCaseConfig + + +def test_gui_config_roundtrip_json(): + cfg = GuiCaseConfig( + topology="two_chain_shared_vest", + N_chain=7, + N_chain_b=5, + profile_kind="step", + step_time_s=0.2, + ) + loaded = GuiCaseConfig.from_json(cfg.to_json()) + assert loaded.topology == "two_chain_shared_vest" + assert loaded.N_chain == 7 + assert loaded.N_chain_b == 5 + assert loaded.profile_kind == "step" + assert loaded.step_time_s == 0.2 + + +def test_gui_config_validation_rejects_bad_enum(): + cfg = GuiCaseConfig(thermo="bad") + with pytest.raises(ValueError): + cfg.validate() + + +def test_gui_config_validation_rejects_nonpositive(): + cfg = GuiCaseConfig(V_cell_m3=0.0) + with pytest.raises(ValueError): + cfg.validate() diff --git a/tests/test_gui_import.py b/tests/test_gui_import.py new file mode 100644 index 0000000..77df70f --- /dev/null +++ b/tests/test_gui_import.py @@ -0,0 +1,6 @@ +def test_gui_modules_import_without_qt_runtime(): + import venting.gui.app as app + import venting.gui.main as main + + assert hasattr(app, "create_main_window") + assert hasattr(main, "main") diff --git a/tests/test_gui_state_layout.py b/tests/test_gui_state_layout.py new file mode 100644 index 0000000..ee7dec0 --- /dev/null +++ b/tests/test_gui_state_layout.py @@ -0,0 +1,31 @@ +import pytest + +from venting.gui.config import GuiCaseConfig +from venting.gui.state_layout import infer_layout, infer_layout_from_modes + + +def test_layout_isothermal(): + cfg = GuiCaseConfig(thermo="isothermal", wall_model="fixed") + layout = infer_layout(cfg, 4) + assert layout.m_slice == slice(0, 4) + assert layout.t_slice is None + assert layout.tw_slice is None + + +def test_layout_intermediate_fixed(): + layout = infer_layout_from_modes("intermediate", "fixed", 3) + assert layout.m_slice == slice(0, 3) + assert layout.t_slice == slice(3, 6) + assert layout.tw_slice is None + + +def test_layout_variable_lumped(): + layout = infer_layout_from_modes("variable", "lumped", 5) + assert layout.m_slice == slice(0, 5) + assert layout.t_slice == slice(5, 10) + assert layout.tw_slice == slice(10, 15) + + +def test_layout_rejects_nonpositive_node_count(): + with pytest.raises(ValueError): + infer_layout_from_modes("intermediate", "fixed", 0) diff --git a/tests/test_jac_sparsity.py b/tests/test_jac_sparsity.py new file mode 100644 index 0000000..1297a1c --- /dev/null +++ b/tests/test_jac_sparsity.py @@ -0,0 +1,35 @@ +from pathlib import Path + +from venting.cases import CaseConfig, NetworkConfig +from venting.graph import build_branching_network +from venting.presets import get_default_panel_preset_v9 +from venting.profiles import make_profile_linear +from venting.solver import solve_case + + +def test_solver_does_not_use_identity_only_jac_sparsity_hint(): + src = Path("src/venting/solver.py").read_text(encoding="utf-8") + assert "jac_sparsity" not in src + + +def test_large_network_solves_with_radau(): + preset = get_default_panel_preset_v9() + cfg = NetworkConfig( + N_chain=20, + 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=2.0, + n_int_per_interface=1, + d_exit_mm=2.0, + n_exit=1, + Cd_int=0.62, + Cd_exit=0.62, + ) + profile = make_profile_linear(101325.0, 20.0) + nodes, edges, bcs = build_branching_network(cfg, profile) + case = CaseConfig("intermediate", 0.0, 300.0, 1.0, 120) + sol = solve_case(nodes, edges, bcs, case) + assert sol.success diff --git a/tests/test_montecarlo.py b/tests/test_montecarlo.py new file mode 100644 index 0000000..10bf9ae --- /dev/null +++ b/tests/test_montecarlo.py @@ -0,0 +1,46 @@ +from venting.cases import CaseConfig, NetworkConfig +from venting.montecarlo import run_mc +from venting.presets import get_default_panel_preset_v9 +from venting.profiles import Profile + + +def _base_net() -> NetworkConfig: + p = get_default_panel_preset_v9() + return NetworkConfig( + N_chain=2, + N_par=1, + V_cell=p.V_cell, + V_vest=p.V_vest, + A_wall_cell=p.A_wall_cell, + A_wall_vest=p.A_wall_vest, + d_int_mm=2.0, + n_int_per_interface=1, + d_exit_mm=2.0, + n_exit=1, + Cd_int=0.62, + Cd_exit=0.62, + ) + + +def _case() -> CaseConfig: + return CaseConfig("isothermal", 0.0, 300.0, 0.05, 40) + + +def test_mc_deterministic(): + prof = Profile("vac", lambda t: 0.0, events=((0.0, "vac"),)) + r1 = run_mc(_base_net(), _case(), prof, (0.5, 0.7), (0.55, 0.65), 5, seed=42) + r2 = run_mc(_base_net(), _case(), prof, (0.5, 0.7), (0.55, 0.65), 5, seed=42) + assert r1["samples"] == r2["samples"] + + +def test_mc_zero_range(): + prof = Profile("vac", lambda t: 0.0, events=((0.0, "vac"),)) + r = run_mc(_base_net(), _case(), prof, (0.62, 0.62), (0.62, 0.62), 5, seed=1) + for edge in r["summary"].values(): + assert abs(edge["std"]) < 1e-12 + + +def test_mc_has_spread(): + prof = Profile("vac", lambda t: 0.0, events=((0.0, "vac"),)) + r = run_mc(_base_net(), _case(), prof, (0.5, 0.7), (0.55, 0.65), 5, seed=1) + assert any(edge["std"] > 0 for edge in r["summary"].values()) diff --git a/tests/test_regression_cli_baseline.py b/tests/test_regression_cli_baseline.py new file mode 100644 index 0000000..1b16c94 --- /dev/null +++ b/tests/test_regression_cli_baseline.py @@ -0,0 +1,11 @@ +from venting.gates import gate_single, gate_two + + +def test_gate_baseline_regression_stable(): + g1 = gate_single() + g2 = gate_two() + assert g1.errP < 1e-6 + assert g1.errT < 1e-6 + assert g1.errMass < 1e-4 + assert g2.errMass < 1e-4 + assert g2.errEnergy is not None and g2.errEnergy < 1e-4 diff --git a/tests/test_short_tube.py b/tests/test_short_tube.py new file mode 100644 index 0000000..f1b65cd --- /dev/null +++ b/tests/test_short_tube.py @@ -0,0 +1,217 @@ +import numpy as np + +from venting.cases import CaseConfig +from venting.constants import T0 +from venting.diagnostics import summarize_result +from venting.flow import mdot_fanno_tube, mdot_orifice_pos, mdot_short_tube_pos +from venting.geometry import circle_area_from_d_mm +from venting.graph import CdConst, ExternalBC, GasNode, ShortTubeEdge +from venting.profiles import Profile +from venting.solver import solve_case + + +def test_short_tube_L0_matches_orifice(): + p_up = 101325.0 + p_dn = 50000.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + + m_or = mdot_orifice_pos(p_up, t_up, p_dn, cd0, area) + m_st = mdot_short_tube_pos(p_up, t_up, p_dn, cd0, area, diam, 0.0, 0.0, 0.0, 0.0) + rel = abs(m_st - m_or) / max(m_or, 1e-20) + assert rel < 0.02 + + +def test_short_tube_mdot_decreases_with_L(): + p_up = 101325.0 + p_dn = 30000.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + + m_05 = mdot_short_tube_pos(p_up, t_up, p_dn, cd0, area, diam, 0.5e-3, 0.0, 0.5, 1.0) + m_10 = mdot_short_tube_pos(p_up, t_up, p_dn, cd0, area, diam, 1.0e-3, 0.0, 0.5, 1.0) + m_20 = mdot_short_tube_pos(p_up, t_up, p_dn, cd0, area, diam, 2.0e-3, 0.0, 0.5, 1.0) + + assert m_05 > m_10 > m_20 + + +def test_short_tube_mdot_decreases_with_roughness(): + p_up = 2.0e5 + p_dn = 1.0e5 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(4.0) + diam = 4.0e-3 + + m_smooth = mdot_short_tube_pos( + p_up, t_up, p_dn, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0 + ) + m_rough = mdot_short_tube_pos( + p_up, t_up, p_dn, cd0, area, diam, 3e-3, 200e-6, 0.5, 1.0 + ) + assert m_smooth > m_rough + + +def test_short_tube_validity_fields_present(): + vol = 131.6e-6 + area = circle_area_from_d_mm(2.0) + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", vol, 181.6e-4)] + edges = [ + ShortTubeEdge( + a=0, + b=-1, + A_total=area, + D=2.0e-3, + L=1.0e-3, + eps=0.0, + K_in=0.5, + K_out=1.0, + Cd_model=CdConst(0.62), + label="exit_short_tube", + ) + ] + bcs = [ExternalBC(0, profile, T_ext=T0)] + case = CaseConfig("intermediate", 0.0, T0, 0.5, 200) + sol = solve_case(nodes, edges, bcs, case) + res = summarize_result(nodes, edges, bcs, case, sol) + flag = res.meta["validity_flags"]["short_tube_flow"] + + assert "Mach_max" in flag + assert "Re_max" in flag + assert "Cd_eff_min" in flag + assert "K_tot_max" in flag + assert "frac_fric" in flag + + +def test_short_tube_in_network_runs(): + vol = 131.6e-6 + area = circle_area_from_d_mm(2.0) + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", vol, 181.6e-4)] + edges = [ + ShortTubeEdge( + a=0, + b=-1, + A_total=area, + D=2.0e-3, + L=1.0e-3, + eps=0.0, + K_in=0.5, + K_out=1.0, + Cd_model=CdConst(0.62), + label="exit_short_tube", + ) + ] + bcs = [ExternalBC(0, profile, T_ext=T0)] + case = CaseConfig("intermediate", 0.0, T0, 0.5, 200) + sol = solve_case(nodes, edges, bcs, case) + assert sol.success + assert np.isfinite(sol.y).all() + + +def test_fanno_choked_less_than_isentropic(): + p_up = 101325.0 + p_dn = 0.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + m_fanno = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0) + m_lossy = mdot_short_tube_pos( + p_up, t_up, p_dn, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0 + ) + assert m_fanno < m_lossy + + +def test_fanno_zero_length_matches_orifice(): + p_up = 101325.0 + p_dn = 50000.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + m_or = mdot_orifice_pos(p_up, t_up, p_dn, cd0, area) + m_f = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, 0.0, 0.0, 0.5, 1.0) + assert abs(m_f - m_or) / max(m_or, 1e-20) < 1e-3 + + +def test_fanno_low_mach_matches_lossy_nozzle(): + p_up = 101325.0 + p_dn = 0.9 * p_up + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + m_fanno = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0) + m_lossy = mdot_short_tube_pos( + p_up, t_up, p_dn, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0 + ) + assert abs(m_fanno - m_lossy) / max(m_lossy, 1e-20) < 0.05 + + +def test_fanno_choking_mach_limit(): + p_up = 101325.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + m1 = mdot_fanno_tube(p_up, t_up, 5000.0, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0) + m2 = mdot_fanno_tube(p_up, t_up, 500.0, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0) + assert m2 <= 1.05 * m1 + + +def test_fanno_sees_K_losses(): + """Fanno mdot must decrease when K_in/K_out increase.""" + p_up = 101325.0 + p_dn = 0.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + L = 3e-3 + + m_k0 = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, L, 0.0, 0.0, 0.0) + m_k15 = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, L, 0.0, 0.5, 1.0) + m_k3 = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, L, 0.0, 1.0, 2.0) + + assert m_k15 < m_k0, f"K=1.5 should give less flow: {m_k15} vs {m_k0}" + assert m_k3 < m_k15, f"K=3.0 should give less flow: {m_k3} vs {m_k15}" + + +def test_fanno_with_K_below_lossy(): + """For typical parameters, Fanno with K losses should be below lossy nozzle.""" + p_up = 101325.0 + p_dn = 0.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + L = 3e-3 + + m_fanno = mdot_fanno_tube(p_up, t_up, p_dn, cd0, area, diam, L, 0.0, 0.5, 1.0) + m_lossy = mdot_short_tube_pos(p_up, t_up, p_dn, cd0, area, diam, L, 0.0, 0.5, 1.0) + assert m_fanno < m_lossy, ( + f"Fanno with K should give LESS flow than lossy nozzle: " + f"fanno={m_fanno:.6f}, lossy={m_lossy:.6f}" + ) + assert m_fanno / m_lossy < 1.0 + + +def test_fanno_monotonic_with_K(): + """Flow should be monotonic non-decreasing as downstream pressure decreases.""" + p_up = 101325.0 + t_up = 300.0 + cd0 = 0.62 + area = circle_area_from_d_mm(2.0) + diam = 2.0e-3 + + prev = 0.0 + for r in [0.99, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]: + m = mdot_fanno_tube(p_up, t_up, r * p_up, cd0, area, diam, 3e-3, 0.0, 0.5, 1.0) + assert m >= 0.98 * prev - 1e-12, f"Non-monotonic at r={r}: {m} < {prev}" + prev = m diff --git a/tests/test_solver_stream.py b/tests/test_solver_stream.py new file mode 100644 index 0000000..00df15a --- /dev/null +++ b/tests/test_solver_stream.py @@ -0,0 +1,92 @@ +import numpy as np + +from venting.cases import CaseConfig +from venting.constants import T0 +from venting.diagnostics import summarize_result +from venting.geometry import circle_area_from_d_mm +from venting.graph import EXT_NODE, CdConst, ExternalBC, GasNode, OrificeEdge +from venting.profiles import Profile +from venting.solver import solve_case, solve_case_stream + +TOL_STREAM = 1e-3 + + +def _canonical_single_node(): + v = 131.6e-6 + a = circle_area_from_d_mm(2.0) + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", v, 181.6e-4)] + edges = [OrificeEdge(0, EXT_NODE, a, CdConst(0.62), label="exit")] + bcs = [ExternalBC(0, profile, T_ext=T0)] + case = CaseConfig("intermediate", 0.0, T0, 0.4, 120) + return nodes, edges, bcs, case + + +def test_streaming_matches_batch_pressure(): + nodes, edges, bcs, case = _canonical_single_node() + sol_batch = solve_case(nodes, edges, bcs, case) + sol_stream = solve_case_stream(nodes, edges, bcs, case) + + res_b = summarize_result(nodes, edges, bcs, case, sol_batch) + res_s = summarize_result(nodes, edges, bcs, case, sol_stream) + + rel = np.max(np.abs(res_b.P - res_s.P) / np.maximum(np.abs(res_b.P), 1.0)) + assert float(rel) < TOL_STREAM + + +def test_streaming_callback_receives_progressive_chunks(): + nodes, edges, bcs, case = _canonical_single_node() + seen = [] + + def cb(payload): + seen.append((payload["progress"], payload["t"].shape[0], payload["t"][-1])) + + solve_case_stream(nodes, edges, bcs, case, callback=cb, n_chunks=8) + assert len(seen) >= 2 + assert seen[-1][0] == 1.0 + assert all(seen[i][1] <= seen[i + 1][1] for i in range(len(seen) - 1)) + assert all(seen[i][2] <= seen[i + 1][2] for i in range(len(seen) - 1)) + assert seen[0][2] < case.duration + + +def test_streaming_can_cancel_early(): + nodes, edges, bcs, case = _canonical_single_node() + stop = {"value": False} + + def cb(_payload): + stop["value"] = True + + sol = solve_case_stream( + nodes, + edges, + bcs, + case, + callback=cb, + n_chunks=10, + should_stop=lambda: stop["value"], + ) + assert sol.t[-1] < case.duration + + +def test_dt_chunk_default(): + nodes, edges, bcs, case = _canonical_single_node() + case = CaseConfig(case.thermo, case.h_conv, case.T_wall, 10.0, 101) + seen = [] + + 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 + + +def test_dt_chunk_small_duration(): + nodes, edges, bcs, case = _canonical_single_node() + case = CaseConfig(case.thermo, case.h_conv, case.T_wall, 0.01, 10) + seen = [] + + def cb(payload): + seen.append(payload["progress"]) + + solve_case_stream(nodes, edges, bcs, case, callback=cb, dt_chunk_s=2.0) + assert len(seen) >= 1 diff --git a/tests/test_state_layout.py b/tests/test_state_layout.py new file mode 100644 index 0000000..b636ac2 --- /dev/null +++ b/tests/test_state_layout.py @@ -0,0 +1,46 @@ +import numpy as np + +from venting.state_layout import infer_layout_from_modes, split_state + + +def _make_y(n_vars: int, n_t: int = 4): + return np.arange(n_vars * n_t, dtype=float).reshape(n_vars, n_t) + + +def test_state_layout_isothermal_n1_n5(): + for n in (1, 5): + y = _make_y(n) + layout = infer_layout_from_modes("isothermal", "fixed", n) + assert layout.m_slice == slice(0, n) + assert layout.t_slice is None + assert layout.tw_slice is None + m, t, tw = split_state(y, n, "isothermal", "fixed") + assert m.shape == (n, y.shape[1]) + assert t is None and tw is None + + +def test_state_layout_intermediate_fixed_n1_n5(): + for n in (1, 5): + y = _make_y(2 * n) + m, t, tw = split_state(y, n, "intermediate", "fixed") + assert m.shape == (n, y.shape[1]) + assert t.shape == (n, y.shape[1]) + assert tw is None + + +def test_state_layout_variable_fixed_n1_n5(): + for n in (1, 5): + y = _make_y(2 * n) + m, t, tw = split_state(y, n, "variable", "fixed") + assert m.shape == (n, y.shape[1]) + assert t.shape == (n, y.shape[1]) + assert tw is None + + +def test_state_layout_intermediate_lumped_n1_n5(): + for n in (1, 5): + y = _make_y(3 * n) + m, t, tw = split_state(y, n, "intermediate", "lumped") + assert m.shape == (n, y.shape[1]) + assert t.shape == (n, y.shape[1]) + assert tw.shape == (n, y.shape[1]) diff --git a/tests/test_thermo_reference.py b/tests/test_thermo_reference.py new file mode 100644 index 0000000..4c6743a --- /dev/null +++ b/tests/test_thermo_reference.py @@ -0,0 +1,21 @@ +import numpy as np + +from venting.thermo import cp_air, gamma_air, h_air + + +def test_cp_at_300K(): + assert abs(cp_air(300.0) - 1006.0) < 5.0 + + +def test_gamma_at_300K(): + assert abs(gamma_air(300.0) - 1.4) < 0.005 + + +def test_cp_monotonic_200_to_600(): + temps = np.linspace(200.0, 600.0, 20) + cps = [cp_air(float(t)) for t in temps] + assert all(cps[i] <= cps[i + 1] for i in range(len(cps) - 1)) + + +def test_h_zero_at_reference(): + assert abs(h_air(298.15)) < 1.0 diff --git a/tests/test_topology_and_scalars.py b/tests/test_topology_and_scalars.py new file mode 100644 index 0000000..842a698 --- /dev/null +++ b/tests/test_topology_and_scalars.py @@ -0,0 +1,70 @@ +import numpy as np + +from venting.cases import CaseConfig, NetworkConfig +from venting.diagnostics import summarize_result +from venting.graph import build_branching_network +from venting.presets import get_default_panel_preset_v9 +from venting.profiles import make_profile_step +from venting.solver import solve_case + + +def _base_cfg() -> NetworkConfig: + p = get_default_panel_preset_v9() + return NetworkConfig( + N_chain=4, + N_par=1, + V_cell=p.V_cell, + V_vest=p.V_vest, + A_wall_cell=p.A_wall_cell, + A_wall_vest=p.A_wall_vest, + d_int_mm=2.0, + n_int_per_interface=1, + d_exit_mm=2.0, + n_exit=1, + Cd_int=0.62, + Cd_exit=0.62, + ) + + +def test_scalar_or_list_expansion_compatible(): + prof = make_profile_step(101325.0, 0.01) + cfg = _base_cfg() + nodes_s, edges_s, _ = build_branching_network(cfg, prof) + + cfg_list = NetworkConfig( + **{ + **cfg.__dict__, + "V_cell": [cfg.V_cell] * cfg.N_chain, + "A_wall_cell": [cfg.A_wall_cell] * cfg.N_chain, + "d_int_mm": [2.0] * cfg.N_chain, + "n_int_per_interface": [1] * cfg.N_chain, + "Cd_int": [0.62] * cfg.N_chain, + "L_int_mm": [0.0] * cfg.N_chain, + "eps_int_um": [0.0] * cfg.N_chain, + "K_in_int": [0.5] * cfg.N_chain, + "K_out_int": [1.0] * cfg.N_chain, + } + ) + nodes_l, edges_l, _ = build_branching_network(cfg_list, prof) + assert len(nodes_s) == len(nodes_l) + assert len(edges_s) == len(edges_l) + + +def test_two_chain_shared_vest_symmetry(): + prof = make_profile_step(101325.0, 0.01) + cfg = _base_cfg() + cfg = NetworkConfig( + **{**cfg.__dict__, "topology": "two_chain_shared_vest", "N_chain_b": 4} + ) + nodes, edges, bcs = build_branching_network(cfg, prof) + case = CaseConfig("intermediate", 0.0, 300.0, 0.2, 100) + sol = solve_case(nodes, edges, bcs, case) + res = summarize_result(nodes, edges, bcs, case, sol) + + # A_cell1 and B_cell1 are symmetric for identical chains + idx_a1 = next(i for i, n in enumerate(sol.nodes_local) if n.name == "A_cell1") + idx_b1 = next(i for i, n in enumerate(sol.nodes_local) if n.name == "B_cell1") + rel = np.max( + np.abs(res.P[idx_a1] - res.P[idx_b1]) / np.maximum(np.abs(res.P[idx_a1]), 1.0) + ) + assert float(rel) < 1e-6 diff --git a/tests/test_v9_features.py b/tests/test_v9_features.py new file mode 100644 index 0000000..f6c9f19 --- /dev/null +++ b/tests/test_v9_features.py @@ -0,0 +1,113 @@ +import numpy as np + +from venting.cases import CaseConfig +from venting.constants import T0 +from venting.diagnostics import summarize_result +from venting.flow import mdot_orifice_pos_props +from venting.geometry import circle_area_from_d_mm +from venting.graph import EXT_NODE, CdConst, ExternalBC, GasNode, OrificeEdge +from venting.profiles import Profile +from venting.solver import solve_case +from venting.thermo import gamma_air + + +def _single_node(): + vol = 131.6e-6 + area = circle_area_from_d_mm(2.0) + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", vol, 181.6e-4)] + edges = [OrificeEdge(0, EXT_NODE, area, CdConst(0.62), label="exit")] + bcs = [ExternalBC(0, profile, T_ext=T0)] + return nodes, edges, bcs + + +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) + r_int = summarize_result( + nodes, edges, bcs, c_int, solve_case(nodes, edges, bcs, c_int) + ) + r_var = summarize_result( + nodes, edges, bcs, c_var, solve_case(nodes, edges, bcs, c_var) + ) + 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 + + +def test_variable_gamma_changes_mdot(): + p_up = 101325.0 + p_dn = 30000.0 + area = np.pi * (0.002**2) / 4.0 + m_cold = mdot_orifice_pos_props( + p_up, 260.0, p_dn, 0.62, area, gamma=gamma_air(260.0) + ) + m_hot = mdot_orifice_pos_props( + p_up, 360.0, p_dn, 0.62, area, gamma=gamma_air(360.0) + ) + assert abs(m_cold - m_hot) > 0.0 + + +def test_lumped_wall_large_capacity_matches_fixed(): + nodes, edges, bcs = _single_node() + fixed = CaseConfig("intermediate", 5.0, 300.0, 0.1, 40, wall_model="fixed") + lumped = CaseConfig( + "intermediate", + 5.0, + 300.0, + 0.1, + 40, + wall_model="lumped", + wall_C_per_area=1e12, + wall_h_out=0.0, + ) + r_fixed = summarize_result( + nodes, edges, bcs, fixed, solve_case(nodes, edges, bcs, fixed) + ) + r_lumped = summarize_result( + nodes, edges, bcs, lumped, solve_case(nodes, edges, bcs, lumped) + ) + rel = np.max(np.abs(r_lumped.P[0] - r_fixed.P[0]) / np.maximum(r_fixed.P[0], 1.0)) + assert float(rel) <= 0.01 + + +def test_h_in_zero_matches_adiabatic(): + nodes, edges, bcs = _single_node() + c_fixed = CaseConfig("intermediate", 0.0, 300.0, 0.1, 40, wall_model="fixed") + c_lumped = CaseConfig( + "intermediate", + 0.0, + 300.0, + 0.1, + 40, + wall_model="lumped", + wall_C_per_area=1e5, + wall_h_out=4.0, + ) + r_fixed = summarize_result( + nodes, edges, bcs, c_fixed, solve_case(nodes, edges, bcs, c_fixed) + ) + r_lumped = summarize_result( + nodes, edges, bcs, c_lumped, solve_case(nodes, edges, bcs, c_lumped) + ) + rel = np.max(np.abs(r_lumped.P[0] - r_fixed.P[0]) / np.maximum(r_fixed.P[0], 1.0)) + assert float(rel) <= 0.01 + + +def test_dynamic_pump_smoke(): + nodes, edges, bcs = _single_node() + case = CaseConfig( + "intermediate", + 0.0, + 300.0, + 0.1, + 40, + external_model="dynamic_pump", + V_ext=0.2, + T_ext=300.0, + pump_speed_m3s=0.01, + P_ult_Pa=10.0, + ) + sol = solve_case(nodes, edges, bcs, case) + assert sol.success + assert np.isfinite(sol.y).all() diff --git a/tests/test_validity_and_units.py b/tests/test_validity_and_units.py new file mode 100644 index 0000000..4ba8508 --- /dev/null +++ b/tests/test_validity_and_units.py @@ -0,0 +1,59 @@ +from pathlib import Path + +import pytest + +from venting.cases import CaseConfig +from venting.constants import T0 +from venting.diagnostics import summarize_result +from venting.geometry import circle_area_from_d_mm +from venting.graph import EXT_NODE, CdConst, ExternalBC, GasNode, OrificeEdge +from venting.profiles import Profile, make_profile_from_table +from venting.solver import solve_case + + +def test_profile_table_rejects_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") + + +def test_profile_table_accepts_mmhg_when_explicit(tmp_path: Path): + p = tmp_path / "profile.csv" + p.write_text("0,760\n1,740\n", encoding="utf-8") + prof = make_profile_from_table("tab", p, pressure_unit="mmHg") + assert prof.P(0.0) > 1.0e5 + + +def test_validity_flags_are_emitted_in_meta(): + V = 131.6e-6 + A = circle_area_from_d_mm(2.0) + profile = Profile("vacuum", lambda t: 0.0, events=((0.0, "vacuum"),)) + nodes = [GasNode("node", V, 181.6e-4)] + edges = [OrificeEdge(0, EXT_NODE, A, CdConst(0.62), label="exit")] + bcs = [ExternalBC(0, profile, T_ext=T0)] + case = CaseConfig("intermediate", 0.0, T0, 0.8, 300) + sol = solve_case(nodes, edges, bcs, case) + res = summarize_result(nodes, edges, bcs, case, sol) + + flags = res.meta.get("validity_flags", {}) + assert "state_integrity" in flags + assert "acoustic_uniformity_0D" in flags + assert "external_pressure_units" in flags + assert "thermo_fit_range" in flags + assert "knudsen_regime" in flags + assert flags["state_integrity"]["status"] in {"ok", "warning", "fail"} + + +def test_knudsen_reference_values(): + import math + + from venting.constants import D_MOL_AIR, K_BOLTZMANN + + def kn(p, t=300.0, d=2e-3): + lam = K_BOLTZMANN * t / (math.sqrt(2.0) * math.pi * (D_MOL_AIR**2) * p) + return lam / d + + assert kn(101325.0) < 0.01 + assert kn(1.0) >= 0.1 + assert 0.01 <= kn(100.0) < 0.1