PythonModels/reporting/amesim_results.py

@@ -47,21 +47,32 @@ _DATA_PATH_RE = re.compile(r"Data_Path=(\S+)")
 _PARAM_ID_RE = re.compile(r"Param_Id=(\d+)")
 
 
-def load_test_mql_amesim_results(archive_path: str | Path) -> AmesimResults:
+def load_test_mql_amesim_results(
+    archive_path: str | Path,
+    *,
+    time_stop_s: float | None = None,
+) -> AmesimResults:
     return load_amesim_results_from_archive(
         archive_path=archive_path,
-        var_member="test_mql_.var",
+        var_member=None,
-        results_member="test_mql_.results",
+        results_member=None,
+        time_stop_s=time_stop_s,
     )
 
 
 def load_amesim_results_from_archive(
     *,
     archive_path: str | Path,
-    var_member: str,
+    var_member: str | None,
-    results_member: str,
+    results_member: str | None,
+    time_stop_s: float | None = None,
 ) -> AmesimResults:
     with tarfile.open(archive_path) as archive:
+        var_member, results_member = _resolve_result_members(
+            archive,
+            var_member=var_member,
+            results_member=results_member,
+        )
         var_file = archive.extractfile(var_member)
         results_file = archive.extractfile(results_member)
         if var_file is None:
