上传PythonModels文件

PythonModels/.gitignore

@@ -0,0 +1,2 @@
+__pycache__/
+*.pyc

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.<branch_name>.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 原始结果，参数大概率还要继续调。

PythonModels/__init__.py

@@ -0,0 +1,2 @@
+"""Python port scaffold for the Modelica-based pressurization system."""
+

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%

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

PythonModels/components/__init__.py

@@ -0,0 +1,2 @@
+"""Component implementations for the Python system model."""
+

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)

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)
+

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)

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)

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)

PythonModels/core/__init__.py

@@ -0,0 +1,2 @@
+"""Core abstractions for the Python system model."""
+

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."""

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)

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)
+

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
+

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,
+    )

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])
+

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",
+]

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'<line x1="{x:.2f}" y1="{top}" x2="{x:.2f}" y2="{top + plot_height}" '
+            'stroke="#d9e2ec" stroke-width="1" />'
+        )
+        x_tick_markup.append(
+            f'<text x="{x:.2f}" y="{height - 30}" text-anchor="middle" '
+            'font-size="14" fill="#102a43">'
+            f"{time_value:.1f}</text>"
+        )
+
+    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'<line x1="{left}" y1="{y:.2f}" x2="{left + plot_width}" y2="{y:.2f}" '
+            'stroke="#d9e2ec" stroke-width="1" />'
+        )
+        y_tick_markup.append(
+            f'<text x="{left - 12}" y="{y + 5:.2f}" text-anchor="end" '
+            'font-size="14" fill="#102a43">'
+            f"{temp_value:.1f}</text>"
+        )
+
+    svg_content = f"""<svg xmlns="http://www.w3.org/2000/svg" width="{width}" height="{height}" viewBox="0 0 {width} {height}">
+  <rect width="{width}" height="{height}" fill="#f7fafc" rx="18" ry="18" />
+  <text x="{width / 2:.0f}" y="32" text-anchor="middle" font-size="24" fill="#102a43">Python Testmodel Tank Temperature</text>
+  <text x="{width / 2:.0f}" y="{height - 8}" text-anchor="middle" font-size="16" fill="#486581">Time (s)</text>
+  <text x="26" y="{height / 2:.0f}" text-anchor="middle" font-size="16" fill="#486581" transform="rotate(-90 26 {height / 2:.0f})">Temperature (K)</text>
+  <rect x="{left}" y="{top}" width="{plot_width}" height="{plot_height}" fill="#ffffff" stroke="#bcccdc" stroke-width="1.5" />
+  {''.join(x_tick_markup)}
+  {''.join(y_tick_markup)}
+  <polyline fill="none" stroke="#d64545" stroke-width="3" stroke-linejoin="round" stroke-linecap="round" points="{points}" />
+</svg>
+"""
+    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,
+    )

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()

PythonModels/systems/__init__.py

@@ -0,0 +1,2 @@
+"""System assembly modules."""
+

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

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,
+        ]