PythonModels/systems/test_mql.py

@@ -3879,6 +3879,33 @@ _ELASTIC_ENDSTOP_FORCE_BINDINGS = (
 )
 
 
+def _default_full_state_data_paths() -> tuple[str, ...]:
+    chamber_paths: list[str] = []
+    piston_paths: list[str] = []
+    for _mass_alias, piston_alias, chamber_field in _PISTON_FORCE_BINDINGS:
+        chamber_alias = _PNCH012_ALIAS_BY_SNAPSHOT_FIELD[chamber_field]
+        chamber_paths.extend([f"press@{chamber_alias}", f"vol@{chamber_alias}"])
+        piston_paths.extend([f"vol1@{piston_alias}", f"vvol1@{piston_alias}"])
+    mass_paths = [
+        f"{signal}@mass_friction_endstops_{index}"
+        for index in range(10, 20)
+        for signal in ("x1", "v1", "acc1")
+    ]
+    return tuple([*chamber_paths, *piston_paths, *mass_paths])
+
+
+_PNCH012_ALIAS_BY_SNAPSHOT_FIELD = {
+    "p4_port3_remote_primary_chamber": "pn_c1_8",
+    "p4_primary_chamber": "pn_c1_9",
+    "p4_port1_remote_primary_chamber": "pn_c1_10",
+    "p4_port1_far_primary_chamber": "pn_c1_11",
+    "p4_port1_next_primary_chamber": "pn_c1_12",
+    "p4_bridge_primary_chamber": "pn_c1_13",
+    "p4_port3_next_primary_chamber": "pn_c1_14",
+    "p4_port3_far_primary_chamber": "pn_c1_15",
+}
+
+
 class TestMqlFullStateClosure:
     def __init__(self, *, pneumatic_closure: object, mechanical_closure: object) -> None:
         self.pneumatic_closure = pneumatic_closure
@@ -3961,6 +3988,116 @@ class TestMqlFullStateClosure:
                 ),
             )
 
