PythonModels/scripts/run_test_mql_full_state_comparison.py

@@ -2,7 +2,7 @@ from __future__ import annotations
 
 import argparse
 from bisect import bisect_left
-from dataclasses import dataclass, field
+from dataclasses import dataclass, field, replace
 from datetime import UTC, datetime
 from math import nextafter, sqrt
 from pathlib import Path
@@ -731,6 +731,30 @@ class TestMqlPnvoEventWindowDiagnostic:
         )
 
 
+@dataclass(frozen=True)
+class TestMqlPnvoEventWindowCandidateMetricComparison:
+    time_s: float
+    metric_name: str
+    default_abs_error: float
+    candidate_abs_error: float
+
+    @property
+    def abs_error_delta(self) -> float:
+        return self.candidate_abs_error - self.default_abs_error
+
+
+@dataclass(frozen=True)
+class TestMqlPnvoEventWindowCandidateComparisonDiagnostic:
+    orifice_alias: str
+    event_time_s: float
+    final_time_s: float
+    default_diagnostic: TestMqlPnvoEventWindowDiagnostic
+    candidate_diagnostic: TestMqlPnvoEventWindowDiagnostic
+    metric_comparisons: tuple[
+        TestMqlPnvoEventWindowCandidateMetricComparison, ...
+    ]
+
+
 @dataclass(frozen=True)
 class TestMqlFullStateComparisonPathConfig:
     archive_path: Path = field(
@@ -3245,12 +3269,14 @@ def run_test_mql_pnvo_event_window_diagnostic(
     config: TestMqlFullStateComparisonScriptConfig | None = None,
     *,
     orifice_alias: str = "pn_morifice_1",
+    final_time_s: float = 0.05,
 ) -> TestMqlPnvoEventWindowDiagnostic:
     config = config or TestMqlFullStateComparisonScriptConfig()
     system = TestMqlSystem(archive_path=config.paths.archive_path)
     control = system.pneumatic_assembly.variable_orifice_controls[orifice_alias]
     event_time_s = control.step.step_time_s
-    final_time_s = 0.05
+    if final_time_s <= event_time_s:
+        raise ValueError("final_time_s must be greater than the PNVO event time")
     spec = system.discover_pneumatic_branch_topology().chamber_segment_specs[0]
     closure_kwargs = {
         "inlet_node_pressure_pa": config.execution.inlet_node_pressure_pa,
@@ -3268,7 +3294,7 @@ def run_test_mql_pnvo_event_window_diagnostic(
     closure = system.full_state_closure_from_spec(spec, **closure_kwargs)
     state_vector = closure.initial_state_vector()
     solver_template = config.execution.solver
-    segments = (
+    segment_templates = (
         SolveIVPConfig(
             t_start=solver_template.t_start,
             t_stop=nextafter(event_time_s, 0.0),
@@ -3320,13 +3346,18 @@ def run_test_mql_pnvo_event_window_diagnostic(
         ),
         SolveIVPConfig(
             t_start=0.048,
-            t_stop=final_time_s,
+            t_stop=0.05,
             method="BDF",
             rtol=1.0e-5,
             atol=1.0e-8,
             max_step=1.0e-5,
         ),
     )
+    segments = tuple(
+        replace(segment, t_stop=min(segment.t_stop, final_time_s))
+        for segment in segment_templates
+        if segment.t_start < final_time_s
+    )
     data_paths = (
         "press@pn_c1_8",
         "temp@pn_c1_8",
@@ -3342,7 +3373,11 @@ def run_test_mql_pnvo_event_window_diagnostic(
         f"xv@{orifice_alias}",
         f"dm2@{orifice_alias}",
     )
-    sample_times = (
+    sample_times = tuple(
+        sorted(
+            {
+                sample_time
+                for sample_time in (
                     nextafter(event_time_s, 0.0),
                     0.0401,
                     0.0402,
@@ -3352,6 +3387,10 @@ def run_test_mql_pnvo_event_window_diagnostic(
                     0.048,
                     final_time_s,
                 )
+                if sample_time <= final_time_s
+            }
+        )
+    )
     state_vector_by_sample_time: dict[float, list[float]] = {}
     segment_diagnostics: list[TestMqlPnvoEventWindowSegmentDiagnostic] = []
     for segment in segments:
@@ -3573,6 +3612,214 @@ def run_test_mql_pnvo_event_window_diagnostic(
     )
 
 
+
+def _event_window_config_with_pneumatic_96_reference_rhs(
+    config: TestMqlFullStateComparisonScriptConfig | None,
+    *,
+    enabled: bool,
+) -> TestMqlFullStateComparisonScriptConfig:
+    base = config or TestMqlFullStateComparisonScriptConfig()
+    return replace(
+        base,
+        execution=replace(
+            base.execution,
+            use_pneumatic_96_reference_rhs=enabled,
+        ),
+    )
+
+
+def _pnvo_event_window_sample_by_time(
+    diagnostic: TestMqlPnvoEventWindowDiagnostic,
+) -> dict[float, TestMqlPnvoEventWindowSampleDiagnostic]:
+    return {sample.time_s: sample for sample in diagnostic.sample_diagnostics}
+
+
+def _append_metric_comparison(
+    comparisons: list[TestMqlPnvoEventWindowCandidateMetricComparison],
+    *,
+    time_s: float,
+    metric_name: str,
+    default_abs_error: float,
+    candidate_abs_error: float,
+) -> None:
+    comparisons.append(
+        TestMqlPnvoEventWindowCandidateMetricComparison(
+            time_s=time_s,
+            metric_name=metric_name,
+            default_abs_error=default_abs_error,
+            candidate_abs_error=candidate_abs_error,
+        )
+    )
+
+
+def _append_pneumatic_96_candidate_metric_comparisons(
+    comparisons: list[TestMqlPnvoEventWindowCandidateMetricComparison],
+    *,
+    time_s: float,
+    default_sample: TestMqlPnvoEventWindowSampleDiagnostic,
+    candidate_sample: TestMqlPnvoEventWindowSampleDiagnostic,
+) -> None:
+    if (
+        not default_sample.pneumatic_96_inlet_line_diagnostics
+        or not candidate_sample.pneumatic_96_inlet_line_diagnostics
+    ):
+        return
+    default_line = default_sample.pneumatic_96_inlet_line_diagnostics[0]
+    candidate_line = candidate_sample.pneumatic_96_inlet_line_diagnostics[0]
+    for metric_name, default_value, candidate_value in (
+        (
+            "p96:p2_abs_error_pa",
+            default_line.python_line_pressure_error_pa,
+            candidate_line.python_line_pressure_error_pa,
+        ),
+        (
+            "p96:t2_abs_error_k",
+            default_line.python_line_temperature_error_k,
+            candidate_line.python_line_temperature_error_k,
+        ),
+        (
+            "p96:storage_dm_abs_error_g_s",
+            default_line.python_storage_mass_derivative_error_g_s,
+            candidate_line.python_storage_mass_derivative_error_g_s,
+        ),
+        (
+            "p96:active_dT_abs_error_k_s",
+            default_line.python_active_temperature_derivative_error_k_s,
+            candidate_line.python_active_temperature_derivative_error_k_s,
+        ),
+        (
+            "p96:active_dP_abs_error_pa_s",
+            default_line.python_active_pressure_derivative_eos_error_pa_s,
+            candidate_line.python_active_pressure_derivative_eos_error_pa_s,
+        ),
+    ):
+        _append_metric_comparison(
+            comparisons,
+            time_s=time_s,
+            metric_name=metric_name,
+            default_abs_error=abs(default_value),
+            candidate_abs_error=abs(candidate_value),
+        )
+
+
+def _append_pneumatic_87_candidate_metric_comparisons(
+    comparisons: list[TestMqlPnvoEventWindowCandidateMetricComparison],
+    *,
+    time_s: float,
+    default_sample: TestMqlPnvoEventWindowSampleDiagnostic,
+    candidate_sample: TestMqlPnvoEventWindowSampleDiagnostic,
+) -> None:
+    if (
+        not default_sample.pnl0003_energy_diagnostics
+        or not candidate_sample.pnl0003_energy_diagnostics
+    ):
+        return
+    default_line = default_sample.pnl0003_energy_diagnostics[0]
+    candidate_line = candidate_sample.pnl0003_energy_diagnostics[0]
+    for metric_name, default_value, candidate_value in (
+        (
+            "p87:node_to_line_dm_abs_error_g_s",
+            default_line.python_pn3_node_to_line_mass_flow_error_g_s,
+            candidate_line.python_pn3_node_to_line_mass_flow_error_g_s,
+        ),
+        (
+            "p87:center_dm_abs_error_g_s",
+            default_line.amesim_center_flow_python_sign_error_g_s,
+            candidate_line.amesim_center_flow_python_sign_error_g_s,
+        ),
+        (
+            "p87:sdm1_abs_error_g_s",
+            default_line.python_port_1_storage_mass_derivative_error_g_s,
+            candidate_line.python_port_1_storage_mass_derivative_error_g_s,
+        ),
+        (
+            "p87:sdm2_abs_error_g_s",
+            default_line.python_port_2_storage_mass_derivative_error_g_s,
+            candidate_line.python_port_2_storage_mass_derivative_error_g_s,
+        ),
+        (
+            "p87:t2_abs_error_k",
+            default_line.python_port_2_temperature_error_to_amesim_k,
+            candidate_line.python_port_2_temperature_error_to_amesim_k,
+        ),
+    ):
+        _append_metric_comparison(
+            comparisons,
+            time_s=time_s,
+            metric_name=metric_name,
+            default_abs_error=abs(default_value),
+            candidate_abs_error=abs(candidate_value),
+        )
+
+
+def _pnvo_event_window_candidate_metric_comparisons(
+    default_diagnostic: TestMqlPnvoEventWindowDiagnostic,
+    candidate_diagnostic: TestMqlPnvoEventWindowDiagnostic,
+) -> tuple[TestMqlPnvoEventWindowCandidateMetricComparison, ...]:
+    candidate_samples = _pnvo_event_window_sample_by_time(candidate_diagnostic)
+    comparisons: list[TestMqlPnvoEventWindowCandidateMetricComparison] = []
+    for default_sample in default_diagnostic.sample_diagnostics:
+        candidate_sample = candidate_samples.get(default_sample.time_s)
+        if candidate_sample is None:
+            continue
+        for data_path in default_sample.data_paths:
+            _append_metric_comparison(
+                comparisons,
+                time_s=default_sample.time_s,
+                metric_name=f"path:{data_path}",
+                default_abs_error=default_sample.abs_error(data_path),
+                candidate_abs_error=candidate_sample.abs_error(data_path),
+            )
+        _append_pneumatic_96_candidate_metric_comparisons(
+            comparisons,
+            time_s=default_sample.time_s,
+            default_sample=default_sample,
+            candidate_sample=candidate_sample,
+        )
+        _append_pneumatic_87_candidate_metric_comparisons(
+            comparisons,
+            time_s=default_sample.time_s,
+            default_sample=default_sample,
+            candidate_sample=candidate_sample,
+        )
+    return tuple(comparisons)
+
+
+def run_test_mql_pnvo_event_window_candidate_comparison(
+    config: TestMqlFullStateComparisonScriptConfig | None = None,
+    *,
+    orifice_alias: str = "pn_morifice_1",
+    final_time_s: float = 0.041,
+) -> TestMqlPnvoEventWindowCandidateComparisonDiagnostic:
+    default_diagnostic = run_test_mql_pnvo_event_window_diagnostic(
+        _event_window_config_with_pneumatic_96_reference_rhs(
+            config,
+            enabled=False,
+        ),
+        orifice_alias=orifice_alias,
+        final_time_s=final_time_s,
+    )
+    candidate_diagnostic = run_test_mql_pnvo_event_window_diagnostic(
+        _event_window_config_with_pneumatic_96_reference_rhs(
+            config,
+            enabled=True,
+        ),
+        orifice_alias=orifice_alias,
+        final_time_s=final_time_s,
+    )
+    return TestMqlPnvoEventWindowCandidateComparisonDiagnostic(
+        orifice_alias=orifice_alias,
+        event_time_s=default_diagnostic.event_time_s,
+        final_time_s=default_diagnostic.final_time_s,
+        default_diagnostic=default_diagnostic,
+        candidate_diagnostic=candidate_diagnostic,
+        metric_comparisons=_pnvo_event_window_candidate_metric_comparisons(
+            default_diagnostic,
+            candidate_diagnostic,
+        ),
+    )
+
+
 def format_test_mql_pnvo_event_window_summary(
     diagnostic: TestMqlPnvoEventWindowDiagnostic,
 ) -> str:
@@ -4094,6 +4341,38 @@ def format_test_mql_pnvo_event_window_summary(
     return "\n".join(lines) + "\n"
 
 
+
+def format_test_mql_pnvo_event_window_candidate_comparison_summary(
+    diagnostic: TestMqlPnvoEventWindowCandidateComparisonDiagnostic,
+) -> str:
+    lines = [
+        "Model: test_mql",
+        f"Mode: PNVO event window candidate comparison ({diagnostic.orifice_alias})",
+        f"Event time: {diagnostic.event_time_s}",
+        f"Final time: {diagnostic.final_time_s}",
+        "Default RHS: current pneumatic_96 internal-energy RHS",
+        "Candidate RHS: pneumatic_96 reference-enthalpy/mass-offset RHS",
+        "Metric comparisons: negative delta means the candidate reduced error",
+    ]
+    for comparison in diagnostic.metric_comparisons:
+        lines.append(
+            f"  t={comparison.time_s} - {comparison.metric_name}: "
+            f"default_abs_error={comparison.default_abs_error}, "
+            f"candidate_abs_error={comparison.candidate_abs_error}, "
+            f"delta={comparison.abs_error_delta}"
+        )
+    lines.append("Final comparison:")
+    for data_path in diagnostic.default_diagnostic.data_paths:
+        lines.append(
+            f"  - {data_path}: "
+            f"default_abs_error={diagnostic.default_diagnostic.abs_error(data_path)}, "
+            f"candidate_abs_error={diagnostic.candidate_diagnostic.abs_error(data_path)}, "
+            f"delta="
+            f"{diagnostic.candidate_diagnostic.abs_error(data_path) - diagnostic.default_diagnostic.abs_error(data_path)}"
+        )
+    return "\n".join(lines) + "\n"
+
+
 def format_test_mql_full_state_comparison_summary(
     run: TestMqlFullStateComparisonRun,
 ) -> str:
@@ -4251,6 +4530,16 @@ def main() -> None:
         action="store_true",
         help="run the segmented PNVO opening window through the t=0.05 s save point",
     )
+    parser.add_argument(
+        "--pnvo-event-window-candidate-comparison",
+        action="store_true",
+        help="compare the PNVO event window with and without the pneumatic_96 RHS candidate",
+    )
+    parser.add_argument(
+        "--pnvo-event-window-final-time",
+        type=float,
+        help="override PNVO event window final time; candidate comparison defaults to 0.041 s",
+    )
     parser.add_argument(
         "--amesim-results-archive",
         type=Path,
@@ -4276,8 +4565,26 @@ def main() -> None:
         diagnostic = run_test_mql_pnvo_event_boundary_diagnostic(config)
         print(format_test_mql_pnvo_event_boundary_summary(diagnostic), end="")
         return
+    if args.pnvo_event_window_candidate_comparison:
+        diagnostic = run_test_mql_pnvo_event_window_candidate_comparison(
+            config,
+            final_time_s=args.pnvo_event_window_final_time or 0.041,
+        )
+        print(
+            format_test_mql_pnvo_event_window_candidate_comparison_summary(
+                diagnostic
+            ),
+            end="",
+        )
+        return
     if args.pnvo_event_window:
+        if args.pnvo_event_window_final_time is None:
             diagnostic = run_test_mql_pnvo_event_window_diagnostic(config)
+        else:
+            diagnostic = run_test_mql_pnvo_event_window_diagnostic(
+                config,
+                final_time_s=args.pnvo_event_window_final_time,
+            )
         print(format_test_mql_pnvo_event_window_summary(diagnostic), end="")
         return
 

tests/test_run_test_mql_full_state_comparison.py

@@ -19,12 +19,15 @@ from PythonModels.scripts.run_test_mql_full_state_comparison import (
     TestMqlPnvoUpstreamEnthalpyDiagnostic,
     TestMqlPnl0001PressureLossCalibrationDiagnostic,
     TestMqlPnvoEventBoundaryDiagnostic,
+    TestMqlPnvoEventWindowCandidateComparisonDiagnostic,
+    TestMqlPnvoEventWindowCandidateMetricComparison,
     TestMqlPnvoEventWindowDiagnostic,
     TestMqlPnvoEventWindowSampleDiagnostic,
     TestMqlPnvoFlowParameterDiagnostic,
     TestMqlPnvoEventWindowSegmentDiagnostic,
     format_test_mql_full_state_comparison_summary,
     format_test_mql_pnvo_event_boundary_summary,
+    format_test_mql_pnvo_event_window_candidate_comparison_summary,
     format_test_mql_pnvo_event_window_summary,
     run_test_mql_full_state_comparison,
     _ideal_pn2vol_reference_dtemp,
@@ -899,6 +902,73 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase):
         self.assertIn("Final comparison:", summary)
         self.assertTrue(summary.endswith("\n"))
 
+    def test_formats_pnvo_event_window_candidate_comparison(self) -> None:
+        data_paths = ("press@pn_c1_8", "dm2@pn_morifice_1")
+        default_diagnostic = TestMqlPnvoEventWindowDiagnostic(
+            orifice_alias="pn_morifice_1",
+            event_time_s=0.04,
+            final_time_s=0.05,
+            data_paths=data_paths,
+            segment_diagnostics=(),
+            sample_diagnostics=(),
+            python_values_by_data_path={
+                "press@pn_c1_8": 1.0,
+                "dm2@pn_morifice_1": 10.0,
+            },
+            amesim_values_by_data_path={
+                "press@pn_c1_8": 4.0,
+                "dm2@pn_morifice_1": 8.0,
+            },
+        )
+        candidate_diagnostic = TestMqlPnvoEventWindowDiagnostic(
+            orifice_alias="pn_morifice_1",
+            event_time_s=0.04,
+            final_time_s=0.05,
+            data_paths=data_paths,
+            segment_diagnostics=(),
+            sample_diagnostics=(),
+            python_values_by_data_path={
+                "press@pn_c1_8": 3.0,
+                "dm2@pn_morifice_1": 13.0,
+            },
+            amesim_values_by_data_path={
+                "press@pn_c1_8": 4.0,
+                "dm2@pn_morifice_1": 8.0,
+            },
+        )
+        metric = TestMqlPnvoEventWindowCandidateMetricComparison(
+            time_s=0.0401,
+            metric_name="p96:p2_abs_error_pa",
+            default_abs_error=32000.0,
+            candidate_abs_error=3600.0,
+        )
+        diagnostic = TestMqlPnvoEventWindowCandidateComparisonDiagnostic(
+            orifice_alias="pn_morifice_1",
+            event_time_s=0.04,
+            final_time_s=0.05,
+            default_diagnostic=default_diagnostic,
+            candidate_diagnostic=candidate_diagnostic,
+            metric_comparisons=(metric,),
+        )
+
+        summary = format_test_mql_pnvo_event_window_candidate_comparison_summary(
+            diagnostic
+        )
+
+        self.assertEqual(metric.abs_error_delta, -28400.0)
+        self.assertIn("PNVO event window candidate comparison", summary)
+        self.assertIn("negative delta means the candidate reduced error", summary)
+        self.assertIn("t=0.0401 - p96:p2_abs_error_pa", summary)
+        self.assertIn("default_abs_error=32000.0", summary)
+        self.assertIn("candidate_abs_error=3600.0", summary)
+        self.assertIn("delta=-28400.0", summary)
+        self.assertIn("Final comparison:", summary)
+        self.assertIn("press@pn_c1_8", summary)
+        self.assertIn("delta=-2.0", summary)
+        self.assertIn("dm2@pn_morifice_1", summary)
+        self.assertIn("delta=3.0", summary)
+        self.assertTrue(summary.endswith("\n"))
+
     def test_full_state_script_can_skip_artifact_files(self) -> None:
         with tempfile.TemporaryDirectory() as tmpdir:
             output_dir = Path(tmpdir)