diff --git a/PythonModels/scripts/run_test_mql_full_state_comparison.py b/PythonModels/scripts/run_test_mql_full_state_comparison.py index 13ededc..dd68ee9 100644 --- a/PythonModels/scripts/run_test_mql_full_state_comparison.py +++ b/PythonModels/scripts/run_test_mql_full_state_comparison.py @@ -53,6 +53,7 @@ class TestMqlFullStateSignalDiagnostic: @dataclass(frozen=True) class TestMqlFullStateComparisonRun: system: TestMqlSystem + closure: object amesim_results: AmesimResults output_schema: TestMqlOutputSchema result: TestMqlSimulationResult @@ -128,6 +129,13 @@ class TestMqlFullStateComparisonRun: diagnostics = self.diagnostics_by_final_abs_error() return diagnostics[0] if diagnostics else None + def chamber_rhs_diagnostic(self, chamber_alias: str, sample_index: int = -1): + state_vector = [row[sample_index] for row in self.result.y] + return self.closure.variable_chamber_rhs_diagnostic( + chamber_alias=chamber_alias, + state_vector=state_vector, + ) + @dataclass(frozen=True) class TestMqlFullStateComparisonPathConfig: @@ -193,6 +201,17 @@ def run_test_mql_full_state_comparison( t_eval=list(config.execution.t_eval) if config.execution.t_eval is not None else None, data_paths=selected_paths, ) + closure = system.full_state_closure_from_spec( + spec, + inlet_node_pressure_pa=config.execution.inlet_node_pressure_pa, + resistance_boundary_pressure_pa=( + config.execution.resistance_boundary_pressure_pa + ), + inlet_node_temperature_k=config.execution.inlet_node_temperature_k, + resistance_boundary_temperature_k=( + config.execution.resistance_boundary_temperature_k + ), + ) series_by_data_path = { data_path: result.series[data_path] for data_path in result.series @@ -213,6 +232,7 @@ def run_test_mql_full_state_comparison( ) run = TestMqlFullStateComparisonRun( system=system, + closure=closure, amesim_results=amesim_results, output_schema=output_schema, result=result, @@ -278,9 +298,38 @@ def format_test_mql_full_state_comparison_summary( f"final_amesim={diagnostic.final_amesim_value}, " f"final_abs_error={diagnostic.final_abs_error}" ) + chamber_diagnostic = _largest_chamber_rhs_diagnostic(run) + if chamber_diagnostic is not None: + lines.append( + "Largest endpoint chamber RHS breakdown: " + f"{chamber_diagnostic.chamber_alias}" + ) + lines.extend( + [ + f" - piston_alias={chamber_diagnostic.piston_alias}", + f" - pressure_pa={chamber_diagnostic.chamber_pressure_pa}", + f" - volume_m3={chamber_diagnostic.chamber_volume_m3}", + f" - volume_rate_m3_s={chamber_diagnostic.chamber_volume_rate_m3_s}", + f" - mass_derivative_kg_s={chamber_diagnostic.mass_derivative_kg_s}", + f" - port_a_energy_flow_w={chamber_diagnostic.port_a_energy_flow_w}", + f" - boundary_work_w={chamber_diagnostic.boundary_work_w}", + f" - energy_derivative_w={chamber_diagnostic.energy_derivative_w}", + ] + ) return "\n".join(lines) + "\n" +def _largest_chamber_rhs_diagnostic(run: TestMqlFullStateComparisonRun): + diagnostic = run.largest_final_abs_error_diagnostic + if diagnostic is None or "@" not in diagnostic.data_path: + return None + _signal, alias = diagnostic.data_path.split("@", 1) + try: + return run.chamber_rhs_diagnostic(alias) + except KeyError: + return None + + def main() -> None: run, output_dir = run_test_mql_full_state_comparison() print(format_test_mql_full_state_comparison_summary(run), end="") diff --git a/PythonModels/systems/test_mql.py b/PythonModels/systems/test_mql.py index 97c7c6f..f88639b 100644 --- a/PythonModels/systems/test_mql.py +++ b/PythonModels/systems/test_mql.py @@ -3856,6 +3856,23 @@ class TestMqlFullStateSnapshot: return self.mechanical.state_count +@dataclass(frozen=True) +class TestMqlVariableChamberRhsDiagnostic: + chamber_alias: str + piston_alias: str + chamber_pressure_pa: float + chamber_temperature_k: float + chamber_volume_m3: float + chamber_volume_rate_m3_s: float + port_a_mass_flow_kg_s: float + port_b_mass_flow_kg_s: float + mass_derivative_kg_s: float + port_a_energy_flow_w: float + port_b_energy_flow_w: float + boundary_work_w: float + energy_derivative_w: float + + _PISTON_FORCE_BINDINGS = ( ("mass_friction_endstops_10", "pn_brp2_8", "p4_port3_remote_primary_chamber"), ("mass_friction_endstops_11", "pn_brp2_9", "p4_primary_chamber"), @@ -3906,6 +3923,50 @@ _PNCH012_ALIAS_BY_SNAPSHOT_FIELD = { } +_PRIMARY_CHAMBER_DIAGNOSTIC_FIELDS = { + "pn_c1_8": ( + "p4_port3_remote_primary_chamber", + "p4_port3_remote_primary_line", + "p4_port3_remote_chamber_to_line_flow", + ), + "pn_c1_9": ( + "p4_primary_chamber", + "p4_primary_line", + "p4_primary_chamber_to_line_flow", + ), + "pn_c1_10": ( + "p4_port1_remote_primary_chamber", + "p4_port1_remote_primary_line", + "p4_port1_remote_chamber_to_line_flow", + ), + "pn_c1_11": ( + "p4_port1_far_primary_chamber", + "p4_port1_far_primary_line", + "p4_port1_far_chamber_to_line_flow", + ), + "pn_c1_12": ( + "p4_port1_next_primary_chamber", + "p4_port1_next_primary_line", + "p4_port1_next_chamber_to_line_flow", + ), + "pn_c1_13": ( + "p4_bridge_primary_chamber", + "p4_bridge_primary_line", + "p4_bridge_chamber_to_line_flow", + ), + "pn_c1_14": ( + "p4_port3_next_primary_chamber", + "p4_port3_next_primary_line", + "p4_port3_next_chamber_to_line_flow", + ), + "pn_c1_15": ( + "p4_port3_far_primary_chamber", + "p4_port3_far_primary_line", + "p4_port3_far_chamber_to_line_flow", + ), +} + + class TestMqlFullStateClosure: def __init__(self, *, pneumatic_closure: object, mechanical_closure: object) -> None: self.pneumatic_closure = pneumatic_closure @@ -3988,6 +4049,60 @@ class TestMqlFullStateClosure: ), ) + def variable_chamber_rhs_diagnostic( + self, + *, + chamber_alias: str, + state_vector: list[float], + ) -> TestMqlVariableChamberRhsDiagnostic: + if chamber_alias not in _PRIMARY_CHAMBER_DIAGNOSTIC_FIELDS: + raise KeyError(chamber_alias) + chamber_field, line_field, flow_field = _PRIMARY_CHAMBER_DIAGNOSTIC_FIELDS[ + chamber_alias + ] + piston_alias = next( + piston_alias + for _mass_alias, piston_alias, bound_chamber_field in _PISTON_FORCE_BINDINGS + if bound_chamber_field == chamber_field + ) + snapshot = self.snapshot(state_vector) + chamber = getattr(self.pneumatic_closure.components, chamber_field) + chamber_properties = getattr(snapshot.pneumatic, chamber_field) + line_properties = getattr(snapshot.pneumatic, line_field) + chamber_to_line_flow = getattr(snapshot.pneumatic, flow_field) + port_a_mass_flow = -chamber_to_line_flow + port_b_mass_flow = 0.0 + port_a_inlet_h = chamber.connection_inlet_enthalpy( + port_m_flow=port_a_mass_flow, + connected_h=line_properties.h, + internal_h=chamber_properties.h, + ) + port_b_inlet_h = chamber.connection_inlet_enthalpy( + port_m_flow=port_b_mass_flow, + connected_h=chamber_properties.h, + internal_h=chamber_properties.h, + ) + port_a_energy_flow = port_a_mass_flow * port_a_inlet_h + port_b_energy_flow = port_b_mass_flow * port_b_inlet_h + boundary_work = -chamber_properties.p * chamber.volume_rate_m3_s() + return TestMqlVariableChamberRhsDiagnostic( + chamber_alias=chamber_alias, + piston_alias=piston_alias, + chamber_pressure_pa=chamber_properties.p, + chamber_temperature_k=chamber_properties.T, + chamber_volume_m3=chamber.volume, + chamber_volume_rate_m3_s=chamber.volume_rate_m3_s(), + port_a_mass_flow_kg_s=port_a_mass_flow, + port_b_mass_flow_kg_s=port_b_mass_flow, + mass_derivative_kg_s=port_a_mass_flow + port_b_mass_flow, + port_a_energy_flow_w=port_a_energy_flow, + port_b_energy_flow_w=port_b_energy_flow, + boundary_work_w=boundary_work, + energy_derivative_w=( + port_a_energy_flow + port_b_energy_flow + boundary_work + ), + ) + def data_path_values( self, *, diff --git a/tests/test_run_test_mql_full_state_comparison.py b/tests/test_run_test_mql_full_state_comparison.py index 55981f3..1d975ff 100644 --- a/tests/test_run_test_mql_full_state_comparison.py +++ b/tests/test_run_test_mql_full_state_comparison.py @@ -59,6 +59,16 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): run.diagnostics_by_final_abs_error()[0].data_path, "press@pn_c1_8", ) + chamber_diagnostic = run.chamber_rhs_diagnostic("pn_c1_8") + self.assertEqual(chamber_diagnostic.chamber_alias, "pn_c1_8") + self.assertEqual(chamber_diagnostic.piston_alias, "pn_brp2_8") + self.assertAlmostEqual(chamber_diagnostic.chamber_volume_m3, 0.015) + self.assertAlmostEqual( + chamber_diagnostic.energy_derivative_w, + chamber_diagnostic.port_a_energy_flow_w + + chamber_diagnostic.port_b_energy_flow_w + + chamber_diagnostic.boundary_work_w, + ) self.assertTrue(csv_summary_exists) self.assertIn("Mode: Python 132 full-state closure comparison", summary) self.assertIn("Compared signals: 2", summary) @@ -66,6 +76,10 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): self.assertIn("Largest final endpoint error: press@pn_c1_8=", summary) self.assertIn("Metrics by max absolute error:", summary) self.assertIn("Endpoint diagnostics by final absolute error:", summary) + self.assertIn("Largest endpoint chamber RHS breakdown: pn_c1_8", summary) + self.assertIn("piston_alias=pn_brp2_8", summary) + self.assertIn("mass_derivative_kg_s=", summary) + self.assertIn("boundary_work_w=", summary) self.assertIn("final_python=", summary) self.assertIn("final_amesim=", summary) self.assertIn("python.press@pn_c1_8", csv_header) @@ -118,6 +132,7 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): self.assertIn("Largest final endpoint error: press@pn_c1_8=", summary) self.assertIn("Metrics by max absolute error:", summary) self.assertIn("Endpoint diagnostics by final absolute error:", summary) + self.assertIn("Largest endpoint chamber RHS breakdown: pn_c1_8", summary) self.assertTrue(summary.endswith("\n")) diff --git a/tests/test_test_mql_pnl0001_segment.py b/tests/test_test_mql_pnl0001_segment.py index 1dded6e..af34639 100644 --- a/tests/test_test_mql_pnl0001_segment.py +++ b/tests/test_test_mql_pnl0001_segment.py @@ -684,6 +684,32 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase): self.assertAlmostEqual(rhs[128], expected_contact_force / support_mass) self.assertAlmostEqual(rhs[129], 0.0) + def test_full_state_closure_diagnoses_variable_chamber_rhs_terms(self) -> None: + closure = self.system.full_state_closure_from_spec( + self.spec, + inlet_node_pressure_pa=15.31e6, + resistance_boundary_pressure_pa=15.29e6, + ) + state = closure.initial_state_vector() + + diagnostic = closure.variable_chamber_rhs_diagnostic( + chamber_alias="pn_c1_8", + state_vector=state, + ) + rhs = closure.rhs(state) + + self.assertEqual(diagnostic.piston_alias, "pn_brp2_8") + self.assertAlmostEqual(diagnostic.chamber_volume_m3, 0.015) + self.assertAlmostEqual(diagnostic.chamber_volume_rate_m3_s, 0.0) + self.assertAlmostEqual(diagnostic.mass_derivative_kg_s, rhs[46]) + self.assertAlmostEqual(diagnostic.energy_derivative_w, rhs[47]) + self.assertAlmostEqual( + diagnostic.energy_derivative_w, + diagnostic.port_a_energy_flow_w + + diagnostic.port_b_energy_flow_w + + diagnostic.boundary_work_w, + ) + def test_simulates_full_state_key_data_paths_for_amesim_comparison(self) -> None: data_paths = ( "press@pn_c1_8",