+    def data_path_values(
+        self,
+        *,
+        time_s: float,
+        state_vector: list[float],
+        data_paths: tuple[str, ...] | list[str] | None = None,
+    ) -> dict[str, float]:
+        from PythonModels.systems.test_mql_pneumatic import pressure_to_amesim_gauge_pa
+
+        selected_paths = tuple(data_paths) if data_paths is not None else _default_full_state_data_paths()
+        snapshot = self.snapshot(state_vector)
+        rhs = self.rhs_at(time_s, state_vector)
+        chamber_properties_by_alias = {}
+        chamber_by_alias = {}
+        for _mass_alias, _piston_alias, chamber_field in _PISTON_FORCE_BINDINGS:
+            chamber = getattr(self.pneumatic_closure.components, chamber_field)
+            chamber_by_alias[chamber.name] = chamber
+            chamber_properties_by_alias[chamber.name] = getattr(snapshot.pneumatic, chamber_field)
+
+        state_by_alias = {state.alias: state for state in snapshot.mechanical.states}
+        acceleration_by_mass_alias = {
+            alias: rhs[self.pneumatic_state_count + 2 * index]
+            for index, alias in enumerate(self.mechanical_closure.mass_aliases)
+        }
+        values_by_data_path = {}
+        for data_path in selected_paths:
+            signal, alias = self._split_data_path(data_path)
+            if alias in chamber_properties_by_alias:
+                chamber_properties = chamber_properties_by_alias[alias]
+                chamber = chamber_by_alias[alias]
+                if signal == "press":
+                    values_by_data_path[data_path] = pressure_to_amesim_gauge_pa(
+                        chamber_properties.p
+                    )
+                elif signal == "temp":
+                    values_by_data_path[data_path] = chamber_properties.T
+                elif signal == "vol":
+                    values_by_data_path[data_path] = chamber.volume_cm3()
+                else:
+                    raise KeyError(data_path)
+            elif alias in self.mechanical_closure.assembly.pistons:
+                piston = self.mechanical_closure.assembly.pistons[alias]
+                kinematics = snapshot.mechanical.piston_kinematics_by_alias[alias]
+                geometry = piston.geometry()
+                if signal == "vol1":
+                    values_by_data_path[data_path] = geometry.chamber_volume_cm3(
+                        kinematics.port_3_displacement_m,
+                        kinematics.port_2_displacement_m,
+                    )
+                elif signal == "vvol1":
+                    values_by_data_path[data_path] = geometry.chamber_volume_rate_l_min(
+                        kinematics.port_3_velocity_m_s,
+                        kinematics.port_2_velocity_m_s,
+                    )
+                elif signal == "length":
+                    values_by_data_path[data_path] = geometry.chamber_length_mm(
+                        kinematics.port_3_displacement_m,
+                        kinematics.port_2_displacement_m,
+                    )
+                elif signal == "x4":
+                    values_by_data_path[data_path] = kinematics.port_3_displacement_m
+                elif signal == "x5":
+                    values_by_data_path[data_path] = kinematics.port_2_displacement_m
+                elif signal == "v4":
+                    values_by_data_path[data_path] = kinematics.port_3_velocity_m_s
+                elif signal == "v5":
+                    values_by_data_path[data_path] = kinematics.port_2_velocity_m_s
+                else:
+                    raise KeyError(data_path)
+            elif alias in state_by_alias:
+                mass_state = state_by_alias[alias]
+                if signal == "x1":
+                    values_by_data_path[data_path] = mass_state.displacement_m
+                elif signal == "v1":
+                    values_by_data_path[data_path] = mass_state.velocity_m_s
+                elif signal == "acc1":
+                    values_by_data_path[data_path] = acceleration_by_mass_alias[alias]
+                else:
+                    raise KeyError(data_path)
+            else:
+                raise KeyError(data_path)
+        return values_by_data_path
+
+    def data_path_series(
+        self,
+        *,
+        times: tuple[float, ...] | list[float],
+        state_rows: list[list[float]],
+        data_paths: tuple[str, ...] | list[str] | None = None,
+    ) -> dict[str, list[float]]:
+        selected_paths = tuple(data_paths) if data_paths is not None else _default_full_state_data_paths()
+        series = {data_path: [] for data_path in selected_paths}
+        for sample_index, time_value in enumerate(times):
+            state_vector = [row[sample_index] for row in state_rows]
+            values = self.data_path_values(
+                time_s=float(time_value),
+                state_vector=state_vector,
+                data_paths=selected_paths,
+            )
+            for data_path in selected_paths:
+                series[data_path].append(values[data_path])
+        return series
+
+    @staticmethod
+    def _split_data_path(data_path: str) -> tuple[str, str]:
+        if "@" not in data_path:
+            raise KeyError(data_path)
+        signal, alias = data_path.split("@", 1)
+        return signal, alias
+
     def _force_by_mass_alias(
         self,
         time_s: float,
@@ -4406,6 +4543,43 @@ class TestMqlSystem:
             t_eval=t_eval,
         )
 
+    def simulate_full_state_series_from_spec(
+        self,
+        spec,
+        *,
+        inlet_node_pressure_pa: float,
+        resistance_boundary_pressure_pa: float,
+        inlet_node_temperature_k: float = 293.15,
+        resistance_boundary_temperature_k: float = 293.15,
+        config: SolveIVPConfig | None = None,
+        t_eval: list[float] | None = None,
+        data_paths: tuple[str, ...] | list[str] | None = None,
+    ) -> TestMqlSimulationResult:
+        closure = self.full_state_closure_from_spec(
+            spec,
+            inlet_node_pressure_pa=inlet_node_pressure_pa,
+            resistance_boundary_pressure_pa=resistance_boundary_pressure_pa,
+            inlet_node_temperature_k=inlet_node_temperature_k,
+            resistance_boundary_temperature_k=resistance_boundary_temperature_k,
+        )
+        solution = integrate_ode(
+            rhs=lambda t, state: closure.rhs_at(t, state),
+            initial_state=closure.initial_state_vector(),
+            config=config or SolveIVPConfig(t_stop=1.0e-4, max_step=1.0e-5),
+            t_eval=t_eval,
+        )
+        times = [float(value) for value in solution.t]
+        state_rows = [list(row) for row in solution.y]
+        series = {
+            "time": times,
+            **closure.data_path_series(
+                times=times,
+                state_rows=state_rows,
+                data_paths=data_paths,
+            ),
+        }
+        return TestMqlSimulationResult(t=times, y=state_rows, series=series)
+
     @staticmethod
     def pneumatic_branch_spec(
         *,

tests/test_test_mql_pnl0001_segment.py

@@ -1,11 +1,19 @@
 from __future__ import annotations
 
 import unittest
+from pathlib import Path
 
 from PythonModels.core.solver import SolveIVPConfig
+from PythonModels.reporting.amesim_results import load_test_mql_amesim_results
+from PythonModels.reporting.test_mql_output_schema import build_test_mql_output_schema
+from PythonModels.reporting.test_mql_output_validation import compare_validated_test_mql_output
 from PythonModels.systems.test_mql import TestMqlSystem
 
 
+REPO_ROOT = Path(__file__).resolve().parents[1]
+TEST_MQL_AME = REPO_ROOT / "AmesimModels" / "test_mql.ame"
+
+
 class TestMqlPnl0001ChamberSegmentTests(unittest.TestCase):
     def setUp(self) -> None:
         self.system = TestMqlSystem()
@@ -247,6 +255,11 @@ class TestMqlPn3NodeChamberSegmentTests(unittest.TestCase):
 
 
 class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls) -> None:
