From 70c91ed019c710c48842ee91697bac28c8a58ed1 Mon Sep 17 00:00:00 2001 From: ljz <425868052@qq.com> Date: Sat, 11 Jul 2026 09:33:25 +0800 Subject: [PATCH] =?UTF-8?q?=E4=B8=8A=E4=BC=A0PythonModels=E6=96=87?= =?UTF-8?q?=E4=BB=B6?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- PythonModels/.gitignore | 2 + PythonModels/README.md | 340 +++++++++ PythonModels/__init__.py | 2 + .../testmodel_modelica_comparison_summary.txt | 4 + .../testmodel/testmodel_primary_series.csv | 202 ++++++ PythonModels/components/__init__.py | 2 + PythonModels/components/cylinder.py | 55 ++ PythonModels/components/orifice.py | 28 + PythonModels/components/pipe.py | 133 ++++ PythonModels/components/tank.py | 55 ++ PythonModels/components/tee.py | 172 +++++ PythonModels/core/__init__.py | 2 + PythonModels/core/base.py | 48 ++ PythonModels/core/medium.py | 96 +++ PythonModels/core/network.py | 78 ++ PythonModels/core/ports.py | 13 + PythonModels/core/solver.py | 102 +++ PythonModels/core/state.py | 21 + PythonModels/reporting/__init__.py | 23 + PythonModels/reporting/testmodel_outputs.py | 393 +++++++++++ PythonModels/scripts/run_testmodel.py | 219 ++++++ PythonModels/systems/__init__.py | 2 + PythonModels/systems/testmodel.py | 303 ++++++++ PythonModels/systems/testmodel_closure.py | 668 ++++++++++++++++++ 24 files changed, 2963 insertions(+) create mode 100644 PythonModels/.gitignore create mode 100644 PythonModels/README.md create mode 100644 PythonModels/__init__.py create mode 100644 PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt create mode 100644 PythonModels/baselines/testmodel/testmodel_primary_series.csv create mode 100644 PythonModels/components/__init__.py create mode 100644 PythonModels/components/cylinder.py create mode 100644 PythonModels/components/orifice.py create mode 100644 PythonModels/components/pipe.py create mode 100644 PythonModels/components/tank.py create mode 100644 PythonModels/components/tee.py create mode 100644 PythonModels/core/__init__.py create mode 100644 PythonModels/core/base.py create mode 100644 PythonModels/core/medium.py create mode 100644 PythonModels/core/network.py create mode 100644 PythonModels/core/ports.py create mode 100644 PythonModels/core/solver.py create mode 100644 PythonModels/core/state.py create mode 100644 PythonModels/reporting/__init__.py create mode 100644 PythonModels/reporting/testmodel_outputs.py create mode 100644 PythonModels/scripts/run_testmodel.py create mode 100644 PythonModels/systems/__init__.py create mode 100644 PythonModels/systems/testmodel.py create mode 100644 PythonModels/systems/testmodel_closure.py diff --git a/PythonModels/.gitignore b/PythonModels/.gitignore new file mode 100644 index 0000000..7a60b85 --- /dev/null +++ b/PythonModels/.gitignore @@ -0,0 +1,2 @@ +__pycache__/ +*.pyc diff --git a/PythonModels/README.md b/PythonModels/README.md new file mode 100644 index 0000000..35ccea7 --- /dev/null +++ b/PythonModels/README.md @@ -0,0 +1,340 @@ +# PythonModels + +`PythonModels` 用于承接 `ModelicaModels` 的 Python 平台移植。 + +目标不是把 `.mo` 文件逐行翻译成 Python，而是建立一个可运行、可对比、可逐步逼近 `OpenModelica` 行为的 Python 仿真框架。 + +当前状态不是“只有骨架”，而是“`Testmodel` 已有一版可运行的 ODE 近似实现，并具备基础结果导出与对比能力”。 + +## 当前目录 + +- `core/`: 通用基础设施 + 包含组件基类、状态与端口数据结构、介质模型、网络装配、积分入口。 +- `components/`: 元件级 Python 实现 + 目前有 `Cylinder`、`Tank`、`Pipe`、`Orifice`、`Tee` 五类元件。 +- `systems/`: 系统级装配与闭合 + 当前只有 `TestModelSystem`，对应 `ModelicaModels/Testmodel.mo`。 +- `reporting/`: 结果导出与对比 + 当前承接主变量 CSV、温度 CSV/SVG、Python 对 OpenModelica 的对比表与误差摘要导出。 +- `scripts/`: 运行脚本 + 当前入口是 `run_testmodel.py`。 +- `baselines/`: 提交进仓库的稳定基线 + 当前承接 Python 主变量基线和 Python 对 Modelica 的误差摘要基线。 +- `runs/`: 每次实际运行的默认输出目录 + 当前脚本默认会在这里创建带时间戳的子目录，用来放这次运行生成的产物。 + +当前关键文件： + +- `core/medium.py`: 理想气体近似介质 `IdealGasMedium` +- `core/medium.py`: 温度相关的空气近似介质 `IdealGasMedium` +- `core/network.py`: `SimulationNetwork`，负责组件注册、连接拓扑和状态向量拼装 +- `core/solver.py`: `integrate_ode()`，优先走 `SciPy solve_ivp`，缺依赖时回退到内置 RK4，并支持 `t_start == t_stop` 的零时长返回 +- `components/pipe.py`: 单阻容管道近似，入口压降 + 出口直连内容腔 +- `components/tee.py`: 三通的最小 stream 混合 helper +- `systems/testmodel.py`: `Testmodel` 的系统装配壳与外部运行入口 +- `systems/testmodel_closure.py`: `Testmodel` 当前专用的闭合、初始化投影、分支求解与端口回写 +- `reporting/testmodel_outputs.py`: `Testmodel` 的 CSV/SVG/对比摘要导出 +- `scripts/run_testmodel.py`: 基线运行与程序化执行入口 +- `tests/test_pythonmodels_regression.py`: 当前 Python 基线回归测试 + +## 当前阶段进度 + +这一阶段原先有 4 件重点工作，现在的状态如下： + +1. `mytee1` 的 stream/焓传播语义：已完成当前阶段收紧 + 现在如果只有一条支路发生倒流，下游来流焓统一按 `tank.h` 处理，不再临时借另一条支路的焓来凑。 +2. 下游初始化/约束处理：已完成当前阶段收口 + 之前是“直接改对象状态再开始积分”，现在已经收成显式的 `consistent_initial_state_vector()` 初始化入口。当前这一步会在不改下游总质量、总内能的前提下，把几段直接相连的体积拉回同一个连接压力。 +3. 自动校验：已完成当前阶段首版 + 已经补了标准库 `unittest` 回归测试，先把初始化投影是否守恒、是否污染原始状态，以及 4 个主变量的提交基线锁住。 +4. 更严格介质模型：已完成当前阶段首版 + 已经从固定 `cp/cv` 的理想气体近似，推进到随温度变化的空气近似，并接上了内能反解和初始化求根。 + +如果只看结果，可以把这一阶段理解成： + +- 连接器语义：首轮收紧已完成 +- 初始化入口：首轮收口已完成 +- 基线验证：首轮保护已完成 +- 介质精化：首轮近似已完成 + +## 当前阶段收口 + +上一轮 `N0-N3` 已全部完成首版，当前可以简单理解为： + +1. `N0`：系统层里最明显的流向/焓判断已经继续下沉到组件 helper。 +2. `N1`：模型参数和运行参数已经收口到配置对象。 +3. `N2`：运行接口已经分成“准备请求”和“执行请求”两层。 +4. `N3`：结果导出和命令行报告格式化已经统一收口到 `reporting/`。 + +这一轮结束后，项目已经不缺“能不能跑”的能力，下一步更重要的是把后续开发最容易卡住的地方先处理掉。 + +## 本次推送更新 + +本次推送已经把上一轮建议里的 `M2-M5` 推进到下面这个状态： + +1. `M2`：已完成当前阶段首版 + - 已把 `Testmodel` 的专用闭合、初始化投影、分支入口流量求解、下游支路出口流量闭合、端口状态回写，从 `systems/testmodel.py` 拆到新的 `systems/testmodel_closure.py` + - `TestModelSystem` 现在主要承担组件装配、网络注册和对闭合器的委托，不再继续堆积系统级手写细节 + +2. `M3`：已完成当前阶段首版 + - 已给两条支路入口流量固定点求解、下游公共压力投影补了显式诊断 + - 诊断内容至少包含 `converged / iterations / residual` + - 已支持严格模式；内部求解不收敛时可以直接抛错，而不是静默返回最后一个近似值 + - `run_testmodel()` 的结构化结果和 `testmodel_run_report.txt` 已能带出最后一次内部闭合求解诊断 + +3. `M4`：已完成当前阶段首版 + - 自动测试已不再只盯最终主变量结果 + - 现在已经覆盖： + - 改支路参数后，初始支路入口流量是否按预期变化 + - 更偏激配置下，初始化和内部闭合是否仍然收敛 + - 有无 Modelica 参考两种运行路径下，程序接口与产物行为是否一致 + +4. `M5`：已启动 + - 当前已经明确选择优先走“更容易扩展”的方向，而不是先追求更贴近 Modelica + - 已完成第一步：把闭合器内部原来大量写死的 `upper/lower` 双支路逻辑，收成可复用的 `BranchClosureComponents / BranchClosureState` 结构 + - 当前已继续推进到 `G1-G5` 的首轮兼容层改造：`snapshot` 已提供通用分支集合，系统层结果生成已拆成“通用键生成 + 旧键别名派生”两层，报告层已开始优先消费通用分支键，旧导出列名仍通过兼容映射保留，兼容测试已显式保护分支顺序和旧导出语义 + +## 下一阶段接手建议 + +如果继续往前推进，建议按下面顺序做，而不是再零散补功能： + +1. `G1`：已完成当前阶段首轮兼容接入 + - `TestModelSnapshot` 已新增 `branches` 集合 + - 每个分支当前至少带 `name / pipe / inlet_flow / outlet_flow / inlet_h / inlet_flow_diagnostics` + - `pipe_upper / pipe_lower / branch_inlet_flows / branch_outlet_flows` 目前仍保留为兼容属性，供旧调用方继续使用 + +2. `G2`：已完成当前阶段首轮内部迁移 + - `evaluate_solution()` 已改成从 `snapshot.branches` 读取数据，再通过显式分支名映射写回当前旧列名 + - `rhs()` 里的分支导数计算已改成通过通用 helper 按分支循环生成，再按当前状态向量顺序拼回 + - 当前外部导出列名仍保持兼容： + - `mypipe.p` + - `mypipe1.p` + - `branch_upper.in/out` + - `branch_lower.in/out` + +3. `G3`：已完成当前阶段首轮兼容测试 + - 当前测试已经显式保护： + - `branches` 顺序是否稳定 + - `snapshot` 新字段和兼容字段是否一致 + - 旧导出列名是否仍映射到正确分支语义 + - 参数变化后 `upper/lower` 的名字和顺序是否不会被打乱 + +4. `G4`：已完成当前阶段首轮兼容拆层 + - `evaluate_solution()` 现在会同时产出： + - 通用分支键：`branch..p/in/out` + - 旧兼容键：`mypipe.p`、`mypipe1.p`、`branch_upper.*`、`branch_lower.*` + - 报告层当前已开始优先读取通用分支键，旧键只作为兼容后备 + - 当前已经把“内部统一表达”和“旧接口兼容导出”拆成两层，但还没有把所有报告/导出逻辑都迁干净 + +5. `G5`：已完成当前阶段首轮兼容收口 + - `evaluate_solution()` 当前会先生成通用分支键，再统一派生旧兼容键 + - 报告层当前已支持“通用键优先、旧键兼容后备” + - 当前已经把系统层和 reporting 层的主要旧专名读取入口收口到少量 helper 上，后续继续迁移不会再到处散改 + +6. `P1`：下一阶段建议从这里接手 + 当前更合适的下一步，不是继续深挖内核通用化，而是切回结果导向主线： + - 定义一份稳定的外部输入参数 schema + - 明确这些结构化参数如何映射到 `TestModelConfig / TestModelRunConfig` + - 建立“结构化参数 -> 仿真执行 -> 结果产物/摘要”的稳定接口 + 这样可以直接服务后续文档解析、网页入口和报告生成，而不是继续在 `Testmodel` 内部做边际收益越来越低的抽象整理 + +7. `P2`：在 `P1` 完成后，再推进文档解析或报告生成链路 + 更现实的顺序应是： + - 先把结构化输入跑通 + - 再把结果摘要/产物组织成更接近最终产品的输出包 + - 最后再接 Word 解析或页面入口 + +如果后续继续推进，这个 README 也要一起更新，不要长期保留已经失效的路线描述。 + +## 当前实现了什么 + +当前代码已经实现： + +1. `m`、`U` 作为动态元件主状态，`p`、`T`、`rho`、`u`、`h` 作为派生量。 +2. `Cylinder`、`Tank`、`Pipe` 的刚性绝热容腔近似。 +3. `Orifice` 的压差开方流量关系。 +4. `Tee` 的简化混合焓处理。 +5. `Testmodel` 的系统级拓扑映射和一版可运行的 `rhs(t, x)`。 +6. 基于 `solve_ivp` 的积分入口，以及 SciPy 不可用时的 RK4 回退。 +7. 温度相关空气近似介质，包括 `cp(T)`、`h(T)`、`u(T)` 以及 `u -> T` 反解。 +8. 显式一致初值入口 `consistent_initial_state_vector()`，以及可迭代初始化器 `initialize_consistent_state()`。 +9. Python 主变量结果导出： + `mytank.p`、`mytank.T`、`mycylinder.p`、`mycylinder.T` +10. 贮箱温度曲线导出： + `testmodel_tank_temperature.csv` + `testmodel_tank_temperature.svg` +11. 基于 `ModelicaModels/Simulation/Testmodel_res.csv` 的逐时刻对比与误差摘要导出。 +12. 基于 `unittest` 的自动回归测试，当前已覆盖初始化守恒、主变量基线、运行接口、内部闭合诊断、通用分支兼容层、通用结果键与旧键别名一致性，以及部分中间闭合过程行为。 + +当前没有实现： + +- 通用 DAE 初始化器 +- `Modelica.Media.Air.SimpleAir` 的严格复刻 +- 面向任意拓扑的通用 connector/stream 求解器 + +## 当前怎么运行 + +最小运行方式： + +```bash +python3 -m PythonModels.scripts.run_testmodel +``` + +如果要改模型参数或运行参数，建议直接改配置对象，而不是改源码里的默认值。例如： + +```python +from PythonModels.core.solver import SolveIVPConfig +from PythonModels.scripts.run_testmodel import ( + TestModelRunConfig, + TestModelSamplingConfig, + run_testmodel, +) +from PythonModels.systems.testmodel import ( + BranchConfig, + CylinderConfig, + OrificeConfig, + PipeConfig, + TankConfig, + TestModelConfig, +) + +run_config = TestModelRunConfig( + model=TestModelConfig( + cylinder=CylinderConfig(p0=30e6), + upper_branch=BranchConfig( + orifice=OrificeConfig(K=8e-6), + pipe=PipeConfig(length=6.0, diameter=0.03), + ), + tank=TankConfig(volume=0.12), + ), + solver=SolveIVPConfig(t_start=0.0, t_stop=10.0, method="BDF"), + sampling=TestModelSamplingConfig(step=0.05), +) + +result = run_testmodel(run_config=run_config) +``` + +如果调用方想先确认“这次运行最后到底会用哪些路径、哪些采样点”，可以先准备请求，再执行： + +```python +from PythonModels.scripts.run_testmodel import ( + prepare_testmodel_run, + run_prepared_testmodel, + TestModelRunConfig, +) + +prepared = prepare_testmodel_run(run_config=TestModelRunConfig()) +print(prepared.output_dir) +print(prepared.t_eval) + +result = run_prepared_testmodel(prepared) +print(result.artifacts.primary_csv_path) +print(result.used_modelica_reference) +``` + +当前脚本会： + +1. 构建 `TestModelSystem` +2. 打印原始初值向量与约束一致后的初值向量 +3. 运行 `0 s -> 20 s` 的仿真，默认采样间隔 `0.1 s` +4. 将结果写入 `PythonModels/runs/` 下本次运行专属的时间戳目录 +5. 若存在 `ModelicaModels/Simulation/Testmodel_res.csv`，自动生成 Python 与 OpenModelica 对比结果 + +当前脚本默认不会再把运行结果直接写到提交基线目录，而是会在 `PythonModels/runs/` 下创建一个带时间戳的子目录，例如： + +- `PythonModels/runs/testmodel_20260512_103000_123456/` + +该目录里通常会包含： + +- `testmodel_primary_series.csv` +- `testmodel_tank_temperature.csv` +- `testmodel_tank_temperature.svg` +- `testmodel_run_report.txt` +- `testmodel_modelica_comparison.csv` +- `testmodel_modelica_comparison_summary.txt` + +## 基线结果 + +当前基线对比摘要来自： +[testmodel_modelica_comparison_summary.txt](/home/lujz/projects/pressurization-transfer-system/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt) + +当前四个主变量的最大误差为： + +- `mytank.p`: `max_abs_error = 134.960857 Pa`, `max_rel_error = 0.006798%` +- `mytank.T`: `max_abs_error = 0.035507 K`, `max_rel_error = 0.009016%` +- `mycylinder.p`: `max_abs_error = 1391.986349 Pa`, `max_rel_error = 0.009447%` +- `mycylinder.T`: `max_abs_error = 0.009069 K`, `max_rel_error = 0.003870%` + +这说明在当前基线工况下，Python 版主变量已经能较好贴近 OpenModelica 结果。 + +## 当前架构判断 + +如果按“组件正确 -> 网络闭合 -> 积分可跑 -> 结果对齐 -> 去近似”来看，当前大致处于： + +- 组件级：已完成首版 +- 系统闭合：已完成首版 +- 积分入口：已完成首版 +- 基线结果对齐：已具备初步能力 +- 去近似：仍在进行中 + +所以当前最准确的说法不是“已完成移植”，而是： + +`Testmodel` 已有一版可运行、可导出、可对比的 Python 近似实现。 + +## 已知限制 + +当前最主要的限制可以直接理解成下面几条： + +- 介质模型已从常 `cp/cv` 推进到温度相关空气近似，但仍不是 `Modelica.Media.Air.SimpleAir` 的严格复刻。 +- 系统整体仍是 ODE 化近似，不是原始 Modelica DAE 的直接复现。 +- `mytee1 -> mytank` 这一段虽然已经去掉早期的“虚拟出口导通系数”，改成了基于压力一致性的下游能量闭合，但本质上仍是工程近似。 +- 当前 `Tee` 的 stream 语义只覆盖了当前 `Testmodel` 需要的最小集合，还不是通用的 `inStream/actualStream` 框架。 +- 当前一致初值仍是 ODE 入口处的约束投影，不等同于真正的 DAE 初始化求解。 +- 当前自动校验主要锁的是 Python 提交基线，还不是稳定的 Modelica 阈值回归。 +- 当前闭合器、系统层和 reporting 层虽然已经开始做“双支路结构化”，但对外结果序列、报告字段和部分导出命名仍然保留 `Testmodel` 专名兼容层，还没有完全转成通用表达。 +- 当前内核已经足够支撑下一阶段“结构化参数 -> 仿真执行 -> 产物输出”的链路开发，但还没有现成的 Word 参数解析入口和正式报告生成链路。 + +所以，当前版本适合： + +- 架构验证 +- 组件接口验证 +- 基线工况对比 +- 结果导出与误差定位 + +但当前版本还不适合： + +- 直接宣称与 OpenModelica 严格等价 +- 作为最终工程结论的唯一依据 +- 直接扩展到更复杂拓扑而不补通用连接器语义 + +## 文件级现状 + +按代码现状逐项看： + +- `core/base.py`: 正常 + 只提供最小抽象层，没有明显冗余。 +- `core/ports.py`: 正常 + `PortState` 目前只保留 `p`、`m_flow`、`h_outflow` 三个必要字段。 +- `core/state.py`: 正常 + `VolumeState` 只负责 `[m, U]` 状态打包。 +- `core/network.py`: 正常 + 负责状态向量拼装和连接摘要，不参与物理求解。 +- `core/solver.py`: 正常 + 已支持 SciPy、RK4 回退和零时长仿真。 +- `components/*.py`: 正常 + 都是当前一版近似模型，没有发现与 README 明显冲突的“未记录能力”。 +- `systems/testmodel.py`: 是当前最重要的技术债集中区 + 这里承载了下游流向切换、焓混合、压力投影等近似逻辑，后续演进应主要落在这里。 +- `scripts/run_testmodel.py`: 正常 + 已不是“最小打印脚本”，而是当前结果导出和对比入口。 +- `baselines/`: 是当前稳定基线，不应该随着日常运行频繁改动。 +- `runs/`: 是当前默认运行产物目录，不是手写源代码，也不应该当作提交基线使用。 + +## 当前主技术债 + +目前最主要的技术债，可以直接理解成下面 4 件事： + +1. 当前初始化虽然已经引入迭代诊断，但本质上仍是 ODE 入口近似，不是真正的 DAE 初始化器。 +2. `systems/testmodel.py` 还是承载了太多系统级闭合和初始化逻辑，只是主要端口的手写 stream 方向判断已经搬到组件 helper 里了，装配参数本身已经基本收口到配置对象。 +3. 自动校验现在主要锁的是 Python 这一版自己的基线，还不是稳定的 Modelica 阈值回归。 +4. 当前空气物性已经完成首轮基线校准，但还不是 `SimpleAir` 的严格复刻。以后如果换工况，或者拿到更多 Modelica 原始结果，参数大概率还要继续调。 diff --git a/PythonModels/__init__.py b/PythonModels/__init__.py new file mode 100644 index 0000000..a392a55 --- /dev/null +++ b/PythonModels/__init__.py @@ -0,0 +1,2 @@ +"""Python port scaffold for the Modelica-based pressurization system.""" + diff --git a/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt b/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt new file mode 100644 index 0000000..b8835bd --- /dev/null +++ b/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt @@ -0,0 +1,4 @@ +mytank.p: max_abs_error=134.960858, max_rel_error=0.006798% +mytank.T: max_abs_error=0.035507, max_rel_error=0.009016% +mycylinder.p: max_abs_error=1391.986349, max_rel_error=0.009447% +mycylinder.T: max_abs_error=0.009069, max_rel_error=0.003870% diff --git a/PythonModels/baselines/testmodel/testmodel_primary_series.csv b/PythonModels/baselines/testmodel/testmodel_primary_series.csv new file mode 100644 index 0000000..f908440 --- /dev/null +++ b/PythonModels/baselines/testmodel/testmodel_primary_series.csv @@ -0,0 +1,202 @@ +time_s,mytank.p,mytank.T,mycylinder.p,mycylinder.T +0.0,99999.9999998181,299.9999999994543,35000000.0,300.0 +0.1,113767.94796194455,310.9637987101763,34857995.191961475,299.65190130981887 +0.2,127490.75529043324,320.1741038867186,34716455.97160761,299.3039342009557 +0.30000000000000004,141168.22832317703,327.9540327439692,34575384.336400226,298.9561065982938 +0.4,154800.2272891848,334.6028385113769,34434781.727959625,298.60842486666263 +0.5,168386.6833862148,340.3440160419025,34294648.85592305,298.26089368263627 +0.6000000000000001,181927.56886228317,345.34597934661645,34154986.006528884,297.9135169662631 +0.7000000000000001,195422.8867284584,349.7379443903134,34015793.14872051,297.5662980942661 +0.8,208872.67226727083,353.62071952544704,33877069.91858829,297.2192381374397 +0.9,222277.06297793236,357.0742683513725,33738814.897943415,296.8723276411647 +1.0,235635.98150118813,360.16219760966084,33601028.88468161,296.5255848498291 +1.1,248949.42783703818,362.9362624751423,33463711.878802836,296.1790180498007 +1.2000000000000002,262217.75812283135,365.43943715392686,33326860.20704978,295.8325994898035 +1.3,275440.80324992316,367.7064864598734,33190475.613635924,295.48635453607636 +1.4000000000000001,288618.55460473493,369.766650344722,33054558.187403098,295.14029246566224 +1.5,301751.3171281971,371.6448887937757,32919104.783141974,294.79439468638526 +1.6,314839.14341954247,373.36206261131525,32784114.858335685,294.44866485903736 +1.7000000000000002,327882.06907680153,374.93594217338176,32649588.04582047,294.10310819940065 +1.8,340880.20077679906,376.3818124575356,32515523.2453146,293.7577259967503 +1.9000000000000001,353833.6270518691,377.7128659999238,32381919.54368144,293.41251959012544 +2.0,366742.42291137256,378.94054912194076,32248776.16726251,293.06749036423287 +2.1,379606.6730674497,380.0748368029417,32116092.24232331,292.72263968200787 +2.2,392426.50707225554,381.12448095768656,31983866.432642274,292.37796667994184 +2.3000000000000003,405202.0671563814,382.0971825271288,31852097.271230437,292.03346983037886 +2.4000000000000004,417933.4155500976,382.9997115435689,31720784.11623488,291.6891515780112 +2.5,430620.6510589185,383.838081019396,31589925.948559783,291.3450125479824 +2.6,443263.883710829,384.61764092425494,31459521.633358993,291.0010529534012 +2.7,455863.2542884696,385.3431706778388,31329569.71857793,290.6572718762472 +2.8000000000000003,468418.84583467664,386.0189391423481,31200069.34769955,290.31367046787994 +2.9000000000000004,480930.76401047717,386.64878296917834,31071019.430919208,289.9702490545257 +3.0,493399.1147767602,387.2361564731079,30942418.8753394,289.627007973335 +3.1,505824.0120818939,387.7841699997202,30814266.50567852,289.283948167411 +3.2,518205.5585633492,388.29565044444143,30686561.263317347,288.9410696919777 +3.3000000000000003,530543.860941083,388.7731587646163,30559302.04752924,288.5983729053581 +3.4000000000000004,542839.025935053,389.21902406251183,30432487.75758757,288.2558581528899 +3.5,555091.1532282452,389.63535331051384,30306117.365346134,287.9135275289302 +3.6,567300.3523937837,390.0240928732544,30180189.74065028,287.5713806480376 +3.7,579466.730902255,390.3870196149787,30054703.77503127,287.2294176336731 +3.8000000000000003,591590.3951155862,390.7257605811333,29929658.371455234,286.88763887333175 +3.9000000000000004,603671.4513957056,391.0418105423231,29805052.432888325,286.5460447413264 +4.0,615710.0061045411,391.3365451815722,29680884.862296656,286.20463559857706 +4.1000000000000005,627706.1656040212,391.61123273483645,29557154.562646367,285.8634117923989 +4.2,639660.0356829006,391.86704273384044,29433860.4428154,285.5223738511199 +4.3,651571.716232906,392.1050410519721,29311001.472504564,285.1815242995298 +4.4,663441.3179459961,392.32624321858174,29188576.510019373,284.840861980108 +4.5,675268.9468771729,392.53157932311984,29066584.46149164,284.5003872946073 +4.6000000000000005,687054.7090814385,392.72191331429804,28945024.2330532,284.16010063101663 +4.7,698798.7106137947,392.8980487868277,28823894.730835862,283.8200023633402 +4.800000000000001,710501.057529244,393.06073417111315,28703194.86097143,283.48009285137505 +4.9,722161.8558827877,393.21066739655134,28582923.529591747,283.1403724404863 +5.0,733781.2088290841,393.3484951327475,28463079.672743235,282.8008422184954 +5.1000000000000005,745359.2175579917,393.47481852250473,28343662.24673756,282.4615037792963 +5.2,756895.9922909128,393.5902150621596,28224670.114733644,282.1223563449153 +5.300000000000001,768391.6386869826,393.6952183545109,28106102.186946325,281.78340031124986 +5.4,779846.2624053361,393.7903288312906,27987957.373590477,281.44463605959197 +5.5,791259.96910511,393.8760162992802,27870234.58488091,281.10606395639473 +5.6000000000000005,802632.8644454395,393.9527222554955,27752932.731032487,280.76768435303745 +5.7,813965.0540854601,394.02086199551263,27636050.72226007,280.4294975855874 +5.800000000000001,825256.6344365194,394.08081904264304,27519587.564161647,280.0915054160418 +5.9,836507.7093909897,394.13296077650165,27403542.18517475,279.75370843660113 +6.0,847718.3891245491,394.17763867488134,27287913.44892988,279.4161062418491 +6.1000000000000005,858888.7786008355,394.2151805066602,27172700.272815377,279.0786992215644 +6.2,870018.9827834871,394.2458959595685,27057901.57421954,278.7414877501898 +6.300000000000001,881109.1066361419,394.27007787647636,26943516.2705307,278.40447218658807 +6.4,892159.2551224378,394.2880033916485,26829543.279137183,278.0676528737956 +6.5,903169.5332060129,394.2999349762812,26715981.5174273,277.73103013877227 +6.6000000000000005,914139.9829839376,394.30609016406663,26602830.551205155,277.3946116451124 +6.7,925070.7482002805,394.30672506364067,26490088.897871558,277.0583932025842 +6.800000000000001,935961.9434794089,394.30206871636136,26377755.375172496,276.7223739769823 +6.9,946813.6729747849,394.2923336718555,26265828.9088526,276.3865543393284 +7.0,957626.040839869,394.2777219546483,26154308.424656466,276.05093464451556 +7.1000000000000005,968399.1512281233,394.2584257130283,26043192.84832872,275.7155152310525 +7.2,979133.1082930088,394.23462782052076,25932481.10561397,275.3802964208052 +7.300000000000001,989828.0161879869,394.2065024339623,25822172.12225682,275.04527851873627 +7.4,1000483.9790665191,394.17421551178546,25712264.824001882,274.71046181264023 +7.5,1011101.1010820667,394.13792529578427,25602758.136593778,274.3758465728772 +7.6000000000000005,1021679.4863880915,394.09778275932865,25493650.985777102,274.04143305210266 +7.7,1032219.239138054,394.0539320247227,25384942.297296483,273.707221484995 +7.800000000000001,1042720.4634854163,394.0065107521557,25276630.996896524,273.37321208797937 +7.9,1053183.2635836392,393.95565050247825,25168716.01032184,273.0394050589494 +8.0,1063607.7435861847,393.9014770758345,25061196.263317037,272.705800576985 +8.1,1073994.0076465139,393.8441108280058,24954070.681626726,272.3723988020679 +8.200000000000001,1084342.0075784405,393.78360959211744,24847339.76225091,272.03921728278755 +8.3,1094651.950514763,393.7201237653114,24741001.368788917,271.7062444960846 +8.4,1104923.9566483083,393.65376548392067,24635054.261552785,271.3734787513138 +8.5,1115158.1282653515,393.5846357175065,24529497.385545585,271.0409203804046 +8.6,1125354.567652169,393.5128313057523,24424329.685770374,270.7085696976142 +8.700000000000001,1135513.3770950353,393.43844517133164,24319550.10723021,270.37642699925055 +8.8,1145634.6588802272,393.3615665197664,24215157.594928175,270.0444925633909 +8.9,1155718.5152940191,393.2822810271899,24111151.093867306,269.712766649598 +9.0,1165765.0486226876,393.2006710168648,24007529.54905069,269.3812494986329 +9.1,1175774.361152508,393.11681562522904,23904291.90548136,269.04994133216485 +9.200000000000001,1185746.555169756,393.03079095819226,23801437.108162407,268.7188423524781 +9.3,1195681.732960707,392.94267023834027,23698964.10209688,268.3879527421748 +9.4,1205579.9968116365,392.8525239436615,23596871.83228785,268.05727266387584 +9.5,1215441.4490088206,392.7604199383583,23495159.243738372,267.72680225991735 +9.600000000000001,1225266.1918385345,392.66642359626576,23393825.281451505,267.39654165204473 +9.700000000000001,1235054.3275870536,392.57059791736026,23292868.89043032,267.0664909411029 +9.8,1244805.7424247153,392.47293675670085,23192291.244732097,266.7366751290975 +9.9,1254520.6963966596,392.3735487251935,23092089.662209835,266.40707625875064 +10.0,1264199.3049931854,392.2724947546991,22992262.951678198,266.07769289844276 +10.100000000000001,1273841.6679760527,392.1698289153906,22892810.0841785,265.74852540655036 +10.200000000000001,1283447.8851070204,392.0656034102887,22793730.030752078,265.41957412301247 +10.3,1293018.0561478487,391.9598686580462,22695021.76244025,265.0908393690368 +10.4,1302552.2808602974,391.8526733713811,22596684.250284337,264.76232144680336 +10.5,1312050.6590061258,391.7440646314241,22498716.465325654,264.4340206391643 +10.600000000000001,1321513.290347094,391.63408795822716,22401117.37860553,264.10593720934105 +10.700000000000001,1330940.2746449616,391.5227873776613,22303885.961165283,263.7780714006174 +10.8,1340331.711661488,391.41020548491883,22207021.184046242,263.4504234360302 +10.9,1349687.7011584332,391.29638350481804,22110522.01828973,263.1229935180553 +11.0,1359008.342897557,391.18136134909884,22014387.434937052,262.79578182829164 +11.100000000000001,1368293.7366406189,391.0651776708804,21918616.405029543,262.4687885271409 +11.200000000000001,1377543.9821493784,390.9478699164463,21823207.89960853,262.1420137534841 +11.3,1386759.1791855951,390.82947437450605,21728160.88971532,261.81545762435434 +11.4,1395939.2416852904,390.70997702085253,21633476.26302751,261.4891421257888 +11.5,1405084.4058089943,390.5894488404711,21539151.583747786,261.1630514204039 +11.600000000000001,1414194.7727817593,390.46792355037144,21445185.807824824,260.83718544108797 +11.700000000000001,1423270.4378832965,390.3454323238288,21351577.952528503,260.5115448028764 +11.8,1432311.4962235806,390.22200536933536,21258327.036879413,260.1861301247683 +11.9,1441318.042742851,390.0976719679259,21165432.0816488,259.86094202976824 +12.0,1450290.1722116095,389.97246050877925,21072892.1093586,259.5359811449295 +12.100000000000001,1459227.9792306225,389.84639852319333,20980706.14428146,259.2112481013979 +12.200000000000001,1468131.5582309205,389.71951271701596,20888873.212440677,258.88674353445447 +12.3,1477001.0034737969,389.59182900161403,20797392.34161027,258.5624680835606 +12.4,1485836.4090508097,389.46337252345796,20706262.561314918,258.23842239240224 +12.5,1494637.8688837802,389.3341676923915,20615482.902829994,257.9146071089349 +12.600000000000001,1503405.4767247946,389.20423820865625,20525052.399181567,257.5910228854296 +12.700000000000001,1512139.3261562001,389.07360708873284,20434970.085146382,257.26767037851835 +12.8,1520839.5105906113,388.94229669006046,20345234.997251876,256.9445502492407 +12.9,1529506.1233409408,388.81032876227204,20255846.173053782,256.62166315000906 +13.0,1538139.2687433015,388.67772882293264,20166802.53641284,256.2990076411302 +13.100000000000001,1546739.030251722,388.5445136335602,20078103.22657088,255.97658615427704 +13.200000000000001,1555305.5007360894,388.4107031827765,19989747.285653118,255.65439933992917 +13.3,1563838.7729005008,388.27631690939495,19901733.757494725,255.332447852161 +13.4,1572338.9392832702,388.14137372147235,19814061.687640797,255.0107323486784 +13.5,1580806.0922569225,388.0058920145713,19726730.123346385,254.68925349085504 +13.600000000000001,1589240.3240281995,387.86988968927403,19639738.113576483,254.36801194376957 +13.700000000000001,1597641.7266380545,387.73338416798015,19553084.709006038,254.04700837624307 +13.8,1606010.3919616563,387.59639241102605,19466768.962019928,253.72624346087662 +13.9,1614346.4117083861,387.458930932153,19380789.926712975,253.40571787408948 +14.0,1622649.8774218408,387.32101581335974,19295146.658889957,253.0854322961577 +14.100000000000001,1630920.8804798292,387.18266271916264,19209838.21606559,252.76538741125262 +14.200000000000001,1639159.5120943757,387.0438869102949,19124863.657464534,252.4455839074805 +14.3,1647365.8633117182,386.90470325686863,19040222.0440214,252.12602247692183 +14.4,1655540.0250123062,386.76512625102305,18955912.438380726,251.8067038156715 +14.5,1663682.087910807,386.6251700190859,18871933.904897016,251.48762862387906 +14.600000000000001,1671792.1425560997,386.4848483332645,18788285.509634707,251.16879760578948 +14.700000000000001,1679870.2793312764,386.3441746228917,18704966.320368182,250.85021146978437 +14.8,1687916.588453644,386.2031619852454,18621975.406581767,250.5318709284234 +14.9,1695931.1599747243,386.06182319595837,18539311.83946974,250.21377669848638 +15.0,1703914.0837802505,385.920170719039,18456974.691936314,249.8959295010154 +15.100000000000001,1711865.4494154856,385.77821670970445,18374963.040397402,249.57833006564294 +15.200000000000001,1719785.3379250832,385.6359727064743,18293276.048945524,249.26097933285033 +15.3,1727673.8459201432,385.4934505781353,18211912.721118417,248.94387786218365 +15.4,1735531.0625176337,385.3506616315768,18130872.137749474,248.6270263958575 +15.5,1743357.0766656764,385.2076169067951,18050153.381413613,248.31042568075102 +15.600000000000001,1751151.9771435438,385.06432718490083,17969755.53642727,247.99407646845765 +15.700000000000001,1758915.8525616606,384.9208029958382,17889677.688848406,247.67797951533504 +15.8,1766648.7913616048,384.77705462583015,17809918.926476505,247.36213558255577 +15.9,1774350.881816107,384.63309212455715,17730478.33885257,247.04654543615806 +16.0,1782022.2120290494,384.4889253120841,17651355.017259125,246.73120984709743 +16.1,1789662.8699354655,384.3445637855432,17572548.05472022,246.41612959129878 +16.2,1797272.9433015438,384.2000169255848,17494056.546001427,246.10130544970855 +16.3,1804852.5197246224,384.05529390260284,17415879.58760983,245.78673820834797 +16.400000000000002,1812401.6866331936,383.9104036827472,17338016.277794052,245.47242865836654 +16.5,1819920.5312869006,383.76535503372884,17260465.716544226,245.15837759609565 +16.6,1827409.1407765402,383.6201565304272,17183227.005592,244.84458582310356 +16.7,1834867.6020240602,383.47481656030743,17106299.24841057,244.53105414625009 +16.8,1842296.0017825624,383.32934332865534,17029681.550214626,244.2177833777421 +16.900000000000002,1849694.4266362996,383.1837448636361,16953373.017960392,243.9047743351897 +17.0,1857062.9630006768,383.03802902118434,16877372.76034562,243.59202784166294 +17.1,1864401.697122252,382.8922034897328,16801679.887809563,243.27954472574828 +17.2,1871710.7150787353,382.7462757947841,16726293.512533028,242.96732582160686 +17.3,1878990.1030524147,382.60025337999576,16651212.745618157,242.65537190477806 +17.400000000000002,1886239.954455942,382.4541455990818,16576436.623591991,242.343682012199 +17.5,1893460.3482686526,382.3079576538116,16501964.331849081,242.03225853722532 +17.6,1900651.3699906773,382.16169648335216,16427794.988527464,241.72110231102408 +17.7,1907813.104956132,382.0153688850827,16353927.713477474,241.41021416968596 +17.8,1914945.6383331183,381.86898151834066,16280361.628261749,241.09959495427398 +17.900000000000002,1922049.0551237254,381.7225409080495,16207095.856155202,240.78924551087363 +18.0,1929123.4401640275,381.57605344823065,16134129.522145053,240.47916669064318 +18.1,1936168.8781240864,381.4295254054061,16061461.752930798,240.16935934986387 +18.2,1943185.4535079484,381.28296292189384,15989091.67692425,239.8598243499916 +18.3,1950173.2506536485,381.1363720190009,15917018.424249483,239.55056255770802 +18.400000000000002,1957132.3537332045,380.98975860011734,15845241.126742886,239.24157484497272 +18.5,1964062.8467526236,380.84312845371574,15773758.917953137,238.93286208907583 +18.6,1970964.8135518986,380.6964872562572,15702570.933141202,238.6244251726908 +18.7,1977838.3378050062,380.54984057500997,15631676.309280336,238.31626498392774 +18.8,1984683.5030199126,380.4031938707821,15561074.185056096,238.00838241638746 +18.900000000000002,1991500.3925385692,380.25655250057133,15490763.700866321,237.7007783692156 +19.0,1998289.0895369116,380.10992172013556,15420743.998821149,237.39345374715774 +19.1,2005049.677024864,379.96330668648653,15351014.222743012,237.08640946061425 +19.200000000000003,2011782.237846337,379.81671246030885,15281573.518166626,236.77964642569634 +19.3,2018486.8546792252,379.67014400830794,15212421.032339014,236.47316556428217 +19.400000000000002,2025163.610035411,379.5236062054879,15143555.914219463,236.16696780407335 +19.5,2031812.5863479478,379.3771038539995,15074977.31358037,235.861054060844 +19.6,2038433.8680519294,379.2306420805889,15006684.359544702,235.55542481191287 +19.700000000000003,2045027.536485746,379.084225365391,14938676.213175347,235.25008113730536 +19.8,2051593.6739729291,378.9378583245209,14870952.025374293,234.9450238874766 +19.900000000000002,2058132.3627231135,378.79154550141357,14803510.948218277,234.64025390703284 +20.0,2064643.6848320398,378.64529136859073,14736352.134958768,234.33577203462033 diff --git a/PythonModels/components/__init__.py b/PythonModels/components/__init__.py new file mode 100644 index 0000000..d8e493c --- /dev/null +++ b/PythonModels/components/__init__.py @@ -0,0 +1,2 @@ +"""Component implementations for the Python system model.""" + diff --git a/PythonModels/components/cylinder.py b/PythonModels/components/cylinder.py new file mode 100644 index 0000000..066fb27 --- /dev/null +++ b/PythonModels/components/cylinder.py @@ -0,0 +1,55 @@ +from __future__ import annotations + +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +class Cylinder(DynamicComponent): + """Python port of ModelicaModels.Mycylinder.""" + + def __init__( + self, + name: str, + medium: IdealGasMedium, + V: float = 0.01, + p0: float = 35e6, + T0: float = 300.0, + ) -> None: + super().__init__(name=name) + self.medium = medium + self.V = V + m0 = p0 * V / (medium.R_gas * T0) + U0 = m0 * medium.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_b = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + props = self.medium.properties_from_mU(self.state.m, self.state.U, self.V) + self.port_b.p = props.p + self.port_b.h_outflow = props.h + return props + + def derivatives_from_connection( + self, + *, + connected_h: float, + port_m_flow: float, + internal_h: float, + ) -> VolumeState: + inlet_h = self.connection_inlet_enthalpy( + port_m_flow=port_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + return self.derivatives(inlet_h, port_m_flow) + + def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState: + return VolumeState(m=m_flow, U=m_flow * inlet_h) diff --git a/PythonModels/components/orifice.py b/PythonModels/components/orifice.py new file mode 100644 index 0000000..9da6282 --- /dev/null +++ b/PythonModels/components/orifice.py @@ -0,0 +1,28 @@ +from __future__ import annotations + +from math import sqrt + +from PythonModels.core.base import AlgebraicComponent +from PythonModels.core.ports import PortState + + +class Orifice(AlgebraicComponent): + """Python port of ModelicaModels.Myorifice.""" + + def __init__(self, name: str, opening: float = 1.0, K: float = 1e-7) -> None: + super().__init__(name=name) + self.opening = opening + self.K = K + self.port_a = PortState() + self.port_b = PortState() + + @property + def K_eff(self) -> float: + return self.K * max(self.opening, 0.001) + + def mass_flow(self, p_a: float, p_b: float) -> float: + dp = p_a - p_b + if dp == 0.0: + return 0.0 + return self.K_eff * sqrt(abs(dp)) * (1.0 if dp > 0.0 else -1.0) + diff --git a/PythonModels/components/pipe.py b/PythonModels/components/pipe.py new file mode 100644 index 0000000..6fb9fe5 --- /dev/null +++ b/PythonModels/components/pipe.py @@ -0,0 +1,133 @@ +from __future__ import annotations + +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +class Pipe(DynamicComponent): + """Python port of ModelicaModels.Mypipe.""" + + def __init__( + self, + name: str, + medium: IdealGasMedium, + L: float = 5.0, + D: float = 0.02, + lambda_darcy: float = 0.02, + p0: float = 1e5, + T0: float = 300.0, + ) -> None: + super().__init__(name=name) + self.medium = medium + self.L = L + self.D = D + self.lambda_darcy = lambda_darcy + self.area = 3.141592653589793 * D * D / 4.0 + self.V = self.area * L + m0 = p0 * self.V / (medium.R_gas * T0) + U0 = m0 * medium.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_a = PortState() + self.port_b = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + props = self.medium.properties_from_mU(self.state.m, self.state.U, self.V) + self.port_b.p = props.p + self.port_a.h_outflow = props.h + self.port_b.h_outflow = props.h + return props + + def inlet_pressure(self, m_flow_a: float, rho: float, core_pressure: float) -> float: + resistance = self.lambda_darcy * (self.L / self.D) + dynamic_term = m_flow_a * abs(m_flow_a) / (2.0 * rho * self.area * self.area) + return core_pressure + resistance * dynamic_term + + def port_a_inlet_enthalpy( + self, + *, + port_a_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + return self.connection_inlet_enthalpy( + port_m_flow=port_a_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + + def port_b_inlet_enthalpy( + self, + *, + port_b_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + return self.connection_inlet_enthalpy( + port_m_flow=port_b_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + + def connection_inlet_enthalpies( + self, + *, + port_a_m_flow: float, + connected_h_a: float, + port_b_m_flow: float, + connected_h_b: float, + internal_h: float, + ) -> tuple[float, float]: + return ( + self.port_a_inlet_enthalpy( + port_a_m_flow=port_a_m_flow, + connected_h=connected_h_a, + internal_h=internal_h, + ), + self.port_b_inlet_enthalpy( + port_b_m_flow=port_b_m_flow, + connected_h=connected_h_b, + internal_h=internal_h, + ), + ) + + def derivatives_from_connections( + self, + *, + port_a_m_flow: float, + connected_h_a: float, + port_b_m_flow: float, + connected_h_b: float, + internal_h: float, + ) -> VolumeState: + inlet_h_a, inlet_h_b = self.connection_inlet_enthalpies( + port_a_m_flow=port_a_m_flow, + connected_h_a=connected_h_a, + port_b_m_flow=port_b_m_flow, + connected_h_b=connected_h_b, + internal_h=internal_h, + ) + return self.derivatives( + inlet_h_a=inlet_h_a, + inlet_h_b=inlet_h_b, + m_flow_a=port_a_m_flow, + m_flow_b=port_b_m_flow, + ) + + def derivatives( + self, + inlet_h_a: float, + inlet_h_b: float, + m_flow_a: float, + m_flow_b: float, + ) -> VolumeState: + dm_dt = m_flow_a + m_flow_b + dU_dt = m_flow_a * inlet_h_a + m_flow_b * inlet_h_b + return VolumeState(m=dm_dt, U=dU_dt) diff --git a/PythonModels/components/tank.py b/PythonModels/components/tank.py new file mode 100644 index 0000000..8f51d69 --- /dev/null +++ b/PythonModels/components/tank.py @@ -0,0 +1,55 @@ +from __future__ import annotations + +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +class Tank(DynamicComponent): + """Python port of ModelicaModels.Mytank.""" + + def __init__( + self, + name: str, + medium: IdealGasMedium, + V: float = 0.1, + p0: float = 1e5, + T0: float = 300.0, + ) -> None: + super().__init__(name=name) + self.medium = medium + self.V = V + m0 = p0 * V / (medium.R_gas * T0) + U0 = m0 * medium.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_a = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + props = self.medium.properties_from_mU(self.state.m, self.state.U, self.V) + self.port_a.p = props.p + self.port_a.h_outflow = props.h + return props + + def derivatives_from_connection( + self, + *, + connected_h: float, + port_m_flow: float, + internal_h: float, + ) -> VolumeState: + inlet_h = self.connection_inlet_enthalpy( + port_m_flow=port_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + return self.derivatives(inlet_h, port_m_flow) + + def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState: + return VolumeState(m=m_flow, U=m_flow * inlet_h) diff --git a/PythonModels/components/tee.py b/PythonModels/components/tee.py new file mode 100644 index 0000000..6602435 --- /dev/null +++ b/PythonModels/components/tee.py @@ -0,0 +1,172 @@ +from __future__ import annotations + +from PythonModels.core.base import AlgebraicComponent +from PythonModels.core.ports import PortState + + +class Tee(AlgebraicComponent): + """Python port of ModelicaModels.Mytee.""" + + def __init__(self, name: str) -> None: + super().__init__(name=name) + self.port_in = PortState() + self.port_out1 = PortState() + self.port_out2 = PortState() + + def mixed_inlet_enthalpy( + self, + branch1_m_flow: float, + branch1_h: float, + branch2_m_flow: float, + branch2_h: float, + fallback_h: float = 0.0, + ) -> float: + positive_1 = max(branch1_m_flow, 0.0) + positive_2 = max(branch2_m_flow, 0.0) + total = positive_1 + positive_2 + if total <= 1e-9: + return fallback_h + return (positive_1 * branch1_h + positive_2 * branch2_h) / total + + def inlet_stream_enthalpy( + self, + branch1_m_flow: float, + branch1_h: float, + branch2_m_flow: float, + branch2_h: float, + fallback_h: float, + ) -> float: + """Approximate `inStream(port_in.h_outflow)` for the current tee topology.""" + + return self.mixed_inlet_enthalpy( + branch1_m_flow, + branch1_h, + branch2_m_flow, + branch2_h, + fallback_h=fallback_h, + ) + + def branch_actual_stream_enthalpy( + self, + branch_m_flow: float, + branch_h: float, + inlet_h: float, + ) -> float: + """Approximate `actualStream(branch.h_outflow)` for a tee branch port.""" + + return inlet_h if branch_m_flow > 0.0 else branch_h + + @staticmethod + def _solve_linear_2x2( + a11: float, + a12: float, + a21: float, + a22: float, + b1: float, + b2: float, + ) -> tuple[float, float] | None: + determinant = a11 * a22 - a12 * a21 + if abs(determinant) <= 1e-12: + return None + x1 = (b1 * a22 - b2 * a12) / determinant + x2 = (a11 * b2 - a21 * b1) / determinant + return x1, x2 + + def solve_branch_outlet_flows_from_energy_balance( + self, + *, + ratio_branch1: float, + ratio_branch2: float, + inlet_h_branch1: float, + inlet_h_branch2: float, + branch1_h: float, + branch2_h: float, + inlet_h: float, + q_in_branch1: float, + q_in_branch2: float, + tolerance: float = 1e-12, + ) -> tuple[float, float]: + """Solve branch outlet flows for the current three-port downstream tee use-case.""" + + rhs_branch1 = q_in_branch1 * inlet_h_branch1 + rhs_branch2 = q_in_branch2 * inlet_h_branch2 + + def solve_both_forward() -> tuple[float, float] | None: + return self._solve_linear_2x2( + (1.0 + ratio_branch1) * branch1_h, + ratio_branch1 * branch2_h, + ratio_branch2 * branch1_h, + (1.0 + ratio_branch2) * branch2_h, + rhs_branch1, + rhs_branch2, + ) + + def solve_one_reverse( + *, + branch1_reverse: bool, + ) -> tuple[float, float] | None: + if branch1_reverse: + return self._solve_linear_2x2( + inlet_h * (1.0 + ratio_branch1), + ratio_branch1 * inlet_h, + ratio_branch2 * inlet_h, + branch2_h + ratio_branch2 * inlet_h, + rhs_branch1, + rhs_branch2, + ) + + return self._solve_linear_2x2( + branch1_h + ratio_branch1 * inlet_h, + ratio_branch1 * inlet_h, + ratio_branch2 * inlet_h, + inlet_h * (1.0 + ratio_branch2), + rhs_branch1, + rhs_branch2, + ) + + def solve_both_reverse() -> tuple[float, float] | None: + return self._solve_linear_2x2( + inlet_h * (1.0 + ratio_branch1), + ratio_branch1 * inlet_h, + ratio_branch2 * inlet_h, + inlet_h * (1.0 + ratio_branch2), + rhs_branch1, + rhs_branch2, + ) + + candidate_solvers = ( + ( + solve_both_forward, + lambda q1, q2: q1 >= -tolerance and q2 >= -tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=True), + lambda q1, q2: q1 < -tolerance and q2 >= -tolerance and q1 + q2 > tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=True), + lambda q1, q2: q1 < -tolerance and q2 >= -tolerance and q1 + q2 <= tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=False), + lambda q1, q2: q2 < -tolerance and q1 >= -tolerance and q1 + q2 > tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=False), + lambda q1, q2: q2 < -tolerance and q1 >= -tolerance and q1 + q2 <= tolerance, + ), + ( + solve_both_reverse, + lambda q1, q2: q1 < -tolerance and q2 < -tolerance, + ), + ) + + for solver, predicate in candidate_solvers: + candidate = solver() + if candidate is None: + continue + q_out_branch1, q_out_branch2 = candidate + if predicate(q_out_branch1, q_out_branch2): + return q_out_branch1, q_out_branch2 + + return solve_both_forward() or (0.0, 0.0) diff --git a/PythonModels/core/__init__.py b/PythonModels/core/__init__.py new file mode 100644 index 0000000..ae12d41 --- /dev/null +++ b/PythonModels/core/__init__.py @@ -0,0 +1,2 @@ +"""Core abstractions for the Python system model.""" + diff --git a/PythonModels/core/base.py b/PythonModels/core/base.py new file mode 100644 index 0000000..bac666e --- /dev/null +++ b/PythonModels/core/base.py @@ -0,0 +1,48 @@ +from __future__ import annotations + +from abc import ABC, abstractmethod + + +class Component(ABC): + def __init__(self, name: str) -> None: + self.name = name + + +class DynamicComponent(Component): + state_size = 2 + + @staticmethod + def actual_stream_enthalpy( + port_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + """Approximate `actualStream(port.h_outflow)` for a mixed control volume port.""" + + return connected_h if port_m_flow > 0.0 else internal_h + + def connection_inlet_enthalpy( + self, + port_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + """Resolve the enthalpy convected into this control volume through one port.""" + + return self.actual_stream_enthalpy( + port_m_flow=port_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + + @abstractmethod + def get_state_vector(self) -> list[float]: + raise NotImplementedError + + @abstractmethod + def set_state_vector(self, values: list[float]) -> None: + raise NotImplementedError + + +class AlgebraicComponent(Component): + """Stateless element described by algebraic constraints only.""" diff --git a/PythonModels/core/medium.py b/PythonModels/core/medium.py new file mode 100644 index 0000000..3410590 --- /dev/null +++ b/PythonModels/core/medium.py @@ -0,0 +1,96 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass(frozen=True) +class ThermodynamicProperties: + p: float + T: float + rho: float + u: float + h: float + + +@dataclass(frozen=True) +class IdealGasMedium: + """Temperature-dependent ideal-gas air approximation. + + This is still not a strict clone of `Modelica.Media.Air.SimpleAir`. + The small linear `cp(T)` term is kept configurable for calibration, but the + current default is calibrated against the committed Testmodel baseline and + therefore falls back to the constant-heat-capacity limit. + """ + + name: str = "SimpleAirApprox" + R_gas: float = 287.0 + cp_ref: float = 1005.0 + T_ref: float = 300.0 + cp_slope: float = 0.0 + + @property + def cv(self) -> float: + return self.cv_at_temperature(self.T_ref) + + @property + def gamma(self) -> float: + return self.cp_at_temperature(self.T_ref) / self.cv + + def cp_at_temperature(self, T: float) -> float: + return self.cp_ref + self.cp_slope * (T - self.T_ref) + + def cv_at_temperature(self, T: float) -> float: + return self.cp_at_temperature(T) - self.R_gas + + def density(self, p: float, T: float) -> float: + return p / (self.R_gas * T) + + def specific_internal_energy(self, T: float) -> float: + delta_T = T - self.T_ref + return ( + self.cv * self.T_ref + + self.cv * delta_T + + 0.5 * self.cp_slope * delta_T * delta_T + ) + + def specific_enthalpy(self, T: float) -> float: + delta_T = T - self.T_ref + return ( + self.cp_ref * self.T_ref + + self.cp_ref * delta_T + + 0.5 * self.cp_slope * delta_T * delta_T + ) + + def temperature_from_internal_energy(self, u: float) -> float: + reference_internal_energy = self.cv * self.T_ref + delta_u = u - reference_internal_energy + + if abs(self.cp_slope) <= 1e-15: + return self.T_ref + delta_u / self.cv + + a = 0.5 * self.cp_slope + b = self.cv + c = -delta_u + discriminant = max(b * b - 4.0 * a * c, 0.0) + positive_root = (-b + discriminant**0.5) / (2.0 * a) + negative_root = (-b - discriminant**0.5) / (2.0 * a) + delta_T = positive_root if abs(positive_root) <= abs(negative_root) else negative_root + return self.T_ref + delta_T + + def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: + if m <= 0.0: + raise ValueError("Mass must stay positive when recovering temperature.") + return self.temperature_from_internal_energy(U / m) + + def pressure(self, m: float, T: float, V: float) -> float: + if V <= 0.0: + raise ValueError("Volume must stay positive.") + return m * self.R_gas * T / V + + def properties_from_mU(self, m: float, U: float, V: float) -> ThermodynamicProperties: + T = self.temperature_from_mass_internal_energy(m, U) + p = self.pressure(m, T, V) + rho = m / V + u = U / m + h = self.specific_enthalpy(T) + return ThermodynamicProperties(p=p, T=T, rho=rho, u=u, h=h) diff --git a/PythonModels/core/network.py b/PythonModels/core/network.py new file mode 100644 index 0000000..e912867 --- /dev/null +++ b/PythonModels/core/network.py @@ -0,0 +1,78 @@ +from __future__ import annotations + +from dataclasses import dataclass + +from PythonModels.core.base import Component, DynamicComponent + + +@dataclass(frozen=True) +class Connection: + source_component: str + source_port: str + target_component: str + target_port: str + + +class SimulationNetwork: + """Container for components, topology, and state-vector bookkeeping.""" + + def __init__(self, name: str) -> None: + self.name = name + self.components: dict[str, Component] = {} + self.connections: list[Connection] = [] + + def add_component(self, component: Component) -> None: + if component.name in self.components: + raise ValueError(f"Duplicate component name: {component.name}") + self.components[component.name] = component + + def connect( + self, + source_component: str, + source_port: str, + target_component: str, + target_port: str, + ) -> None: + self.connections.append( + Connection( + source_component=source_component, + source_port=source_port, + target_component=target_component, + target_port=target_port, + ) + ) + + def dynamic_components(self) -> list[DynamicComponent]: + return [ + component + for component in self.components.values() + if isinstance(component, DynamicComponent) + ] + + def initial_state_vector(self) -> list[float]: + values: list[float] = [] + for component in self.dynamic_components(): + values.extend(component.get_state_vector()) + return values + + def apply_state_vector(self, values: list[float]) -> None: + cursor = 0 + for component in self.dynamic_components(): + next_cursor = cursor + component.state_size + component.set_state_vector(values[cursor:next_cursor]) + cursor = next_cursor + if cursor != len(values): + raise ValueError("State vector length does not match dynamic components.") + + def summary(self) -> str: + lines = [f"Network: {self.name}", "Components:"] + for name, component in self.components.items(): + lines.append(f" - {name}: {component.__class__.__name__}") + lines.append("Connections:") + for conn in self.connections: + lines.append( + f" - {conn.source_component}.{conn.source_port}" + f" -> {conn.target_component}.{conn.target_port}" + ) + return "\n".join(lines) + diff --git a/PythonModels/core/ports.py b/PythonModels/core/ports.py new file mode 100644 index 0000000..d340db2 --- /dev/null +++ b/PythonModels/core/ports.py @@ -0,0 +1,13 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass +class PortState: + """Python-side analogue of a Modelica fluid port.""" + + p: float = 0.0 + m_flow: float = 0.0 + h_outflow: float = 0.0 + diff --git a/PythonModels/core/solver.py b/PythonModels/core/solver.py new file mode 100644 index 0000000..3e101e1 --- /dev/null +++ b/PythonModels/core/solver.py @@ -0,0 +1,102 @@ +from __future__ import annotations + +from dataclasses import dataclass +from typing import Callable + + +@dataclass(frozen=True) +class SolveIVPConfig: + t_start: float = 0.0 + t_stop: float = 20.0 + method: str = "BDF" + rtol: float = 1e-6 + atol: float = 1e-8 + max_step: float = 1e-3 + + +@dataclass(frozen=True) +class ODESolution: + t: list[float] + y: list[list[float]] + success: bool + message: str + + +def _vector_add(a: list[float], b: list[float], scale: float = 1.0) -> list[float]: + return [x + scale * y for x, y in zip(a, b)] + + +def _runge_kutta_4( + rhs: Callable[[float, list[float]], list[float]], + initial_state: list[float], + config: SolveIVPConfig, + t_eval: list[float] | None, +) -> ODESolution: + if t_eval is None: + point_count = max( + 2, + int((config.t_stop - config.t_start) / max(config.max_step, 1e-6)) + 1, + ) + step = (config.t_stop - config.t_start) / (point_count - 1) + t_eval = [config.t_start + index * step for index in range(point_count)] + + state = list(initial_state) + states = [[value] for value in state] + times = [float(t_eval[0])] + current_time = float(t_eval[0]) + + for target_time in t_eval[1:]: + while current_time < target_time - 1e-15: + dt = min(config.max_step, target_time - current_time) + k1 = rhs(current_time, state) + k2 = rhs(current_time + 0.5 * dt, _vector_add(state, k1, 0.5 * dt)) + k3 = rhs(current_time + 0.5 * dt, _vector_add(state, k2, 0.5 * dt)) + k4 = rhs(current_time + dt, _vector_add(state, k3, dt)) + state = [ + value + (dt / 6.0) * (a + 2.0 * b + 2.0 * c + d) + for value, a, b, c, d in zip(state, k1, k2, k3, k4) + ] + current_time += dt + + times.append(float(target_time)) + for index, value in enumerate(state): + states[index].append(value) + + return ODESolution( + t=times, + y=states, + success=True, + message="Integrated with built-in RK4 fallback because SciPy is unavailable.", + ) + + +def integrate_ode( + rhs: Callable[[float, list[float]], list[float]], + initial_state: list[float], + config: SolveIVPConfig, + t_eval: list[float] | None = None, +): + """Thin wrapper around scipy.integrate.solve_ivp with a pure-Python fallback.""" + + if abs(config.t_stop - config.t_start) <= 1e-15: + return ODESolution( + t=[float(config.t_start)], + y=[[value] for value in initial_state], + success=True, + message="Skipped integration because t_start equals t_stop.", + ) + + try: + from scipy.integrate import solve_ivp + except ImportError: + return _runge_kutta_4(rhs, initial_state, config, t_eval) + + return solve_ivp( + fun=rhs, + t_span=(config.t_start, config.t_stop), + y0=initial_state, + method=config.method, + rtol=config.rtol, + atol=config.atol, + t_eval=t_eval, + ) diff --git a/PythonModels/core/state.py b/PythonModels/core/state.py new file mode 100644 index 0000000..e4cf3b3 --- /dev/null +++ b/PythonModels/core/state.py @@ -0,0 +1,21 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass +class VolumeState: + """Primary dynamic state for rigid adiabatic control volumes.""" + + m: float + U: float + + def as_vector(self) -> list[float]: + return [self.m, self.U] + + @classmethod + def from_vector(cls, values: list[float]) -> "VolumeState": + if len(values) != 2: + raise ValueError("VolumeState requires exactly two values: [m, U].") + return cls(m=values[0], U=values[1]) + diff --git a/PythonModels/reporting/__init__.py b/PythonModels/reporting/__init__.py new file mode 100644 index 0000000..9ca6880 --- /dev/null +++ b/PythonModels/reporting/__init__.py @@ -0,0 +1,23 @@ +from PythonModels.reporting.testmodel_outputs import ( + COMPARISON_KEYS, + MODELICA_COMPARISON_COLUMNS, + PRIMARY_KEYS, + TestModelArtifacts, + export_testmodel_artifacts, + format_testmodel_run_report, + load_modelica_series, + write_testmodel_run_report, + write_modelica_comparison, +) + +__all__ = [ + "COMPARISON_KEYS", + "MODELICA_COMPARISON_COLUMNS", + "PRIMARY_KEYS", + "TestModelArtifacts", + "export_testmodel_artifacts", + "format_testmodel_run_report", + "load_modelica_series", + "write_testmodel_run_report", + "write_modelica_comparison", +] diff --git a/PythonModels/reporting/testmodel_outputs.py b/PythonModels/reporting/testmodel_outputs.py new file mode 100644 index 0000000..9338b1c --- /dev/null +++ b/PythonModels/reporting/testmodel_outputs.py @@ -0,0 +1,393 @@ +from __future__ import annotations + +from bisect import bisect_left +import csv +from dataclasses import dataclass +from pathlib import Path +from typing import Any + + +PRIMARY_KEYS = ( + "mytank.p", + "mytank.T", + "mycylinder.p", + "mycylinder.T", +) + +MODELICA_COMPARISON_COLUMNS = { + "mytank.p": "mytank.p", + "mytank.T": "mytank.T", + "mycylinder.p": "mycylinder.p", + "mycylinder.T": "mycylinder.T", + "branch.upper_branch.p": "mypipe.p", + "branch.upper_branch.in": "myorifice.port_a.m_flow", + "branch.upper_branch.out": "mytee1.port_out2.m_flow", + "branch.lower_branch.p": "mypipe1.p", + "branch.lower_branch.in": "myorifice1.port_a.m_flow", + "branch.lower_branch.out": "mytee1.port_out1.m_flow", +} + +COMPARISON_KEYS = tuple(MODELICA_COMPARISON_COLUMNS.keys()) + + +def _branch_series_values( + series: dict[str, list[float]], + branch_name: str, + legacy_key: str, +) -> list[float]: + generic_key = f"branch.{branch_name}.{legacy_key.split('.')[-1]}" + if generic_key in series: + return series[generic_key] + return series[legacy_key] + + +@dataclass(frozen=True) +class TestModelArtifacts: + primary_csv_path: Path + temperature_csv_path: Path + temperature_svg_path: Path + run_report_path: Path + comparison_csv_path: Path | None = None + comparison_summary_path: Path | None = None + + +def format_testmodel_run_report( + *, + network_summary: str, + initialization: Any, + raw_initial_state: tuple[float, ...], + consistent_initial_state: tuple[float, ...], + solution: Any, + series: dict[str, list[float]], + solve_diagnostics: Any, + artifacts: TestModelArtifacts, + comparison_summary: dict[str, tuple[float, float]] | None, +) -> str: + lines = [ + network_summary, + "", + f"Initialization converged: {initialization.converged}", + f"Initialization iterations: {initialization.iterations}", + f"Initialization max state delta: {initialization.max_state_delta:.6e}", + f"Initialization max flow delta: {initialization.max_flow_delta:.6e}", + f"Initialization max enthalpy delta: {initialization.max_enthalpy_delta:.6e}", + ( + "Initialization downstream pressure spread: " + f"{initialization.downstream_pressure_spread:.6e}" + ), + "", + "Raw initial state vector:", + str(list(raw_initial_state)), + "", + "Constraint-consistent initial state vector:", + str(list(consistent_initial_state)), + "", + f"Solver success: {solution.success}", + f"Solver message: {solution.message}", + f"Final time: {solution.t[-1]:.2f} s", + f"Final tank pressure: {series['mytank.p'][-1]:.3f} Pa", + f"Final tank temperature: {series['mytank.T'][-1]:.3f} K", + f"Final cylinder pressure: {series['mycylinder.p'][-1]:.3f} Pa", + ( + "Final branch inflow: " + f"{_branch_series_values(series, 'upper_branch', 'branch_upper.in')[-1] + _branch_series_values(series, 'lower_branch', 'branch_lower.in')[-1]:.6f} kg/s" + ), + ] + + if solve_diagnostics is not None: + lines.extend( + [ + "", + "Final closure solve diagnostics:", + ( + "Upper branch inlet solve: " + f"converged={solve_diagnostics.upper_branch_inlet.converged}, " + f"iterations={solve_diagnostics.upper_branch_inlet.iterations}, " + f"residual={solve_diagnostics.upper_branch_inlet.residual:.6e}" + ), + ( + "Lower branch inlet solve: " + f"converged={solve_diagnostics.lower_branch_inlet.converged}, " + f"iterations={solve_diagnostics.lower_branch_inlet.iterations}, " + f"residual={solve_diagnostics.lower_branch_inlet.residual:.6e}" + ), + ] + ) + if solve_diagnostics.downstream_pressure_projection is not None: + lines.append( + "Downstream pressure projection: " + f"converged={solve_diagnostics.downstream_pressure_projection.converged}, " + f"iterations={solve_diagnostics.downstream_pressure_projection.iterations}, " + f"residual={solve_diagnostics.downstream_pressure_projection.residual:.6e}" + ) + + lines.extend( + [ + f"Primary series CSV: {artifacts.primary_csv_path}", + f"Temperature CSV: {artifacts.temperature_csv_path}", + f"Temperature plot: {artifacts.temperature_svg_path}", + f"Run report TXT: {artifacts.run_report_path}", + ] + ) + + if ( + artifacts.comparison_csv_path is not None + and artifacts.comparison_summary_path is not None + ): + lines.extend( + [ + f"Modelica comparison CSV: {artifacts.comparison_csv_path}", + f"Modelica comparison summary: {artifacts.comparison_summary_path}", + ] + ) + if comparison_summary is not None: + for key, (max_abs_error, max_rel_error) in comparison_summary.items(): + lines.append( + f"{key} max abs error: {max_abs_error:.6f}, " + f"max rel error: {max_rel_error:.6%}" + ) + + return "\n".join(lines) + "\n" + + +def write_testmodel_run_report(output_dir: Path, report_text: str) -> Path: + report_path = output_dir / "testmodel_run_report.txt" + report_path.write_text(report_text, encoding="utf-8") + return report_path + + +def _write_primary_series_csv(output_dir: Path, series: dict[str, list[float]]) -> Path: + csv_path = output_dir / "testmodel_primary_series.csv" + with csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle) + writer.writerow(["time_s", *PRIMARY_KEYS]) + for index, time_value in enumerate(series["time"]): + writer.writerow([time_value, *(series[key][index] for key in PRIMARY_KEYS)]) + return csv_path + + +def _write_temperature_csv(output_dir: Path, time_values: list[float], temperatures: list[float]) -> Path: + csv_path = output_dir / "testmodel_tank_temperature.csv" + with csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle) + writer.writerow(["time_s", "mytank_T_K"]) + writer.writerows(zip(time_values, temperatures)) + return csv_path + + +def _write_temperature_svg(output_dir: Path, time_values: list[float], temperatures: list[float]) -> Path: + svg_path = output_dir / "testmodel_tank_temperature.svg" + + width = 900 + height = 520 + left = 90 + right = 40 + top = 60 + bottom = 70 + plot_width = width - left - right + plot_height = height - top - bottom + + min_time = min(time_values) + max_time = max(time_values) + min_temp = min(temperatures) + max_temp = max(temperatures) + temp_padding = max(1.0, (max_temp - min_temp) * 0.08) + min_temp -= temp_padding + max_temp += temp_padding + + def scale_x(value: float) -> float: + return left + (value - min_time) / max(max_time - min_time, 1e-12) * plot_width + + def scale_y(value: float) -> float: + return top + (max_temp - value) / max(max_temp - min_temp, 1e-12) * plot_height + + points = " ".join( + f"{scale_x(time_value):.2f},{scale_y(temperature):.2f}" + for time_value, temperature in zip(time_values, temperatures) + ) + + x_ticks = 5 + y_ticks = 5 + x_tick_markup = [] + y_tick_markup = [] + + for index in range(x_ticks + 1): + fraction = index / x_ticks + time_value = min_time + fraction * (max_time - min_time) + x = left + fraction * plot_width + x_tick_markup.append( + f'' + ) + x_tick_markup.append( + f'' + f"{time_value:.1f}" + ) + + for index in range(y_ticks + 1): + fraction = index / y_ticks + temp_value = min_temp + fraction * (max_temp - min_temp) + y = top + plot_height - fraction * plot_height + y_tick_markup.append( + f'' + ) + y_tick_markup.append( + f'' + f"{temp_value:.1f}" + ) + + svg_content = f""" + + Python Testmodel Tank Temperature + Time (s) + Temperature (K) + + {''.join(x_tick_markup)} + {''.join(y_tick_markup)} + + +""" + svg_path.write_text(svg_content, encoding="utf-8") + return svg_path + + +def load_modelica_series(csv_path: Path, variable_names: tuple[str, ...]) -> dict[str, list[float]]: + series = {"time": []} + for variable_name in variable_names: + series[variable_name] = [] + + with csv_path.open("r", newline="", encoding="utf-8") as handle: + reader = csv.DictReader(handle) + available_variable_names = tuple( + variable_name + for variable_name in variable_names + if MODELICA_COMPARISON_COLUMNS.get(variable_name, variable_name) in (reader.fieldnames or ()) + ) + for row in reader: + series["time"].append(float(row["time"])) + for variable_name in available_variable_names: + modelica_column = MODELICA_COMPARISON_COLUMNS.get(variable_name, variable_name) + series[variable_name].append(float(row[modelica_column])) + + return series + + +def _interpolate_series_value(time_values: list[float], values: list[float], target_time: float) -> float: + if target_time <= time_values[0]: + return values[0] + if target_time >= time_values[-1]: + return values[-1] + + right_index = bisect_left(time_values, target_time) + if right_index < len(time_values) and abs(time_values[right_index] - target_time) <= 1e-12: + return values[right_index] + + left_index = right_index - 1 + left_time = time_values[left_index] + right_time = time_values[right_index] + fraction = (target_time - left_time) / (right_time - left_time) + return values[left_index] + fraction * (values[right_index] - values[left_index]) + + +def write_modelica_comparison( + output_dir: Path, + python_series: dict[str, list[float]], + modelica_series: dict[str, list[float]], +) -> tuple[Path, Path, dict[str, tuple[float, float]]]: + comparison_csv_path = output_dir / "testmodel_modelica_comparison.csv" + summary_path = output_dir / "testmodel_modelica_comparison_summary.txt" + summary: dict[str, tuple[float, float]] = {} + + with comparison_csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle) + header = ["time_s"] + comparison_keys = tuple( + key + for key in COMPARISON_KEYS + if key in python_series and key in modelica_series and modelica_series[key] + ) + for key in comparison_keys: + header.extend( + [ + f"python.{key}", + f"modelica.{key}", + f"abs_error.{key}", + f"rel_error.{key}", + ] + ) + writer.writerow(header) + + max_abs_errors = {key: 0.0 for key in comparison_keys} + max_rel_errors = {key: 0.0 for key in comparison_keys} + + for index, time_value in enumerate(python_series["time"]): + row = [time_value] + for key in comparison_keys: + python_value = python_series[key][index] + modelica_value = _interpolate_series_value( + modelica_series["time"], + modelica_series[key], + time_value, + ) + abs_error = abs(python_value - modelica_value) + rel_error = abs_error / max(abs(modelica_value), 1e-9) + max_abs_errors[key] = max(max_abs_errors[key], abs_error) + max_rel_errors[key] = max(max_rel_errors[key], rel_error) + row.extend([python_value, modelica_value, abs_error, rel_error]) + writer.writerow(row) + + summary_lines = [] + for key in comparison_keys: + summary[key] = (max_abs_errors[key], max_rel_errors[key]) + summary_lines.append( + f"{key}: max_abs_error={max_abs_errors[key]:.6f}, " + f"max_rel_error={max_rel_errors[key]:.6%}" + ) + summary_path.write_text("\n".join(summary_lines) + "\n", encoding="utf-8") + return comparison_csv_path, summary_path, summary + + +def export_testmodel_artifacts( + *, + output_dir: Path, + series: dict[str, list[float]], + modelica_series: dict[str, list[float]] | None = None, +) -> tuple[TestModelArtifacts, dict[str, tuple[float, float]] | None]: + output_dir.mkdir(parents=True, exist_ok=True) + + primary_csv_path = _write_primary_series_csv(output_dir, series) + temperature_csv_path = _write_temperature_csv( + output_dir, + series["time"], + series["mytank.T"], + ) + temperature_svg_path = _write_temperature_svg( + output_dir, + series["time"], + series["mytank.T"], + ) + + comparison_csv_path = None + comparison_summary_path = None + comparison_summary = None + if modelica_series is not None: + ( + comparison_csv_path, + comparison_summary_path, + comparison_summary, + ) = write_modelica_comparison(output_dir, series, modelica_series) + + return ( + TestModelArtifacts( + primary_csv_path=primary_csv_path, + temperature_csv_path=temperature_csv_path, + temperature_svg_path=temperature_svg_path, + run_report_path=output_dir / "testmodel_run_report.txt", + comparison_csv_path=comparison_csv_path, + comparison_summary_path=comparison_summary_path, + ), + comparison_summary, + ) diff --git a/PythonModels/scripts/run_testmodel.py b/PythonModels/scripts/run_testmodel.py new file mode 100644 index 0000000..4834677 --- /dev/null +++ b/PythonModels/scripts/run_testmodel.py @@ -0,0 +1,219 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from datetime import UTC, datetime +from pathlib import Path + +from PythonModels.reporting import ( + COMPARISON_KEYS, + PRIMARY_KEYS, + TestModelArtifacts, + export_testmodel_artifacts, + format_testmodel_run_report, + load_modelica_series, + write_testmodel_run_report, +) +from PythonModels.core.solver import SolveIVPConfig +from PythonModels.systems.testmodel import ( + InitializationDiagnostics, + TestModelConfig, + TestModelSystem, +) +from PythonModels.systems.testmodel_closure import TestModelSolveDiagnostics + + +@dataclass(frozen=True) +class TestModelSamplingConfig: + step: float = 0.1 + + +@dataclass(frozen=True) +class TestModelPathConfig: + output_dir: Path | None = None + modelica_result_path: Path | None = None + + +@dataclass(frozen=True) +class TestModelExecutionConfig: + use_modelica_reference_if_available: bool = True + + +@dataclass(frozen=True) +class TestModelRunConfig: + model: TestModelConfig = field(default_factory=TestModelConfig) + solver: SolveIVPConfig = field(default_factory=SolveIVPConfig) + sampling: TestModelSamplingConfig = field(default_factory=TestModelSamplingConfig) + paths: TestModelPathConfig = field(default_factory=TestModelPathConfig) + execution: TestModelExecutionConfig = field(default_factory=TestModelExecutionConfig) + + @property + def sample_step(self) -> float: + return self.sampling.step + + def sample_times(self) -> list[float]: + return _sample_times( + self.solver.t_start, + self.solver.t_stop, + step=self.sampling.step, + ) + + +@dataclass(frozen=True) +class PreparedTestModelRun: + run_config: TestModelRunConfig + repo_root: Path + output_dir: Path + modelica_result_path: Path + t_eval: tuple[float, ...] + use_modelica_reference_if_available: bool + modelica_reference_exists: bool + + +@dataclass(frozen=True) +class TestModelRunResult: + run_config: TestModelRunConfig + prepared_run: PreparedTestModelRun + system: TestModelSystem + initialization: InitializationDiagnostics + raw_initial_state: tuple[float, ...] + consistent_initial_state: tuple[float, ...] + solution: object + series: dict[str, list[float]] + solve_diagnostics: TestModelSolveDiagnostics | None + artifacts: TestModelArtifacts + comparison_summary: dict[str, tuple[float, float]] | None + used_modelica_reference: bool + + +def _sample_times(t_start: float, t_stop: float, step: float) -> list[float]: + point_count = int(round((t_stop - t_start) / step)) + return [t_start + index * step for index in range(point_count + 1)] + + +def _default_run_output_dir(pythonmodels_root: Path) -> Path: + timestamp = datetime.now(UTC).strftime("testmodel_%Y%m%d_%H%M%S_%f") + return pythonmodels_root / "runs" / timestamp + + +def prepare_testmodel_run( + *, + run_config: TestModelRunConfig | None = None, + output_dir: Path | None = None, + modelica_result_path: Path | None = None, +) -> PreparedTestModelRun: + run_config = run_config or TestModelRunConfig() + repo_root = Path(__file__).resolve().parents[2] + pythonmodels_root = Path(__file__).resolve().parents[1] + resolved_output_dir = ( + output_dir + or run_config.paths.output_dir + or _default_run_output_dir(pythonmodels_root) + ) + resolved_modelica_result_path = ( + modelica_result_path + or run_config.paths.modelica_result_path + or repo_root / "ModelicaModels" / "Simulation" / "Testmodel_res.csv" + ) + t_eval = tuple(run_config.sample_times()) + return PreparedTestModelRun( + run_config=run_config, + repo_root=repo_root, + output_dir=resolved_output_dir, + modelica_result_path=resolved_modelica_result_path, + t_eval=t_eval, + use_modelica_reference_if_available=run_config.execution.use_modelica_reference_if_available, + modelica_reference_exists=resolved_modelica_result_path.exists(), + ) + + +def run_prepared_testmodel(prepared_run: PreparedTestModelRun) -> TestModelRunResult: + run_config = prepared_run.run_config + system = TestModelSystem(config=run_config.model) + raw_initial_state = tuple(system.initial_state_vector()) + initialization = system.initialize_consistent_state() + consistent_initial_state = tuple(initialization.state_vector) + solution = system.simulate(config=run_config.solver, t_eval=list(prepared_run.t_eval)) + series = system.evaluate_solution(solution) + solve_diagnostics = system.last_solve_diagnostics + + modelica_series = None + used_modelica_reference = False + if ( + prepared_run.use_modelica_reference_if_available + and prepared_run.modelica_reference_exists + ): + modelica_series = load_modelica_series( + prepared_run.modelica_result_path, + COMPARISON_KEYS, + ) + used_modelica_reference = True + + artifacts, comparison_summary = export_testmodel_artifacts( + output_dir=prepared_run.output_dir, + series=series, + modelica_series=modelica_series, + ) + report_text = format_testmodel_run_report( + network_summary=system.network.summary(), + initialization=initialization, + raw_initial_state=raw_initial_state, + consistent_initial_state=consistent_initial_state, + solution=solution, + series=series, + solve_diagnostics=solve_diagnostics, + artifacts=artifacts, + comparison_summary=comparison_summary, + ) + write_testmodel_run_report(prepared_run.output_dir, report_text) + + return TestModelRunResult( + run_config=run_config, + prepared_run=prepared_run, + system=system, + initialization=initialization, + raw_initial_state=raw_initial_state, + consistent_initial_state=consistent_initial_state, + solution=solution, + series=series, + solve_diagnostics=solve_diagnostics, + artifacts=artifacts, + comparison_summary=comparison_summary, + used_modelica_reference=used_modelica_reference, + ) + + +def run_testmodel( + *, + run_config: TestModelRunConfig | None = None, + output_dir: Path | None = None, + modelica_result_path: Path | None = None, +) -> TestModelRunResult: + prepared_run = prepare_testmodel_run( + run_config=run_config, + output_dir=output_dir, + modelica_result_path=modelica_result_path, + ) + return run_prepared_testmodel(prepared_run) + + +def main() -> None: + run_config = TestModelRunConfig() + result = run_testmodel(run_config=run_config) + print( + format_testmodel_run_report( + network_summary=result.system.network.summary(), + initialization=result.initialization, + raw_initial_state=result.raw_initial_state, + consistent_initial_state=result.consistent_initial_state, + solution=result.solution, + series=result.series, + solve_diagnostics=result.solve_diagnostics, + artifacts=result.artifacts, + comparison_summary=result.comparison_summary, + ), + end="", + ) + + +if __name__ == "__main__": + main() diff --git a/PythonModels/systems/__init__.py b/PythonModels/systems/__init__.py new file mode 100644 index 0000000..187800c --- /dev/null +++ b/PythonModels/systems/__init__.py @@ -0,0 +1,2 @@ +"""System assembly modules.""" + diff --git a/PythonModels/systems/testmodel.py b/PythonModels/systems/testmodel.py new file mode 100644 index 0000000..2652292 --- /dev/null +++ b/PythonModels/systems/testmodel.py @@ -0,0 +1,303 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Any + +from PythonModels.components.cylinder import Cylinder +from PythonModels.components.orifice import Orifice +from PythonModels.components.pipe import Pipe +from PythonModels.components.tank import Tank +from PythonModels.components.tee import Tee +from PythonModels.core.medium import IdealGasMedium +from PythonModels.core.network import SimulationNetwork +from PythonModels.core.solver import SolveIVPConfig, integrate_ode +from PythonModels.systems.testmodel_closure import ( + BranchClosureComponents, + InitializationDiagnostics, + TestModelClosure, + TestModelClosureComponents, + TestModelSnapshot, +) + + +@dataclass(frozen=True) +class CylinderConfig: + volume: float = 0.01 + p0: float = 35e6 + T0: float = 300.0 + + +@dataclass(frozen=True) +class OrificeConfig: + K: float = 1e-5 + + +@dataclass(frozen=True) +class TankConfig: + volume: float = 0.1 + p0: float = 1e5 + T0: float = 300.0 + + +@dataclass(frozen=True) +class PipeConfig: + length: float = 5.0 + diameter: float = 0.02 + lambda_darcy: float = 0.02 + p0: float = 1e5 + T0: float = 300.0 + + +@dataclass(frozen=True) +class BranchConfig: + orifice: OrificeConfig = field(default_factory=OrificeConfig) + pipe: PipeConfig = field(default_factory=PipeConfig) + + +@dataclass(frozen=True) +class TestModelConfig: + cylinder: CylinderConfig = field(default_factory=CylinderConfig) + upper_branch: BranchConfig = field(default_factory=BranchConfig) + lower_branch: BranchConfig = field(default_factory=BranchConfig) + tank: TankConfig = field(default_factory=TankConfig) + + +class TestModelSystem: + """Runnable first-pass Python system for the current Testmodel topology. + + This version keeps the component split from the Modelica model while keeping + the downstream tee-tank pressure coupling in the ODE framework. The original + Modelica system is a tighter DAE because both pipe outlets discharge into an + ideal lossless junction directly connected to the tank. Here the branch + outlet flows are solved from a pressure-consistent energy balance so the + outlet is no longer driven by an arbitrary conductance parameter. + """ + + def __init__( + self, + medium: IdealGasMedium | None = None, + config: TestModelConfig | None = None, + ) -> None: + self.medium = medium or IdealGasMedium() + self.config = config or TestModelConfig() + + self.mycylinder = Cylinder( + name="mycylinder", + medium=self.medium, + V=self.config.cylinder.volume, + p0=self.config.cylinder.p0, + T0=self.config.cylinder.T0, + ) + self.mytee = Tee(name="mytee") + self.myorifice = Orifice(name="myorifice", K=self.config.upper_branch.orifice.K) + self.mypipe = Pipe( + name="mypipe", + medium=self.medium, + L=self.config.upper_branch.pipe.length, + D=self.config.upper_branch.pipe.diameter, + lambda_darcy=self.config.upper_branch.pipe.lambda_darcy, + p0=self.config.upper_branch.pipe.p0, + T0=self.config.upper_branch.pipe.T0, + ) + self.myorifice1 = Orifice(name="myorifice1", K=self.config.lower_branch.orifice.K) + self.mypipe1 = Pipe( + name="mypipe1", + medium=self.medium, + L=self.config.lower_branch.pipe.length, + D=self.config.lower_branch.pipe.diameter, + lambda_darcy=self.config.lower_branch.pipe.lambda_darcy, + p0=self.config.lower_branch.pipe.p0, + T0=self.config.lower_branch.pipe.T0, + ) + self.mytee1 = Tee(name="mytee1") + self.mytank = Tank( + name="mytank", + medium=self.medium, + V=self.config.tank.volume, + p0=self.config.tank.p0, + T0=self.config.tank.T0, + ) + + self.network = SimulationNetwork(name="Testmodel") + for component in ( + self.mycylinder, + self.mytee, + self.myorifice, + self.mypipe, + self.myorifice1, + self.mypipe1, + self.mytee1, + self.mytank, + ): + self.network.add_component(component) + + self.network.connect("mycylinder", "port_b", "mytee", "port_in") + self.network.connect("mytee", "port_out1", "myorifice", "port_a") + self.network.connect("myorifice", "port_b", "mypipe", "port_a") + self.network.connect("mypipe", "port_b", "mytee1", "port_out2") + self.network.connect("mytee", "port_out2", "myorifice1", "port_a") + self.network.connect("myorifice1", "port_b", "mypipe1", "port_a") + self.network.connect("mypipe1", "port_b", "mytee1", "port_out1") + self.network.connect("mytee1", "port_in", "mytank", "port_a") + + self.closure = TestModelClosure( + medium=self.medium, + components=TestModelClosureComponents( + cylinder=self.mycylinder, + upstream_tee=self.mytee, + upper_branch=BranchClosureComponents( + name="upper_branch", + orifice=self.myorifice, + pipe=self.mypipe, + ), + lower_branch=BranchClosureComponents( + name="lower_branch", + orifice=self.myorifice1, + pipe=self.mypipe1, + ), + downstream_tee=self.mytee1, + tank=self.mytank, + ), + initial_state_vector=self.initial_state_vector, + apply_state_vector=self.apply_state_vector, + ) + + def initial_state_vector(self) -> list[float]: + return self.network.initial_state_vector() + + def apply_state_vector(self, values: list[float]) -> None: + self.network.apply_state_vector(values) + + def consistent_initial_state_vector(self) -> list[float]: + return self.closure.consistent_initial_state_vector() + + @property + def last_solve_diagnostics(self): + return self.closure.last_solve_diagnostics + + def initialize_consistent_state( + self, + max_iterations: int = 12, + state_tolerance: float = 1e-9, + flow_tolerance: float = 1e-9, + enthalpy_tolerance: float = 1e-6, + pressure_tolerance: float = 1e-6, + strict_internal_solvers: bool = False, + ) -> InitializationDiagnostics: + return self.closure.initialize_consistent_state( + max_iterations=max_iterations, + state_tolerance=state_tolerance, + flow_tolerance=flow_tolerance, + enthalpy_tolerance=enthalpy_tolerance, + pressure_tolerance=pressure_tolerance, + strict_internal_solvers=strict_internal_solvers, + ) + + def project_downstream_pressure_constraints(self, *, strict: bool = False) -> None: + self.closure.project_downstream_pressure_constraints(strict=strict) + + def snapshot( + self, + state_vector: list[float] | None = None, + *, + strict: bool = False, + ) -> TestModelSnapshot: + return self.closure.snapshot(state_vector, strict=strict) + + def rhs(self, _t: float, state_vector: list[float]) -> list[float]: + return self.closure.rhs(state_vector) + + @staticmethod + def _legacy_branch_series_key_map() -> tuple[tuple[str, str, str], tuple[str, str, str]]: + return ( + ("upper_branch", "branch_upper.in", "branch_upper.out"), + ("lower_branch", "branch_lower.in", "branch_lower.out"), + ) + + @classmethod + def _legacy_branch_series_keys_by_name(cls) -> dict[str, tuple[str, str]]: + return { + branch_name: (inlet_key, outlet_key) + for branch_name, inlet_key, outlet_key in cls._legacy_branch_series_key_map() + } + + @staticmethod + def _generic_branch_series_keys(branch_name: str) -> tuple[str, str, str]: + return ( + f"branch.{branch_name}.p", + f"branch.{branch_name}.in", + f"branch.{branch_name}.out", + ) + + @staticmethod + def _legacy_branch_pressure_keys_by_name() -> dict[str, str]: + return { + "upper_branch": "mypipe.p", + "lower_branch": "mypipe1.p", + } + + @classmethod + def _append_legacy_branch_series_aliases( + cls, + series: dict[str, list[float]], + ) -> dict[str, list[float]]: + legacy_branch_series_keys = cls._legacy_branch_series_keys_by_name() + legacy_branch_pressure_keys = cls._legacy_branch_pressure_keys_by_name() + for branch_name, (legacy_inlet_key, legacy_outlet_key) in legacy_branch_series_keys.items(): + pressure_key, generic_inlet_key, generic_outlet_key = cls._generic_branch_series_keys( + branch_name + ) + series[legacy_branch_pressure_keys[branch_name]] = list(series[pressure_key]) + series[legacy_inlet_key] = list(series[generic_inlet_key]) + series[legacy_outlet_key] = list(series[generic_outlet_key]) + return series + + def simulate( + self, + config: SolveIVPConfig | None = None, + t_eval: list[float] | None = None, + ) -> Any: + return integrate_ode( + rhs=self.rhs, + initial_state=self.consistent_initial_state_vector(), + config=config or SolveIVPConfig(), + t_eval=t_eval, + ) + + def evaluate_solution(self, solution: Any) -> dict[str, list[float]]: + series = { + "time": [], + "mycylinder.p": [], + "mycylinder.T": [], + "mytank.p": [], + "mytank.T": [], + } + for branch_name, _, _ in self._legacy_branch_series_key_map(): + pressure_key, inlet_key, outlet_key = self._generic_branch_series_keys(branch_name) + series[pressure_key] = [] + series[inlet_key] = [] + series[outlet_key] = [] + + for index, time_value in enumerate(solution.t): + state_vector = [row[index] for row in solution.y] + snapshot = self.snapshot(state_vector) + series["time"].append(float(time_value)) + series["mycylinder.p"].append(snapshot.cylinder.p) + series["mycylinder.T"].append(snapshot.cylinder.T) + series["mytank.p"].append(snapshot.tank.p) + series["mytank.T"].append(snapshot.tank.T) + for branch in snapshot.branches: + pressure_key, generic_inlet_key, generic_outlet_key = self._generic_branch_series_keys( + branch.name + ) + series[pressure_key].append(branch.pipe.p) + series[generic_inlet_key].append(branch.inlet_flow) + series[generic_outlet_key].append(branch.outlet_flow) + + return self._append_legacy_branch_series_aliases(series) + + +def build_testmodel() -> SimulationNetwork: + """Compatibility helper for callers that only need the topology.""" + + return TestModelSystem().network diff --git a/PythonModels/systems/testmodel_closure.py b/PythonModels/systems/testmodel_closure.py new file mode 100644 index 0000000..cbcaff8 --- /dev/null +++ b/PythonModels/systems/testmodel_closure.py @@ -0,0 +1,668 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Callable + +from PythonModels.components.cylinder import Cylinder +from PythonModels.components.orifice import Orifice +from PythonModels.components.pipe import Pipe +from PythonModels.components.tank import Tank +from PythonModels.components.tee import Tee +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.state import VolumeState + + +@dataclass(frozen=True) +class BranchInletFlowDiagnostics: + converged: bool + iterations: int + residual: float + m_flow: float + inlet_pressure: float + + +@dataclass(frozen=True) +class DownstreamPressureDiagnostics: + converged: bool + iterations: int + residual: float + pressure: float + target_total_internal_energy: float + + +@dataclass(frozen=True) +class TestModelSolveDiagnostics: + upper_branch_inlet: BranchInletFlowDiagnostics + lower_branch_inlet: BranchInletFlowDiagnostics + downstream_pressure_projection: DownstreamPressureDiagnostics | None + + +@dataclass(frozen=True) +class BranchClosureComponents: + name: str + orifice: Orifice + pipe: Pipe + + +@dataclass(frozen=True) +class BranchClosureState: + name: str + pipe: ThermodynamicProperties + inlet_flow: float + outlet_flow: float + inlet_h: float + inlet_flow_diagnostics: BranchInletFlowDiagnostics + + +@dataclass(frozen=True) +class BranchSnapshot: + name: str + pipe: ThermodynamicProperties + inlet_flow: float + outlet_flow: float + inlet_h: float + inlet_flow_diagnostics: BranchInletFlowDiagnostics + + +@dataclass(frozen=True) +class TestModelSnapshot: + cylinder: ThermodynamicProperties + tank: ThermodynamicProperties + tee_upstream_h: float + tee_downstream_h: float + branches: tuple[BranchSnapshot, ...] = field(default_factory=tuple) + solve_diagnostics: TestModelSolveDiagnostics | None = None + + @property + def pipe_upper(self) -> ThermodynamicProperties: + return self.branches[0].pipe + + @property + def pipe_lower(self) -> ThermodynamicProperties: + return self.branches[1].pipe + + @property + def branch_inlet_flows(self) -> tuple[float, ...]: + return tuple(branch.inlet_flow for branch in self.branches) + + @property + def branch_outlet_flows(self) -> tuple[float, ...]: + return tuple(branch.outlet_flow for branch in self.branches) + + +@dataclass(frozen=True) +class InitializationDiagnostics: + converged: bool + iterations: int + max_state_delta: float + max_flow_delta: float + max_enthalpy_delta: float + downstream_pressure_spread: float + state_vector: tuple[float, ...] + + +@dataclass(frozen=True) +class TestModelClosureComponents: + cylinder: Cylinder + upstream_tee: Tee + upper_branch: BranchClosureComponents + lower_branch: BranchClosureComponents + downstream_tee: Tee + tank: Tank + + def branches(self) -> tuple[BranchClosureComponents, BranchClosureComponents]: + return (self.upper_branch, self.lower_branch) + + +class TestModelClosure: + """Owns Testmodel-specific closure, projection and port-writeback logic.""" + + def __init__( + self, + *, + medium: IdealGasMedium, + components: TestModelClosureComponents, + initial_state_vector: Callable[[], list[float]], + apply_state_vector: Callable[[list[float]], None], + ) -> None: + self.medium = medium + self.components = components + self._initial_state_vector = initial_state_vector + self._apply_state_vector = apply_state_vector + self.last_solve_diagnostics: TestModelSolveDiagnostics | None = None + self.last_downstream_pressure_diagnostics: DownstreamPressureDiagnostics | None = None + + @staticmethod + def _downstream_pressure_spread(snapshot: TestModelSnapshot) -> float: + downstream_pressures = tuple(branch.pipe.p for branch in snapshot.branches) + ( + snapshot.tank.p, + ) + return max(downstream_pressures) - min(downstream_pressures) + + @staticmethod + def _initialization_flow_delta( + previous_snapshot: TestModelSnapshot | None, + current_snapshot: TestModelSnapshot, + ) -> float: + if previous_snapshot is None: + return max(abs(branch.outlet_flow) for branch in current_snapshot.branches) + return max( + abs(curr - prev) + for curr, prev in zip( + current_snapshot.branch_outlet_flows, + previous_snapshot.branch_outlet_flows, + ) + ) + + @staticmethod + def _initialization_enthalpy_delta( + previous_snapshot: TestModelSnapshot | None, + current_snapshot: TestModelSnapshot, + ) -> float: + if previous_snapshot is None: + return abs(current_snapshot.tee_downstream_h - current_snapshot.tank.h) + return max( + abs(current_snapshot.tee_upstream_h - previous_snapshot.tee_upstream_h), + abs(current_snapshot.tee_downstream_h - previous_snapshot.tee_downstream_h), + ) + + def consistent_initial_state_vector(self) -> list[float]: + return list(self.initialize_consistent_state().state_vector) + + def initialize_consistent_state( + self, + max_iterations: int = 12, + state_tolerance: float = 1e-9, + flow_tolerance: float = 1e-9, + enthalpy_tolerance: float = 1e-6, + pressure_tolerance: float = 1e-6, + strict_internal_solvers: bool = False, + ) -> InitializationDiagnostics: + raw_state = self._initial_state_vector() + previous_snapshot: TestModelSnapshot | None = None + diagnostics: InitializationDiagnostics | None = None + + for iteration in range(1, max_iterations + 1): + state_before_projection = self._initial_state_vector() + self.snapshot(state_before_projection, strict=strict_internal_solvers) + + self.project_downstream_pressure_constraints(strict=strict_internal_solvers) + + state_after_projection = self._initial_state_vector() + snapshot_after_projection = self.snapshot( + state_after_projection, + strict=strict_internal_solvers, + ) + + max_state_delta = max( + abs(after - before) + for before, after in zip(state_before_projection, state_after_projection) + ) + max_flow_delta = self._initialization_flow_delta( + previous_snapshot, + snapshot_after_projection, + ) + max_enthalpy_delta = self._initialization_enthalpy_delta( + previous_snapshot, + snapshot_after_projection, + ) + downstream_pressure_spread = self._downstream_pressure_spread( + snapshot_after_projection, + ) + + diagnostics = InitializationDiagnostics( + converged=( + max_state_delta <= state_tolerance + and max_flow_delta <= flow_tolerance + and max_enthalpy_delta <= enthalpy_tolerance + and downstream_pressure_spread <= pressure_tolerance + ), + iterations=iteration, + max_state_delta=max_state_delta, + max_flow_delta=max_flow_delta, + max_enthalpy_delta=max_enthalpy_delta, + downstream_pressure_spread=downstream_pressure_spread, + state_vector=tuple(state_after_projection), + ) + previous_snapshot = snapshot_after_projection + + if diagnostics.converged: + self._apply_state_vector(raw_state) + return diagnostics + + assert diagnostics is not None + self._apply_state_vector(raw_state) + return diagnostics + + def _solve_branch_inlet_flow( + self, + orifice: Orifice, + pipe: Pipe, + p_upstream: float, + pipe_props: ThermodynamicProperties, + *, + strict: bool = False, + ) -> tuple[float, BranchInletFlowDiagnostics]: + m_flow = orifice.mass_flow(p_upstream, pipe_props.p) + rho = max(pipe_props.rho, 1e-9) + p_inlet = pipe.inlet_pressure(m_flow, rho, pipe_props.p) + residual = abs(orifice.mass_flow(p_upstream, p_inlet) - m_flow) + converged = False + iterations = 0 + for iteration in range(1, 9): + p_inlet = pipe.inlet_pressure(m_flow, rho, pipe_props.p) + next_m_flow = orifice.mass_flow(p_upstream, p_inlet) + residual = abs(next_m_flow - m_flow) + iterations = iteration + if residual <= 1e-9 * max(1.0, abs(next_m_flow)): + m_flow = next_m_flow + converged = True + break + m_flow = next_m_flow + diagnostics = BranchInletFlowDiagnostics( + converged=converged, + iterations=iterations, + residual=residual, + m_flow=m_flow, + inlet_pressure=p_inlet, + ) + if strict and not diagnostics.converged: + raise RuntimeError( + f"Branch inlet flow solve did not converge for {pipe.name}: residual={residual:.6e}" + ) + return m_flow, diagnostics + + def _solve_downstream_branch_flows( + self, + cylinder: ThermodynamicProperties, + tank: ThermodynamicProperties, + branch_states: tuple[BranchClosureState, BranchClosureState], + ) -> tuple[float, float]: + return self._solve_downstream_branch_flows_from_state( + inlet_h_upper=branch_states[0].inlet_h, + inlet_h_lower=branch_states[1].inlet_h, + pipe_upper_h=max(branch_states[0].pipe.h, 1e-9), + pipe_lower_h=max(branch_states[1].pipe.h, 1e-9), + tank_h=max(tank.h, 1e-9), + q_in_upper=branch_states[0].inlet_flow, + q_in_lower=branch_states[1].inlet_flow, + ) + + def _project_volume_energy_to_pressure( + self, + component: Pipe | Tank, + target_pressure: float, + ) -> None: + target_temperature = target_pressure * component.V / ( + max(component.state.m, 1e-12) * self.medium.R_gas + ) + target_internal_energy = ( + component.state.m * self.medium.specific_internal_energy(target_temperature) + ) + component.state = VolumeState(m=component.state.m, U=target_internal_energy) + + def _downstream_total_internal_energy_for_pressure( + self, + target_pressure: float, + downstream_components: tuple[Pipe | Tank, ...], + ) -> float: + total_internal_energy = 0.0 + for component in downstream_components: + target_temperature = target_pressure * component.V / ( + max(component.state.m, 1e-12) * self.medium.R_gas + ) + total_internal_energy += ( + component.state.m * self.medium.specific_internal_energy(target_temperature) + ) + return total_internal_energy + + def _solve_downstream_common_pressure( + self, + downstream_components: tuple[Pipe | Tank, ...], + target_total_internal_energy: float, + *, + strict: bool = False, + ) -> tuple[float, DownstreamPressureDiagnostics]: + lower_pressure = 1.0 + upper_pressure = max(component.properties().p for component in downstream_components) + upper_pressure = max(upper_pressure, 1e5) + + def residual(pressure: float) -> float: + return ( + self._downstream_total_internal_energy_for_pressure( + pressure, + downstream_components, + ) + - target_total_internal_energy + ) + + upper_residual = residual(upper_pressure) + iteration_count = 0 + + while upper_residual < 0.0: + upper_pressure *= 2.0 + upper_residual = residual(upper_pressure) + + final_pressure = 0.5 * (lower_pressure + upper_pressure) + final_residual = residual(final_pressure) + converged = False + for iteration in range(1, 81): + middle_pressure = 0.5 * (lower_pressure + upper_pressure) + middle_residual = residual(middle_pressure) + iteration_count = iteration + final_pressure = middle_pressure + final_residual = middle_residual + if abs(middle_residual) <= 1e-12 * max(1.0, target_total_internal_energy): + converged = True + break + if middle_residual > 0.0: + upper_pressure = middle_pressure + else: + lower_pressure = middle_pressure + + diagnostics = DownstreamPressureDiagnostics( + converged=converged, + iterations=iteration_count, + residual=final_residual, + pressure=final_pressure, + target_total_internal_energy=target_total_internal_energy, + ) + if strict and not diagnostics.converged: + raise RuntimeError( + "Downstream common-pressure solve did not converge: " + f"residual={final_residual:.6e}" + ) + return final_pressure, diagnostics + + def project_downstream_pressure_constraints(self, *, strict: bool = False) -> None: + downstream_components = ( + self.components.upper_branch.pipe, + self.components.lower_branch.pipe, + self.components.tank, + ) + total_internal_energy = sum(component.state.U for component in downstream_components) + common_pressure, diagnostics = self._solve_downstream_common_pressure( + downstream_components, + total_internal_energy, + strict=strict, + ) + self.last_downstream_pressure_diagnostics = diagnostics + + for component in downstream_components: + self._project_volume_energy_to_pressure(component, common_pressure) + + def _downstream_connection_enthalpy( + self, + q_out_upper: float, + q_out_lower: float, + pipe_upper_h: float, + pipe_lower_h: float, + tank_h: float, + ) -> float: + return self.components.downstream_tee.inlet_stream_enthalpy( + q_out_lower, + pipe_lower_h, + q_out_upper, + pipe_upper_h, + fallback_h=tank_h, + ) + + def _solve_downstream_branch_flows_from_state( + self, + *, + inlet_h_upper: float, + inlet_h_lower: float, + pipe_upper_h: float, + pipe_lower_h: float, + tank_h: float, + q_in_upper: float, + q_in_lower: float, + ) -> tuple[float, float]: + return self.components.downstream_tee.solve_branch_outlet_flows_from_energy_balance( + ratio_branch1=self.components.upper_branch.pipe.V / self.components.tank.V, + ratio_branch2=self.components.lower_branch.pipe.V / self.components.tank.V, + inlet_h_branch1=inlet_h_upper, + inlet_h_branch2=inlet_h_lower, + branch1_h=pipe_upper_h, + branch2_h=pipe_lower_h, + inlet_h=tank_h, + q_in_branch1=q_in_upper, + q_in_branch2=q_in_lower, + ) + + def _evaluate_branch_states( + self, + cylinder: ThermodynamicProperties, + ) -> tuple[BranchClosureState, BranchClosureState]: + states: list[BranchClosureState] = [] + for branch in self.components.branches(): + pipe_properties = branch.pipe.properties() + inlet_flow, inlet_flow_diagnostics = self._solve_branch_inlet_flow( + branch.orifice, + branch.pipe, + cylinder.p, + pipe_properties, + ) + inlet_h = branch.pipe.port_a_inlet_enthalpy( + port_a_m_flow=inlet_flow, + connected_h=cylinder.h, + internal_h=pipe_properties.h, + ) + states.append( + BranchClosureState( + name=branch.name, + pipe=pipe_properties, + inlet_flow=inlet_flow, + outlet_flow=0.0, + inlet_h=inlet_h, + inlet_flow_diagnostics=inlet_flow_diagnostics, + ) + ) + return (states[0], states[1]) + + @staticmethod + def _with_branch_outlet_flows( + branch_states: tuple[BranchClosureState, BranchClosureState], + outlet_flows: tuple[float, float], + ) -> tuple[BranchClosureState, BranchClosureState]: + return ( + BranchClosureState( + name=branch_states[0].name, + pipe=branch_states[0].pipe, + inlet_flow=branch_states[0].inlet_flow, + outlet_flow=outlet_flows[0], + inlet_h=branch_states[0].inlet_h, + inlet_flow_diagnostics=branch_states[0].inlet_flow_diagnostics, + ), + BranchClosureState( + name=branch_states[1].name, + pipe=branch_states[1].pipe, + inlet_flow=branch_states[1].inlet_flow, + outlet_flow=outlet_flows[1], + inlet_h=branch_states[1].inlet_h, + inlet_flow_diagnostics=branch_states[1].inlet_flow_diagnostics, + ), + ) + + @staticmethod + def _branch_snapshots( + branch_states: tuple[BranchClosureState, BranchClosureState], + ) -> tuple[BranchSnapshot, BranchSnapshot]: + return ( + BranchSnapshot( + name=branch_states[0].name, + pipe=branch_states[0].pipe, + inlet_flow=branch_states[0].inlet_flow, + outlet_flow=branch_states[0].outlet_flow, + inlet_h=branch_states[0].inlet_h, + inlet_flow_diagnostics=branch_states[0].inlet_flow_diagnostics, + ), + BranchSnapshot( + name=branch_states[1].name, + pipe=branch_states[1].pipe, + inlet_flow=branch_states[1].inlet_flow, + outlet_flow=branch_states[1].outlet_flow, + inlet_h=branch_states[1].inlet_h, + inlet_flow_diagnostics=branch_states[1].inlet_flow_diagnostics, + ), + ) + + def snapshot( + self, + state_vector: list[float] | None = None, + *, + strict: bool = False, + ) -> TestModelSnapshot: + if state_vector is not None: + self._apply_state_vector(list(state_vector)) + + cylinder = self.components.cylinder.properties() + tank = self.components.tank.properties() + branch_states = self._evaluate_branch_states(cylinder) + if strict: + for branch_state in branch_states: + if not branch_state.inlet_flow_diagnostics.converged: + raise RuntimeError( + "Branch inlet flow solve did not converge for " + f"{branch_state.name}: residual=" + f"{branch_state.inlet_flow_diagnostics.residual:.6e}" + ) + outlet_flows = self._solve_downstream_branch_flows(cylinder, tank, branch_states) + branch_states = self._with_branch_outlet_flows(branch_states, outlet_flows) + + tee_upstream_h = self.components.upstream_tee.inlet_stream_enthalpy( + -branch_states[0].inlet_flow, + branch_states[0].pipe.h, + -branch_states[1].inlet_flow, + branch_states[1].pipe.h, + fallback_h=cylinder.h, + ) + tee_downstream_h = self._downstream_connection_enthalpy( + branch_states[0].outlet_flow, + branch_states[1].outlet_flow, + branch_states[0].pipe.h, + branch_states[1].pipe.h, + tank.h, + ) + + self._write_port_states( + cylinder, + tank, + branch_states, + tee_upstream_h, + tee_downstream_h, + ) + + solve_diagnostics = TestModelSolveDiagnostics( + upper_branch_inlet=branch_states[0].inlet_flow_diagnostics, + lower_branch_inlet=branch_states[1].inlet_flow_diagnostics, + downstream_pressure_projection=self.last_downstream_pressure_diagnostics, + ) + self.last_solve_diagnostics = solve_diagnostics + branch_snapshots = self._branch_snapshots(branch_states) + + return TestModelSnapshot( + cylinder=cylinder, + tank=tank, + tee_upstream_h=tee_upstream_h, + tee_downstream_h=tee_downstream_h, + branches=branch_snapshots, + solve_diagnostics=solve_diagnostics, + ) + + def _write_port_states( + self, + cylinder: ThermodynamicProperties, + tank: ThermodynamicProperties, + branch_states: tuple[BranchClosureState, BranchClosureState], + tee_upstream_h: float, + tee_downstream_h: float, + ) -> None: + cylinder_m_flow = -sum(branch_state.inlet_flow for branch_state in branch_states) + tank_m_flow = sum(branch_state.outlet_flow for branch_state in branch_states) + + self.components.cylinder.port_b.m_flow = cylinder_m_flow + + self.components.upstream_tee.port_in.p = cylinder.p + self.components.upstream_tee.port_out1.p = cylinder.p + self.components.upstream_tee.port_out2.p = cylinder.p + self.components.upstream_tee.port_in.m_flow = cylinder_m_flow + self.components.upstream_tee.port_in.h_outflow = tee_upstream_h + self.components.upstream_tee.port_out1.h_outflow = cylinder.h + self.components.upstream_tee.port_out2.h_outflow = cylinder.h + self.components.upstream_tee.port_out1.m_flow = -branch_states[0].inlet_flow + self.components.upstream_tee.port_out2.m_flow = -branch_states[1].inlet_flow + + for branch_components, branch_state in zip(self.components.branches(), branch_states): + branch_components.orifice.port_a.p = cylinder.p + branch_components.orifice.port_b.p = branch_components.pipe.inlet_pressure( + branch_state.inlet_flow, + max(branch_state.pipe.rho, 1e-9), + branch_state.pipe.p, + ) + branch_components.orifice.port_a.m_flow = branch_state.inlet_flow + branch_components.orifice.port_b.m_flow = -branch_state.inlet_flow + branch_components.orifice.port_a.h_outflow = cylinder.h + branch_components.orifice.port_b.h_outflow = branch_state.pipe.h + + branch_components.pipe.port_a.p = branch_components.orifice.port_b.p + branch_components.pipe.port_a.m_flow = branch_state.inlet_flow + branch_components.pipe.port_b.m_flow = -branch_state.outlet_flow + branch_components.pipe.port_b.p = branch_state.pipe.p + + self.components.downstream_tee.port_in.p = tank.p + self.components.downstream_tee.port_out1.p = tank.p + self.components.downstream_tee.port_out2.p = tank.p + self.components.downstream_tee.port_in.m_flow = -tank_m_flow + self.components.downstream_tee.port_out1.m_flow = branch_states[1].outlet_flow + self.components.downstream_tee.port_out2.m_flow = branch_states[0].outlet_flow + self.components.downstream_tee.port_in.h_outflow = tee_downstream_h + self.components.downstream_tee.port_out1.h_outflow = tank.h + self.components.downstream_tee.port_out2.h_outflow = tank.h + + self.components.tank.port_a.m_flow = tank_m_flow + + def _branch_derivative_states( + self, + snapshot: TestModelSnapshot, + ) -> tuple[VolumeState, VolumeState]: + derivative_states: list[VolumeState] = [] + for branch_components, branch_snapshot in zip(self.components.branches(), snapshot.branches): + derivative_states.append( + branch_components.pipe.derivatives_from_connections( + port_a_m_flow=branch_snapshot.inlet_flow, + connected_h_a=snapshot.cylinder.h, + port_b_m_flow=-branch_snapshot.outlet_flow, + connected_h_b=snapshot.tank.h, + internal_h=branch_snapshot.pipe.h, + ) + ) + return (derivative_states[0], derivative_states[1]) + + def rhs(self, state_vector: list[float]) -> list[float]: + snapshot = self.snapshot(state_vector) + + cylinder_m_flow = -sum(branch.inlet_flow for branch in snapshot.branches) + tank_m_flow = sum(branch.outlet_flow for branch in snapshot.branches) + d_cylinder = self.components.cylinder.derivatives_from_connection( + connected_h=snapshot.tee_upstream_h, + port_m_flow=cylinder_m_flow, + internal_h=snapshot.cylinder.h, + ) + branch_derivatives = self._branch_derivative_states(snapshot) + d_tank = self.components.tank.derivatives_from_connection( + connected_h=snapshot.tee_downstream_h, + port_m_flow=tank_m_flow, + internal_h=snapshot.tank.h, + ) + + return [ + d_cylinder.m, + d_cylinder.U, + branch_derivatives[0].m, + branch_derivatives[0].U, + branch_derivatives[1].m, + branch_derivatives[1].U, + d_tank.m, + d_tank.U, + ]