@@ -69,11 +80,113 @@ def load_amesim_results_from_archive(
         if results_file is None:
             raise AmesimResultsError(f"Missing AMESim results member: {results_member}")
         var_lines = var_file.read().decode("latin1").splitlines()
+        variables = tuple(
+            _parse_variable_line(index, line) for index, line in enumerate(var_lines)
+        )
+        if time_stop_s is not None:
+            return _parse_amesim_results_window(
+                results_file,
+                variables,
+                time_stop_s=time_stop_s,
+            )
         results_data = results_file.read()
-    variables = tuple(_parse_variable_line(index, line) for index, line in enumerate(var_lines))
     return parse_amesim_results_bytes(results_data, variables)
 
 
+def _resolve_result_members(
+    archive: tarfile.TarFile,
+    *,
+    var_member: str | None,
+    results_member: str | None,
+) -> tuple[str, str]:
+    member_names = set(archive.getnames())
+    if var_member is not None or results_member is not None:
+        if var_member is None or results_member is None:
+            raise AmesimResultsError(
+                "var_member and results_member must either both be set or both be omitted."
+            )
+        return var_member, results_member
+
+    preferred = ("test_mql_.var", "test_mql_.results")
+    if preferred[0] in member_names and preferred[1] in member_names:
+        return preferred
+
+    pairs = sorted(
+        (name, f"{name[:-4]}.results")
+        for name in member_names
+        if name.endswith(".var") and f"{name[:-4]}.results" in member_names
+    )
+    if len(pairs) != 1:
+        raise AmesimResultsError(
+            "Unable to identify a unique AMESim .var/.results member pair."
+        )
+    return pairs[0]
+
+
+def _parse_amesim_results_window(
+    results_file,
+    variables: tuple[AmesimVariable, ...],
+    *,
+    time_stop_s: float,
+) -> AmesimResults:
+    header = results_file.read(8)
+    if len(header) < 8:
+        raise AmesimResultsError("AMESim results data is too small.")
+    point_count, encoded_saved_variable_count = struct.unpack("<2i", header)
+    saved_variable_count = abs(encoded_saved_variable_count)
+    if point_count <= 0 or saved_variable_count <= 0:
+        raise AmesimResultsError("Invalid AMESim results header.")
+
+    mapping_data = results_file.read(saved_variable_count * 4)
+    if len(mapping_data) != saved_variable_count * 4:
+        raise AmesimResultsError("AMESim results variable mapping is truncated.")
+    saved_variable_indices = struct.unpack(
+        f"<{saved_variable_count}i",
+        mapping_data,
+    )
+    if any(index < 0 or index >= len(variables) for index in saved_variable_indices):
+        raise AmesimResultsError(
+            "AMESim results variable mapping references unknown .var rows."
+        )
+
+    row_length = 1 + saved_variable_count
+    row_byte_count = row_length * 8
+    times: list[float] = []
+    series_lists: dict[str, list[float]] = {}
+    saved_paths: list[tuple[int, str]] = []
+    for column, variable_index in enumerate(saved_variable_indices, start=1):
+        data_path = variables[variable_index].data_path
+        if data_path is None:
+            continue
+        series_lists[data_path] = []
+        saved_paths.append((column, data_path))
+
+    for _row_index in range(point_count):
+        row = results_file.read(row_byte_count)
+        if len(row) != row_byte_count:
+            raise AmesimResultsError("AMESim results matrix is truncated.")
+        time_s = struct.unpack_from("<d", row, 0)[0]
+        times.append(time_s)
+        for column, data_path in saved_paths:
+            series_lists[data_path].append(
+                struct.unpack_from("<d", row, column * 8)[0]
+            )
+        # Keep one real sample after the requested stop so endpoint finite
+        # differences do not silently fall back to a backward-only slope.
+        if time_s > time_stop_s + 1.0e-12:
+            break
+
+    return AmesimResults(
+        times=tuple(times),
+        variables=variables,
+        saved_variable_indices=tuple(saved_variable_indices),
+        series_by_data_path={
+            data_path: tuple(values) for data_path, values in series_lists.items()
+        },
+        final_values_by_data_path={},
+    )
+
+
 def parse_amesim_results_bytes(
     results_data: bytes,
     variables: tuple[AmesimVariable, ...],

PythonModels/scripts/run_test_mql_full_state_comparison.py

@@ -556,8 +556,13 @@ class TestMqlFullStateComparisonPathConfig:
         / "AmesimModels"
         / "test_mql.ame"
     )
+    amesim_results_archive_path: Path | None = None
     output_dir: Path | None = None
 
+    @property
+    def resolved_amesim_results_archive_path(self) -> Path:
+        return self.amesim_results_archive_path or self.archive_path
+
 
 @dataclass(frozen=True)
 class TestMqlFullStateComparisonExecutionConfig:
@@ -595,7 +600,10 @@ def run_test_mql_full_state_comparison(
 ) -> tuple[TestMqlFullStateComparisonRun, Path]:
     config = config or TestMqlFullStateComparisonScriptConfig()
     system = TestMqlSystem(archive_path=config.paths.archive_path)
-    amesim_results = load_test_mql_amesim_results(config.paths.archive_path)
+    amesim_results = load_test_mql_amesim_results(
+        config.paths.resolved_amesim_results_archive_path,
+        time_stop_s=config.execution.solver.t_stop,
+    )
     output_schema = build_test_mql_output_schema(amesim_results)
     selected_paths = config.execution.data_paths
     spec = system.discover_pneumatic_branch_topology().chamber_segment_specs[0]
@@ -719,7 +727,10 @@ def run_test_mql_pnvo_event_boundary_diagnostic(
         state_vector=state_vector,
         data_paths=data_paths,
     )
-    amesim_results = load_test_mql_amesim_results(config.paths.archive_path)
+    amesim_results = load_test_mql_amesim_results(
+        config.paths.resolved_amesim_results_archive_path,
+        time_stop_s=event_time_s,
+    )
     amesim_values = {
         data_path: interpolate_series_value(
             amesim_results.times,
@@ -884,6 +895,24 @@ def _pnl0001_linear_conductance(
     return abs(dm1_g_s) * 1.0e-3 * sqrt(temperature_k) / pressure_drop_pa
 
 
+def _ideal_pn2vol_reference_dtemp(
+    *,
+    gas: AmesimPneumaticGas,
+    mass_kg: float,
+    temperature_k: float,
+    mass_derivative_kg_s: float,
+    enthalpy_flow_w: float,
+    heat_flow_w: float,
+    reference_temperature_k: float = 298.15,
+) -> float:
+    reference_offset_flow_w = (
+        gas.cp * reference_temperature_k - gas.cv * temperature_k
+    ) * mass_derivative_kg_s
+    return (
+        enthalpy_flow_w + reference_offset_flow_w + heat_flow_w
+    ) / (mass_kg * gas.cv)
+
+
 def _series_finite_difference_at(
     *,
     times: tuple[float, ...] | list[float],
@@ -1094,21 +1123,6 @@ def _pnl0001_energy_flow_diagnostic(
     )
     python_reference_sum = python_reference_port_1 + python_reference_node_port_2
 
-    def pn2vol2_reference_dtemp(
-        *,
-        mass_kg: float,
-        temperature_k: float,
-        mass_derivative_kg_s: float,
-        enthalpy_flow_w: float,
-        heat_flow_w: float,
-    ) -> float:
-        storage_reference_h = reference_h(temperature_k)
-        return (
-            enthalpy_flow_w
-            - storage_reference_h * mass_derivative_kg_s
-            + heat_flow_w
-        ) / (mass_kg * line.gas.cv)
-
     amesim_mgas_kg = amesim_value(f"mgas@{line_alias}") * 1.0e-3
     amesim_line_temperature = amesim_value(f"t2@{line_alias}")
     amesim_mass_derivative = (amesim_node_dm2 - amesim_dm1) * 1.0e-3
@@ -1129,28 +1143,32 @@ def _pnl0001_energy_flow_diagnostic(
     candidate_pn2pipefr_dm2i_g_s = -amesim_dm1
     candidate_pn2pipefr_dh2i = -reference_enthalpy_estimate
     candidate_storage_sdh = amesim_node_dh2 + candidate_pn2pipefr_dh2i
-    candidate_pn2pipefr_dtemp = pn2vol2_reference_dtemp(
+    candidate_pn2pipefr_dtemp = _ideal_pn2vol_reference_dtemp(
+        gas=line.gas,
         mass_kg=amesim_mgas_kg,
         temperature_k=amesim_line_temperature,
         mass_derivative_kg_s=amesim_mass_derivative,
         enthalpy_flow_w=candidate_storage_sdh,
         heat_flow_w=amesim_heat_flow,
     )
-    amesim_node_plus_dh1_dtemp = pn2vol2_reference_dtemp(
+    amesim_node_plus_dh1_dtemp = _ideal_pn2vol_reference_dtemp(
+        gas=line.gas,
         mass_kg=amesim_mgas_kg,
         temperature_k=amesim_line_temperature,
         mass_derivative_kg_s=amesim_mass_derivative,
         enthalpy_flow_w=amesim_node_dh2 + amesim_dh1,
         heat_flow_w=amesim_heat_flow,
     )
-    amesim_node_minus_dh1_dtemp = pn2vol2_reference_dtemp(
+    amesim_node_minus_dh1_dtemp = _ideal_pn2vol_reference_dtemp(
+        gas=line.gas,
         mass_kg=amesim_mgas_kg,
         temperature_k=amesim_line_temperature,
         mass_derivative_kg_s=amesim_mass_derivative,
         enthalpy_flow_w=amesim_node_dh2 - amesim_dh1,
         heat_flow_w=amesim_heat_flow,
     )
-    python_reference_node_minus_dh1_dtemp = pn2vol2_reference_dtemp(
+    python_reference_node_minus_dh1_dtemp = _ideal_pn2vol_reference_dtemp(
+        gas=line.gas,
         mass_kg=line.state.m,
         temperature_k=line_properties.T,
         mass_derivative_kg_s=rhs_diagnostic.mass_derivative_kg_s,
@@ -2034,7 +2052,16 @@ def run_test_mql_pnvo_event_window_diagnostic(
         f"xv@{orifice_alias}",
         f"dm2@{orifice_alias}",
     )
-    sample_times = (event_time_s, 0.041, 0.042, 0.045, 0.048, final_time_s)
+    sample_times = (
+        nextafter(event_time_s, 0.0),
+        0.0401,
+        0.0402,
+        0.041,
+        0.042,
+        0.045,
+        0.048,
+        final_time_s,
+    )
     state_vector_by_sample_time: dict[float, list[float]] = {}
     segment_diagnostics: list[TestMqlPnvoEventWindowSegmentDiagnostic] = []
     for segment in segments:
@@ -2081,7 +2108,10 @@ def run_test_mql_pnvo_event_window_diagnostic(
         state_vector = [row[-1] for row in solution.y]
         if not solution.success:
             break
-    amesim_results = load_test_mql_amesim_results(config.paths.archive_path)
+    amesim_results = load_test_mql_amesim_results(
+        config.paths.resolved_amesim_results_archive_path,
+        time_stop_s=final_time_s,
+    )
     sample_diagnostics = []
     for sample_time in sample_times:
         if sample_time not in state_vector_by_sample_time:
@@ -2717,17 +2747,29 @@ def main() -> None:
         action="store_true",
         help="run the segmented PNVO opening window through the t=0.05 s save point",
     )
+    parser.add_argument(
+        "--amesim-results-archive",
+        type=Path,
+        help="read AMESim observations from a separate .ame results package",
+    )
     args = parser.parse_args()
+    config = None
+    if args.amesim_results_archive is not None:
+        config = TestMqlFullStateComparisonScriptConfig(
+            paths=TestMqlFullStateComparisonPathConfig(
+                amesim_results_archive_path=args.amesim_results_archive,
+            )
+        )
     if args.pnvo_event_boundary:
-        diagnostic = run_test_mql_pnvo_event_boundary_diagnostic()
+        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:
-        diagnostic = run_test_mql_pnvo_event_window_diagnostic()
+        diagnostic = run_test_mql_pnvo_event_window_diagnostic(config)
         print(format_test_mql_pnvo_event_window_summary(diagnostic), end="")
         return
 
-    run, output_dir = run_test_mql_full_state_comparison()
+    run, output_dir = run_test_mql_full_state_comparison(config)
     print(format_test_mql_full_state_comparison_summary(run), end="")
     print(f"Output directory: {output_dir}")
 

tests/test_amesim_results.py

@@ -2,8 +2,12 @@ from __future__ import annotations
 
 import unittest
 from pathlib import Path
+from unittest.mock import Mock
 
-from PythonModels.reporting.amesim_results import load_test_mql_amesim_results
+from PythonModels.reporting.amesim_results import (
+    _resolve_result_members,
+    load_test_mql_amesim_results,
+)
 
 
 REPO_ROOT = Path(__file__).resolve().parents[1]
@@ -41,6 +45,37 @@ class AmesimResultsTests(unittest.TestCase):
         self.assertAlmostEqual(temp3[-1], self.results.final_value("temp3@pn_c1_8"))
         self.assertAlmostEqual(press3[-1], self.results.final_value("press3@pn_c1_8"))
 
+    def test_window_loader_keeps_one_sample_after_stop(self) -> None:
+        results = load_test_mql_amesim_results(TEST_MQL_AME, time_stop_s=0.04)
+
+        self.assertEqual(results.point_count, 6)
+        self.assertAlmostEqual(results.times[-2], 0.04)
+        self.assertAlmostEqual(results.times[-1], 0.05)
+        self.assertEqual(
+            len(results.series("temp@pn_c1_8")),
+            results.point_count,
+        )
+        self.assertEqual(results.final_values_by_data_path, {})
+
+    def test_auto_detects_renamed_result_members(self) -> None:
+        archive = Mock()
+        archive.getnames.return_value = (
+            "renamed_model_.var",
+            "renamed_model_.results",
+            "renamed_model_.sim",
+        )
+
+        pair = _resolve_result_members(
+            archive,
+            var_member=None,
+            results_member=None,
+        )
+
+        self.assertEqual(
+            pair,
+            ("renamed_model_.var", "renamed_model_.results"),
+        )
+
 
 if __name__ == "__main__":
     unittest.main()

tests/test_run_test_mql_full_state_comparison.py

@@ -26,6 +26,7 @@ from PythonModels.scripts.run_test_mql_full_state_comparison import (
     format_test_mql_pnvo_event_boundary_summary,
     format_test_mql_pnvo_event_window_summary,
     run_test_mql_full_state_comparison,
+    _ideal_pn2vol_reference_dtemp,
     _series_sample_bracket,
 )
 
@@ -35,6 +36,28 @@ TEST_MQL_AME = REPO_ROOT / "AmesimModels" / "test_mql.ame"
 
 
 class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase):
+    def test_ideal_pn2vol_reference_dtemp_uses_internal_energy_mass_offset(
+        self,
+    ) -> None:
+        from PythonModels.components.amesim_pneumatic import HELIUM_PNEUMATIC_GAS
+
+        gas = HELIUM_PNEUMATIC_GAS
+        result = _ideal_pn2vol_reference_dtemp(
+            gas=gas,
+            mass_kg=0.002,
+            temperature_k=320.0,
+            mass_derivative_kg_s=0.1,
+            enthalpy_flow_w=-100.0,
+            heat_flow_w=5.0,
+        )
+
+        expected = (
+            -100.0
+            + (gas.cp * 298.15 - gas.cv * 320.0) * 0.1
+            + 5.0
+        ) / (0.002 * gas.cv)
+        self.assertAlmostEqual(result, expected)
+
     def test_series_sample_bracket_distinguishes_exact_and_interpolated_times(self) -> None:
         times = (0.04000000000000009, 0.050000000000000114)