+        cls.amesim_results = load_test_mql_amesim_results(TEST_MQL_AME)
+        cls.output_schema = build_test_mql_output_schema(cls.amesim_results)
+
     def setUp(self) -> None:
         self.system = TestMqlSystem()
         self.spec = self.system.discover_pneumatic_branch_topology().chamber_segment_specs[0]
@@ -671,6 +684,50 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
         self.assertAlmostEqual(rhs[128], expected_contact_force / support_mass)
         self.assertAlmostEqual(rhs[129], 0.0)
 
+    def test_simulates_full_state_key_data_paths_for_amesim_comparison(self) -> None:
+        data_paths = (
+            "press@pn_c1_8",
+            "vol@pn_c1_8",
+            "vol1@pn_brp2_8",
+            "vvol1@pn_brp2_8",
+            "x1@mass_friction_endstops_10",
+            "v1@mass_friction_endstops_10",
+            "acc1@mass_friction_endstops_10",
+            "x1@mass_friction_endstops_18",
+            "v1@mass_friction_endstops_18",
+            "acc1@mass_friction_endstops_18",
+        )
+
+        result = self.system.simulate_full_state_series_from_spec(
+            self.spec,
+            inlet_node_pressure_pa=15.31e6,
+            resistance_boundary_pressure_pa=15.29e6,
+            config=SolveIVPConfig(t_start=0.0, t_stop=0.0),
+            data_paths=data_paths,
+        )
+        comparison = compare_validated_test_mql_output(
+            times=result.t,
+            series_by_data_path={data_path: result.series[data_path] for data_path in data_paths},
+            schema=self.output_schema,
+            amesim_results=self.amesim_results,
+            data_paths=data_paths,
+        )
+
+        self.assertEqual(result.series["time"], [0.0])
+        self.assertEqual(len(result.y), 132)
+        self.assertEqual(len(comparison.metrics), len(data_paths))
+        self.assertAlmostEqual(result.series["press@pn_c1_8"][0], -1300.0)
+        self.assertAlmostEqual(result.series["vol@pn_c1_8"][0], 15000.0)
+        self.assertAlmostEqual(result.series["vol1@pn_brp2_8"][0], 0.0)
+        self.assertAlmostEqual(result.series["vvol1@pn_brp2_8"][0], 0.0)
+        self.assertAlmostEqual(result.series["x1@mass_friction_endstops_10"][0], 0.0)
+        self.assertAlmostEqual(result.series["v1@mass_friction_endstops_10"][0], 0.0)
+        self.assertAlmostEqual(
+            result.series["acc1@mass_friction_endstops_10"][0],
+            -0.8167936695810917,
+        )
+        self.assertAlmostEqual(result.series["acc1@mass_friction_endstops_18"][0], 0.0)
+
     def test_simulates_full_132_state_closure(self) -> None:
         solution = self.system.simulate_full_state_from_spec(
             self.spec,