diff --git a/PythonModels/systems/test_mql_closure.py b/PythonModels/systems/test_mql_closure.py new file mode 100644 index 0000000..ece137d --- /dev/null +++ b/PythonModels/systems/test_mql_closure.py @@ -0,0 +1,157 @@ +from __future__ import annotations + +from dataclasses import dataclass +from typing import Callable + +from PythonModels.components.amesim_pneumatic import ( + AmesimPneumaticOrifice, + AmesimPneumaticVolume, +) +from PythonModels.core.medium import ThermodynamicProperties +from PythonModels.core.state import VolumeState + + +@dataclass(frozen=True) +class TestMqlPneumaticBranchComponents: + name: str + upstream_volume: AmesimPneumaticVolume + orifice: AmesimPneumaticOrifice + downstream_volume: AmesimPneumaticVolume + + +@dataclass(frozen=True) +class TestMqlPneumaticBranchState: + name: str + upstream: ThermodynamicProperties + downstream: ThermodynamicProperties + flow: float + upstream_inlet_h: float + downstream_inlet_h: float + + +@dataclass(frozen=True) +class TestMqlPneumaticSnapshot: + branch: TestMqlPneumaticBranchState + + @property + def flow(self) -> float: + return self.branch.flow + + @property + def upstream(self) -> ThermodynamicProperties: + return self.branch.upstream + + @property + def downstream(self) -> ThermodynamicProperties: + return self.branch.downstream + + +class TestMqlPneumaticClosure: + """Minimal test_mql pneumatic closure following the existing Testmodel pattern. + + The canonical branch flow is positive from ``upstream_volume`` through the + orifice port_a/port_b into ``downstream_volume``. Port ``m_flow`` values are + written with Modelica-style signs: positive means flow into that component. + """ + + def __init__( + self, + *, + components: TestMqlPneumaticBranchComponents, + initial_state_vector: Callable[[], list[float]], + apply_state_vector: Callable[[list[float]], None], + ) -> None: + self.components = components + self._initial_state_vector = initial_state_vector + self._apply_state_vector = apply_state_vector + + def initial_state_vector(self) -> list[float]: + return self._initial_state_vector() + + def apply_state_vector(self, values: list[float]) -> None: + self._apply_state_vector(values) + + def snapshot( + self, + state_vector: list[float] | None = None, + ) -> TestMqlPneumaticSnapshot: + if state_vector is not None: + self._apply_state_vector(state_vector) + upstream = self.components.upstream_volume.properties() + downstream = self.components.downstream_volume.properties() + flow = self._solve_branch_flow(upstream, downstream) + upstream_inlet_h = self.components.upstream_volume.connection_inlet_enthalpy( + port_m_flow=-flow, + connected_h=downstream.h, + internal_h=upstream.h, + ) + downstream_inlet_h = self.components.downstream_volume.connection_inlet_enthalpy( + port_m_flow=flow, + connected_h=upstream.h, + internal_h=downstream.h, + ) + branch = TestMqlPneumaticBranchState( + name=self.components.name, + upstream=upstream, + downstream=downstream, + flow=flow, + upstream_inlet_h=upstream_inlet_h, + downstream_inlet_h=downstream_inlet_h, + ) + self._write_port_states(branch) + return TestMqlPneumaticSnapshot(branch=branch) + + def _solve_branch_flow( + self, + upstream: ThermodynamicProperties, + downstream: ThermodynamicProperties, + ) -> float: + upstream_temperature = upstream.T if upstream.p >= downstream.p else downstream.T + return self.components.orifice.mass_flow( + upstream.p, + downstream.p, + upstream_temperature, + ) + + def _write_port_states(self, branch: TestMqlPneumaticBranchState) -> None: + upstream_volume = self.components.upstream_volume + downstream_volume = self.components.downstream_volume + orifice = self.components.orifice + + upstream_volume.port_a.p = branch.upstream.p + upstream_volume.port_a.m_flow = -branch.flow + upstream_volume.port_a.h_outflow = branch.upstream.h + + orifice.port_a.p = branch.upstream.p + orifice.port_a.m_flow = branch.flow + orifice.port_a.h_outflow = branch.upstream.h + orifice.port_b.p = branch.downstream.p + orifice.port_b.m_flow = -branch.flow + orifice.port_b.h_outflow = branch.downstream.h + + downstream_volume.port_a.p = branch.downstream.p + downstream_volume.port_a.m_flow = branch.flow + downstream_volume.port_a.h_outflow = branch.downstream.h + + def branch_derivatives( + self, + snapshot: TestMqlPneumaticSnapshot, + ) -> tuple[VolumeState, VolumeState]: + flow = snapshot.flow + upstream_derivative = self.components.upstream_volume.derivatives( + inlet_h=snapshot.branch.upstream_inlet_h, + m_flow=-flow, + ) + downstream_derivative = self.components.downstream_volume.derivatives( + inlet_h=snapshot.branch.downstream_inlet_h, + m_flow=flow, + ) + return upstream_derivative, downstream_derivative + + def rhs(self, state_vector: list[float]) -> list[float]: + snapshot = self.snapshot(state_vector) + upstream_derivative, downstream_derivative = self.branch_derivatives(snapshot) + return [ + *upstream_derivative.as_vector(), + *downstream_derivative.as_vector(), + ] diff --git a/tests/test_test_mql_closure.py b/tests/test_test_mql_closure.py new file mode 100644 index 0000000..dc698fe --- /dev/null +++ b/tests/test_test_mql_closure.py @@ -0,0 +1,138 @@ +from __future__ import annotations + +import unittest + +from PythonModels.components.amesim_pneumatic import ( + AmesimPneumaticOrifice, + AmesimPneumaticVolume, +) +from PythonModels.systems.test_mql_closure import ( + TestMqlPneumaticBranchComponents, + TestMqlPneumaticClosure, +) + + +class TestMqlPneumaticClosureTests(unittest.TestCase): + def _make_closure( + self, + *, + upstream_pressure: float = 15.3e6, + downstream_pressure: float = 1.0e5, + ) -> TestMqlPneumaticClosure: + upstream = AmesimPneumaticVolume.from_liters( + name="source_chamber", + volume_liters=57.0, + p0=upstream_pressure, + T0=293.15, + ) + downstream = AmesimPneumaticVolume.from_liters( + name="target_chamber", + volume_liters=15.0, + p0=downstream_pressure, + T0=293.15, + ) + orifice = AmesimPneumaticOrifice.from_mm2( + name="branch_orifice", + area_mm2=78.5, + flow_coefficient=0.9, + ) + + def initial_state_vector() -> list[float]: + return [*upstream.get_state_vector(), *downstream.get_state_vector()] + + def apply_state_vector(values: list[float]) -> None: + if len(values) != 4: + raise ValueError("two-volume closure state vector requires four values") + upstream.set_state_vector(values[:2]) + downstream.set_state_vector(values[2:]) + + return TestMqlPneumaticClosure( + components=TestMqlPneumaticBranchComponents( + name="source_to_target", + upstream_volume=upstream, + orifice=orifice, + downstream_volume=downstream, + ), + initial_state_vector=initial_state_vector, + apply_state_vector=apply_state_vector, + ) + + def test_snapshot_writes_modelica_style_port_flow_signs(self) -> None: + closure = self._make_closure() + + snapshot = closure.snapshot() + + self.assertGreater(snapshot.flow, 0.0) + self.assertAlmostEqual( + closure.components.upstream_volume.port_a.m_flow, + -snapshot.flow, + ) + self.assertAlmostEqual( + closure.components.orifice.port_a.m_flow, + snapshot.flow, + ) + self.assertAlmostEqual( + closure.components.orifice.port_b.m_flow, + -snapshot.flow, + ) + self.assertAlmostEqual( + closure.components.downstream_volume.port_a.m_flow, + snapshot.flow, + ) + self.assertAlmostEqual( + closure.components.orifice.port_a.p, + snapshot.upstream.p, + ) + self.assertAlmostEqual( + closure.components.orifice.port_b.p, + snapshot.downstream.p, + ) + + def test_rhs_uses_canonical_flow_for_volume_mass_balance(self) -> None: + closure = self._make_closure() + state = closure.initial_state_vector() + + snapshot = closure.snapshot(state) + rhs = closure.rhs(state) + + self.assertEqual(len(rhs), 4) + self.assertAlmostEqual(rhs[0], -snapshot.flow) + self.assertAlmostEqual(rhs[2], snapshot.flow) + self.assertAlmostEqual(rhs[0] + rhs[2], 0.0) + self.assertAlmostEqual(rhs[1], -snapshot.flow * snapshot.upstream.h) + self.assertAlmostEqual(rhs[3], snapshot.flow * snapshot.upstream.h) + + def test_reverse_pressure_reverses_canonical_flow_and_port_signs(self) -> None: + closure = self._make_closure( + upstream_pressure=1.0e5, + downstream_pressure=15.3e6, + ) + + snapshot = closure.snapshot() + rhs = closure.rhs(closure.initial_state_vector()) + + self.assertLess(snapshot.flow, 0.0) + self.assertGreater(closure.components.upstream_volume.port_a.m_flow, 0.0) + self.assertLess(closure.components.downstream_volume.port_a.m_flow, 0.0) + self.assertAlmostEqual(rhs[0], -snapshot.flow) + self.assertAlmostEqual(rhs[2], snapshot.flow) + self.assertAlmostEqual(rhs[0] + rhs[2], 0.0) + self.assertAlmostEqual(rhs[1], -snapshot.flow * snapshot.downstream.h) + self.assertAlmostEqual(rhs[3], snapshot.flow * snapshot.downstream.h) + + def test_snapshot_can_apply_explicit_state_vector(self) -> None: + closure = self._make_closure() + state = closure.initial_state_vector() + state[0] *= 0.99 + + snapshot = closure.snapshot(state) + + self.assertAlmostEqual( + closure.components.upstream_volume.state.m, + state[0], + ) + self.assertGreater(snapshot.upstream.p, snapshot.downstream.p) + + +if __name__ == "__main__": + unittest.main()