diff --git a/.gitignore b/.gitignore index 4daa310..e79f529 100644 --- a/.gitignore +++ b/.gitignore @@ -33,3 +33,6 @@ logs/ # OS files .DS_Store Thumbs.db + +# test files +test/ diff --git a/PythonModels/.gitignore b/PythonModels/.gitignore new file mode 100644 index 0000000..7a60b85 --- /dev/null +++ b/PythonModels/.gitignore @@ -0,0 +1,2 @@ +__pycache__/ +*.pyc diff --git a/PythonModels/README.md b/PythonModels/README.md new file mode 100644 index 0000000..35ccea7 --- /dev/null +++ b/PythonModels/README.md @@ -0,0 +1,340 @@ +# PythonModels + +`PythonModels` 用于承接 `ModelicaModels` 的 Python 平台移植。 + +目标不是把 `.mo` 文件逐行翻译成 Python，而是建立一个可运行、可对比、可逐步逼近 `OpenModelica` 行为的 Python 仿真框架。 + +当前状态不是“只有骨架”，而是“`Testmodel` 已有一版可运行的 ODE 近似实现，并具备基础结果导出与对比能力”。 + +## 当前目录 + +- `core/`: 通用基础设施 + 包含组件基类、状态与端口数据结构、介质模型、网络装配、积分入口。 +- `components/`: 元件级 Python 实现 + 目前有 `Cylinder`、`Tank`、`Pipe`、`Orifice`、`Tee` 五类元件。 +- `systems/`: 系统级装配与闭合 + 当前只有 `TestModelSystem`，对应 `ModelicaModels/Testmodel.mo`。 +- `reporting/`: 结果导出与对比 + 当前承接主变量 CSV、温度 CSV/SVG、Python 对 OpenModelica 的对比表与误差摘要导出。 +- `scripts/`: 运行脚本 + 当前入口是 `run_testmodel.py`。 +- `baselines/`: 提交进仓库的稳定基线 + 当前承接 Python 主变量基线和 Python 对 Modelica 的误差摘要基线。 +- `runs/`: 每次实际运行的默认输出目录 + 当前脚本默认会在这里创建带时间戳的子目录，用来放这次运行生成的产物。 + +当前关键文件： + +- `core/medium.py`: 理想气体近似介质 `IdealGasMedium` +- `core/medium.py`: 温度相关的空气近似介质 `IdealGasMedium` +- `core/network.py`: `SimulationNetwork`，负责组件注册、连接拓扑和状态向量拼装 +- `core/solver.py`: `integrate_ode()`，优先走 `SciPy solve_ivp`，缺依赖时回退到内置 RK4，并支持 `t_start == t_stop` 的零时长返回 +- `components/pipe.py`: 单阻容管道近似，入口压降 + 出口直连内容腔 +- `components/tee.py`: 三通的最小 stream 混合 helper +- `systems/testmodel.py`: `Testmodel` 的系统装配壳与外部运行入口 +- `systems/testmodel_closure.py`: `Testmodel` 当前专用的闭合、初始化投影、分支求解与端口回写 +- `reporting/testmodel_outputs.py`: `Testmodel` 的 CSV/SVG/对比摘要导出 +- `scripts/run_testmodel.py`: 基线运行与程序化执行入口 +- `tests/test_pythonmodels_regression.py`: 当前 Python 基线回归测试 + +## 当前阶段进度 + +这一阶段原先有 4 件重点工作，现在的状态如下： + +1. `mytee1` 的 stream/焓传播语义：已完成当前阶段收紧 + 现在如果只有一条支路发生倒流，下游来流焓统一按 `tank.h` 处理，不再临时借另一条支路的焓来凑。 +2. 下游初始化/约束处理：已完成当前阶段收口 + 之前是“直接改对象状态再开始积分”，现在已经收成显式的 `consistent_initial_state_vector()` 初始化入口。当前这一步会在不改下游总质量、总内能的前提下，把几段直接相连的体积拉回同一个连接压力。 +3. 自动校验：已完成当前阶段首版 + 已经补了标准库 `unittest` 回归测试，先把初始化投影是否守恒、是否污染原始状态，以及 4 个主变量的提交基线锁住。 +4. 更严格介质模型：已完成当前阶段首版 + 已经从固定 `cp/cv` 的理想气体近似，推进到随温度变化的空气近似，并接上了内能反解和初始化求根。 + +如果只看结果，可以把这一阶段理解成： + +- 连接器语义：首轮收紧已完成 +- 初始化入口：首轮收口已完成 +- 基线验证：首轮保护已完成 +- 介质精化：首轮近似已完成 + +## 当前阶段收口 + +上一轮 `N0-N3` 已全部完成首版，当前可以简单理解为： + +1. `N0`：系统层里最明显的流向/焓判断已经继续下沉到组件 helper。 +2. `N1`：模型参数和运行参数已经收口到配置对象。 +3. `N2`：运行接口已经分成“准备请求”和“执行请求”两层。 +4. `N3`：结果导出和命令行报告格式化已经统一收口到 `reporting/`。 + +这一轮结束后，项目已经不缺“能不能跑”的能力，下一步更重要的是把后续开发最容易卡住的地方先处理掉。 + +## 本次推送更新 + +本次推送已经把上一轮建议里的 `M2-M5` 推进到下面这个状态： + +1. `M2`：已完成当前阶段首版 + - 已把 `Testmodel` 的专用闭合、初始化投影、分支入口流量求解、下游支路出口流量闭合、端口状态回写，从 `systems/testmodel.py` 拆到新的 `systems/testmodel_closure.py` + - `TestModelSystem` 现在主要承担组件装配、网络注册和对闭合器的委托，不再继续堆积系统级手写细节 + +2. `M3`：已完成当前阶段首版 + - 已给两条支路入口流量固定点求解、下游公共压力投影补了显式诊断 + - 诊断内容至少包含 `converged / iterations / residual` + - 已支持严格模式；内部求解不收敛时可以直接抛错，而不是静默返回最后一个近似值 + - `run_testmodel()` 的结构化结果和 `testmodel_run_report.txt` 已能带出最后一次内部闭合求解诊断 + +3. `M4`：已完成当前阶段首版 + - 自动测试已不再只盯最终主变量结果 + - 现在已经覆盖： + - 改支路参数后，初始支路入口流量是否按预期变化 + - 更偏激配置下，初始化和内部闭合是否仍然收敛 + - 有无 Modelica 参考两种运行路径下，程序接口与产物行为是否一致 + +4. `M5`：已启动 + - 当前已经明确选择优先走“更容易扩展”的方向，而不是先追求更贴近 Modelica + - 已完成第一步：把闭合器内部原来大量写死的 `upper/lower` 双支路逻辑，收成可复用的 `BranchClosureComponents / BranchClosureState` 结构 + - 当前已继续推进到 `G1-G5` 的首轮兼容层改造：`snapshot` 已提供通用分支集合，系统层结果生成已拆成“通用键生成 + 旧键别名派生”两层，报告层已开始优先消费通用分支键，旧导出列名仍通过兼容映射保留，兼容测试已显式保护分支顺序和旧导出语义 + +## 下一阶段接手建议 + +如果继续往前推进，建议按下面顺序做，而不是再零散补功能： + +1. `G1`：已完成当前阶段首轮兼容接入 + - `TestModelSnapshot` 已新增 `branches` 集合 + - 每个分支当前至少带 `name / pipe / inlet_flow / outlet_flow / inlet_h / inlet_flow_diagnostics` + - `pipe_upper / pipe_lower / branch_inlet_flows / branch_outlet_flows` 目前仍保留为兼容属性，供旧调用方继续使用 + +2. `G2`：已完成当前阶段首轮内部迁移 + - `evaluate_solution()` 已改成从 `snapshot.branches` 读取数据，再通过显式分支名映射写回当前旧列名 + - `rhs()` 里的分支导数计算已改成通过通用 helper 按分支循环生成，再按当前状态向量顺序拼回 + - 当前外部导出列名仍保持兼容： + - `mypipe.p` + - `mypipe1.p` + - `branch_upper.in/out` + - `branch_lower.in/out` + +3. `G3`：已完成当前阶段首轮兼容测试 + - 当前测试已经显式保护： + - `branches` 顺序是否稳定 + - `snapshot` 新字段和兼容字段是否一致 + - 旧导出列名是否仍映射到正确分支语义 + - 参数变化后 `upper/lower` 的名字和顺序是否不会被打乱 + +4. `G4`：已完成当前阶段首轮兼容拆层 + - `evaluate_solution()` 现在会同时产出： + - 通用分支键：`branch..p/in/out` + - 旧兼容键：`mypipe.p`、`mypipe1.p`、`branch_upper.*`、`branch_lower.*` + - 报告层当前已开始优先读取通用分支键，旧键只作为兼容后备 + - 当前已经把“内部统一表达”和“旧接口兼容导出”拆成两层，但还没有把所有报告/导出逻辑都迁干净 + +5. `G5`：已完成当前阶段首轮兼容收口 + - `evaluate_solution()` 当前会先生成通用分支键，再统一派生旧兼容键 + - 报告层当前已支持“通用键优先、旧键兼容后备” + - 当前已经把系统层和 reporting 层的主要旧专名读取入口收口到少量 helper 上，后续继续迁移不会再到处散改 + +6. `P1`：下一阶段建议从这里接手 + 当前更合适的下一步，不是继续深挖内核通用化，而是切回结果导向主线： + - 定义一份稳定的外部输入参数 schema + - 明确这些结构化参数如何映射到 `TestModelConfig / TestModelRunConfig` + - 建立“结构化参数 -> 仿真执行 -> 结果产物/摘要”的稳定接口 + 这样可以直接服务后续文档解析、网页入口和报告生成，而不是继续在 `Testmodel` 内部做边际收益越来越低的抽象整理 + +7. `P2`：在 `P1` 完成后，再推进文档解析或报告生成链路 + 更现实的顺序应是： + - 先把结构化输入跑通 + - 再把结果摘要/产物组织成更接近最终产品的输出包 + - 最后再接 Word 解析或页面入口 + +如果后续继续推进，这个 README 也要一起更新，不要长期保留已经失效的路线描述。 + +## 当前实现了什么 + +当前代码已经实现： + +1. `m`、`U` 作为动态元件主状态，`p`、`T`、`rho`、`u`、`h` 作为派生量。 +2. `Cylinder`、`Tank`、`Pipe` 的刚性绝热容腔近似。 +3. `Orifice` 的压差开方流量关系。 +4. `Tee` 的简化混合焓处理。 +5. `Testmodel` 的系统级拓扑映射和一版可运行的 `rhs(t, x)`。 +6. 基于 `solve_ivp` 的积分入口，以及 SciPy 不可用时的 RK4 回退。 +7. 温度相关空气近似介质，包括 `cp(T)`、`h(T)`、`u(T)` 以及 `u -> T` 反解。 +8. 显式一致初值入口 `consistent_initial_state_vector()`，以及可迭代初始化器 `initialize_consistent_state()`。 +9. Python 主变量结果导出： + `mytank.p`、`mytank.T`、`mycylinder.p`、`mycylinder.T` +10. 贮箱温度曲线导出： + `testmodel_tank_temperature.csv` + `testmodel_tank_temperature.svg` +11. 基于 `ModelicaModels/Simulation/Testmodel_res.csv` 的逐时刻对比与误差摘要导出。 +12. 基于 `unittest` 的自动回归测试，当前已覆盖初始化守恒、主变量基线、运行接口、内部闭合诊断、通用分支兼容层、通用结果键与旧键别名一致性，以及部分中间闭合过程行为。 + +当前没有实现： + +- 通用 DAE 初始化器 +- `Modelica.Media.Air.SimpleAir` 的严格复刻 +- 面向任意拓扑的通用 connector/stream 求解器 + +## 当前怎么运行 + +最小运行方式： + +```bash +python3 -m PythonModels.scripts.run_testmodel +``` + +如果要改模型参数或运行参数，建议直接改配置对象，而不是改源码里的默认值。例如： + +```python +from PythonModels.core.solver import SolveIVPConfig +from PythonModels.scripts.run_testmodel import ( + TestModelRunConfig, + TestModelSamplingConfig, + run_testmodel, +) +from PythonModels.systems.testmodel import ( + BranchConfig, + CylinderConfig, + OrificeConfig, + PipeConfig, + TankConfig, + TestModelConfig, +) + +run_config = TestModelRunConfig( + model=TestModelConfig( + cylinder=CylinderConfig(p0=30e6), + upper_branch=BranchConfig( + orifice=OrificeConfig(K=8e-6), + pipe=PipeConfig(length=6.0, diameter=0.03), + ), + tank=TankConfig(volume=0.12), + ), + solver=SolveIVPConfig(t_start=0.0, t_stop=10.0, method="BDF"), + sampling=TestModelSamplingConfig(step=0.05), +) + +result = run_testmodel(run_config=run_config) +``` + +如果调用方想先确认“这次运行最后到底会用哪些路径、哪些采样点”，可以先准备请求，再执行： + +```python +from PythonModels.scripts.run_testmodel import ( + prepare_testmodel_run, + run_prepared_testmodel, + TestModelRunConfig, +) + +prepared = prepare_testmodel_run(run_config=TestModelRunConfig()) +print(prepared.output_dir) +print(prepared.t_eval) + +result = run_prepared_testmodel(prepared) +print(result.artifacts.primary_csv_path) +print(result.used_modelica_reference) +``` + +当前脚本会： + +1. 构建 `TestModelSystem` +2. 打印原始初值向量与约束一致后的初值向量 +3. 运行 `0 s -> 20 s` 的仿真，默认采样间隔 `0.1 s` +4. 将结果写入 `PythonModels/runs/` 下本次运行专属的时间戳目录 +5. 若存在 `ModelicaModels/Simulation/Testmodel_res.csv`，自动生成 Python 与 OpenModelica 对比结果 + +当前脚本默认不会再把运行结果直接写到提交基线目录，而是会在 `PythonModels/runs/` 下创建一个带时间戳的子目录，例如： + +- `PythonModels/runs/testmodel_20260512_103000_123456/` + +该目录里通常会包含： + +- `testmodel_primary_series.csv` +- `testmodel_tank_temperature.csv` +- `testmodel_tank_temperature.svg` +- `testmodel_run_report.txt` +- `testmodel_modelica_comparison.csv` +- `testmodel_modelica_comparison_summary.txt` + +## 基线结果 + +当前基线对比摘要来自： +[testmodel_modelica_comparison_summary.txt](/home/lujz/projects/pressurization-transfer-system/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt) + +当前四个主变量的最大误差为： + +- `mytank.p`: `max_abs_error = 134.960857 Pa`, `max_rel_error = 0.006798%` +- `mytank.T`: `max_abs_error = 0.035507 K`, `max_rel_error = 0.009016%` +- `mycylinder.p`: `max_abs_error = 1391.986349 Pa`, `max_rel_error = 0.009447%` +- `mycylinder.T`: `max_abs_error = 0.009069 K`, `max_rel_error = 0.003870%` + +这说明在当前基线工况下，Python 版主变量已经能较好贴近 OpenModelica 结果。 + +## 当前架构判断 + +如果按“组件正确 -> 网络闭合 -> 积分可跑 -> 结果对齐 -> 去近似”来看，当前大致处于： + +- 组件级：已完成首版 +- 系统闭合：已完成首版 +- 积分入口：已完成首版 +- 基线结果对齐：已具备初步能力 +- 去近似：仍在进行中 + +所以当前最准确的说法不是“已完成移植”，而是： + +`Testmodel` 已有一版可运行、可导出、可对比的 Python 近似实现。 + +## 已知限制 + +当前最主要的限制可以直接理解成下面几条： + +- 介质模型已从常 `cp/cv` 推进到温度相关空气近似，但仍不是 `Modelica.Media.Air.SimpleAir` 的严格复刻。 +- 系统整体仍是 ODE 化近似，不是原始 Modelica DAE 的直接复现。 +- `mytee1 -> mytank` 这一段虽然已经去掉早期的“虚拟出口导通系数”，改成了基于压力一致性的下游能量闭合，但本质上仍是工程近似。 +- 当前 `Tee` 的 stream 语义只覆盖了当前 `Testmodel` 需要的最小集合，还不是通用的 `inStream/actualStream` 框架。 +- 当前一致初值仍是 ODE 入口处的约束投影，不等同于真正的 DAE 初始化求解。 +- 当前自动校验主要锁的是 Python 提交基线，还不是稳定的 Modelica 阈值回归。 +- 当前闭合器、系统层和 reporting 层虽然已经开始做“双支路结构化”，但对外结果序列、报告字段和部分导出命名仍然保留 `Testmodel` 专名兼容层，还没有完全转成通用表达。 +- 当前内核已经足够支撑下一阶段“结构化参数 -> 仿真执行 -> 产物输出”的链路开发，但还没有现成的 Word 参数解析入口和正式报告生成链路。 + +所以，当前版本适合： + +- 架构验证 +- 组件接口验证 +- 基线工况对比 +- 结果导出与误差定位 + +但当前版本还不适合： + +- 直接宣称与 OpenModelica 严格等价 +- 作为最终工程结论的唯一依据 +- 直接扩展到更复杂拓扑而不补通用连接器语义 + +## 文件级现状 + +按代码现状逐项看： + +- `core/base.py`: 正常 + 只提供最小抽象层，没有明显冗余。 +- `core/ports.py`: 正常 + `PortState` 目前只保留 `p`、`m_flow`、`h_outflow` 三个必要字段。 +- `core/state.py`: 正常 + `VolumeState` 只负责 `[m, U]` 状态打包。 +- `core/network.py`: 正常 + 负责状态向量拼装和连接摘要，不参与物理求解。 +- `core/solver.py`: 正常 + 已支持 SciPy、RK4 回退和零时长仿真。 +- `components/*.py`: 正常 + 都是当前一版近似模型，没有发现与 README 明显冲突的“未记录能力”。 +- `systems/testmodel.py`: 是当前最重要的技术债集中区 + 这里承载了下游流向切换、焓混合、压力投影等近似逻辑，后续演进应主要落在这里。 +- `scripts/run_testmodel.py`: 正常 + 已不是“最小打印脚本”，而是当前结果导出和对比入口。 +- `baselines/`: 是当前稳定基线，不应该随着日常运行频繁改动。 +- `runs/`: 是当前默认运行产物目录，不是手写源代码，也不应该当作提交基线使用。 + +## 当前主技术债 + +目前最主要的技术债，可以直接理解成下面 4 件事： + +1. 当前初始化虽然已经引入迭代诊断，但本质上仍是 ODE 入口近似，不是真正的 DAE 初始化器。 +2. `systems/testmodel.py` 还是承载了太多系统级闭合和初始化逻辑，只是主要端口的手写 stream 方向判断已经搬到组件 helper 里了，装配参数本身已经基本收口到配置对象。 +3. 自动校验现在主要锁的是 Python 这一版自己的基线，还不是稳定的 Modelica 阈值回归。 +4. 当前空气物性已经完成首轮基线校准，但还不是 `SimpleAir` 的严格复刻。以后如果换工况，或者拿到更多 Modelica 原始结果，参数大概率还要继续调。 diff --git a/PythonModels/__init__.py b/PythonModels/__init__.py new file mode 100644 index 0000000..a392a55 --- /dev/null +++ b/PythonModels/__init__.py @@ -0,0 +1,2 @@ +"""Python port scaffold for the Modelica-based pressurization system.""" + diff --git a/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt b/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt new file mode 100644 index 0000000..b8835bd --- /dev/null +++ b/PythonModels/baselines/testmodel/testmodel_modelica_comparison_summary.txt @@ -0,0 +1,4 @@ +mytank.p: max_abs_error=134.960858, max_rel_error=0.006798% +mytank.T: max_abs_error=0.035507, max_rel_error=0.009016% +mycylinder.p: max_abs_error=1391.986349, max_rel_error=0.009447% +mycylinder.T: max_abs_error=0.009069, max_rel_error=0.003870% diff --git a/PythonModels/baselines/testmodel/testmodel_primary_series.csv b/PythonModels/baselines/testmodel/testmodel_primary_series.csv new file mode 100644 index 0000000..f908440 --- /dev/null +++ b/PythonModels/baselines/testmodel/testmodel_primary_series.csv @@ -0,0 +1,202 @@ +time_s,mytank.p,mytank.T,mycylinder.p,mycylinder.T +0.0,99999.9999998181,299.9999999994543,35000000.0,300.0 +0.1,113767.94796194455,310.9637987101763,34857995.191961475,299.65190130981887 +0.2,127490.75529043324,320.1741038867186,34716455.97160761,299.3039342009557 +0.30000000000000004,141168.22832317703,327.9540327439692,34575384.336400226,298.9561065982938 +0.4,154800.2272891848,334.6028385113769,34434781.727959625,298.60842486666263 +0.5,168386.6833862148,340.3440160419025,34294648.85592305,298.26089368263627 +0.6000000000000001,181927.56886228317,345.34597934661645,34154986.006528884,297.9135169662631 +0.7000000000000001,195422.8867284584,349.7379443903134,34015793.14872051,297.5662980942661 +0.8,208872.67226727083,353.62071952544704,33877069.91858829,297.2192381374397 +0.9,222277.06297793236,357.0742683513725,33738814.897943415,296.8723276411647 +1.0,235635.98150118813,360.16219760966084,33601028.88468161,296.5255848498291 +1.1,248949.42783703818,362.9362624751423,33463711.878802836,296.1790180498007 +1.2000000000000002,262217.75812283135,365.43943715392686,33326860.20704978,295.8325994898035 +1.3,275440.80324992316,367.7064864598734,33190475.613635924,295.48635453607636 +1.4000000000000001,288618.55460473493,369.766650344722,33054558.187403098,295.14029246566224 +1.5,301751.3171281971,371.6448887937757,32919104.783141974,294.79439468638526 +1.6,314839.14341954247,373.36206261131525,32784114.858335685,294.44866485903736 +1.7000000000000002,327882.06907680153,374.93594217338176,32649588.04582047,294.10310819940065 +1.8,340880.20077679906,376.3818124575356,32515523.2453146,293.7577259967503 +1.9000000000000001,353833.6270518691,377.7128659999238,32381919.54368144,293.41251959012544 +2.0,366742.42291137256,378.94054912194076,32248776.16726251,293.06749036423287 +2.1,379606.6730674497,380.0748368029417,32116092.24232331,292.72263968200787 +2.2,392426.50707225554,381.12448095768656,31983866.432642274,292.37796667994184 +2.3000000000000003,405202.0671563814,382.0971825271288,31852097.271230437,292.03346983037886 +2.4000000000000004,417933.4155500976,382.9997115435689,31720784.11623488,291.6891515780112 +2.5,430620.6510589185,383.838081019396,31589925.948559783,291.3450125479824 +2.6,443263.883710829,384.61764092425494,31459521.633358993,291.0010529534012 +2.7,455863.2542884696,385.3431706778388,31329569.71857793,290.6572718762472 +2.8000000000000003,468418.84583467664,386.0189391423481,31200069.34769955,290.31367046787994 +2.9000000000000004,480930.76401047717,386.64878296917834,31071019.430919208,289.9702490545257 +3.0,493399.1147767602,387.2361564731079,30942418.8753394,289.627007973335 +3.1,505824.0120818939,387.7841699997202,30814266.50567852,289.283948167411 +3.2,518205.5585633492,388.29565044444143,30686561.263317347,288.9410696919777 +3.3000000000000003,530543.860941083,388.7731587646163,30559302.04752924,288.5983729053581 +3.4000000000000004,542839.025935053,389.21902406251183,30432487.75758757,288.2558581528899 +3.5,555091.1532282452,389.63535331051384,30306117.365346134,287.9135275289302 +3.6,567300.3523937837,390.0240928732544,30180189.74065028,287.5713806480376 +3.7,579466.730902255,390.3870196149787,30054703.77503127,287.2294176336731 +3.8000000000000003,591590.3951155862,390.7257605811333,29929658.371455234,286.88763887333175 +3.9000000000000004,603671.4513957056,391.0418105423231,29805052.432888325,286.5460447413264 +4.0,615710.0061045411,391.3365451815722,29680884.862296656,286.20463559857706 +4.1000000000000005,627706.1656040212,391.61123273483645,29557154.562646367,285.8634117923989 +4.2,639660.0356829006,391.86704273384044,29433860.4428154,285.5223738511199 +4.3,651571.716232906,392.1050410519721,29311001.472504564,285.1815242995298 +4.4,663441.3179459961,392.32624321858174,29188576.510019373,284.840861980108 +4.5,675268.9468771729,392.53157932311984,29066584.46149164,284.5003872946073 +4.6000000000000005,687054.7090814385,392.72191331429804,28945024.2330532,284.16010063101663 +4.7,698798.7106137947,392.8980487868277,28823894.730835862,283.8200023633402 +4.800000000000001,710501.057529244,393.06073417111315,28703194.86097143,283.48009285137505 +4.9,722161.8558827877,393.21066739655134,28582923.529591747,283.1403724404863 +5.0,733781.2088290841,393.3484951327475,28463079.672743235,282.8008422184954 +5.1000000000000005,745359.2175579917,393.47481852250473,28343662.24673756,282.4615037792963 +5.2,756895.9922909128,393.5902150621596,28224670.114733644,282.1223563449153 +5.300000000000001,768391.6386869826,393.6952183545109,28106102.186946325,281.78340031124986 +5.4,779846.2624053361,393.7903288312906,27987957.373590477,281.44463605959197 +5.5,791259.96910511,393.8760162992802,27870234.58488091,281.10606395639473 +5.6000000000000005,802632.8644454395,393.9527222554955,27752932.731032487,280.76768435303745 +5.7,813965.0540854601,394.02086199551263,27636050.72226007,280.4294975855874 +5.800000000000001,825256.6344365194,394.08081904264304,27519587.564161647,280.0915054160418 +5.9,836507.7093909897,394.13296077650165,27403542.18517475,279.75370843660113 +6.0,847718.3891245491,394.17763867488134,27287913.44892988,279.4161062418491 +6.1000000000000005,858888.7786008355,394.2151805066602,27172700.272815377,279.0786992215644 +6.2,870018.9827834871,394.2458959595685,27057901.57421954,278.7414877501898 +6.300000000000001,881109.1066361419,394.27007787647636,26943516.2705307,278.40447218658807 +6.4,892159.2551224378,394.2880033916485,26829543.279137183,278.0676528737956 +6.5,903169.5332060129,394.2999349762812,26715981.5174273,277.73103013877227 +6.6000000000000005,914139.9829839376,394.30609016406663,26602830.551205155,277.3946116451124 +6.7,925070.7482002805,394.30672506364067,26490088.897871558,277.0583932025842 +6.800000000000001,935961.9434794089,394.30206871636136,26377755.375172496,276.7223739769823 +6.9,946813.6729747849,394.2923336718555,26265828.9088526,276.3865543393284 +7.0,957626.040839869,394.2777219546483,26154308.424656466,276.05093464451556 +7.1000000000000005,968399.1512281233,394.2584257130283,26043192.84832872,275.7155152310525 +7.2,979133.1082930088,394.23462782052076,25932481.10561397,275.3802964208052 +7.300000000000001,989828.0161879869,394.2065024339623,25822172.12225682,275.04527851873627 +7.4,1000483.9790665191,394.17421551178546,25712264.824001882,274.71046181264023 +7.5,1011101.1010820667,394.13792529578427,25602758.136593778,274.3758465728772 +7.6000000000000005,1021679.4863880915,394.09778275932865,25493650.985777102,274.04143305210266 +7.7,1032219.239138054,394.0539320247227,25384942.297296483,273.707221484995 +7.800000000000001,1042720.4634854163,394.0065107521557,25276630.996896524,273.37321208797937 +7.9,1053183.2635836392,393.95565050247825,25168716.01032184,273.0394050589494 +8.0,1063607.7435861847,393.9014770758345,25061196.263317037,272.705800576985 +8.1,1073994.0076465139,393.8441108280058,24954070.681626726,272.3723988020679 +8.200000000000001,1084342.0075784405,393.78360959211744,24847339.76225091,272.03921728278755 +8.3,1094651.950514763,393.7201237653114,24741001.368788917,271.7062444960846 +8.4,1104923.9566483083,393.65376548392067,24635054.261552785,271.3734787513138 +8.5,1115158.1282653515,393.5846357175065,24529497.385545585,271.0409203804046 +8.6,1125354.567652169,393.5128313057523,24424329.685770374,270.7085696976142 +8.700000000000001,1135513.3770950353,393.43844517133164,24319550.10723021,270.37642699925055 +8.8,1145634.6588802272,393.3615665197664,24215157.594928175,270.0444925633909 +8.9,1155718.5152940191,393.2822810271899,24111151.093867306,269.712766649598 +9.0,1165765.0486226876,393.2006710168648,24007529.54905069,269.3812494986329 +9.1,1175774.361152508,393.11681562522904,23904291.90548136,269.04994133216485 +9.200000000000001,1185746.555169756,393.03079095819226,23801437.108162407,268.7188423524781 +9.3,1195681.732960707,392.94267023834027,23698964.10209688,268.3879527421748 +9.4,1205579.9968116365,392.8525239436615,23596871.83228785,268.05727266387584 +9.5,1215441.4490088206,392.7604199383583,23495159.243738372,267.72680225991735 +9.600000000000001,1225266.1918385345,392.66642359626576,23393825.281451505,267.39654165204473 +9.700000000000001,1235054.3275870536,392.57059791736026,23292868.89043032,267.0664909411029 +9.8,1244805.7424247153,392.47293675670085,23192291.244732097,266.7366751290975 +9.9,1254520.6963966596,392.3735487251935,23092089.662209835,266.40707625875064 +10.0,1264199.3049931854,392.2724947546991,22992262.951678198,266.07769289844276 +10.100000000000001,1273841.6679760527,392.1698289153906,22892810.0841785,265.74852540655036 +10.200000000000001,1283447.8851070204,392.0656034102887,22793730.030752078,265.41957412301247 +10.3,1293018.0561478487,391.9598686580462,22695021.76244025,265.0908393690368 +10.4,1302552.2808602974,391.8526733713811,22596684.250284337,264.76232144680336 +10.5,1312050.6590061258,391.7440646314241,22498716.465325654,264.4340206391643 +10.600000000000001,1321513.290347094,391.63408795822716,22401117.37860553,264.10593720934105 +10.700000000000001,1330940.2746449616,391.5227873776613,22303885.961165283,263.7780714006174 +10.8,1340331.711661488,391.41020548491883,22207021.184046242,263.4504234360302 +10.9,1349687.7011584332,391.29638350481804,22110522.01828973,263.1229935180553 +11.0,1359008.342897557,391.18136134909884,22014387.434937052,262.79578182829164 +11.100000000000001,1368293.7366406189,391.0651776708804,21918616.405029543,262.4687885271409 +11.200000000000001,1377543.9821493784,390.9478699164463,21823207.89960853,262.1420137534841 +11.3,1386759.1791855951,390.82947437450605,21728160.88971532,261.81545762435434 +11.4,1395939.2416852904,390.70997702085253,21633476.26302751,261.4891421257888 +11.5,1405084.4058089943,390.5894488404711,21539151.583747786,261.1630514204039 +11.600000000000001,1414194.7727817593,390.46792355037144,21445185.807824824,260.83718544108797 +11.700000000000001,1423270.4378832965,390.3454323238288,21351577.952528503,260.5115448028764 +11.8,1432311.4962235806,390.22200536933536,21258327.036879413,260.1861301247683 +11.9,1441318.042742851,390.0976719679259,21165432.0816488,259.86094202976824 +12.0,1450290.1722116095,389.97246050877925,21072892.1093586,259.5359811449295 +12.100000000000001,1459227.9792306225,389.84639852319333,20980706.14428146,259.2112481013979 +12.200000000000001,1468131.5582309205,389.71951271701596,20888873.212440677,258.88674353445447 +12.3,1477001.0034737969,389.59182900161403,20797392.34161027,258.5624680835606 +12.4,1485836.4090508097,389.46337252345796,20706262.561314918,258.23842239240224 +12.5,1494637.8688837802,389.3341676923915,20615482.902829994,257.9146071089349 +12.600000000000001,1503405.4767247946,389.20423820865625,20525052.399181567,257.5910228854296 +12.700000000000001,1512139.3261562001,389.07360708873284,20434970.085146382,257.26767037851835 +12.8,1520839.5105906113,388.94229669006046,20345234.997251876,256.9445502492407 +12.9,1529506.1233409408,388.81032876227204,20255846.173053782,256.62166315000906 +13.0,1538139.2687433015,388.67772882293264,20166802.53641284,256.2990076411302 +13.100000000000001,1546739.030251722,388.5445136335602,20078103.22657088,255.97658615427704 +13.200000000000001,1555305.5007360894,388.4107031827765,19989747.285653118,255.65439933992917 +13.3,1563838.7729005008,388.27631690939495,19901733.757494725,255.332447852161 +13.4,1572338.9392832702,388.14137372147235,19814061.687640797,255.0107323486784 +13.5,1580806.0922569225,388.0058920145713,19726730.123346385,254.68925349085504 +13.600000000000001,1589240.3240281995,387.86988968927403,19639738.113576483,254.36801194376957 +13.700000000000001,1597641.7266380545,387.73338416798015,19553084.709006038,254.04700837624307 +13.8,1606010.3919616563,387.59639241102605,19466768.962019928,253.72624346087662 +13.9,1614346.4117083861,387.458930932153,19380789.926712975,253.40571787408948 +14.0,1622649.8774218408,387.32101581335974,19295146.658889957,253.0854322961577 +14.100000000000001,1630920.8804798292,387.18266271916264,19209838.21606559,252.76538741125262 +14.200000000000001,1639159.5120943757,387.0438869102949,19124863.657464534,252.4455839074805 +14.3,1647365.8633117182,386.90470325686863,19040222.0440214,252.12602247692183 +14.4,1655540.0250123062,386.76512625102305,18955912.438380726,251.8067038156715 +14.5,1663682.087910807,386.6251700190859,18871933.904897016,251.48762862387906 +14.600000000000001,1671792.1425560997,386.4848483332645,18788285.509634707,251.16879760578948 +14.700000000000001,1679870.2793312764,386.3441746228917,18704966.320368182,250.85021146978437 +14.8,1687916.588453644,386.2031619852454,18621975.406581767,250.5318709284234 +14.9,1695931.1599747243,386.06182319595837,18539311.83946974,250.21377669848638 +15.0,1703914.0837802505,385.920170719039,18456974.691936314,249.8959295010154 +15.100000000000001,1711865.4494154856,385.77821670970445,18374963.040397402,249.57833006564294 +15.200000000000001,1719785.3379250832,385.6359727064743,18293276.048945524,249.26097933285033 +15.3,1727673.8459201432,385.4934505781353,18211912.721118417,248.94387786218365 +15.4,1735531.0625176337,385.3506616315768,18130872.137749474,248.6270263958575 +15.5,1743357.0766656764,385.2076169067951,18050153.381413613,248.31042568075102 +15.600000000000001,1751151.9771435438,385.06432718490083,17969755.53642727,247.99407646845765 +15.700000000000001,1758915.8525616606,384.9208029958382,17889677.688848406,247.67797951533504 +15.8,1766648.7913616048,384.77705462583015,17809918.926476505,247.36213558255577 +15.9,1774350.881816107,384.63309212455715,17730478.33885257,247.04654543615806 +16.0,1782022.2120290494,384.4889253120841,17651355.017259125,246.73120984709743 +16.1,1789662.8699354655,384.3445637855432,17572548.05472022,246.41612959129878 +16.2,1797272.9433015438,384.2000169255848,17494056.546001427,246.10130544970855 +16.3,1804852.5197246224,384.05529390260284,17415879.58760983,245.78673820834797 +16.400000000000002,1812401.6866331936,383.9104036827472,17338016.277794052,245.47242865836654 +16.5,1819920.5312869006,383.76535503372884,17260465.716544226,245.15837759609565 +16.6,1827409.1407765402,383.6201565304272,17183227.005592,244.84458582310356 +16.7,1834867.6020240602,383.47481656030743,17106299.24841057,244.53105414625009 +16.8,1842296.0017825624,383.32934332865534,17029681.550214626,244.2177833777421 +16.900000000000002,1849694.4266362996,383.1837448636361,16953373.017960392,243.9047743351897 +17.0,1857062.9630006768,383.03802902118434,16877372.76034562,243.59202784166294 +17.1,1864401.697122252,382.8922034897328,16801679.887809563,243.27954472574828 +17.2,1871710.7150787353,382.7462757947841,16726293.512533028,242.96732582160686 +17.3,1878990.1030524147,382.60025337999576,16651212.745618157,242.65537190477806 +17.400000000000002,1886239.954455942,382.4541455990818,16576436.623591991,242.343682012199 +17.5,1893460.3482686526,382.3079576538116,16501964.331849081,242.03225853722532 +17.6,1900651.3699906773,382.16169648335216,16427794.988527464,241.72110231102408 +17.7,1907813.104956132,382.0153688850827,16353927.713477474,241.41021416968596 +17.8,1914945.6383331183,381.86898151834066,16280361.628261749,241.09959495427398 +17.900000000000002,1922049.0551237254,381.7225409080495,16207095.856155202,240.78924551087363 +18.0,1929123.4401640275,381.57605344823065,16134129.522145053,240.47916669064318 +18.1,1936168.8781240864,381.4295254054061,16061461.752930798,240.16935934986387 +18.2,1943185.4535079484,381.28296292189384,15989091.67692425,239.8598243499916 +18.3,1950173.2506536485,381.1363720190009,15917018.424249483,239.55056255770802 +18.400000000000002,1957132.3537332045,380.98975860011734,15845241.126742886,239.24157484497272 +18.5,1964062.8467526236,380.84312845371574,15773758.917953137,238.93286208907583 +18.6,1970964.8135518986,380.6964872562572,15702570.933141202,238.6244251726908 +18.7,1977838.3378050062,380.54984057500997,15631676.309280336,238.31626498392774 +18.8,1984683.5030199126,380.4031938707821,15561074.185056096,238.00838241638746 +18.900000000000002,1991500.3925385692,380.25655250057133,15490763.700866321,237.7007783692156 +19.0,1998289.0895369116,380.10992172013556,15420743.998821149,237.39345374715774 +19.1,2005049.677024864,379.96330668648653,15351014.222743012,237.08640946061425 +19.200000000000003,2011782.237846337,379.81671246030885,15281573.518166626,236.77964642569634 +19.3,2018486.8546792252,379.67014400830794,15212421.032339014,236.47316556428217 +19.400000000000002,2025163.610035411,379.5236062054879,15143555.914219463,236.16696780407335 +19.5,2031812.5863479478,379.3771038539995,15074977.31358037,235.861054060844 +19.6,2038433.8680519294,379.2306420805889,15006684.359544702,235.55542481191287 +19.700000000000003,2045027.536485746,379.084225365391,14938676.213175347,235.25008113730536 +19.8,2051593.6739729291,378.9378583245209,14870952.025374293,234.9450238874766 +19.900000000000002,2058132.3627231135,378.79154550141357,14803510.948218277,234.64025390703284 +20.0,2064643.6848320398,378.64529136859073,14736352.134958768,234.33577203462033 diff --git a/PythonModels/components/__init__.py b/PythonModels/components/__init__.py new file mode 100644 index 0000000..d8e493c --- /dev/null +++ b/PythonModels/components/__init__.py @@ -0,0 +1,2 @@ +"""Component implementations for the Python system model.""" + diff --git a/PythonModels/components/cylinder.py b/PythonModels/components/cylinder.py new file mode 100644 index 0000000..066fb27 --- /dev/null +++ b/PythonModels/components/cylinder.py @@ -0,0 +1,55 @@ +from __future__ import annotations + +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +class Cylinder(DynamicComponent): + """Python port of ModelicaModels.Mycylinder.""" + + def __init__( + self, + name: str, + medium: IdealGasMedium, + V: float = 0.01, + p0: float = 35e6, + T0: float = 300.0, + ) -> None: + super().__init__(name=name) + self.medium = medium + self.V = V + m0 = p0 * V / (medium.R_gas * T0) + U0 = m0 * medium.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_b = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + props = self.medium.properties_from_mU(self.state.m, self.state.U, self.V) + self.port_b.p = props.p + self.port_b.h_outflow = props.h + return props + + def derivatives_from_connection( + self, + *, + connected_h: float, + port_m_flow: float, + internal_h: float, + ) -> VolumeState: + inlet_h = self.connection_inlet_enthalpy( + port_m_flow=port_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + return self.derivatives(inlet_h, port_m_flow) + + def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState: + return VolumeState(m=m_flow, U=m_flow * inlet_h) diff --git a/PythonModels/components/orifice.py b/PythonModels/components/orifice.py new file mode 100644 index 0000000..9da6282 --- /dev/null +++ b/PythonModels/components/orifice.py @@ -0,0 +1,28 @@ +from __future__ import annotations + +from math import sqrt + +from PythonModels.core.base import AlgebraicComponent +from PythonModels.core.ports import PortState + + +class Orifice(AlgebraicComponent): + """Python port of ModelicaModels.Myorifice.""" + + def __init__(self, name: str, opening: float = 1.0, K: float = 1e-7) -> None: + super().__init__(name=name) + self.opening = opening + self.K = K + self.port_a = PortState() + self.port_b = PortState() + + @property + def K_eff(self) -> float: + return self.K * max(self.opening, 0.001) + + def mass_flow(self, p_a: float, p_b: float) -> float: + dp = p_a - p_b + if dp == 0.0: + return 0.0 + return self.K_eff * sqrt(abs(dp)) * (1.0 if dp > 0.0 else -1.0) + diff --git a/PythonModels/components/pipe.py b/PythonModels/components/pipe.py new file mode 100644 index 0000000..6fb9fe5 --- /dev/null +++ b/PythonModels/components/pipe.py @@ -0,0 +1,133 @@ +from __future__ import annotations + +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +class Pipe(DynamicComponent): + """Python port of ModelicaModels.Mypipe.""" + + def __init__( + self, + name: str, + medium: IdealGasMedium, + L: float = 5.0, + D: float = 0.02, + lambda_darcy: float = 0.02, + p0: float = 1e5, + T0: float = 300.0, + ) -> None: + super().__init__(name=name) + self.medium = medium + self.L = L + self.D = D + self.lambda_darcy = lambda_darcy + self.area = 3.141592653589793 * D * D / 4.0 + self.V = self.area * L + m0 = p0 * self.V / (medium.R_gas * T0) + U0 = m0 * medium.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_a = PortState() + self.port_b = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + props = self.medium.properties_from_mU(self.state.m, self.state.U, self.V) + self.port_b.p = props.p + self.port_a.h_outflow = props.h + self.port_b.h_outflow = props.h + return props + + def inlet_pressure(self, m_flow_a: float, rho: float, core_pressure: float) -> float: + resistance = self.lambda_darcy * (self.L / self.D) + dynamic_term = m_flow_a * abs(m_flow_a) / (2.0 * rho * self.area * self.area) + return core_pressure + resistance * dynamic_term + + def port_a_inlet_enthalpy( + self, + *, + port_a_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + return self.connection_inlet_enthalpy( + port_m_flow=port_a_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + + def port_b_inlet_enthalpy( + self, + *, + port_b_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + return self.connection_inlet_enthalpy( + port_m_flow=port_b_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + + def connection_inlet_enthalpies( + self, + *, + port_a_m_flow: float, + connected_h_a: float, + port_b_m_flow: float, + connected_h_b: float, + internal_h: float, + ) -> tuple[float, float]: + return ( + self.port_a_inlet_enthalpy( + port_a_m_flow=port_a_m_flow, + connected_h=connected_h_a, + internal_h=internal_h, + ), + self.port_b_inlet_enthalpy( + port_b_m_flow=port_b_m_flow, + connected_h=connected_h_b, + internal_h=internal_h, + ), + ) + + def derivatives_from_connections( + self, + *, + port_a_m_flow: float, + connected_h_a: float, + port_b_m_flow: float, + connected_h_b: float, + internal_h: float, + ) -> VolumeState: + inlet_h_a, inlet_h_b = self.connection_inlet_enthalpies( + port_a_m_flow=port_a_m_flow, + connected_h_a=connected_h_a, + port_b_m_flow=port_b_m_flow, + connected_h_b=connected_h_b, + internal_h=internal_h, + ) + return self.derivatives( + inlet_h_a=inlet_h_a, + inlet_h_b=inlet_h_b, + m_flow_a=port_a_m_flow, + m_flow_b=port_b_m_flow, + ) + + def derivatives( + self, + inlet_h_a: float, + inlet_h_b: float, + m_flow_a: float, + m_flow_b: float, + ) -> VolumeState: + dm_dt = m_flow_a + m_flow_b + dU_dt = m_flow_a * inlet_h_a + m_flow_b * inlet_h_b + return VolumeState(m=dm_dt, U=dU_dt) diff --git a/PythonModels/components/tank.py b/PythonModels/components/tank.py new file mode 100644 index 0000000..8f51d69 --- /dev/null +++ b/PythonModels/components/tank.py @@ -0,0 +1,55 @@ +from __future__ import annotations + +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +class Tank(DynamicComponent): + """Python port of ModelicaModels.Mytank.""" + + def __init__( + self, + name: str, + medium: IdealGasMedium, + V: float = 0.1, + p0: float = 1e5, + T0: float = 300.0, + ) -> None: + super().__init__(name=name) + self.medium = medium + self.V = V + m0 = p0 * V / (medium.R_gas * T0) + U0 = m0 * medium.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_a = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + props = self.medium.properties_from_mU(self.state.m, self.state.U, self.V) + self.port_a.p = props.p + self.port_a.h_outflow = props.h + return props + + def derivatives_from_connection( + self, + *, + connected_h: float, + port_m_flow: float, + internal_h: float, + ) -> VolumeState: + inlet_h = self.connection_inlet_enthalpy( + port_m_flow=port_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + return self.derivatives(inlet_h, port_m_flow) + + def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState: + return VolumeState(m=m_flow, U=m_flow * inlet_h) diff --git a/PythonModels/components/tee.py b/PythonModels/components/tee.py new file mode 100644 index 0000000..6602435 --- /dev/null +++ b/PythonModels/components/tee.py @@ -0,0 +1,172 @@ +from __future__ import annotations + +from PythonModels.core.base import AlgebraicComponent +from PythonModels.core.ports import PortState + + +class Tee(AlgebraicComponent): + """Python port of ModelicaModels.Mytee.""" + + def __init__(self, name: str) -> None: + super().__init__(name=name) + self.port_in = PortState() + self.port_out1 = PortState() + self.port_out2 = PortState() + + def mixed_inlet_enthalpy( + self, + branch1_m_flow: float, + branch1_h: float, + branch2_m_flow: float, + branch2_h: float, + fallback_h: float = 0.0, + ) -> float: + positive_1 = max(branch1_m_flow, 0.0) + positive_2 = max(branch2_m_flow, 0.0) + total = positive_1 + positive_2 + if total <= 1e-9: + return fallback_h + return (positive_1 * branch1_h + positive_2 * branch2_h) / total + + def inlet_stream_enthalpy( + self, + branch1_m_flow: float, + branch1_h: float, + branch2_m_flow: float, + branch2_h: float, + fallback_h: float, + ) -> float: + """Approximate `inStream(port_in.h_outflow)` for the current tee topology.""" + + return self.mixed_inlet_enthalpy( + branch1_m_flow, + branch1_h, + branch2_m_flow, + branch2_h, + fallback_h=fallback_h, + ) + + def branch_actual_stream_enthalpy( + self, + branch_m_flow: float, + branch_h: float, + inlet_h: float, + ) -> float: + """Approximate `actualStream(branch.h_outflow)` for a tee branch port.""" + + return inlet_h if branch_m_flow > 0.0 else branch_h + + @staticmethod + def _solve_linear_2x2( + a11: float, + a12: float, + a21: float, + a22: float, + b1: float, + b2: float, + ) -> tuple[float, float] | None: + determinant = a11 * a22 - a12 * a21 + if abs(determinant) <= 1e-12: + return None + x1 = (b1 * a22 - b2 * a12) / determinant + x2 = (a11 * b2 - a21 * b1) / determinant + return x1, x2 + + def solve_branch_outlet_flows_from_energy_balance( + self, + *, + ratio_branch1: float, + ratio_branch2: float, + inlet_h_branch1: float, + inlet_h_branch2: float, + branch1_h: float, + branch2_h: float, + inlet_h: float, + q_in_branch1: float, + q_in_branch2: float, + tolerance: float = 1e-12, + ) -> tuple[float, float]: + """Solve branch outlet flows for the current three-port downstream tee use-case.""" + + rhs_branch1 = q_in_branch1 * inlet_h_branch1 + rhs_branch2 = q_in_branch2 * inlet_h_branch2 + + def solve_both_forward() -> tuple[float, float] | None: + return self._solve_linear_2x2( + (1.0 + ratio_branch1) * branch1_h, + ratio_branch1 * branch2_h, + ratio_branch2 * branch1_h, + (1.0 + ratio_branch2) * branch2_h, + rhs_branch1, + rhs_branch2, + ) + + def solve_one_reverse( + *, + branch1_reverse: bool, + ) -> tuple[float, float] | None: + if branch1_reverse: + return self._solve_linear_2x2( + inlet_h * (1.0 + ratio_branch1), + ratio_branch1 * inlet_h, + ratio_branch2 * inlet_h, + branch2_h + ratio_branch2 * inlet_h, + rhs_branch1, + rhs_branch2, + ) + + return self._solve_linear_2x2( + branch1_h + ratio_branch1 * inlet_h, + ratio_branch1 * inlet_h, + ratio_branch2 * inlet_h, + inlet_h * (1.0 + ratio_branch2), + rhs_branch1, + rhs_branch2, + ) + + def solve_both_reverse() -> tuple[float, float] | None: + return self._solve_linear_2x2( + inlet_h * (1.0 + ratio_branch1), + ratio_branch1 * inlet_h, + ratio_branch2 * inlet_h, + inlet_h * (1.0 + ratio_branch2), + rhs_branch1, + rhs_branch2, + ) + + candidate_solvers = ( + ( + solve_both_forward, + lambda q1, q2: q1 >= -tolerance and q2 >= -tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=True), + lambda q1, q2: q1 < -tolerance and q2 >= -tolerance and q1 + q2 > tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=True), + lambda q1, q2: q1 < -tolerance and q2 >= -tolerance and q1 + q2 <= tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=False), + lambda q1, q2: q2 < -tolerance and q1 >= -tolerance and q1 + q2 > tolerance, + ), + ( + lambda: solve_one_reverse(branch1_reverse=False), + lambda q1, q2: q2 < -tolerance and q1 >= -tolerance and q1 + q2 <= tolerance, + ), + ( + solve_both_reverse, + lambda q1, q2: q1 < -tolerance and q2 < -tolerance, + ), + ) + + for solver, predicate in candidate_solvers: + candidate = solver() + if candidate is None: + continue + q_out_branch1, q_out_branch2 = candidate + if predicate(q_out_branch1, q_out_branch2): + return q_out_branch1, q_out_branch2 + + return solve_both_forward() or (0.0, 0.0) diff --git a/PythonModels/core/__init__.py b/PythonModels/core/__init__.py new file mode 100644 index 0000000..ae12d41 --- /dev/null +++ b/PythonModels/core/__init__.py @@ -0,0 +1,2 @@ +"""Core abstractions for the Python system model.""" + diff --git a/PythonModels/core/base.py b/PythonModels/core/base.py new file mode 100644 index 0000000..bac666e --- /dev/null +++ b/PythonModels/core/base.py @@ -0,0 +1,48 @@ +from __future__ import annotations + +from abc import ABC, abstractmethod + + +class Component(ABC): + def __init__(self, name: str) -> None: + self.name = name + + +class DynamicComponent(Component): + state_size = 2 + + @staticmethod + def actual_stream_enthalpy( + port_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + """Approximate `actualStream(port.h_outflow)` for a mixed control volume port.""" + + return connected_h if port_m_flow > 0.0 else internal_h + + def connection_inlet_enthalpy( + self, + port_m_flow: float, + connected_h: float, + internal_h: float, + ) -> float: + """Resolve the enthalpy convected into this control volume through one port.""" + + return self.actual_stream_enthalpy( + port_m_flow=port_m_flow, + connected_h=connected_h, + internal_h=internal_h, + ) + + @abstractmethod + def get_state_vector(self) -> list[float]: + raise NotImplementedError + + @abstractmethod + def set_state_vector(self, values: list[float]) -> None: + raise NotImplementedError + + +class AlgebraicComponent(Component): + """Stateless element described by algebraic constraints only.""" diff --git a/PythonModels/core/medium.py b/PythonModels/core/medium.py new file mode 100644 index 0000000..3410590 --- /dev/null +++ b/PythonModels/core/medium.py @@ -0,0 +1,96 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass(frozen=True) +class ThermodynamicProperties: + p: float + T: float + rho: float + u: float + h: float + + +@dataclass(frozen=True) +class IdealGasMedium: + """Temperature-dependent ideal-gas air approximation. + + This is still not a strict clone of `Modelica.Media.Air.SimpleAir`. + The small linear `cp(T)` term is kept configurable for calibration, but the + current default is calibrated against the committed Testmodel baseline and + therefore falls back to the constant-heat-capacity limit. + """ + + name: str = "SimpleAirApprox" + R_gas: float = 287.0 + cp_ref: float = 1005.0 + T_ref: float = 300.0 + cp_slope: float = 0.0 + + @property + def cv(self) -> float: + return self.cv_at_temperature(self.T_ref) + + @property + def gamma(self) -> float: + return self.cp_at_temperature(self.T_ref) / self.cv + + def cp_at_temperature(self, T: float) -> float: + return self.cp_ref + self.cp_slope * (T - self.T_ref) + + def cv_at_temperature(self, T: float) -> float: + return self.cp_at_temperature(T) - self.R_gas + + def density(self, p: float, T: float) -> float: + return p / (self.R_gas * T) + + def specific_internal_energy(self, T: float) -> float: + delta_T = T - self.T_ref + return ( + self.cv * self.T_ref + + self.cv * delta_T + + 0.5 * self.cp_slope * delta_T * delta_T + ) + + def specific_enthalpy(self, T: float) -> float: + delta_T = T - self.T_ref + return ( + self.cp_ref * self.T_ref + + self.cp_ref * delta_T + + 0.5 * self.cp_slope * delta_T * delta_T + ) + + def temperature_from_internal_energy(self, u: float) -> float: + reference_internal_energy = self.cv * self.T_ref + delta_u = u - reference_internal_energy + + if abs(self.cp_slope) <= 1e-15: + return self.T_ref + delta_u / self.cv + + a = 0.5 * self.cp_slope + b = self.cv + c = -delta_u + discriminant = max(b * b - 4.0 * a * c, 0.0) + positive_root = (-b + discriminant**0.5) / (2.0 * a) + negative_root = (-b - discriminant**0.5) / (2.0 * a) + delta_T = positive_root if abs(positive_root) <= abs(negative_root) else negative_root + return self.T_ref + delta_T + + def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: + if m <= 0.0: + raise ValueError("Mass must stay positive when recovering temperature.") + return self.temperature_from_internal_energy(U / m) + + def pressure(self, m: float, T: float, V: float) -> float: + if V <= 0.0: + raise ValueError("Volume must stay positive.") + return m * self.R_gas * T / V + + def properties_from_mU(self, m: float, U: float, V: float) -> ThermodynamicProperties: + T = self.temperature_from_mass_internal_energy(m, U) + p = self.pressure(m, T, V) + rho = m / V + u = U / m + h = self.specific_enthalpy(T) + return ThermodynamicProperties(p=p, T=T, rho=rho, u=u, h=h) diff --git a/PythonModels/core/network.py b/PythonModels/core/network.py new file mode 100644 index 0000000..e912867 --- /dev/null +++ b/PythonModels/core/network.py @@ -0,0 +1,78 @@ +from __future__ import annotations + +from dataclasses import dataclass + +from PythonModels.core.base import Component, DynamicComponent + + +@dataclass(frozen=True) +class Connection: + source_component: str + source_port: str + target_component: str + target_port: str + + +class SimulationNetwork: + """Container for components, topology, and state-vector bookkeeping.""" + + def __init__(self, name: str) -> None: + self.name = name + self.components: dict[str, Component] = {} + self.connections: list[Connection] = [] + + def add_component(self, component: Component) -> None: + if component.name in self.components: + raise ValueError(f"Duplicate component name: {component.name}") + self.components[component.name] = component + + def connect( + self, + source_component: str, + source_port: str, + target_component: str, + target_port: str, + ) -> None: + self.connections.append( + Connection( + source_component=source_component, + source_port=source_port, + target_component=target_component, + target_port=target_port, + ) + ) + + def dynamic_components(self) -> list[DynamicComponent]: + return [ + component + for component in self.components.values() + if isinstance(component, DynamicComponent) + ] + + def initial_state_vector(self) -> list[float]: + values: list[float] = [] + for component in self.dynamic_components(): + values.extend(component.get_state_vector()) + return values + + def apply_state_vector(self, values: list[float]) -> None: + cursor = 0 + for component in self.dynamic_components(): + next_cursor = cursor + component.state_size + component.set_state_vector(values[cursor:next_cursor]) + cursor = next_cursor + if cursor != len(values): + raise ValueError("State vector length does not match dynamic components.") + + def summary(self) -> str: + lines = [f"Network: {self.name}", "Components:"] + for name, component in self.components.items(): + lines.append(f" - {name}: {component.__class__.__name__}") + lines.append("Connections:") + for conn in self.connections: + lines.append( + f" - {conn.source_component}.{conn.source_port}" + f" -> {conn.target_component}.{conn.target_port}" + ) + return "\n".join(lines) + diff --git a/PythonModels/core/ports.py b/PythonModels/core/ports.py new file mode 100644 index 0000000..d340db2 --- /dev/null +++ b/PythonModels/core/ports.py @@ -0,0 +1,13 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass +class PortState: + """Python-side analogue of a Modelica fluid port.""" + + p: float = 0.0 + m_flow: float = 0.0 + h_outflow: float = 0.0 + diff --git a/PythonModels/core/solver.py b/PythonModels/core/solver.py new file mode 100644 index 0000000..3e101e1 --- /dev/null +++ b/PythonModels/core/solver.py @@ -0,0 +1,102 @@ +from __future__ import annotations + +from dataclasses import dataclass +from typing import Callable + + +@dataclass(frozen=True) +class SolveIVPConfig: + t_start: float = 0.0 + t_stop: float = 20.0 + method: str = "BDF" + rtol: float = 1e-6 + atol: float = 1e-8 + max_step: float = 1e-3 + + +@dataclass(frozen=True) +class ODESolution: + t: list[float] + y: list[list[float]] + success: bool + message: str + + +def _vector_add(a: list[float], b: list[float], scale: float = 1.0) -> list[float]: + return [x + scale * y for x, y in zip(a, b)] + + +def _runge_kutta_4( + rhs: Callable[[float, list[float]], list[float]], + initial_state: list[float], + config: SolveIVPConfig, + t_eval: list[float] | None, +) -> ODESolution: + if t_eval is None: + point_count = max( + 2, + int((config.t_stop - config.t_start) / max(config.max_step, 1e-6)) + 1, + ) + step = (config.t_stop - config.t_start) / (point_count - 1) + t_eval = [config.t_start + index * step for index in range(point_count)] + + state = list(initial_state) + states = [[value] for value in state] + times = [float(t_eval[0])] + current_time = float(t_eval[0]) + + for target_time in t_eval[1:]: + while current_time < target_time - 1e-15: + dt = min(config.max_step, target_time - current_time) + k1 = rhs(current_time, state) + k2 = rhs(current_time + 0.5 * dt, _vector_add(state, k1, 0.5 * dt)) + k3 = rhs(current_time + 0.5 * dt, _vector_add(state, k2, 0.5 * dt)) + k4 = rhs(current_time + dt, _vector_add(state, k3, dt)) + state = [ + value + (dt / 6.0) * (a + 2.0 * b + 2.0 * c + d) + for value, a, b, c, d in zip(state, k1, k2, k3, k4) + ] + current_time += dt + + times.append(float(target_time)) + for index, value in enumerate(state): + states[index].append(value) + + return ODESolution( + t=times, + y=states, + success=True, + message="Integrated with built-in RK4 fallback because SciPy is unavailable.", + ) + + +def integrate_ode( + rhs: Callable[[float, list[float]], list[float]], + initial_state: list[float], + config: SolveIVPConfig, + t_eval: list[float] | None = None, +): + """Thin wrapper around scipy.integrate.solve_ivp with a pure-Python fallback.""" + + if abs(config.t_stop - config.t_start) <= 1e-15: + return ODESolution( + t=[float(config.t_start)], + y=[[value] for value in initial_state], + success=True, + message="Skipped integration because t_start equals t_stop.", + ) + + try: + from scipy.integrate import solve_ivp + except ImportError: + return _runge_kutta_4(rhs, initial_state, config, t_eval) + + return solve_ivp( + fun=rhs, + t_span=(config.t_start, config.t_stop), + y0=initial_state, + method=config.method, + rtol=config.rtol, + atol=config.atol, + t_eval=t_eval, + ) diff --git a/PythonModels/core/state.py b/PythonModels/core/state.py new file mode 100644 index 0000000..e4cf3b3 --- /dev/null +++ b/PythonModels/core/state.py @@ -0,0 +1,21 @@ +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass +class VolumeState: + """Primary dynamic state for rigid adiabatic control volumes.""" + + m: float + U: float + + def as_vector(self) -> list[float]: + return [self.m, self.U] + + @classmethod + def from_vector(cls, values: list[float]) -> "VolumeState": + if len(values) != 2: + raise ValueError("VolumeState requires exactly two values: [m, U].") + return cls(m=values[0], U=values[1]) + diff --git a/PythonModels/reporting/__init__.py b/PythonModels/reporting/__init__.py new file mode 100644 index 0000000..9ca6880 --- /dev/null +++ b/PythonModels/reporting/__init__.py @@ -0,0 +1,23 @@ +from PythonModels.reporting.testmodel_outputs import ( + COMPARISON_KEYS, + MODELICA_COMPARISON_COLUMNS, + PRIMARY_KEYS, + TestModelArtifacts, + export_testmodel_artifacts, + format_testmodel_run_report, + load_modelica_series, + write_testmodel_run_report, + write_modelica_comparison, +) + +__all__ = [ + "COMPARISON_KEYS", + "MODELICA_COMPARISON_COLUMNS", + "PRIMARY_KEYS", + "TestModelArtifacts", + "export_testmodel_artifacts", + "format_testmodel_run_report", + "load_modelica_series", + "write_testmodel_run_report", + "write_modelica_comparison", +] diff --git a/PythonModels/reporting/testmodel_outputs.py b/PythonModels/reporting/testmodel_outputs.py new file mode 100644 index 0000000..9338b1c --- /dev/null +++ b/PythonModels/reporting/testmodel_outputs.py @@ -0,0 +1,393 @@ +from __future__ import annotations + +from bisect import bisect_left +import csv +from dataclasses import dataclass +from pathlib import Path +from typing import Any + + +PRIMARY_KEYS = ( + "mytank.p", + "mytank.T", + "mycylinder.p", + "mycylinder.T", +) + +MODELICA_COMPARISON_COLUMNS = { + "mytank.p": "mytank.p", + "mytank.T": "mytank.T", + "mycylinder.p": "mycylinder.p", + "mycylinder.T": "mycylinder.T", + "branch.upper_branch.p": "mypipe.p", + "branch.upper_branch.in": "myorifice.port_a.m_flow", + "branch.upper_branch.out": "mytee1.port_out2.m_flow", + "branch.lower_branch.p": "mypipe1.p", + "branch.lower_branch.in": "myorifice1.port_a.m_flow", + "branch.lower_branch.out": "mytee1.port_out1.m_flow", +} + +COMPARISON_KEYS = tuple(MODELICA_COMPARISON_COLUMNS.keys()) + + +def _branch_series_values( + series: dict[str, list[float]], + branch_name: str, + legacy_key: str, +) -> list[float]: + generic_key = f"branch.{branch_name}.{legacy_key.split('.')[-1]}" + if generic_key in series: + return series[generic_key] + return series[legacy_key] + + +@dataclass(frozen=True) +class TestModelArtifacts: + primary_csv_path: Path + temperature_csv_path: Path + temperature_svg_path: Path + run_report_path: Path + comparison_csv_path: Path | None = None + comparison_summary_path: Path | None = None + + +def format_testmodel_run_report( + *, + network_summary: str, + initialization: Any, + raw_initial_state: tuple[float, ...], + consistent_initial_state: tuple[float, ...], + solution: Any, + series: dict[str, list[float]], + solve_diagnostics: Any, + artifacts: TestModelArtifacts, + comparison_summary: dict[str, tuple[float, float]] | None, +) -> str: + lines = [ + network_summary, + "", + f"Initialization converged: {initialization.converged}", + f"Initialization iterations: {initialization.iterations}", + f"Initialization max state delta: {initialization.max_state_delta:.6e}", + f"Initialization max flow delta: {initialization.max_flow_delta:.6e}", + f"Initialization max enthalpy delta: {initialization.max_enthalpy_delta:.6e}", + ( + "Initialization downstream pressure spread: " + f"{initialization.downstream_pressure_spread:.6e}" + ), + "", + "Raw initial state vector:", + str(list(raw_initial_state)), + "", + "Constraint-consistent initial state vector:", + str(list(consistent_initial_state)), + "", + f"Solver success: {solution.success}", + f"Solver message: {solution.message}", + f"Final time: {solution.t[-1]:.2f} s", + f"Final tank pressure: {series['mytank.p'][-1]:.3f} Pa", + f"Final tank temperature: {series['mytank.T'][-1]:.3f} K", + f"Final cylinder pressure: {series['mycylinder.p'][-1]:.3f} Pa", + ( + "Final branch inflow: " + f"{_branch_series_values(series, 'upper_branch', 'branch_upper.in')[-1] + _branch_series_values(series, 'lower_branch', 'branch_lower.in')[-1]:.6f} kg/s" + ), + ] + + if solve_diagnostics is not None: + lines.extend( + [ + "", + "Final closure solve diagnostics:", + ( + "Upper branch inlet solve: " + f"converged={solve_diagnostics.upper_branch_inlet.converged}, " + f"iterations={solve_diagnostics.upper_branch_inlet.iterations}, " + f"residual={solve_diagnostics.upper_branch_inlet.residual:.6e}" + ), + ( + "Lower branch inlet solve: " + f"converged={solve_diagnostics.lower_branch_inlet.converged}, " + f"iterations={solve_diagnostics.lower_branch_inlet.iterations}, " + f"residual={solve_diagnostics.lower_branch_inlet.residual:.6e}" + ), + ] + ) + if solve_diagnostics.downstream_pressure_projection is not None: + lines.append( + "Downstream pressure projection: " + f"converged={solve_diagnostics.downstream_pressure_projection.converged}, " + f"iterations={solve_diagnostics.downstream_pressure_projection.iterations}, " + f"residual={solve_diagnostics.downstream_pressure_projection.residual:.6e}" + ) + + lines.extend( + [ + f"Primary series CSV: {artifacts.primary_csv_path}", + f"Temperature CSV: {artifacts.temperature_csv_path}", + f"Temperature plot: {artifacts.temperature_svg_path}", + f"Run report TXT: {artifacts.run_report_path}", + ] + ) + + if ( + artifacts.comparison_csv_path is not None + and artifacts.comparison_summary_path is not None + ): + lines.extend( + [ + f"Modelica comparison CSV: {artifacts.comparison_csv_path}", + f"Modelica comparison summary: {artifacts.comparison_summary_path}", + ] + ) + if comparison_summary is not None: + for key, (max_abs_error, max_rel_error) in comparison_summary.items(): + lines.append( + f"{key} max abs error: {max_abs_error:.6f}, " + f"max rel error: {max_rel_error:.6%}" + ) + + return "\n".join(lines) + "\n" + + +def write_testmodel_run_report(output_dir: Path, report_text: str) -> Path: + report_path = output_dir / "testmodel_run_report.txt" + report_path.write_text(report_text, encoding="utf-8") + return report_path + + +def _write_primary_series_csv(output_dir: Path, series: dict[str, list[float]]) -> Path: + csv_path = output_dir / "testmodel_primary_series.csv" + with csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle) + writer.writerow(["time_s", *PRIMARY_KEYS]) + for index, time_value in enumerate(series["time"]): + writer.writerow([time_value, *(series[key][index] for key in PRIMARY_KEYS)]) + return csv_path + + +def _write_temperature_csv(output_dir: Path, time_values: list[float], temperatures: list[float]) -> Path: + csv_path = output_dir / "testmodel_tank_temperature.csv" + with csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle) + writer.writerow(["time_s", "mytank_T_K"]) + writer.writerows(zip(time_values, temperatures)) + return csv_path + + +def _write_temperature_svg(output_dir: Path, time_values: list[float], temperatures: list[float]) -> Path: + svg_path = output_dir / "testmodel_tank_temperature.svg" + + width = 900 + height = 520 + left = 90 + right = 40 + top = 60 + bottom = 70 + plot_width = width - left - right + plot_height = height - top - bottom + + min_time = min(time_values) + max_time = max(time_values) + min_temp = min(temperatures) + max_temp = max(temperatures) + temp_padding = max(1.0, (max_temp - min_temp) * 0.08) + min_temp -= temp_padding + max_temp += temp_padding + + def scale_x(value: float) -> float: + return left + (value - min_time) / max(max_time - min_time, 1e-12) * plot_width + + def scale_y(value: float) -> float: + return top + (max_temp - value) / max(max_temp - min_temp, 1e-12) * plot_height + + points = " ".join( + f"{scale_x(time_value):.2f},{scale_y(temperature):.2f}" + for time_value, temperature in zip(time_values, temperatures) + ) + + x_ticks = 5 + y_ticks = 5 + x_tick_markup = [] + y_tick_markup = [] + + for index in range(x_ticks + 1): + fraction = index / x_ticks + time_value = min_time + fraction * (max_time - min_time) + x = left + fraction * plot_width + x_tick_markup.append( + f'' + ) + x_tick_markup.append( + f'' + f"{time_value:.1f}" + ) + + for index in range(y_ticks + 1): + fraction = index / y_ticks + temp_value = min_temp + fraction * (max_temp - min_temp) + y = top + plot_height - fraction * plot_height + y_tick_markup.append( + f'' + ) + y_tick_markup.append( + f'' + f"{temp_value:.1f}" + ) + + svg_content = f""" + + Python Testmodel Tank Temperature + Time (s) + Temperature (K) + + {''.join(x_tick_markup)} + {''.join(y_tick_markup)} + + +""" + svg_path.write_text(svg_content, encoding="utf-8") + return svg_path + + +def load_modelica_series(csv_path: Path, variable_names: tuple[str, ...]) -> dict[str, list[float]]: + series = {"time": []} + for variable_name in variable_names: + series[variable_name] = [] + + with csv_path.open("r", newline="", encoding="utf-8") as handle: + reader = csv.DictReader(handle) + available_variable_names = tuple( + variable_name + for variable_name in variable_names + if MODELICA_COMPARISON_COLUMNS.get(variable_name, variable_name) in (reader.fieldnames or ()) + ) + for row in reader: + series["time"].append(float(row["time"])) + for variable_name in available_variable_names: + modelica_column = MODELICA_COMPARISON_COLUMNS.get(variable_name, variable_name) + series[variable_name].append(float(row[modelica_column])) + + return series + + +def _interpolate_series_value(time_values: list[float], values: list[float], target_time: float) -> float: + if target_time <= time_values[0]: + return values[0] + if target_time >= time_values[-1]: + return values[-1] + + right_index = bisect_left(time_values, target_time) + if right_index < len(time_values) and abs(time_values[right_index] - target_time) <= 1e-12: + return values[right_index] + + left_index = right_index - 1 + left_time = time_values[left_index] + right_time = time_values[right_index] + fraction = (target_time - left_time) / (right_time - left_time) + return values[left_index] + fraction * (values[right_index] - values[left_index]) + + +def write_modelica_comparison( + output_dir: Path, + python_series: dict[str, list[float]], + modelica_series: dict[str, list[float]], +) -> tuple[Path, Path, dict[str, tuple[float, float]]]: + comparison_csv_path = output_dir / "testmodel_modelica_comparison.csv" + summary_path = output_dir / "testmodel_modelica_comparison_summary.txt" + summary: dict[str, tuple[float, float]] = {} + + with comparison_csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle) + header = ["time_s"] + comparison_keys = tuple( + key + for key in COMPARISON_KEYS + if key in python_series and key in modelica_series and modelica_series[key] + ) + for key in comparison_keys: + header.extend( + [ + f"python.{key}", + f"modelica.{key}", + f"abs_error.{key}", + f"rel_error.{key}", + ] + ) + writer.writerow(header) + + max_abs_errors = {key: 0.0 for key in comparison_keys} + max_rel_errors = {key: 0.0 for key in comparison_keys} + + for index, time_value in enumerate(python_series["time"]): + row = [time_value] + for key in comparison_keys: + python_value = python_series[key][index] + modelica_value = _interpolate_series_value( + modelica_series["time"], + modelica_series[key], + time_value, + ) + abs_error = abs(python_value - modelica_value) + rel_error = abs_error / max(abs(modelica_value), 1e-9) + max_abs_errors[key] = max(max_abs_errors[key], abs_error) + max_rel_errors[key] = max(max_rel_errors[key], rel_error) + row.extend([python_value, modelica_value, abs_error, rel_error]) + writer.writerow(row) + + summary_lines = [] + for key in comparison_keys: + summary[key] = (max_abs_errors[key], max_rel_errors[key]) + summary_lines.append( + f"{key}: max_abs_error={max_abs_errors[key]:.6f}, " + f"max_rel_error={max_rel_errors[key]:.6%}" + ) + summary_path.write_text("\n".join(summary_lines) + "\n", encoding="utf-8") + return comparison_csv_path, summary_path, summary + + +def export_testmodel_artifacts( + *, + output_dir: Path, + series: dict[str, list[float]], + modelica_series: dict[str, list[float]] | None = None, +) -> tuple[TestModelArtifacts, dict[str, tuple[float, float]] | None]: + output_dir.mkdir(parents=True, exist_ok=True) + + primary_csv_path = _write_primary_series_csv(output_dir, series) + temperature_csv_path = _write_temperature_csv( + output_dir, + series["time"], + series["mytank.T"], + ) + temperature_svg_path = _write_temperature_svg( + output_dir, + series["time"], + series["mytank.T"], + ) + + comparison_csv_path = None + comparison_summary_path = None + comparison_summary = None + if modelica_series is not None: + ( + comparison_csv_path, + comparison_summary_path, + comparison_summary, + ) = write_modelica_comparison(output_dir, series, modelica_series) + + return ( + TestModelArtifacts( + primary_csv_path=primary_csv_path, + temperature_csv_path=temperature_csv_path, + temperature_svg_path=temperature_svg_path, + run_report_path=output_dir / "testmodel_run_report.txt", + comparison_csv_path=comparison_csv_path, + comparison_summary_path=comparison_summary_path, + ), + comparison_summary, + ) diff --git a/PythonModels/scripts/run_testmodel.py b/PythonModels/scripts/run_testmodel.py new file mode 100644 index 0000000..4834677 --- /dev/null +++ b/PythonModels/scripts/run_testmodel.py @@ -0,0 +1,219 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from datetime import UTC, datetime +from pathlib import Path + +from PythonModels.reporting import ( + COMPARISON_KEYS, + PRIMARY_KEYS, + TestModelArtifacts, + export_testmodel_artifacts, + format_testmodel_run_report, + load_modelica_series, + write_testmodel_run_report, +) +from PythonModels.core.solver import SolveIVPConfig +from PythonModels.systems.testmodel import ( + InitializationDiagnostics, + TestModelConfig, + TestModelSystem, +) +from PythonModels.systems.testmodel_closure import TestModelSolveDiagnostics + + +@dataclass(frozen=True) +class TestModelSamplingConfig: + step: float = 0.1 + + +@dataclass(frozen=True) +class TestModelPathConfig: + output_dir: Path | None = None + modelica_result_path: Path | None = None + + +@dataclass(frozen=True) +class TestModelExecutionConfig: + use_modelica_reference_if_available: bool = True + + +@dataclass(frozen=True) +class TestModelRunConfig: + model: TestModelConfig = field(default_factory=TestModelConfig) + solver: SolveIVPConfig = field(default_factory=SolveIVPConfig) + sampling: TestModelSamplingConfig = field(default_factory=TestModelSamplingConfig) + paths: TestModelPathConfig = field(default_factory=TestModelPathConfig) + execution: TestModelExecutionConfig = field(default_factory=TestModelExecutionConfig) + + @property + def sample_step(self) -> float: + return self.sampling.step + + def sample_times(self) -> list[float]: + return _sample_times( + self.solver.t_start, + self.solver.t_stop, + step=self.sampling.step, + ) + + +@dataclass(frozen=True) +class PreparedTestModelRun: + run_config: TestModelRunConfig + repo_root: Path + output_dir: Path + modelica_result_path: Path + t_eval: tuple[float, ...] + use_modelica_reference_if_available: bool + modelica_reference_exists: bool + + +@dataclass(frozen=True) +class TestModelRunResult: + run_config: TestModelRunConfig + prepared_run: PreparedTestModelRun + system: TestModelSystem + initialization: InitializationDiagnostics + raw_initial_state: tuple[float, ...] + consistent_initial_state: tuple[float, ...] + solution: object + series: dict[str, list[float]] + solve_diagnostics: TestModelSolveDiagnostics | None + artifacts: TestModelArtifacts + comparison_summary: dict[str, tuple[float, float]] | None + used_modelica_reference: bool + + +def _sample_times(t_start: float, t_stop: float, step: float) -> list[float]: + point_count = int(round((t_stop - t_start) / step)) + return [t_start + index * step for index in range(point_count + 1)] + + +def _default_run_output_dir(pythonmodels_root: Path) -> Path: + timestamp = datetime.now(UTC).strftime("testmodel_%Y%m%d_%H%M%S_%f") + return pythonmodels_root / "runs" / timestamp + + +def prepare_testmodel_run( + *, + run_config: TestModelRunConfig | None = None, + output_dir: Path | None = None, + modelica_result_path: Path | None = None, +) -> PreparedTestModelRun: + run_config = run_config or TestModelRunConfig() + repo_root = Path(__file__).resolve().parents[2] + pythonmodels_root = Path(__file__).resolve().parents[1] + resolved_output_dir = ( + output_dir + or run_config.paths.output_dir + or _default_run_output_dir(pythonmodels_root) + ) + resolved_modelica_result_path = ( + modelica_result_path + or run_config.paths.modelica_result_path + or repo_root / "ModelicaModels" / "Simulation" / "Testmodel_res.csv" + ) + t_eval = tuple(run_config.sample_times()) + return PreparedTestModelRun( + run_config=run_config, + repo_root=repo_root, + output_dir=resolved_output_dir, + modelica_result_path=resolved_modelica_result_path, + t_eval=t_eval, + use_modelica_reference_if_available=run_config.execution.use_modelica_reference_if_available, + modelica_reference_exists=resolved_modelica_result_path.exists(), + ) + + +def run_prepared_testmodel(prepared_run: PreparedTestModelRun) -> TestModelRunResult: + run_config = prepared_run.run_config + system = TestModelSystem(config=run_config.model) + raw_initial_state = tuple(system.initial_state_vector()) + initialization = system.initialize_consistent_state() + consistent_initial_state = tuple(initialization.state_vector) + solution = system.simulate(config=run_config.solver, t_eval=list(prepared_run.t_eval)) + series = system.evaluate_solution(solution) + solve_diagnostics = system.last_solve_diagnostics + + modelica_series = None + used_modelica_reference = False + if ( + prepared_run.use_modelica_reference_if_available + and prepared_run.modelica_reference_exists + ): + modelica_series = load_modelica_series( + prepared_run.modelica_result_path, + COMPARISON_KEYS, + ) + used_modelica_reference = True + + artifacts, comparison_summary = export_testmodel_artifacts( + output_dir=prepared_run.output_dir, + series=series, + modelica_series=modelica_series, + ) + report_text = format_testmodel_run_report( + network_summary=system.network.summary(), + initialization=initialization, + raw_initial_state=raw_initial_state, + consistent_initial_state=consistent_initial_state, + solution=solution, + series=series, + solve_diagnostics=solve_diagnostics, + artifacts=artifacts, + comparison_summary=comparison_summary, + ) + write_testmodel_run_report(prepared_run.output_dir, report_text) + + return TestModelRunResult( + run_config=run_config, + prepared_run=prepared_run, + system=system, + initialization=initialization, + raw_initial_state=raw_initial_state, + consistent_initial_state=consistent_initial_state, + solution=solution, + series=series, + solve_diagnostics=solve_diagnostics, + artifacts=artifacts, + comparison_summary=comparison_summary, + used_modelica_reference=used_modelica_reference, + ) + + +def run_testmodel( + *, + run_config: TestModelRunConfig | None = None, + output_dir: Path | None = None, + modelica_result_path: Path | None = None, +) -> TestModelRunResult: + prepared_run = prepare_testmodel_run( + run_config=run_config, + output_dir=output_dir, + modelica_result_path=modelica_result_path, + ) + return run_prepared_testmodel(prepared_run) + + +def main() -> None: + run_config = TestModelRunConfig() + result = run_testmodel(run_config=run_config) + print( + format_testmodel_run_report( + network_summary=result.system.network.summary(), + initialization=result.initialization, + raw_initial_state=result.raw_initial_state, + consistent_initial_state=result.consistent_initial_state, + solution=result.solution, + series=result.series, + solve_diagnostics=result.solve_diagnostics, + artifacts=result.artifacts, + comparison_summary=result.comparison_summary, + ), + end="", + ) + + +if __name__ == "__main__": + main() diff --git a/PythonModels/systems/__init__.py b/PythonModels/systems/__init__.py new file mode 100644 index 0000000..187800c --- /dev/null +++ b/PythonModels/systems/__init__.py @@ -0,0 +1,2 @@ +"""System assembly modules.""" + diff --git a/PythonModels/systems/testmodel.py b/PythonModels/systems/testmodel.py new file mode 100644 index 0000000..2652292 --- /dev/null +++ b/PythonModels/systems/testmodel.py @@ -0,0 +1,303 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Any + +from PythonModels.components.cylinder import Cylinder +from PythonModels.components.orifice import Orifice +from PythonModels.components.pipe import Pipe +from PythonModels.components.tank import Tank +from PythonModels.components.tee import Tee +from PythonModels.core.medium import IdealGasMedium +from PythonModels.core.network import SimulationNetwork +from PythonModels.core.solver import SolveIVPConfig, integrate_ode +from PythonModels.systems.testmodel_closure import ( + BranchClosureComponents, + InitializationDiagnostics, + TestModelClosure, + TestModelClosureComponents, + TestModelSnapshot, +) + + +@dataclass(frozen=True) +class CylinderConfig: + volume: float = 0.01 + p0: float = 35e6 + T0: float = 300.0 + + +@dataclass(frozen=True) +class OrificeConfig: + K: float = 1e-5 + + +@dataclass(frozen=True) +class TankConfig: + volume: float = 0.1 + p0: float = 1e5 + T0: float = 300.0 + + +@dataclass(frozen=True) +class PipeConfig: + length: float = 5.0 + diameter: float = 0.02 + lambda_darcy: float = 0.02 + p0: float = 1e5 + T0: float = 300.0 + + +@dataclass(frozen=True) +class BranchConfig: + orifice: OrificeConfig = field(default_factory=OrificeConfig) + pipe: PipeConfig = field(default_factory=PipeConfig) + + +@dataclass(frozen=True) +class TestModelConfig: + cylinder: CylinderConfig = field(default_factory=CylinderConfig) + upper_branch: BranchConfig = field(default_factory=BranchConfig) + lower_branch: BranchConfig = field(default_factory=BranchConfig) + tank: TankConfig = field(default_factory=TankConfig) + + +class TestModelSystem: + """Runnable first-pass Python system for the current Testmodel topology. + + This version keeps the component split from the Modelica model while keeping + the downstream tee-tank pressure coupling in the ODE framework. The original + Modelica system is a tighter DAE because both pipe outlets discharge into an + ideal lossless junction directly connected to the tank. Here the branch + outlet flows are solved from a pressure-consistent energy balance so the + outlet is no longer driven by an arbitrary conductance parameter. + """ + + def __init__( + self, + medium: IdealGasMedium | None = None, + config: TestModelConfig | None = None, + ) -> None: + self.medium = medium or IdealGasMedium() + self.config = config or TestModelConfig() + + self.mycylinder = Cylinder( + name="mycylinder", + medium=self.medium, + V=self.config.cylinder.volume, + p0=self.config.cylinder.p0, + T0=self.config.cylinder.T0, + ) + self.mytee = Tee(name="mytee") + self.myorifice = Orifice(name="myorifice", K=self.config.upper_branch.orifice.K) + self.mypipe = Pipe( + name="mypipe", + medium=self.medium, + L=self.config.upper_branch.pipe.length, + D=self.config.upper_branch.pipe.diameter, + lambda_darcy=self.config.upper_branch.pipe.lambda_darcy, + p0=self.config.upper_branch.pipe.p0, + T0=self.config.upper_branch.pipe.T0, + ) + self.myorifice1 = Orifice(name="myorifice1", K=self.config.lower_branch.orifice.K) + self.mypipe1 = Pipe( + name="mypipe1", + medium=self.medium, + L=self.config.lower_branch.pipe.length, + D=self.config.lower_branch.pipe.diameter, + lambda_darcy=self.config.lower_branch.pipe.lambda_darcy, + p0=self.config.lower_branch.pipe.p0, + T0=self.config.lower_branch.pipe.T0, + ) + self.mytee1 = Tee(name="mytee1") + self.mytank = Tank( + name="mytank", + medium=self.medium, + V=self.config.tank.volume, + p0=self.config.tank.p0, + T0=self.config.tank.T0, + ) + + self.network = SimulationNetwork(name="Testmodel") + for component in ( + self.mycylinder, + self.mytee, + self.myorifice, + self.mypipe, + self.myorifice1, + self.mypipe1, + self.mytee1, + self.mytank, + ): + self.network.add_component(component) + + self.network.connect("mycylinder", "port_b", "mytee", "port_in") + self.network.connect("mytee", "port_out1", "myorifice", "port_a") + self.network.connect("myorifice", "port_b", "mypipe", "port_a") + self.network.connect("mypipe", "port_b", "mytee1", "port_out2") + self.network.connect("mytee", "port_out2", "myorifice1", "port_a") + self.network.connect("myorifice1", "port_b", "mypipe1", "port_a") + self.network.connect("mypipe1", "port_b", "mytee1", "port_out1") + self.network.connect("mytee1", "port_in", "mytank", "port_a") + + self.closure = TestModelClosure( + medium=self.medium, + components=TestModelClosureComponents( + cylinder=self.mycylinder, + upstream_tee=self.mytee, + upper_branch=BranchClosureComponents( + name="upper_branch", + orifice=self.myorifice, + pipe=self.mypipe, + ), + lower_branch=BranchClosureComponents( + name="lower_branch", + orifice=self.myorifice1, + pipe=self.mypipe1, + ), + downstream_tee=self.mytee1, + tank=self.mytank, + ), + initial_state_vector=self.initial_state_vector, + apply_state_vector=self.apply_state_vector, + ) + + def initial_state_vector(self) -> list[float]: + return self.network.initial_state_vector() + + def apply_state_vector(self, values: list[float]) -> None: + self.network.apply_state_vector(values) + + def consistent_initial_state_vector(self) -> list[float]: + return self.closure.consistent_initial_state_vector() + + @property + def last_solve_diagnostics(self): + return self.closure.last_solve_diagnostics + + def initialize_consistent_state( + self, + max_iterations: int = 12, + state_tolerance: float = 1e-9, + flow_tolerance: float = 1e-9, + enthalpy_tolerance: float = 1e-6, + pressure_tolerance: float = 1e-6, + strict_internal_solvers: bool = False, + ) -> InitializationDiagnostics: + return self.closure.initialize_consistent_state( + max_iterations=max_iterations, + state_tolerance=state_tolerance, + flow_tolerance=flow_tolerance, + enthalpy_tolerance=enthalpy_tolerance, + pressure_tolerance=pressure_tolerance, + strict_internal_solvers=strict_internal_solvers, + ) + + def project_downstream_pressure_constraints(self, *, strict: bool = False) -> None: + self.closure.project_downstream_pressure_constraints(strict=strict) + + def snapshot( + self, + state_vector: list[float] | None = None, + *, + strict: bool = False, + ) -> TestModelSnapshot: + return self.closure.snapshot(state_vector, strict=strict) + + def rhs(self, _t: float, state_vector: list[float]) -> list[float]: + return self.closure.rhs(state_vector) + + @staticmethod + def _legacy_branch_series_key_map() -> tuple[tuple[str, str, str], tuple[str, str, str]]: + return ( + ("upper_branch", "branch_upper.in", "branch_upper.out"), + ("lower_branch", "branch_lower.in", "branch_lower.out"), + ) + + @classmethod + def _legacy_branch_series_keys_by_name(cls) -> dict[str, tuple[str, str]]: + return { + branch_name: (inlet_key, outlet_key) + for branch_name, inlet_key, outlet_key in cls._legacy_branch_series_key_map() + } + + @staticmethod + def _generic_branch_series_keys(branch_name: str) -> tuple[str, str, str]: + return ( + f"branch.{branch_name}.p", + f"branch.{branch_name}.in", + f"branch.{branch_name}.out", + ) + + @staticmethod + def _legacy_branch_pressure_keys_by_name() -> dict[str, str]: + return { + "upper_branch": "mypipe.p", + "lower_branch": "mypipe1.p", + } + + @classmethod + def _append_legacy_branch_series_aliases( + cls, + series: dict[str, list[float]], + ) -> dict[str, list[float]]: + legacy_branch_series_keys = cls._legacy_branch_series_keys_by_name() + legacy_branch_pressure_keys = cls._legacy_branch_pressure_keys_by_name() + for branch_name, (legacy_inlet_key, legacy_outlet_key) in legacy_branch_series_keys.items(): + pressure_key, generic_inlet_key, generic_outlet_key = cls._generic_branch_series_keys( + branch_name + ) + series[legacy_branch_pressure_keys[branch_name]] = list(series[pressure_key]) + series[legacy_inlet_key] = list(series[generic_inlet_key]) + series[legacy_outlet_key] = list(series[generic_outlet_key]) + return series + + def simulate( + self, + config: SolveIVPConfig | None = None, + t_eval: list[float] | None = None, + ) -> Any: + return integrate_ode( + rhs=self.rhs, + initial_state=self.consistent_initial_state_vector(), + config=config or SolveIVPConfig(), + t_eval=t_eval, + ) + + def evaluate_solution(self, solution: Any) -> dict[str, list[float]]: + series = { + "time": [], + "mycylinder.p": [], + "mycylinder.T": [], + "mytank.p": [], + "mytank.T": [], + } + for branch_name, _, _ in self._legacy_branch_series_key_map(): + pressure_key, inlet_key, outlet_key = self._generic_branch_series_keys(branch_name) + series[pressure_key] = [] + series[inlet_key] = [] + series[outlet_key] = [] + + for index, time_value in enumerate(solution.t): + state_vector = [row[index] for row in solution.y] + snapshot = self.snapshot(state_vector) + series["time"].append(float(time_value)) + series["mycylinder.p"].append(snapshot.cylinder.p) + series["mycylinder.T"].append(snapshot.cylinder.T) + series["mytank.p"].append(snapshot.tank.p) + series["mytank.T"].append(snapshot.tank.T) + for branch in snapshot.branches: + pressure_key, generic_inlet_key, generic_outlet_key = self._generic_branch_series_keys( + branch.name + ) + series[pressure_key].append(branch.pipe.p) + series[generic_inlet_key].append(branch.inlet_flow) + series[generic_outlet_key].append(branch.outlet_flow) + + return self._append_legacy_branch_series_aliases(series) + + +def build_testmodel() -> SimulationNetwork: + """Compatibility helper for callers that only need the topology.""" + + return TestModelSystem().network diff --git a/PythonModels/systems/testmodel_closure.py b/PythonModels/systems/testmodel_closure.py new file mode 100644 index 0000000..cbcaff8 --- /dev/null +++ b/PythonModels/systems/testmodel_closure.py @@ -0,0 +1,668 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Callable + +from PythonModels.components.cylinder import Cylinder +from PythonModels.components.orifice import Orifice +from PythonModels.components.pipe import Pipe +from PythonModels.components.tank import Tank +from PythonModels.components.tee import Tee +from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties +from PythonModels.core.state import VolumeState + + +@dataclass(frozen=True) +class BranchInletFlowDiagnostics: + converged: bool + iterations: int + residual: float + m_flow: float + inlet_pressure: float + + +@dataclass(frozen=True) +class DownstreamPressureDiagnostics: + converged: bool + iterations: int + residual: float + pressure: float + target_total_internal_energy: float + + +@dataclass(frozen=True) +class TestModelSolveDiagnostics: + upper_branch_inlet: BranchInletFlowDiagnostics + lower_branch_inlet: BranchInletFlowDiagnostics + downstream_pressure_projection: DownstreamPressureDiagnostics | None + + +@dataclass(frozen=True) +class BranchClosureComponents: + name: str + orifice: Orifice + pipe: Pipe + + +@dataclass(frozen=True) +class BranchClosureState: + name: str + pipe: ThermodynamicProperties + inlet_flow: float + outlet_flow: float + inlet_h: float + inlet_flow_diagnostics: BranchInletFlowDiagnostics + + +@dataclass(frozen=True) +class BranchSnapshot: + name: str + pipe: ThermodynamicProperties + inlet_flow: float + outlet_flow: float + inlet_h: float + inlet_flow_diagnostics: BranchInletFlowDiagnostics + + +@dataclass(frozen=True) +class TestModelSnapshot: + cylinder: ThermodynamicProperties + tank: ThermodynamicProperties + tee_upstream_h: float + tee_downstream_h: float + branches: tuple[BranchSnapshot, ...] = field(default_factory=tuple) + solve_diagnostics: TestModelSolveDiagnostics | None = None + + @property + def pipe_upper(self) -> ThermodynamicProperties: + return self.branches[0].pipe + + @property + def pipe_lower(self) -> ThermodynamicProperties: + return self.branches[1].pipe + + @property + def branch_inlet_flows(self) -> tuple[float, ...]: + return tuple(branch.inlet_flow for branch in self.branches) + + @property + def branch_outlet_flows(self) -> tuple[float, ...]: + return tuple(branch.outlet_flow for branch in self.branches) + + +@dataclass(frozen=True) +class InitializationDiagnostics: + converged: bool + iterations: int + max_state_delta: float + max_flow_delta: float + max_enthalpy_delta: float + downstream_pressure_spread: float + state_vector: tuple[float, ...] + + +@dataclass(frozen=True) +class TestModelClosureComponents: + cylinder: Cylinder + upstream_tee: Tee + upper_branch: BranchClosureComponents + lower_branch: BranchClosureComponents + downstream_tee: Tee + tank: Tank + + def branches(self) -> tuple[BranchClosureComponents, BranchClosureComponents]: + return (self.upper_branch, self.lower_branch) + + +class TestModelClosure: + """Owns Testmodel-specific closure, projection and port-writeback logic.""" + + def __init__( + self, + *, + medium: IdealGasMedium, + components: TestModelClosureComponents, + initial_state_vector: Callable[[], list[float]], + apply_state_vector: Callable[[list[float]], None], + ) -> None: + self.medium = medium + self.components = components + self._initial_state_vector = initial_state_vector + self._apply_state_vector = apply_state_vector + self.last_solve_diagnostics: TestModelSolveDiagnostics | None = None + self.last_downstream_pressure_diagnostics: DownstreamPressureDiagnostics | None = None + + @staticmethod + def _downstream_pressure_spread(snapshot: TestModelSnapshot) -> float: + downstream_pressures = tuple(branch.pipe.p for branch in snapshot.branches) + ( + snapshot.tank.p, + ) + return max(downstream_pressures) - min(downstream_pressures) + + @staticmethod + def _initialization_flow_delta( + previous_snapshot: TestModelSnapshot | None, + current_snapshot: TestModelSnapshot, + ) -> float: + if previous_snapshot is None: + return max(abs(branch.outlet_flow) for branch in current_snapshot.branches) + return max( + abs(curr - prev) + for curr, prev in zip( + current_snapshot.branch_outlet_flows, + previous_snapshot.branch_outlet_flows, + ) + ) + + @staticmethod + def _initialization_enthalpy_delta( + previous_snapshot: TestModelSnapshot | None, + current_snapshot: TestModelSnapshot, + ) -> float: + if previous_snapshot is None: + return abs(current_snapshot.tee_downstream_h - current_snapshot.tank.h) + return max( + abs(current_snapshot.tee_upstream_h - previous_snapshot.tee_upstream_h), + abs(current_snapshot.tee_downstream_h - previous_snapshot.tee_downstream_h), + ) + + def consistent_initial_state_vector(self) -> list[float]: + return list(self.initialize_consistent_state().state_vector) + + def initialize_consistent_state( + self, + max_iterations: int = 12, + state_tolerance: float = 1e-9, + flow_tolerance: float = 1e-9, + enthalpy_tolerance: float = 1e-6, + pressure_tolerance: float = 1e-6, + strict_internal_solvers: bool = False, + ) -> InitializationDiagnostics: + raw_state = self._initial_state_vector() + previous_snapshot: TestModelSnapshot | None = None + diagnostics: InitializationDiagnostics | None = None + + for iteration in range(1, max_iterations + 1): + state_before_projection = self._initial_state_vector() + self.snapshot(state_before_projection, strict=strict_internal_solvers) + + self.project_downstream_pressure_constraints(strict=strict_internal_solvers) + + state_after_projection = self._initial_state_vector() + snapshot_after_projection = self.snapshot( + state_after_projection, + strict=strict_internal_solvers, + ) + + max_state_delta = max( + abs(after - before) + for before, after in zip(state_before_projection, state_after_projection) + ) + max_flow_delta = self._initialization_flow_delta( + previous_snapshot, + snapshot_after_projection, + ) + max_enthalpy_delta = self._initialization_enthalpy_delta( + previous_snapshot, + snapshot_after_projection, + ) + downstream_pressure_spread = self._downstream_pressure_spread( + snapshot_after_projection, + ) + + diagnostics = InitializationDiagnostics( + converged=( + max_state_delta <= state_tolerance + and max_flow_delta <= flow_tolerance + and max_enthalpy_delta <= enthalpy_tolerance + and downstream_pressure_spread <= pressure_tolerance + ), + iterations=iteration, + max_state_delta=max_state_delta, + max_flow_delta=max_flow_delta, + max_enthalpy_delta=max_enthalpy_delta, + downstream_pressure_spread=downstream_pressure_spread, + state_vector=tuple(state_after_projection), + ) + previous_snapshot = snapshot_after_projection + + if diagnostics.converged: + self._apply_state_vector(raw_state) + return diagnostics + + assert diagnostics is not None + self._apply_state_vector(raw_state) + return diagnostics + + def _solve_branch_inlet_flow( + self, + orifice: Orifice, + pipe: Pipe, + p_upstream: float, + pipe_props: ThermodynamicProperties, + *, + strict: bool = False, + ) -> tuple[float, BranchInletFlowDiagnostics]: + m_flow = orifice.mass_flow(p_upstream, pipe_props.p) + rho = max(pipe_props.rho, 1e-9) + p_inlet = pipe.inlet_pressure(m_flow, rho, pipe_props.p) + residual = abs(orifice.mass_flow(p_upstream, p_inlet) - m_flow) + converged = False + iterations = 0 + for iteration in range(1, 9): + p_inlet = pipe.inlet_pressure(m_flow, rho, pipe_props.p) + next_m_flow = orifice.mass_flow(p_upstream, p_inlet) + residual = abs(next_m_flow - m_flow) + iterations = iteration + if residual <= 1e-9 * max(1.0, abs(next_m_flow)): + m_flow = next_m_flow + converged = True + break + m_flow = next_m_flow + diagnostics = BranchInletFlowDiagnostics( + converged=converged, + iterations=iterations, + residual=residual, + m_flow=m_flow, + inlet_pressure=p_inlet, + ) + if strict and not diagnostics.converged: + raise RuntimeError( + f"Branch inlet flow solve did not converge for {pipe.name}: residual={residual:.6e}" + ) + return m_flow, diagnostics + + def _solve_downstream_branch_flows( + self, + cylinder: ThermodynamicProperties, + tank: ThermodynamicProperties, + branch_states: tuple[BranchClosureState, BranchClosureState], + ) -> tuple[float, float]: + return self._solve_downstream_branch_flows_from_state( + inlet_h_upper=branch_states[0].inlet_h, + inlet_h_lower=branch_states[1].inlet_h, + pipe_upper_h=max(branch_states[0].pipe.h, 1e-9), + pipe_lower_h=max(branch_states[1].pipe.h, 1e-9), + tank_h=max(tank.h, 1e-9), + q_in_upper=branch_states[0].inlet_flow, + q_in_lower=branch_states[1].inlet_flow, + ) + + def _project_volume_energy_to_pressure( + self, + component: Pipe | Tank, + target_pressure: float, + ) -> None: + target_temperature = target_pressure * component.V / ( + max(component.state.m, 1e-12) * self.medium.R_gas + ) + target_internal_energy = ( + component.state.m * self.medium.specific_internal_energy(target_temperature) + ) + component.state = VolumeState(m=component.state.m, U=target_internal_energy) + + def _downstream_total_internal_energy_for_pressure( + self, + target_pressure: float, + downstream_components: tuple[Pipe | Tank, ...], + ) -> float: + total_internal_energy = 0.0 + for component in downstream_components: + target_temperature = target_pressure * component.V / ( + max(component.state.m, 1e-12) * self.medium.R_gas + ) + total_internal_energy += ( + component.state.m * self.medium.specific_internal_energy(target_temperature) + ) + return total_internal_energy + + def _solve_downstream_common_pressure( + self, + downstream_components: tuple[Pipe | Tank, ...], + target_total_internal_energy: float, + *, + strict: bool = False, + ) -> tuple[float, DownstreamPressureDiagnostics]: + lower_pressure = 1.0 + upper_pressure = max(component.properties().p for component in downstream_components) + upper_pressure = max(upper_pressure, 1e5) + + def residual(pressure: float) -> float: + return ( + self._downstream_total_internal_energy_for_pressure( + pressure, + downstream_components, + ) + - target_total_internal_energy + ) + + upper_residual = residual(upper_pressure) + iteration_count = 0 + + while upper_residual < 0.0: + upper_pressure *= 2.0 + upper_residual = residual(upper_pressure) + + final_pressure = 0.5 * (lower_pressure + upper_pressure) + final_residual = residual(final_pressure) + converged = False + for iteration in range(1, 81): + middle_pressure = 0.5 * (lower_pressure + upper_pressure) + middle_residual = residual(middle_pressure) + iteration_count = iteration + final_pressure = middle_pressure + final_residual = middle_residual + if abs(middle_residual) <= 1e-12 * max(1.0, target_total_internal_energy): + converged = True + break + if middle_residual > 0.0: + upper_pressure = middle_pressure + else: + lower_pressure = middle_pressure + + diagnostics = DownstreamPressureDiagnostics( + converged=converged, + iterations=iteration_count, + residual=final_residual, + pressure=final_pressure, + target_total_internal_energy=target_total_internal_energy, + ) + if strict and not diagnostics.converged: + raise RuntimeError( + "Downstream common-pressure solve did not converge: " + f"residual={final_residual:.6e}" + ) + return final_pressure, diagnostics + + def project_downstream_pressure_constraints(self, *, strict: bool = False) -> None: + downstream_components = ( + self.components.upper_branch.pipe, + self.components.lower_branch.pipe, + self.components.tank, + ) + total_internal_energy = sum(component.state.U for component in downstream_components) + common_pressure, diagnostics = self._solve_downstream_common_pressure( + downstream_components, + total_internal_energy, + strict=strict, + ) + self.last_downstream_pressure_diagnostics = diagnostics + + for component in downstream_components: + self._project_volume_energy_to_pressure(component, common_pressure) + + def _downstream_connection_enthalpy( + self, + q_out_upper: float, + q_out_lower: float, + pipe_upper_h: float, + pipe_lower_h: float, + tank_h: float, + ) -> float: + return self.components.downstream_tee.inlet_stream_enthalpy( + q_out_lower, + pipe_lower_h, + q_out_upper, + pipe_upper_h, + fallback_h=tank_h, + ) + + def _solve_downstream_branch_flows_from_state( + self, + *, + inlet_h_upper: float, + inlet_h_lower: float, + pipe_upper_h: float, + pipe_lower_h: float, + tank_h: float, + q_in_upper: float, + q_in_lower: float, + ) -> tuple[float, float]: + return self.components.downstream_tee.solve_branch_outlet_flows_from_energy_balance( + ratio_branch1=self.components.upper_branch.pipe.V / self.components.tank.V, + ratio_branch2=self.components.lower_branch.pipe.V / self.components.tank.V, + inlet_h_branch1=inlet_h_upper, + inlet_h_branch2=inlet_h_lower, + branch1_h=pipe_upper_h, + branch2_h=pipe_lower_h, + inlet_h=tank_h, + q_in_branch1=q_in_upper, + q_in_branch2=q_in_lower, + ) + + def _evaluate_branch_states( + self, + cylinder: ThermodynamicProperties, + ) -> tuple[BranchClosureState, BranchClosureState]: + states: list[BranchClosureState] = [] + for branch in self.components.branches(): + pipe_properties = branch.pipe.properties() + inlet_flow, inlet_flow_diagnostics = self._solve_branch_inlet_flow( + branch.orifice, + branch.pipe, + cylinder.p, + pipe_properties, + ) + inlet_h = branch.pipe.port_a_inlet_enthalpy( + port_a_m_flow=inlet_flow, + connected_h=cylinder.h, + internal_h=pipe_properties.h, + ) + states.append( + BranchClosureState( + name=branch.name, + pipe=pipe_properties, + inlet_flow=inlet_flow, + outlet_flow=0.0, + inlet_h=inlet_h, + inlet_flow_diagnostics=inlet_flow_diagnostics, + ) + ) + return (states[0], states[1]) + + @staticmethod + def _with_branch_outlet_flows( + branch_states: tuple[BranchClosureState, BranchClosureState], + outlet_flows: tuple[float, float], + ) -> tuple[BranchClosureState, BranchClosureState]: + return ( + BranchClosureState( + name=branch_states[0].name, + pipe=branch_states[0].pipe, + inlet_flow=branch_states[0].inlet_flow, + outlet_flow=outlet_flows[0], + inlet_h=branch_states[0].inlet_h, + inlet_flow_diagnostics=branch_states[0].inlet_flow_diagnostics, + ), + BranchClosureState( + name=branch_states[1].name, + pipe=branch_states[1].pipe, + inlet_flow=branch_states[1].inlet_flow, + outlet_flow=outlet_flows[1], + inlet_h=branch_states[1].inlet_h, + inlet_flow_diagnostics=branch_states[1].inlet_flow_diagnostics, + ), + ) + + @staticmethod + def _branch_snapshots( + branch_states: tuple[BranchClosureState, BranchClosureState], + ) -> tuple[BranchSnapshot, BranchSnapshot]: + return ( + BranchSnapshot( + name=branch_states[0].name, + pipe=branch_states[0].pipe, + inlet_flow=branch_states[0].inlet_flow, + outlet_flow=branch_states[0].outlet_flow, + inlet_h=branch_states[0].inlet_h, + inlet_flow_diagnostics=branch_states[0].inlet_flow_diagnostics, + ), + BranchSnapshot( + name=branch_states[1].name, + pipe=branch_states[1].pipe, + inlet_flow=branch_states[1].inlet_flow, + outlet_flow=branch_states[1].outlet_flow, + inlet_h=branch_states[1].inlet_h, + inlet_flow_diagnostics=branch_states[1].inlet_flow_diagnostics, + ), + ) + + def snapshot( + self, + state_vector: list[float] | None = None, + *, + strict: bool = False, + ) -> TestModelSnapshot: + if state_vector is not None: + self._apply_state_vector(list(state_vector)) + + cylinder = self.components.cylinder.properties() + tank = self.components.tank.properties() + branch_states = self._evaluate_branch_states(cylinder) + if strict: + for branch_state in branch_states: + if not branch_state.inlet_flow_diagnostics.converged: + raise RuntimeError( + "Branch inlet flow solve did not converge for " + f"{branch_state.name}: residual=" + f"{branch_state.inlet_flow_diagnostics.residual:.6e}" + ) + outlet_flows = self._solve_downstream_branch_flows(cylinder, tank, branch_states) + branch_states = self._with_branch_outlet_flows(branch_states, outlet_flows) + + tee_upstream_h = self.components.upstream_tee.inlet_stream_enthalpy( + -branch_states[0].inlet_flow, + branch_states[0].pipe.h, + -branch_states[1].inlet_flow, + branch_states[1].pipe.h, + fallback_h=cylinder.h, + ) + tee_downstream_h = self._downstream_connection_enthalpy( + branch_states[0].outlet_flow, + branch_states[1].outlet_flow, + branch_states[0].pipe.h, + branch_states[1].pipe.h, + tank.h, + ) + + self._write_port_states( + cylinder, + tank, + branch_states, + tee_upstream_h, + tee_downstream_h, + ) + + solve_diagnostics = TestModelSolveDiagnostics( + upper_branch_inlet=branch_states[0].inlet_flow_diagnostics, + lower_branch_inlet=branch_states[1].inlet_flow_diagnostics, + downstream_pressure_projection=self.last_downstream_pressure_diagnostics, + ) + self.last_solve_diagnostics = solve_diagnostics + branch_snapshots = self._branch_snapshots(branch_states) + + return TestModelSnapshot( + cylinder=cylinder, + tank=tank, + tee_upstream_h=tee_upstream_h, + tee_downstream_h=tee_downstream_h, + branches=branch_snapshots, + solve_diagnostics=solve_diagnostics, + ) + + def _write_port_states( + self, + cylinder: ThermodynamicProperties, + tank: ThermodynamicProperties, + branch_states: tuple[BranchClosureState, BranchClosureState], + tee_upstream_h: float, + tee_downstream_h: float, + ) -> None: + cylinder_m_flow = -sum(branch_state.inlet_flow for branch_state in branch_states) + tank_m_flow = sum(branch_state.outlet_flow for branch_state in branch_states) + + self.components.cylinder.port_b.m_flow = cylinder_m_flow + + self.components.upstream_tee.port_in.p = cylinder.p + self.components.upstream_tee.port_out1.p = cylinder.p + self.components.upstream_tee.port_out2.p = cylinder.p + self.components.upstream_tee.port_in.m_flow = cylinder_m_flow + self.components.upstream_tee.port_in.h_outflow = tee_upstream_h + self.components.upstream_tee.port_out1.h_outflow = cylinder.h + self.components.upstream_tee.port_out2.h_outflow = cylinder.h + self.components.upstream_tee.port_out1.m_flow = -branch_states[0].inlet_flow + self.components.upstream_tee.port_out2.m_flow = -branch_states[1].inlet_flow + + for branch_components, branch_state in zip(self.components.branches(), branch_states): + branch_components.orifice.port_a.p = cylinder.p + branch_components.orifice.port_b.p = branch_components.pipe.inlet_pressure( + branch_state.inlet_flow, + max(branch_state.pipe.rho, 1e-9), + branch_state.pipe.p, + ) + branch_components.orifice.port_a.m_flow = branch_state.inlet_flow + branch_components.orifice.port_b.m_flow = -branch_state.inlet_flow + branch_components.orifice.port_a.h_outflow = cylinder.h + branch_components.orifice.port_b.h_outflow = branch_state.pipe.h + + branch_components.pipe.port_a.p = branch_components.orifice.port_b.p + branch_components.pipe.port_a.m_flow = branch_state.inlet_flow + branch_components.pipe.port_b.m_flow = -branch_state.outlet_flow + branch_components.pipe.port_b.p = branch_state.pipe.p + + self.components.downstream_tee.port_in.p = tank.p + self.components.downstream_tee.port_out1.p = tank.p + self.components.downstream_tee.port_out2.p = tank.p + self.components.downstream_tee.port_in.m_flow = -tank_m_flow + self.components.downstream_tee.port_out1.m_flow = branch_states[1].outlet_flow + self.components.downstream_tee.port_out2.m_flow = branch_states[0].outlet_flow + self.components.downstream_tee.port_in.h_outflow = tee_downstream_h + self.components.downstream_tee.port_out1.h_outflow = tank.h + self.components.downstream_tee.port_out2.h_outflow = tank.h + + self.components.tank.port_a.m_flow = tank_m_flow + + def _branch_derivative_states( + self, + snapshot: TestModelSnapshot, + ) -> tuple[VolumeState, VolumeState]: + derivative_states: list[VolumeState] = [] + for branch_components, branch_snapshot in zip(self.components.branches(), snapshot.branches): + derivative_states.append( + branch_components.pipe.derivatives_from_connections( + port_a_m_flow=branch_snapshot.inlet_flow, + connected_h_a=snapshot.cylinder.h, + port_b_m_flow=-branch_snapshot.outlet_flow, + connected_h_b=snapshot.tank.h, + internal_h=branch_snapshot.pipe.h, + ) + ) + return (derivative_states[0], derivative_states[1]) + + def rhs(self, state_vector: list[float]) -> list[float]: + snapshot = self.snapshot(state_vector) + + cylinder_m_flow = -sum(branch.inlet_flow for branch in snapshot.branches) + tank_m_flow = sum(branch.outlet_flow for branch in snapshot.branches) + d_cylinder = self.components.cylinder.derivatives_from_connection( + connected_h=snapshot.tee_upstream_h, + port_m_flow=cylinder_m_flow, + internal_h=snapshot.cylinder.h, + ) + branch_derivatives = self._branch_derivative_states(snapshot) + d_tank = self.components.tank.derivatives_from_connection( + connected_h=snapshot.tee_downstream_h, + port_m_flow=tank_m_flow, + internal_h=snapshot.tank.h, + ) + + return [ + d_cylinder.m, + d_cylinder.U, + branch_derivatives[0].m, + branch_derivatives[0].U, + branch_derivatives[1].m, + branch_derivatives[1].U, + d_tank.m, + d_tank.U, + ] diff --git a/README.md b/README.md index bcf5f6e..c1945dd 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,448 @@ # SystemSimulationApp +SystemSimulationApp 是一个面向系统仿真的可视化建模原型。当前项目已经具备两条基础能力： + +- 使用 FastAPI 嵌入 draw.io，提供可视化拖拽建模界面。 +- 使用 `PythonModels` 中已有的 Python 仿真模型，运行当前固定拓扑的 `TestModelSystem`。 + +下一阶段的目标是把这两条能力接起来，先完成一个小 MVP： + +```text +draw.io 画图 + -> 解析 draw.io XML + -> 生成 SystemModel + -> 映射到 PythonModels.TestModelConfig + -> 调用 run_testmodel() + -> 返回并展示仿真结果 +``` + +## 当前状态 + +### 已有前端和接口基础 + +- `app/main.py`: FastAPI 主入口。 +- `app/static/index.html`: 当前 draw.io iframe 页面。 +- `/api/component-library.drawiolib`: 当前 draw.io 自定义组件库接口。 +- `/api/system-xml`: 当前 draw.io XML 到简单系统 XML 的导出接口。 +- `/drawio/`: FastAPI 到 draw.io Docker 的同源代理入口。 + +### 已有 Python 仿真模型 + +`PythonModels` 目录里已经有一套 Python 仿真骨架： + +```text +PythonModels/components/ + cylinder.py -> Cylinder，端口 port_b + tank.py -> Tank，端口 port_a + pipe.py -> Pipe，端口 port_a / port_b + orifice.py -> Orifice，端口 port_a / port_b + tee.py -> Tee，端口 port_in / port_out1 / port_out2 + +PythonModels/core/ + network.py -> SimulationNetwork / Connection + ports.py -> PortState + solver.py -> integrate_ode / SolveIVPConfig + +PythonModels/systems/ + testmodel.py -> 当前唯一可运行系统 TestModelSystem +``` + +当前 `TestModelSystem` 支持的是固定拓扑： + +```text +mycylinder + -> mytee + -> myorifice -> mypipe + -> myorifice1 -> mypipe1 + -> mytee1 + -> mytank +``` + +因此第一版不要做任意拓扑仿真，而是先把 draw.io 中画出的这个固定拓扑映射到现有 `TestModelSystem`。 + +## 核心原则 + +不要让仿真逻辑直接依赖 draw.io XML。 + +推荐分层： + +```text +draw.io XML + -> SystemModel + -> PythonModels adapter + -> PythonModels runner + -> SimulationResult +``` + +其中： + +- draw.io XML 只是可视化画布格式。 +- `SystemModel` 是系统仿真的核心结构化数据。 +- PythonModels adapter 负责把 `SystemModel` 转成 `TestModelConfig / TestModelRunConfig`。 +- 仿真器只认结构化模型，不直接解析画布。 + +这样后续即使替换 draw.io，仿真层也不用推倒重来。 + +## MVP 范围 + +第一版只支持当前 `TestModelSystem` 对应的拓扑： + +```text +cylinder -> tee -> valve/orifice -> pipe -> tee -> tank + -> valve/orifice -> pipe -> +``` + +第一版支持的组件： + +| draw.io 类型 | PythonModels 类型 | 配置对象 | 端口 | +| --- | --- | --- | --- | +| `cylinder` | `Cylinder` | `CylinderConfig` | `port_b` | +| `tank` | `Tank` | `TankConfig` | `port_a` | +| `pipe` | `Pipe` | `PipeConfig` | `port_a`, `port_b` | +| `valve` / `orifice` | `Orifice` | `OrificeConfig` | `port_a`, `port_b` | +| `tee` | `Tee` | 暂无独立配置 | `port_in`, `port_out1`, `port_out2` | + +前端可以继续显示“阀门”，但底层建议使用： + +```xml +modelType="orifice" +``` + +因为 PythonModels 里实际对应的是 `Orifice`。 + +## SystemModel 设计 + +建议新增统一结构： + +```json +{ + "name": "testmodel_from_drawio", + "components": [ + { + "id": "mycylinder", + "type": "cylinder", + "role": "cylinder", + "parameters": { + "volume": 0.01, + "p0": 35000000, + "T0": 300 + } + } + ], + "connections": [ + { + "source": { + "component": "mycylinder", + "port": "port_b" + }, + "target": { + "component": "mytee", + "port": "port_in" + } + } + ], + "simulation": { + "t_start": 0, + "t_stop": 20, + "step": 0.1, + "method": "BDF" + } +} +``` + +## draw.io XML 扩展约定 + +组件节点需要补充仿真语义： + +```xml + + ... + +``` + +管段示例： + +```xml + + ... + +``` + +孔板 / 阀门示例： + +```xml + + ... + +``` + +连接需要从“组件级”升级到“端口级”： + +```xml + + + +``` + +如果第一版来不及做端口选择，可以后端按组件类型和拓扑位置推断默认端口，但这只能作为 MVP 兜底。 + +## 参数映射 + +从 draw.io / SystemModel 解析出的参数映射到 PythonModels 配置对象： + +```text +cylinder.volume -> CylinderConfig.volume +cylinder.p0 -> CylinderConfig.p0 +cylinder.T0 -> CylinderConfig.T0 + +tank.volume -> TankConfig.volume +tank.p0 -> TankConfig.p0 +tank.T0 -> TankConfig.T0 + +pipe.length -> PipeConfig.length +pipe.diameter -> PipeConfig.diameter +pipe.lambda_darcy -> PipeConfig.lambda_darcy +pipe.p0 -> PipeConfig.p0 +pipe.T0 -> PipeConfig.T0 + +orifice.K -> OrificeConfig.K +orifice.opening -> OrificeConfig.opening +``` + +注意：当前 `Orifice` 类支持 `opening`，但 `OrificeConfig` 暂时只有 `K`。建议后续给 `OrificeConfig` 增加 `opening` 字段，这样前端“阀门开度”才有实际落点。 + +## 建议新增后端模块 + +```text +app/models/system.py +app/services/drawio_parser.py +app/services/model_validator.py +app/services/pythonmodels_adapter.py +app/services/simulation_runner.py +``` + +职责： + +```text +app/models/system.py + 定义 SystemModel、ComponentInstance、Connection、SimulationConfig。 + +app/services/drawio_parser.py + 负责 draw.io XML -> SystemModel。 + +app/services/model_validator.py + 检查组件、端口、参数、拓扑是否合法。 + +app/services/pythonmodels_adapter.py + 负责 SystemModel -> TestModelConfig / TestModelRunConfig。 + +app/services/simulation_runner.py + 调用 PythonModels，并把结果整理成前端可消费的 JSON。 +``` + +## 建议新增接口 + +```text +POST /api/models/from-drawio + 输入 draw.io XML + 输出 SystemModel JSON + +POST /api/models/validate + 输入 SystemModel + 输出校验结果 + +POST /api/simulations/testmodel + 输入 SystemModel + 仿真配置 + 输出仿真结果 + +POST /api/simulations/testmodel/from-drawio + 输入 draw.io XML + 一步完成解析、校验、仿真 +``` + +第一版可以优先实现： + +```text +POST /api/simulations/testmodel/from-drawio +``` + +这样页面点击“运行仿真”时可以直接把 draw.io XML 发给后端。 + +## 第一版拓扑校验规则 + +第一版只支持固定 `TestModelSystem` 拓扑，因此校验可以严格一点： + +- 必须有 1 个 `cylinder`。 +- 必须有 1 个 `tank`。 +- 必须有 2 个 `tee`。 +- 必须有 2 个 `orifice` 或 `valve`。 +- 必须有 2 个 `pipe`。 +- 连接必须形成两条支路。 +- 每条支路必须是 `orifice -> pipe`。 +- 上游必须从 `cylinder` 到第一个 `tee`。 +- 下游必须汇入第二个 `tee` 再到 `tank`。 + +第一版不支持任意拓扑。 + +## 实现步骤和进度 + +### A. draw.io 侧建模增强 + +- [ ] 扩展组件库，增加 `cylinder / tee / orifice / pipe / tank`。 +- [ ] 给每个组件写入 `modelType`、`role`、默认参数。 +- [ ] 给组件补充端口语义。 +- [ ] 给连接 edge 补充 `sourcePort / targetPort`。 +- [ ] 前端增加“运行仿真”按钮。 + +### B. SystemModel 数据层 + +- [ ] 新增 `app/models/system.py`。 +- [ ] 定义 `SystemModel`。 +- [ ] 定义 `ComponentInstance`。 +- [ ] 定义端口级 `Connection`。 +- [ ] 定义 `SimulationConfig`。 + +### C. draw.io 解析 + +- [ ] 新增 `app/services/drawio_parser.py`。 +- [ ] 支持 `mxGraphModel` 输入。 +- [ ] 支持 `mxfile` 输入。 +- [ ] 解析 `` 组件。 +- [ ] 解析 `` 连接。 +- [ ] 把 draw.io 参数转换成正确的 Python 类型。 +- [ ] 输出标准 `SystemModel`。 + +### D. 模型校验 + +- [ ] 新增 `app/services/model_validator.py`。 +- [ ] 校验组件类型是否合法。 +- [ ] 校验必填参数是否存在。 +- [ ] 校验数值参数范围。 +- [ ] 校验端口名是否合法。 +- [ ] 校验连接 source / target 是否存在。 +- [ ] 校验是否符合固定 TestModel 拓扑。 +- [ ] 返回结构化错误列表。 + +### E. PythonModels 适配层 + +- [ ] 新增 `app/services/pythonmodels_adapter.py`。 +- [ ] 实现 `SystemModel -> CylinderConfig`。 +- [ ] 实现 `SystemModel -> TankConfig`。 +- [ ] 实现 `SystemModel -> PipeConfig`。 +- [ ] 实现 `SystemModel -> OrificeConfig`。 +- [ ] 实现 `SystemModel -> TestModelConfig`。 +- [ ] 实现 `SimulationConfig -> SolveIVPConfig`。 +- [ ] 如需要，给 `OrificeConfig` 增加 `opening` 字段。 + +### F. 仿真运行 + +- [ ] 新增 `app/services/simulation_runner.py`。 +- [ ] 调用 `PythonModels.scripts.run_testmodel.run_testmodel()`。 +- [ ] 整理返回的 `series`。 +- [ ] 返回 `mytank.p / mytank.T / mycylinder.p / mycylinder.T`。 +- [ ] 返回 `branch..p / in / out`。 +- [ ] 保留初始化诊断和求解诊断。 + +### G. FastAPI 接口 + +- [ ] 新增 `POST /api/models/from-drawio`。 +- [ ] 新增 `POST /api/models/validate`。 +- [ ] 新增 `POST /api/simulations/testmodel`。 +- [ ] 新增 `POST /api/simulations/testmodel/from-drawio`。 + +### H. 前端结果展示 + +- [ ] 点击“运行仿真”后导出当前 draw.io XML。 +- [ ] 调用 `/api/simulations/testmodel/from-drawio`。 +- [ ] 展示校验错误。 +- [ ] 展示仿真成功 / 失败状态。 +- [ ] 展示结果表格。 +- [ ] 展示压力、温度、流量曲线。 +- [ ] 支持下载结果 JSON / CSV。 + +## MVP 验收标准 + +第一版完成后，应能做到： + +- [ ] 打开 FastAPI 页面。 +- [ ] 在 draw.io 画布中拖入气瓶、三通、阀门/孔板、管段、贮箱。 +- [ ] 按固定 TestModel 拓扑连接组件。 +- [ ] 配置主要参数。 +- [ ] 点击运行仿真。 +- [ ] 后端解析 draw.io XML。 +- [ ] 后端生成 `SystemModel`。 +- [ ] 后端映射到 `TestModelConfig`。 +- [ ] 后端调用 `run_testmodel()`。 +- [ ] 前端展示 `mytank.p`、`mytank.T`、`mycylinder.p`、`mycylinder.T`。 +- [ ] 前端展示两条支路的压力和流量。 + +## 后续阶段：通用拓扑仿真 + +固定 TestModel 跑通后，再考虑通用拓扑： + +```text +SystemModel + -> instantiate components + -> SimulationNetwork.add_component() + -> SimulationNetwork.connect() + -> 自动生成 rhs + -> integrate_ode() +``` + +这一步需要新增： + +- 通用 connector / stream 求解。 +- 通用代数闭合。 +- 任意 `Tee / Pipe / Orifice` 网络流量求解。 +- 自动状态向量组装。 +- 自动 `rhs` 生成。 + +当前不建议立刻做这一步，因为 `PythonModels/systems/testmodel_closure.py` 仍然是为 `TestModelSystem` 专门写的闭合器，不是任意拓扑求解器。 + +## 当前优先级 + +下一步优先做： + +```text +draw.io XML + -> SystemModel + -> TestModelConfig + -> run_testmodel() +``` + +暂时不要直接做： + +```text +draw.io XML + -> PythonModels 任意网络 +``` + +这样最贴合现有代码，成本最低，技术债也最可控。 diff --git a/app/main.py b/app/main.py index c9718fd..40cd7f9 100644 --- a/app/main.py +++ b/app/main.py @@ -2,6 +2,7 @@ from __future__ import annotations import base64 import json +import math import os import urllib.error import urllib.request @@ -23,55 +24,938 @@ DRAWIO_UPSTREAM_URL = ( or os.getenv("DRAWIO_BASE_URL") or "http://127.0.0.1:8081/" ) +DRAWIO_APP_SCRIPT_PATH = "js/app.min.js" +GROUP_PROTECTION_STYLE = ( + "editable=0;resizable=0;rotatable=0;movable=1;deletable=1;connectable=0;" +) +CHILD_PROTECTION_STYLE = ( + "part=1;pointerEvents=0;editable=0;movable=0;resizable=0;rotatable=0;" + "deletable=0;connectable=0;" +) +PORT_PROTECTION_STYLE = ( + "editable=0;movable=0;resizable=0;rotatable=0;deletable=0;connectable=1;" +) + + +def library_cell( + cell_id: str, + style: str, + x: float, + y: float, + width: float, + height: float, + value: str = "", +) -> dict[str, str]: + return { + "id": cell_id, + "value": value, + "style": style, + "x": f"{x:g}", + "y": f"{y:g}", + "width": f"{width:g}", + "height": f"{height:g}", + } + + +def component_library_entry( + *, + title: str, + width: int, + height: int, + component_id: str, + label: str, + component_type: str, + model_type: str, + ports: str, + parameters: dict[str, object], + cells: list[dict[str, str]], +) -> dict[str, object]: + attributes = { + "id": component_id, + "label": label, + "componentType": component_type, + "modelType": model_type, + "ports": ports, + } + attributes.update({key: str(value) for key, value in parameters.items()}) + + graph = ET.Element("mxGraphModel") + root = ET.SubElement(graph, "root") + ET.SubElement(root, "mxCell", {"id": "0"}) + ET.SubElement(root, "mxCell", {"id": "1", "parent": "0"}) + + component = ET.SubElement(root, "object", attributes) + group_cell = ET.SubElement( + component, + "mxCell", + { + "style": ( + "group;html=1;whiteSpace=wrap;container=1;recursiveResize=0;" + "collapsible=0;labelPosition=center;verticalLabelPosition=bottom;" + "align=center;verticalAlign=top;spacingTop=4;fontSize=12;" + f"{GROUP_PROTECTION_STYLE}" + ), + "vertex": "1", + "connectable": "0", + "parent": "1", + }, + ) + ET.SubElement( + group_cell, + "mxGeometry", + {"width": str(width), "height": str(height), "as": "geometry"}, + ) + + for cell in cells: + cell_id = f"{component_id}_{cell['id']}" + is_port = cell["id"].startswith("port_") + + if is_port: + cell_parent = ET.SubElement( + root, + "object", + { + "id": cell_id, + "label": cell["value"], + "portName": cell["id"], + "portType": "fluid", + }, + ) + style = f"{cell['style']}{PORT_PROTECTION_STYLE}" + connectable = "1" + else: + cell_parent = root + style = f"{cell['style']}{CHILD_PROTECTION_STYLE}" + connectable = "0" + + cell_attributes = { + "value": cell["value"], + "style": style, + "vertex": "1", + "connectable": connectable, + "parent": component_id, + } + if not is_port: + cell_attributes["id"] = cell_id + + child = ET.SubElement(cell_parent, "mxCell", cell_attributes) + ET.SubElement( + child, + "mxGeometry", + { + "x": cell["x"], + "y": cell["y"], + "width": cell["width"], + "height": cell["height"], + "as": "geometry", + }, + ) + + return { + "title": title, + "w": width, + "h": height, + "xml": ET.tostring(graph, encoding="unicode"), + } + + +STROKE_DARK = "strokeColor=#3d4752;strokeWidth=2;" +PORT_STYLE = "ellipse;whiteSpace=wrap;html=1;fillColor=#ffffff;strokeColor=#3d4752;strokeWidth=2;" +LABEL_STYLE = ( + "text;html=1;strokeColor=none;fillColor=none;align=center;" + "verticalAlign=middle;whiteSpace=wrap;fontSize=12;fontStyle=1;" +) + +COMPONENT_SCHEMAS = { + "cylinder": { + "label": "\u6c14\u74f6", + "parameters": { + "volume": { + "label": "\u5bb9\u79ef", + "type": "float", + "unit": "m3", + "default": 0.01, + "min": 0, + }, + "p0": { + "label": "\u521d\u59cb\u538b\u529b", + "type": "float", + "unit": "Pa", + "default": 35000000, + "min": 0, + }, + "T0": { + "label": "\u521d\u59cb\u6e29\u5ea6", + "type": "float", + "unit": "K", + "default": 300, + "min": 0, + }, + }, + }, + "tank": { + "label": "\u8d2e\u7bb1", + "parameters": { + "volume": { + "label": "\u5bb9\u79ef", + "type": "float", + "unit": "m3", + "default": 0.1, + "min": 0, + }, + "p0": { + "label": "\u521d\u59cb\u538b\u529b", + "type": "float", + "unit": "Pa", + "default": 100000, + "min": 0, + }, + "T0": { + "label": "\u521d\u59cb\u6e29\u5ea6", + "type": "float", + "unit": "K", + "default": 300, + "min": 0, + }, + }, + }, + "pipe": { + "label": "\u7ba1\u6bb5", + "parameters": { + "length": { + "label": "\u957f\u5ea6", + "type": "float", + "unit": "m", + "default": 5.0, + "min": 0, + }, + "diameter": { + "label": "\u76f4\u5f84", + "type": "float", + "unit": "m", + "default": 0.02, + "min": 0, + }, + "lambda_darcy": { + "label": "\u6469\u963b\u7cfb\u6570", + "type": "float", + "default": 0.02, + "min": 0, + }, + "p0": { + "label": "\u521d\u59cb\u538b\u529b", + "type": "float", + "unit": "Pa", + "default": 100000, + "min": 0, + }, + "T0": { + "label": "\u521d\u59cb\u6e29\u5ea6", + "type": "float", + "unit": "K", + "default": 300, + "min": 0, + }, + }, + }, + "orifice": { + "label": "\u5b54\u677f/\u9600\u95e8", + "parameters": { + "K": { + "label": "\u6d41\u91cf\u7cfb\u6570", + "type": "float", + "default": 0.00001, + "min": 0, + }, + "opening": { + "label": "\u5f00\u5ea6", + "type": "float", + "default": 1.0, + "min": 0, + "max": 1, + }, + }, + }, + "tee": { + "label": "\u4e09\u901a", + "parameters": {}, + }, +} + + COMPONENT_LIBRARY = [ - { - "title": "\u8d2e\u7bb1", - "w": 80, - "h": 80, - "xml": ( - '' - '' - "" - ), - }, - { - "title": "\u9600\u95e8", - "w": 72, - "h": 56, - "xml": ( - '' - '' - ), - }, - { - "title": "\u7ba1\u6bb5", - "w": 96, - "h": 36, - "xml": ( - '' - '' - "" - ), - }, + component_library_entry( + title="\u6c14\u74f6", + width=90, + height=120, + component_id="cylinder", + label="\u6c14\u74f6", + component_type="cylinder", + model_type="cylinder", + ports="port_b", + parameters={"volume": 0.01, "p0": 35000000, "T0": 300}, + cells=[ + library_cell( + "cap", + "rounded=1;whiteSpace=wrap;html=1;arcSize=35;fillColor=#e5e7eb;" + f"{STROKE_DARK}", + 31, + 2, + 28, + 10, + ), + library_cell( + "neck", + f"rounded=0;whiteSpace=wrap;html=1;fillColor=#9ca3af;{STROKE_DARK}", + 37, + 11, + 16, + 10, + ), + library_cell( + "body", + "rounded=1;whiteSpace=wrap;html=1;arcSize=28;fillColor=#dae8fc;" + "strokeColor=#2f5597;strokeWidth=2;", + 24, + 18, + 42, + 80, + ), + library_cell( + "shoulder", + "ellipse;whiteSpace=wrap;html=1;fillColor=none;" + "strokeColor=#2f5597;strokeWidth=2;", + 29, + 18, + 32, + 16, + ), + library_cell( + "foot", + f"rounded=0;whiteSpace=wrap;html=1;fillColor=#9ca3af;{STROKE_DARK}", + 28, + 98, + 34, + 6, + ), + library_cell("port_b", PORT_STYLE, 67, 50, 10, 10), + library_cell("label", LABEL_STYLE, 14, 104, 62, 16, "\u6c14\u74f6"), + ], + ), + component_library_entry( + title="\u8d2e\u7bb1", + width=100, + height=112, + component_id="tank", + label="\u8d2e\u7bb1", + component_type="tank", + model_type="tank", + ports="port_a", + parameters={"volume": 0.1, "p0": 100000, "T0": 300}, + cells=[ + library_cell( + "vessel", + "shape=cylinder;whiteSpace=wrap;html=1;boundedLbl=1;" + "backgroundOutline=1;size=14;fillColor=#d5e8d4;" + "strokeColor=#4f7f45;strokeWidth=2;", + 18, + 8, + 64, + 72, + ), + library_cell( + "level", + "rounded=0;whiteSpace=wrap;html=1;fillColor=#ffffff;" + "strokeColor=#4f7f45;strokeWidth=1;dashed=1;", + 27, + 43, + 46, + 1, + ), + library_cell( + "leg_left", + f"rounded=0;whiteSpace=wrap;html=1;fillColor=#9ca3af;{STROKE_DARK}", + 30, + 80, + 8, + 16, + ), + library_cell( + "leg_right", + f"rounded=0;whiteSpace=wrap;html=1;fillColor=#9ca3af;{STROKE_DARK}", + 62, + 80, + 8, + 16, + ), + library_cell("port_a", PORT_STYLE, 84, 39, 10, 10), + library_cell("label", LABEL_STYLE, 18, 94, 64, 16, "\u8d2e\u7bb1"), + ], + ), + component_library_entry( + title="\u7ba1\u6bb5", + width=122, + height=52, + component_id="pipe", + label="\u7ba1\u6bb5", + component_type="pipe", + model_type="pipe", + ports="port_a,port_b", + parameters={ + "length": 5.0, + "diameter": 0.02, + "lambda_darcy": 0.02, + "p0": 100000, + "T0": 300, + }, + cells=[ + library_cell( + "body", + "rounded=1;whiteSpace=wrap;html=1;arcSize=12;fillColor=#f5f5f5;" + "strokeColor=#6b7280;strokeWidth=2;", + 14, + 17, + 94, + 16, + ), + library_cell( + "centerline", + "rounded=0;whiteSpace=wrap;html=1;fillColor=#6b7280;strokeColor=none;", + 20, + 24, + 82, + 2, + ), + library_cell( + "flange_a", + f"rounded=0;whiteSpace=wrap;html=1;fillColor=#d1d5db;{STROKE_DARK}", + 7, + 12, + 8, + 26, + ), + library_cell( + "flange_b", + f"rounded=0;whiteSpace=wrap;html=1;fillColor=#d1d5db;{STROKE_DARK}", + 107, + 12, + 8, + 26, + ), + library_cell("port_a", PORT_STYLE, 0, 20, 10, 10), + library_cell("port_b", PORT_STYLE, 112, 20, 10, 10), + library_cell("label", LABEL_STYLE, 29, 36, 64, 14, "\u7ba1\u6bb5"), + ], + ), + component_library_entry( + title="\u5b54\u677f/\u9600\u95e8", + width=96, + height=76, + component_id="orifice", + label="\u5b54\u677f/\u9600\u95e8", + component_type="valve", + model_type="orifice", + ports="port_a,port_b", + parameters={"K": 0.00001, "opening": 1.0}, + cells=[ + library_cell( + "line_a", + "rounded=0;whiteSpace=wrap;html=1;fillColor=#6b7280;strokeColor=none;", + 0, + 35, + 14, + 4, + ), + library_cell( + "tri_left", + "shape=triangle;direction=east;whiteSpace=wrap;html=1;" + "fillColor=#ffe6cc;strokeColor=#d79b00;strokeWidth=2;", + 13, + 22, + 32, + 30, + ), + library_cell( + "tri_right", + "shape=triangle;direction=west;whiteSpace=wrap;html=1;" + "fillColor=#ffe6cc;strokeColor=#d79b00;strokeWidth=2;", + 51, + 22, + 32, + 30, + ), + library_cell( + "plate", + "rounded=0;whiteSpace=wrap;html=1;fillColor=#d79b00;" + "strokeColor=#a86f00;strokeWidth=1;", + 46, + 16, + 4, + 42, + ), + library_cell( + "stem", + "rounded=0;whiteSpace=wrap;html=1;fillColor=#6b7280;strokeColor=none;", + 47, + 7, + 2, + 10, + ), + library_cell( + "handle", + f"rounded=1;whiteSpace=wrap;html=1;arcSize=30;fillColor=#f3f4f6;{STROKE_DARK}", + 38, + 0, + 20, + 8, + ), + library_cell( + "line_b", + "rounded=0;whiteSpace=wrap;html=1;fillColor=#6b7280;strokeColor=none;", + 82, + 35, + 14, + 4, + ), + library_cell("port_a", PORT_STYLE, 0, 32, 10, 10), + library_cell("port_b", PORT_STYLE, 86, 32, 10, 10), + library_cell( + "label", + LABEL_STYLE, + 18, + 58, + 60, + 16, + "\u5b54\u677f/\u9600\u95e8", + ), + ], + ), + component_library_entry( + title="\u4e09\u901a", + width=100, + height=96, + component_id="tee", + label="\u4e09\u901a", + component_type="tee", + model_type="tee", + ports="port_in,port_out1,port_out2", + parameters={}, + cells=[ + library_cell( + "horizontal", + "rounded=1;whiteSpace=wrap;html=1;arcSize=12;fillColor=#f5f5f5;" + "strokeColor=#6b7280;strokeWidth=2;", + 13, + 40, + 74, + 16, + ), + library_cell( + "vertical", + "rounded=1;whiteSpace=wrap;html=1;arcSize=12;fillColor=#f5f5f5;" + "strokeColor=#6b7280;strokeWidth=2;", + 42, + 10, + 16, + 62, + ), + library_cell( + "junction", + "ellipse;whiteSpace=wrap;html=1;fillColor=#ffffff;" + "strokeColor=#6b7280;strokeWidth=2;", + 38, + 36, + 24, + 24, + ), + library_cell("port_out1", PORT_STYLE, 90, 43, 10, 10), + library_cell("port_out2", PORT_STYLE, 0, 43, 10, 10), + library_cell("port_in", PORT_STYLE, 45, 0, 10, 10), + library_cell("label", LABEL_STYLE, 24, 74, 52, 16, "\u4e09\u901a"), + ], + ), ] +DRAWIO_SELECTION_PLUGIN = r""" +(function() { + function postPluginStatus(status, extra) { + try { + var payload = { + event: "systemSimulationPluginStatus", + status: status, + }; + + if (extra) { + for (var key in extra) { + payload[key] = extra[key]; + } + } + + window.parent.postMessage(JSON.stringify(payload), "*"); + } catch (error) { + if (window.console) { + console.warn("System simulation plugin status failed", error); + } + } + } + + function bridgeCellId(cell) { + if (!cell) { + return null; + } + if (cell.id != null) { + return String(cell.id); + } + if (typeof cell.getId === "function") { + return String(cell.getId()); + } + return null; + } + + function bridgeParentOf(cell) { + if (!cell) { + return null; + } + if (typeof cell.getParent === "function") { + return cell.getParent(); + } + return cell.parent || null; + } + + function bridgeCellAndAncestorIds(cell) { + var ids = []; + var seen = {}; + var current = cell; + var depth = 0; + + while (current && depth < 10) { + var id = bridgeCellId(current); + if (id != null && !seen[id]) { + seen[id] = true; + ids.push(id); + } + current = bridgeParentOf(current); + depth += 1; + } + + return ids; + } + + function postSelectionForGraph(graph, source) { + try { + if (!graph || typeof graph.getSelectionCells !== "function") { + return; + } + + var cells = graph.getSelectionCells() || []; + var ids = []; + var details = []; + + for (var i = 0; i < cells.length; i++) { + var cell = cells[i]; + var id = bridgeCellId(cell); + var parent = bridgeParentOf(cell); + var parentId = bridgeCellId(parent); + var ancestorIds = bridgeCellAndAncestorIds(cell); + + for (var j = 0; j < ancestorIds.length; j++) { + ids.push(ancestorIds[j]); + } + details.push({ + id: id, + parentId: parentId, + ancestorIds: ancestorIds, + }); + } + + window.parent.postMessage(JSON.stringify({ + event: "systemSimulationSelection", + source: source || "mxgraph-prototype", + ids: ids, + cells: details, + }), "*"); + } catch (error) { + if (window.console) { + console.warn("System simulation prototype bridge failed", error); + } + } + } + + function scheduleSelectionForGraph(graph, source) { + window.setTimeout(function() { + postSelectionForGraph(graph, source); + }, 0); + } + + function wrapPrototypeMethod(prototype, methodName, graphFromContext) { + var original = prototype && prototype[methodName]; + if ( + typeof original !== "function" || + original.__systemSimulationSelectionWrapped + ) { + return false; + } + + var wrapped = function() { + var result = original.apply(this, arguments); + var graph = graphFromContext ? graphFromContext(this) : this; + scheduleSelectionForGraph(graph, "prototype:" + methodName); + return result; + }; + wrapped.__systemSimulationSelectionWrapped = true; + prototype[methodName] = wrapped; + return true; + } + + function installPrototypeSelectionBridge() { + if (window.__systemSimulationPrototypeSelectionBridgeInstalled) { + return true; + } + + var patched = false; + var graphPrototype = + window.mxGraph && window.mxGraph.prototype ? window.mxGraph.prototype : null; + if (graphPrototype) { + [ + "setSelectionCell", + "setSelectionCells", + "addSelectionCell", + "addSelectionCells", + "removeSelectionCell", + "removeSelectionCells", + "clearSelection", + ].forEach(function(methodName) { + patched = wrapPrototypeMethod(graphPrototype, methodName) || patched; + }); + } + + var selectionPrototype = + window.mxGraphSelectionModel && window.mxGraphSelectionModel.prototype + ? window.mxGraphSelectionModel.prototype + : null; + if (selectionPrototype) { + patched = + wrapPrototypeMethod(selectionPrototype, "changeSelection", function(model) { + return model.graph || null; + }) || patched; + } + + if (patched) { + window.__systemSimulationPrototypeSelectionBridgeInstalled = true; + postPluginStatus("prototype-ready"); + return true; + } + + return false; + } + + function registerSelectionPlugin() { + if (window.__systemSimulationSelectionPluginRegistered) { + return true; + } + + if (!window.Draw || typeof window.Draw.loadPlugin !== "function") { + return false; + } + + window.__systemSimulationSelectionPluginRegistered = true; + window.Draw.loadPlugin(function(ui) { + var graph = ui && ui.editor ? ui.editor.graph : null; + if (!graph) { + postPluginStatus("graph-missing", { + hasUi: !!ui, + hasEditor: !!(ui && ui.editor), + }); + return; + } + + function cellId(cell) { + if (!cell) { + return null; + } + if (cell.id != null) { + return String(cell.id); + } + if (typeof cell.getId === "function") { + return String(cell.getId()); + } + return null; + } + + function parentOf(cell) { + if (!cell) { + return null; + } + if (typeof cell.getParent === "function") { + return cell.getParent(); + } + return cell.parent || null; + } + + function eventName(mxName, fallback) { + var mx = window.mxEvent || (typeof mxEvent !== "undefined" ? mxEvent : null); + return mx && mx[mxName] ? mx[mxName] : fallback; + } + + function collectCellAndAncestorIds(cell) { + var ids = []; + var seen = {}; + var current = cell; + var depth = 0; + + while (current && depth < 10) { + var id = cellId(current); + if (id != null && !seen[id]) { + seen[id] = true; + ids.push(id); + } + current = parentOf(current); + depth += 1; + } + + return ids; + } + + function postSelection() { + try { + var cells = graph.getSelectionCells ? graph.getSelectionCells() : []; + var ids = []; + var details = []; + + for (var i = 0; i < cells.length; i++) { + var cell = cells[i]; + var id = cellId(cell); + var parent = parentOf(cell); + var parentId = cellId(parent); + var ancestorIds = collectCellAndAncestorIds(cell); + + for (var j = 0; j < ancestorIds.length; j++) { + ids.push(ancestorIds[j]); + } + details.push({ + id: id, + parentId: parentId, + ancestorIds: ancestorIds, + }); + } + + window.parent.postMessage(JSON.stringify({ + event: "systemSimulationSelection", + ids: ids, + cells: details, + }), "*"); + } catch (error) { + if (window.console) { + console.warn("System simulation selection plugin failed", error); + } + } + } + + function schedulePostSelection() { + window.setTimeout(postSelection, 0); + } + + function wrapSelectionMethod(methodName) { + var original = graph[methodName]; + if ( + typeof original !== "function" || + original.__systemSimulationSelectionWrapped + ) { + return; + } + + var wrapped = function() { + var result = original.apply(this, arguments); + schedulePostSelection(); + return result; + }; + wrapped.__systemSimulationSelectionWrapped = true; + graph[methodName] = wrapped; + } + + var selectionModel = graph.getSelectionModel ? graph.getSelectionModel() : null; + if (selectionModel && typeof selectionModel.addListener === "function") { + selectionModel.addListener(eventName("CHANGE", "change"), function() { + schedulePostSelection(); + }); + } + + if (typeof graph.addListener === "function") { + graph.addListener(eventName("CLICK", "click"), function() { + schedulePostSelection(); + }); + graph.addListener(eventName("DOUBLE_CLICK", "doubleClick"), function() { + schedulePostSelection(); + }); + } + + [ + "setSelectionCell", + "setSelectionCells", + "addSelectionCell", + "addSelectionCells", + "removeSelectionCell", + "removeSelectionCells", + "clearSelection", + ].forEach(wrapSelectionMethod); + + if (graph.container && typeof graph.container.addEventListener === "function") { + graph.container.addEventListener("mouseup", schedulePostSelection, true); + graph.container.addEventListener("pointerup", schedulePostSelection, true); + graph.container.addEventListener("keyup", schedulePostSelection, true); + } + + window.parent.postMessage(JSON.stringify({ + event: "systemSimulationPluginReady", + }), "*"); + postPluginStatus("ready"); + schedulePostSelection(); + }); + + return true; + } + + if (!installPrototypeSelectionBridge() && !registerSelectionPlugin()) { + var attempts = 0; + postPluginStatus("waiting-for-draw-loadplugin", { attempts: attempts }); + var timer = window.setInterval(function() { + attempts += 1; + var prototypeReady = installPrototypeSelectionBridge(); + var pluginReady = registerSelectionPlugin(); + + if (pluginReady) { + postPluginStatus("registered", { attempts: attempts }); + } + + if (prototypeReady || pluginReady) { + window.clearInterval(timer); + } else if ( + attempts === 10 || + attempts === 50 || + attempts === 100 || + attempts === 300 || + attempts % 600 === 0 + ) { + postPluginStatus("waiting-for-draw-loadplugin", { attempts: attempts }); + } + }, 100); + } +})(); +""" + + class DiagramPayload(BaseModel): xml: str +class UpdateComponentParametersPayload(BaseModel): + xml: str + component_id: str + name: str | None = None + parameters: dict[str, object] + + @app.get("/") def index(request: Request) -> HTMLResponse: html = INDEX_TEMPLATE.read_text(encoding="utf-8") @@ -94,11 +978,38 @@ def component_library() -> Response: ) +@app.get("/api/component-schemas") +def component_schemas() -> dict[str, object]: + return COMPONENT_SCHEMAS + + +@app.get("/api/drawio-selection-plugin.js") +def drawio_selection_plugin() -> Response: + return Response( + content=DRAWIO_SELECTION_PLUGIN, + media_type="application/javascript", + headers={"Cache-Control": "no-store"}, + ) + + +@app.post("/api/diagram-components") +def diagram_components(payload: DiagramPayload) -> dict[str, object]: + try: + model = read_graph_model(payload.xml) + root = get_graph_root(model) + components = list(read_components(root)) + except ValueError as exc: + raise HTTPException(status_code=400, detail=str(exc)) from exc + + return {"components": components} + + @app.api_route("/drawio", methods=["GET", "POST", "HEAD", "OPTIONS"]) @app.api_route("/drawio/{path:path}", methods=["GET", "POST", "HEAD", "OPTIONS"]) async def drawio_proxy(request: Request, path: str = "") -> Response: upstream_url = build_drawio_upstream_url(path, request.url.query) body = await request.body() + normalized_path = normalize_drawio_proxy_path(path) headers = { key: value for key, value in request.headers.items() @@ -108,6 +1019,9 @@ async def drawio_proxy(request: Request, path: str = "") -> Response: "connection", "content-length", "accept-encoding", + "if-modified-since", + "if-none-match", + "if-range", } } @@ -121,10 +1035,25 @@ async def drawio_proxy(request: Request, path: str = "") -> Response: try: with urllib.request.urlopen(upstream_request, timeout=30) as upstream_response: content = upstream_response.read() + content = inject_drawio_selection_plugin( + content, + upstream_response.headers, + request.url.query, + normalized_path, + ) + response_headers = filter_proxy_headers(upstream_response.headers) + if is_drawio_injection_target( + normalized_path, + upstream_response.headers, + request.url.query, + ): + response_headers["Cache-Control"] = "no-store" + response_headers.pop("ETag", None) + response_headers.pop("Last-Modified", None) return Response( content=content, status_code=upstream_response.status, - headers=filter_proxy_headers(upstream_response.headers), + headers=response_headers, ) except urllib.error.HTTPError as exc: return Response( @@ -150,6 +1079,27 @@ def export_system_xml(payload: DiagramPayload) -> Response: return Response(content=system_xml, media_type="application/xml") +@app.post("/api/update-component-parameters") +def update_component_parameters_endpoint( + payload: UpdateComponentParametersPayload, +) -> Response: + try: + model = read_graph_model(payload.xml) + update_component_parameters( + model, + payload.component_id, + payload.name, + payload.parameters, + ) + except ValueError as exc: + raise HTTPException(status_code=400, detail=str(exc)) from exc + + return Response( + content=ET.tostring(model, encoding="utf-8"), + media_type="application/xml", + ) + + def default_drawio_base_url(request: Request) -> str: hostname = request.url.hostname or "127.0.0.1" if ":" in hostname and not hostname.startswith("["): @@ -166,6 +1116,10 @@ def build_drawio_upstream_url(path: str, query: str) -> str: return target +def normalize_drawio_proxy_path(path: str) -> str: + return path.replace("\\", "/").lstrip("/") + + def filter_proxy_headers(headers) -> dict[str, str]: excluded = { "connection", @@ -186,6 +1140,64 @@ def filter_proxy_headers(headers) -> dict[str, str]: } +def is_drawio_injection_target(path: str, headers, query: str) -> bool: + content_type = headers.get("Content-Type", "") + if "text/html" in content_type and "embed=1" in query: + return True + return path == DRAWIO_APP_SCRIPT_PATH and is_javascript_content(content_type) + + +def is_javascript_content(content_type: str) -> bool: + lowered = content_type.lower() + return "javascript" in lowered or "ecmascript" in lowered + + +def inject_drawio_selection_plugin( + content: bytes, + headers, + query: str, + path: str = "", +) -> bytes: + content_type = headers.get("Content-Type", "") + + if path == DRAWIO_APP_SCRIPT_PATH and is_javascript_content(content_type): + return append_drawio_selection_plugin_to_app_script(content) + + if "text/html" not in content_type or "embed=1" not in query: + return content + + try: + html = content.decode("utf-8") + except UnicodeDecodeError: + return content + + plugin_tag = '' + if plugin_tag in html: + return content + + main_script = '' + if main_script in html: + html = html.replace(main_script, f"{plugin_tag}\n{main_script}", 1) + elif "" in html: + html = html.replace("", f"{plugin_tag}\n", 1) + else: + return content + + return html.encode("utf-8") + + +def append_drawio_selection_plugin_to_app_script(content: bytes) -> bytes: + marker = b"systemSimulationSelectionPluginRegistered" + if marker in content: + return content + + return ( + content + + b"\n;\n/* SystemSimulationApp selection bridge */\n" + + DRAWIO_SELECTION_PLUGIN.encode("utf-8") + ) + + def read_graph_model(xml_text: str) -> ET.Element: xml_text = xml_text.strip() if not xml_text: @@ -250,11 +1262,15 @@ def parse_graph_model(xml_text: str) -> ET.Element: return model -def build_system_xml(model: ET.Element) -> bytes: +def get_graph_root(model: ET.Element) -> ET.Element: root = next((child for child in model if tag_name(child) == "root"), None) if root is None: raise ValueError("mxGraphModel does not contain a root node.") + return root + +def build_system_xml(model: ET.Element) -> bytes: + root = get_graph_root(model) system = ET.Element("System") components_node = ET.SubElement(system, "Components") connections_node = ET.SubElement(system, "Connections") @@ -283,10 +1299,135 @@ def build_system_xml(model: ET.Element) -> bytes: return ET.tostring(system, encoding="utf-8", xml_declaration=True) +def update_component_parameters( + model: ET.Element, + component_id: str, + name: str | None, + parameters: dict[str, object], +) -> None: + root = get_graph_root(model) + component = find_component_element(root, component_id) + if component is None: + raise ValueError(f"Component '{component_id}' was not found.") + + if name is not None: + normalized_name = normalize_component_name(name) + component.set("name", normalized_name) + component.set("label", normalized_name) + update_component_label_cells(root, component_id, normalized_name) + + schema = schema_for_component(component) + allowed_parameters = set(schema.get("parameters", {})) + if not allowed_parameters: + if parameters: + raise ValueError("This component does not have configurable parameters.") + return + + for name, value in parameters.items(): + if name not in allowed_parameters: + raise ValueError(f"Parameter '{name}' is not supported by this component.") + + definition = schema["parameters"][name] + component.set(name, normalize_parameter_value(name, value, definition)) + + +def normalize_component_name(name: str) -> str: + normalized = name.strip() + if not normalized: + raise ValueError("Component name cannot be empty.") + return normalized + + +def update_component_label_cells( + root: ET.Element, + component_id: str, + label: str, +) -> None: + for cell in root.iter(): + if tag_name(cell) != "mxCell": + continue + if cell.attrib.get("parent") != component_id: + continue + style = cell.attrib.get("style", "") + if style.startswith("text") or "text;" in style: + cell.set("value", label) + + +def find_component_element(root: ET.Element, component_id: str) -> ET.Element | None: + for element in root: + if tag_name(element) != "object": + continue + if "portName" in element.attrib: + continue + if element.attrib.get("id") == component_id: + return element + return None + + +def schema_for_component(component: ET.Element) -> dict[str, object]: + model_type = component.attrib.get("modelType") + component_type = component.attrib.get("componentType") + if model_type in COMPONENT_SCHEMAS: + return COMPONENT_SCHEMAS[model_type] + if component_type in COMPONENT_SCHEMAS: + return COMPONENT_SCHEMAS[component_type] + return {"label": component_type or "component", "parameters": {}} + + +def normalize_parameter_value( + name: str, + value: object, + definition: dict[str, object], +) -> str: + value_type = definition.get("type", "string") + if value is None or value == "": + raise ValueError(f"Parameter '{name}' cannot be empty.") + + if value_type == "float": + try: + number = float(value) + except (TypeError, ValueError) as exc: + raise ValueError(f"Parameter '{name}' must be a number.") from exc + if not math.isfinite(number): + raise ValueError(f"Parameter '{name}' must be finite.") + validate_number_range(name, number, definition) + return f"{number:g}" + + if value_type == "int": + try: + number = int(float(value)) + except (TypeError, ValueError) as exc: + raise ValueError(f"Parameter '{name}' must be an integer.") from exc + validate_number_range(name, number, definition) + return str(number) + + if value_type == "bool": + if isinstance(value, bool): + return "true" if value else "false" + return "true" if str(value).lower() in {"1", "true", "yes"} else "false" + + return str(value) + + +def validate_number_range( + name: str, + value: float, + definition: dict[str, object], +) -> None: + min_value = definition.get("min") + max_value = definition.get("max") + if min_value is not None and value < float(min_value): + raise ValueError(f"Parameter '{name}' must be >= {min_value}.") + if max_value is not None and value > float(max_value): + raise ValueError(f"Parameter '{name}' must be <= {max_value}.") + + def read_components(root: ET.Element) -> Iterable[dict]: for element in root: if tag_name(element) != "object": continue + if "portName" in element.attrib: + continue cell = next((child for child in element if tag_name(child) == "mxCell"), None) if cell is None or cell.attrib.get("vertex") != "1": @@ -311,20 +1452,182 @@ def read_components(root: ET.Element) -> Iterable[dict]: def read_connections(root: ET.Element) -> Iterable[dict[str, str]]: + endpoints = read_connection_endpoints(root) + for cell in root.iter(): if tag_name(cell) != "mxCell" or cell.attrib.get("edge") != "1": continue - source = cell.attrib.get("source") - target = cell.attrib.get("target") - if not source or not target: + source_cell = cell.attrib.get("source") + target_cell = cell.attrib.get("target") + if not source_cell or not target_cell: continue - edge = {"source": source, "target": target} + source = endpoints.get(source_cell, {"component": source_cell}) + target = endpoints.get(target_cell, {"component": target_cell}) + edge = { + "source": source["component"], + "target": target["component"], + } + if "port" in source: + edge["sourcePort"] = source["port"] + edge["sourceCell"] = source_cell + if "port" in target: + edge["targetPort"] = target["port"] + edge["targetCell"] = target_cell if "id" in cell.attrib: edge["id"] = cell.attrib["id"] yield edge +def read_connection_endpoints(root: ET.Element) -> dict[str, dict[str, str]]: + endpoints: dict[str, dict[str, str]] = {} + component_ports: dict[str, list[str]] = {} + inferred_port_cells: dict[str, list[dict[str, object]]] = {} + + for element in root.iter(): + if tag_name(element) != "object": + continue + + component_id = element.attrib.get("id") + ports = element.attrib.get("ports") + if component_id and ports and "portName" not in element.attrib: + component_ports[component_id] = [ + port.strip() + for port in ports.split(",") + if port.strip() + ] + + endpoint_id = element.attrib.get("id") + port_name = element.attrib.get("portName") + if not endpoint_id or not port_name: + continue + + cell = next((child for child in element if tag_name(child) == "mxCell"), None) + if cell is None: + continue + + component_id = cell.attrib.get("parent") + if not component_id: + continue + + endpoint = { + "component": component_id, + "port": port_name, + } + endpoints[endpoint_id] = endpoint + + cell_id = cell.attrib.get("id") + if cell_id: + endpoints[cell_id] = endpoint + + for cell in root.iter(): + if tag_name(cell) != "mxCell": + continue + if cell.attrib.get("vertex") != "1" or not is_connectable(cell): + continue + + cell_id = cell.attrib.get("id") + component_id = cell.attrib.get("parent") + if not cell_id or not component_id or cell_id in endpoints: + continue + if component_id not in component_ports: + continue + + geometry = next( + (child for child in cell if tag_name(child) == "mxGeometry"), + None, + ) + inferred_port_cells.setdefault(component_id, []).append( + { + "id": cell_id, + "x": read_float_attribute(geometry, "x"), + "y": read_float_attribute(geometry, "y"), + "width": read_float_attribute(geometry, "width"), + "height": read_float_attribute(geometry, "height"), + } + ) + + for component_id, cells in inferred_port_cells.items(): + for cell_id, port_name in infer_port_cell_names( + cells, + component_ports.get(component_id, []), + ).items(): + endpoints.setdefault( + cell_id, + { + "component": component_id, + "port": port_name, + }, + ) + + return endpoints + + +def is_connectable(cell: ET.Element) -> bool: + connectable = cell.attrib.get("connectable") + if connectable is not None: + return connectable == "1" + + style = cell.attrib.get("style", "") + style_flags = { + item.split("=", 1)[0]: item.split("=", 1)[1] + for item in style.split(";") + if "=" in item + } + return style_flags.get("connectable") == "1" + + +def infer_port_cell_names( + cells: list[dict[str, object]], + port_names: list[str], +) -> dict[str, str]: + if len(cells) != len(port_names): + return {} + + if len(port_names) == 1: + return {str(cells[0]["id"]): port_names[0]} + + if set(port_names) == {"port_a", "port_b"} and len(cells) == 2: + sorted_cells = sorted(cells, key=lambda cell: cell_center(cell)[0]) + return { + str(sorted_cells[0]["id"]): "port_a", + str(sorted_cells[1]["id"]): "port_b", + } + + if set(port_names) == {"port_in", "port_out1", "port_out2"} and len(cells) == 3: + sorted_by_y = sorted(cells, key=lambda cell: cell_center(cell)[1]) + top_cell = sorted_by_y[0] + remaining = sorted(sorted_by_y[1:], key=lambda cell: cell_center(cell)[0]) + return { + str(top_cell["id"]): "port_in", + str(remaining[0]["id"]): "port_out2", + str(remaining[1]["id"]): "port_out1", + } + + sorted_cells = sorted(cells, key=lambda cell: (cell_center(cell)[1], cell_center(cell)[0])) + return { + str(cell["id"]): port_names[index] + for index, cell in enumerate(sorted_cells) + } + + +def cell_center(cell: dict[str, object]) -> tuple[float, float]: + return ( + float(cell["x"]) + float(cell["width"]) / 2.0, + float(cell["y"]) + float(cell["height"]) / 2.0, + ) + + +def read_float_attribute(element: ET.Element | None, name: str) -> float: + if element is None: + return 0.0 + + try: + return float(element.attrib.get(name, "0")) + except ValueError: + return 0.0 + + def tag_name(element: ET.Element) -> str: return element.tag.rsplit("}", 1)[-1] diff --git a/app/static/index.html b/app/static/index.html index 45942ec..808d41e 100644 --- a/app/static/index.html +++ b/app/static/index.html @@ -12,14 +12,14 @@ body { margin: 0; height: 100vh; - color: #18202a; + color: #1f2933; font-family: system-ui, -apple-system, BlinkMacSystemFont, "Segoe UI", sans-serif; - background: #f6f8fb; + background: #f4f7fb; } .app { @@ -34,7 +34,7 @@ justify-content: space-between; gap: 16px; padding: 0 14px; - border-bottom: 1px solid #d8dee8; + border-bottom: 1px solid #d6dde8; background: #ffffff; } @@ -51,11 +51,11 @@ .actions { display: flex; align-items: center; - gap: 10px; + gap: 8px; } .status { - max-width: 220px; + max-width: 260px; overflow: hidden; color: #5d6877; font-size: 13px; @@ -66,27 +66,218 @@ button { height: 32px; padding: 0 12px; - border: 1px solid #1b6fb8; + border: 1px solid #1d6fb8; border-radius: 6px; color: #ffffff; font: inherit; font-size: 14px; - background: #1b6fb8; + background: #1d6fb8; cursor: pointer; } + button.secondary { + border-color: #c7d0dc; + color: #273646; + background: #ffffff; + } + button:disabled { - border-color: #98a2af; - background: #98a2af; + border-color: #a7b0bd; + color: #ffffff; + background: #a7b0bd; cursor: wait; } + .workbench { + display: grid; + grid-template-columns: minmax(0, 1fr) 340px; + min-height: 0; + } + + .canvas { + min-width: 0; + min-height: 0; + } + iframe { display: block; width: 100%; height: 100%; border: 0; } + + .side-panel { + display: grid; + grid-template-rows: auto minmax(140px, 32%) minmax(0, 1fr); + min-height: 0; + border-left: 1px solid #d6dde8; + background: #ffffff; + } + + .panel-header, + .panel-section { + padding: 12px; + } + + .panel-header { + display: flex; + align-items: center; + justify-content: space-between; + gap: 8px; + border-bottom: 1px solid #e3e8ef; + } + + .panel-title { + font-size: 14px; + font-weight: 650; + } + + .component-list { + min-height: 0; + overflow: auto; + border-bottom: 1px solid #e3e8ef; + } + + .component-row { + display: block; + width: 100%; + height: auto; + min-height: 54px; + padding: 9px 12px; + border: 0; + border-bottom: 1px solid #edf1f5; + border-radius: 0; + color: #1f2933; + text-align: left; + background: #ffffff; + } + + .component-row:hover, + .component-row.active { + background: #eef6ff; + } + + .component-name { + overflow: hidden; + font-size: 14px; + font-weight: 620; + text-overflow: ellipsis; + white-space: nowrap; + } + + .component-meta { + margin-top: 3px; + overflow: hidden; + color: #667387; + font-size: 12px; + text-overflow: ellipsis; + white-space: nowrap; + } + + .empty-state { + padding: 18px 12px; + color: #667387; + font-size: 13px; + } + + .parameter-panel { + min-height: 0; + overflow: auto; + } + + .selected-heading { + margin: 0; + overflow-wrap: anywhere; + font-size: 15px; + font-weight: 650; + } + + .selected-meta { + margin-top: 4px; + overflow-wrap: anywhere; + color: #667387; + font-size: 12px; + } + + .field-grid { + display: grid; + gap: 12px; + margin-top: 14px; + } + + .field { + display: grid; + gap: 5px; + } + + .name-field { + margin-top: 14px; + } + + .field label { + display: flex; + align-items: baseline; + justify-content: space-between; + gap: 8px; + color: #334155; + font-size: 13px; + font-weight: 590; + } + + .field-key { + color: #7b8797; + font-size: 12px; + font-weight: 400; + } + + .input-row { + display: grid; + grid-template-columns: minmax(0, 1fr) auto; + align-items: center; + gap: 8px; + } + + input { + width: 100%; + height: 34px; + padding: 0 9px; + border: 1px solid #c7d0dc; + border-radius: 6px; + color: #1f2933; + font: inherit; + font-size: 14px; + background: #ffffff; + } + + input:focus { + border-color: #1d6fb8; + outline: 2px solid rgba(29, 111, 184, 0.16); + } + + .unit { + min-width: 26px; + color: #667387; + font-size: 12px; + } + + .form-actions { + display: flex; + justify-content: flex-end; + margin-top: 16px; + } + + @media (max-width: 980px) { + .workbench { + grid-template-columns: minmax(0, 1fr); + grid-template-rows: minmax(0, 1fr) 360px; + } + + .side-panel { + grid-template-rows: auto minmax(120px, 42%) minmax(0, 1fr); + border-top: 1px solid #d6dde8; + border-left: 0; + } + } @@ -95,23 +286,74 @@
System Simulation App
加载中 +
- +
+
+ +
+ +