PythonModels/scripts/run_test_mql_full_state_comparison.py

@@ -26,6 +26,7 @@ from PythonModels.systems.test_mql import (
     TestMqlPnl0001LineRhsDiagnostic,
     TestMqlSimulationResult,
     TestMqlSystem,
+    TestMqlVariableChamberRhsDiagnostic,
 )
 from PythonModels.systems.test_mql_pneumatic import AMESIM_REFERENCE_PRESSURE_PA
 
@@ -307,6 +308,76 @@ class TestMqlPnvoFlowParameterDiagnostic:
     python_to_amesim_cm_ratio: float
 
 
+@dataclass(frozen=True)
+class TestMqlPnl0001EnergyFlowDiagnostic:
+    line_alias: str
+    chamber_alias: str
+    node_alias: str
+    time_s: float
+    amesim_port_1_enthalpy_flow_w: float
+    amesim_port_1_mass_flow_g_s: float
+    amesim_node_port_2_enthalpy_flow_w: float
+    amesim_node_port_2_mass_flow_g_s: float
+    amesim_observed_port_enthalpy_sum_w: float
+    python_port_1_internal_energy_flow_w: float
+    python_port_2_internal_energy_flow_w: float
+    python_current_energy_derivative_w: float
+    python_direct_enthalpy_port_1_flow_w: float
+    python_direct_enthalpy_port_2_flow_w: float
+    python_direct_enthalpy_energy_derivative_w: float
+    python_direct_minus_current_energy_derivative_w: float
+    python_p4_node_port_2_enthalpy_flow_w: float
+    python_p4_node_port_2_mass_flow_g_s: float
+    python_p4_node_port_2_connected_h_j_kg: float
+    python_p4_node_port_2_connected_temperature_k: float
+    python_p4_node_port_2_connected_h_delta_to_line_h_j_kg: float
+    python_p4_node_port_2_counterfactual_dtemp_k_s: float
+    reference_temperature_k: float
+    amesim_port_1_reference_enthalpy_estimate_w: float
+    amesim_port_1_reference_enthalpy_error_w: float
+    amesim_candidate_pn2pipefr_dm2i_g_s: float
+    amesim_candidate_pn2pipefr_dh2i_w: float
+    amesim_candidate_storage_enthalpy_sum_w: float
+    amesim_candidate_pn2pipefr_dtemp_k_s: float
+    amesim_candidate_pn2pipefr_dtemp_residual_k_s: float
+    python_reference_port_1_enthalpy_flow_w: float
+    python_reference_node_port_2_enthalpy_flow_w: float
+    python_reference_observed_port_enthalpy_sum_w: float
+    python_reference_port_1_to_amesim_error_w: float
+    python_reference_node_port_2_to_amesim_error_w: float
+    python_reference_sum_to_amesim_error_w: float
+    python_current_line_temperature_derivative_k_s: float
+    amesim_line_temperature_derivative_fd_k_s: float
+    amesim_pn2vol2_dtemp_node_plus_dh1_k_s: float
+    amesim_pn2vol2_dtemp_node_minus_dh1_k_s: float
+    python_reference_pn2vol2_dtemp_node_minus_dh1_k_s: float
+
+
+@dataclass(frozen=True)
+class TestMqlPnch012EnergyEquationDiagnostic:
+    chamber_alias: str
+    piston_alias: str
+    line_alias: str
+    time_s: float
+    amesim_port_1_enthalpy_flow_w: float
+    amesim_port_1_mass_flow_g_s: float
+    amesim_chamber_mass_derivative_fd_g_s: float
+    amesim_chamber_mass_derivative_residual_g_s: float
+    amesim_chamber_temperature_derivative_fd_k_s: float
+    amesim_chamber_volume_rate_fd_m3_s: float
+    amesim_piston_volume_rate_m3_s: float
+    amesim_heat_flow_w: float
+    amesim_boundary_work_fd_volume_w: float
+    amesim_pn2vol_reference_dtemp_fd_volume_k_s: float
+    amesim_pn2vol_reference_dtemp_fd_volume_residual_k_s: float
+    amesim_pn2vol_reference_dtemp_piston_volume_k_s: float
+    amesim_pn2vol_reference_dtemp_piston_volume_residual_k_s: float
+    python_port_a_mass_flow_kg_s: float
+    python_current_energy_derivative_w: float
+    python_current_chamber_temperature_derivative_k_s: float
+    reference_temperature_k: float
+
+
 @dataclass(frozen=True)
 class TestMqlPnvoEventWindowSampleDiagnostic:
     time_s: float
@@ -322,6 +393,12 @@ class TestMqlPnvoEventWindowSampleDiagnostic:
     pnvo_flow_parameter_diagnostics: tuple[
         TestMqlPnvoFlowParameterDiagnostic, ...
     ] = field(default_factory=tuple)
+    pnl0001_energy_flow_diagnostics: tuple[
+        TestMqlPnl0001EnergyFlowDiagnostic, ...
+    ] = field(default_factory=tuple)
+    pnch012_energy_equation_diagnostics: tuple[
+        TestMqlPnch012EnergyEquationDiagnostic, ...
+    ] = field(default_factory=tuple)
 
     def abs_error(self, data_path: str) -> float:
         return abs(
@@ -683,6 +760,394 @@ def _pnl0001_linear_conductance(
     return abs(dm1_g_s) * 1.0e-3 * sqrt(temperature_k) / pressure_drop_pa
 
 
+def _series_finite_difference_at(
+    *,
+    times: tuple[float, ...] | list[float],
+    values: tuple[float, ...] | list[float],
+    time_s: float,
+) -> float:
+    if len(times) != len(values):
+        raise ValueError("times and values must have equal length")
+    if len(times) < 2:
+        raise ValueError("at least two samples are required")
+    index = min(range(len(times)), key=lambda idx: abs(times[idx] - time_s))
+    if index == 0:
+        left = 0
+        right = 1
+    elif index == len(times) - 1:
+        left = len(times) - 2
+        right = len(times) - 1
+    else:
+        left = index - 1
+        right = index + 1
+    dt = times[right] - times[left]
+    if dt == 0.0:
+        raise ValueError("finite difference time interval must be non-zero")
+    return (values[right] - values[left]) / dt
+
+
+def _pnl0001_energy_flow_diagnostic(
+    *,
+    closure: object,
+    amesim_results: AmesimResults,
+    state_vector: list[float],
+    rhs_diagnostic: TestMqlPnl0001LineRhsDiagnostic,
+    line_alias: str,
+    chamber_alias: str,
+    node_alias: str,
+    time_s: float,
+) -> TestMqlPnl0001EnergyFlowDiagnostic:
+    if line_alias != "pneumatic_69" or chamber_alias != "pn_c1_8":
+        raise KeyError(f"Unsupported PNL0001 energy diagnostic: {line_alias}")
+    if node_alias != "pnnode4_16":
+        raise KeyError(f"Unsupported PNL0001 node diagnostic: {node_alias}")
+
+    snapshot = closure.snapshot_at(time_s, state_vector).pneumatic
+    line = closure.pneumatic_closure.components.p4_port3_remote_primary_line
+    line_properties = snapshot.p4_port3_remote_primary_line
+    chamber_properties = snapshot.p4_port3_remote_primary_chamber
+
+    def amesim_value(data_path: str) -> float:
+        return interpolate_series_value(
+            amesim_results.times,
+            amesim_results.series(data_path),
+            time_s,
+        )
+
+    amesim_dh1 = amesim_value(f"dh1@{line_alias}")
+    amesim_dm1 = amesim_value(f"dm1@{line_alias}")
+    amesim_node_dh2 = amesim_value(f"dh2@{node_alias}")
+    amesim_node_dm2 = amesim_value(f"dm2@{node_alias}")
+    direct_port_1 = rhs_diagnostic.chamber_to_line_flow_kg_s * (
+        chamber_properties.h
+        if rhs_diagnostic.chamber_to_line_flow_kg_s > 0.0
+        else line_properties.h
+    )
+    direct_port_2 = rhs_diagnostic.node_to_line_flow_kg_s * line_properties.h
+    direct_energy_derivative = (
+        direct_port_1 + direct_port_2 + rhs_diagnostic.thermal_energy_flow_w
+    )
+    p4_node_port_2_mass_flow_kg_s = (
+        snapshot.p4_port3_remote_balance.port_2_mass_flow_g_s * 1.0e-3
+    )
+    if abs(p4_node_port_2_mass_flow_kg_s) <= 1.0e-12:
+        p4_node_port_2_connected_h = line_properties.h
+    else:
+        p4_node_port_2_connected_h = (
+            snapshot.p4_port3_remote_balance.port_2_enthalpy_flow_w
+            / p4_node_port_2_mass_flow_kg_s
+        )
+    p4_node_port_2_inlet_u = (
+        p4_node_port_2_connected_h / line.gas.gamma
+        if rhs_diagnostic.node_to_line_flow_kg_s > 0.0
+        else line_properties.u
+    )
+    p4_node_port_2_counterfactual_energy_derivative = (
+        rhs_diagnostic.port_1_energy_flow_w
+        + rhs_diagnostic.node_to_line_flow_kg_s * p4_node_port_2_inlet_u
+        + rhs_diagnostic.thermal_energy_flow_w
+    )
+    p4_node_port_2_counterfactual_dtemp = (
+        p4_node_port_2_counterfactual_energy_derivative
+        - line_properties.u * rhs_diagnostic.mass_derivative_kg_s
+    ) / (line.state.m * line.gas.cv)
+    reference_temperature_k = 298.15
+
+    def reference_h(temperature_k: float) -> float:
+        return line.gas.cp * (temperature_k - reference_temperature_k)
+
+    def reference_enthalpy_flow(
+        *,
+        flow_kg_s: float,
+        node_temperature_k: float,
+        connected_temperature_k: float,
+    ) -> float:
+        stream_temperature_k = (
+            node_temperature_k if flow_kg_s >= 0.0 else connected_temperature_k
+        )
+        return flow_kg_s * reference_h(stream_temperature_k)
+
+    amesim_stream_temperature = (
+        amesim_value(f"t2@{line_alias}")
+        if amesim_dm1 >= 0.0
+        else amesim_value(f"temp@{chamber_alias}")
+    )
+    reference_enthalpy_estimate = (
+        amesim_dm1 * 1.0e-3 * reference_h(amesim_stream_temperature)
+    )
+    python_dm1_kg_s = -rhs_diagnostic.chamber_to_line_flow_kg_s
+    python_port_1_stream_temperature_k = (
+        line_properties.T if python_dm1_kg_s >= 0.0 else chamber_properties.T
+    )
+    python_reference_port_1 = python_dm1_kg_s * reference_h(
+        python_port_1_stream_temperature_k
+    )
+    python_reference_node_port_2 = sum(
+        (
+            reference_enthalpy_flow(
+                flow_kg_s=snapshot.p4_port3_remote_node_to_line_flow,
+                node_temperature_k=snapshot.p4_port3_remote_primary_line.T,
+                connected_temperature_k=snapshot.p4_port3_line.T,
+            ),
+            reference_enthalpy_flow(
+                flow_kg_s=snapshot.p4_port3_remote_to_port3_line_flow,
+                node_temperature_k=snapshot.p4_port3_remote_primary_line.T,
+                connected_temperature_k=snapshot.p4_port3_remote_port3_line.T,
+            ),
+            reference_enthalpy_flow(
+                flow_kg_s=snapshot.p4_port3_remote_orifice_to_node_flow,
+                node_temperature_k=snapshot.p4_port3_remote_orifice_line_port_2.T,
+                connected_temperature_k=snapshot.p4_port3_remote_primary_line.T,
+            ),
+        )
+    )
+    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
+    amesim_heat_flow = (
+        line.heat_transfer_coefficient
+        * line.heat_transfer_area
+        * (line.external_temperature - amesim_line_temperature)
+    )
+    current_line_dtemp = (
+        rhs_diagnostic.energy_derivative_w
+        - line_properties.u * rhs_diagnostic.mass_derivative_kg_s
+    ) / (line.state.m * line.gas.cv)
+    amesim_line_dtemp_fd = _series_finite_difference_at(
+        times=amesim_results.times,
+        values=amesim_results.series(f"t2@{line_alias}"),
+        time_s=time_s,
+    )
+    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(
+        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(
+        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(
+        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(
+        mass_kg=line.state.m,
+        temperature_k=line_properties.T,
+        mass_derivative_kg_s=rhs_diagnostic.mass_derivative_kg_s,
+        enthalpy_flow_w=python_reference_node_port_2 - python_reference_port_1,
+        heat_flow_w=rhs_diagnostic.thermal_energy_flow_w,
+    )
+    return TestMqlPnl0001EnergyFlowDiagnostic(
+        line_alias=line_alias,
+        chamber_alias=chamber_alias,
+        node_alias=node_alias,
+        time_s=time_s,
+        amesim_port_1_enthalpy_flow_w=amesim_dh1,
+        amesim_port_1_mass_flow_g_s=amesim_dm1,
+        amesim_node_port_2_enthalpy_flow_w=amesim_node_dh2,
+        amesim_node_port_2_mass_flow_g_s=amesim_node_dm2,
+        amesim_observed_port_enthalpy_sum_w=amesim_dh1 + amesim_node_dh2,
+        python_port_1_internal_energy_flow_w=rhs_diagnostic.port_1_energy_flow_w,
+        python_port_2_internal_energy_flow_w=rhs_diagnostic.port_2_energy_flow_w,
+        python_current_energy_derivative_w=rhs_diagnostic.energy_derivative_w,
+        python_direct_enthalpy_port_1_flow_w=direct_port_1,
+        python_direct_enthalpy_port_2_flow_w=direct_port_2,
+        python_direct_enthalpy_energy_derivative_w=direct_energy_derivative,
+        python_direct_minus_current_energy_derivative_w=(
+            direct_energy_derivative - rhs_diagnostic.energy_derivative_w
+        ),
+        python_p4_node_port_2_enthalpy_flow_w=(
+            snapshot.p4_port3_remote_balance.port_2_enthalpy_flow_w
+        ),
+        python_p4_node_port_2_mass_flow_g_s=(
+            snapshot.p4_port3_remote_balance.port_2_mass_flow_g_s
+        ),
+        python_p4_node_port_2_connected_h_j_kg=p4_node_port_2_connected_h,
+        python_p4_node_port_2_connected_temperature_k=(
+            p4_node_port_2_connected_h / line.gas.cp
+        ),
+        python_p4_node_port_2_connected_h_delta_to_line_h_j_kg=(
+            p4_node_port_2_connected_h - line_properties.h
+        ),
+        python_p4_node_port_2_counterfactual_dtemp_k_s=(
+            p4_node_port_2_counterfactual_dtemp
+        ),
+        reference_temperature_k=reference_temperature_k,
+        amesim_port_1_reference_enthalpy_estimate_w=reference_enthalpy_estimate,
+        amesim_port_1_reference_enthalpy_error_w=(
+            reference_enthalpy_estimate - amesim_dh1
+        ),
+        amesim_candidate_pn2pipefr_dm2i_g_s=candidate_pn2pipefr_dm2i_g_s,
+        amesim_candidate_pn2pipefr_dh2i_w=candidate_pn2pipefr_dh2i,
+        amesim_candidate_storage_enthalpy_sum_w=candidate_storage_sdh,
+        amesim_candidate_pn2pipefr_dtemp_k_s=candidate_pn2pipefr_dtemp,
+        amesim_candidate_pn2pipefr_dtemp_residual_k_s=(
+            candidate_pn2pipefr_dtemp - amesim_line_dtemp_fd
+        ),
+        python_reference_port_1_enthalpy_flow_w=python_reference_port_1,
+        python_reference_node_port_2_enthalpy_flow_w=python_reference_node_port_2,
+        python_reference_observed_port_enthalpy_sum_w=python_reference_sum,
+        python_reference_port_1_to_amesim_error_w=(
+            python_reference_port_1 - amesim_dh1
+        ),
+        python_reference_node_port_2_to_amesim_error_w=(
+            python_reference_node_port_2 - amesim_node_dh2
+        ),
+        python_reference_sum_to_amesim_error_w=(
+            python_reference_sum - (amesim_dh1 + amesim_node_dh2)
+        ),
+        python_current_line_temperature_derivative_k_s=current_line_dtemp,
+        amesim_line_temperature_derivative_fd_k_s=amesim_line_dtemp_fd,
+        amesim_pn2vol2_dtemp_node_plus_dh1_k_s=amesim_node_plus_dh1_dtemp,
+        amesim_pn2vol2_dtemp_node_minus_dh1_k_s=amesim_node_minus_dh1_dtemp,
+        python_reference_pn2vol2_dtemp_node_minus_dh1_k_s=(
+            python_reference_node_minus_dh1_dtemp
+        ),
+    )
+
+
+
+def _pnch012_energy_equation_diagnostic(
+    *,
+    closure: object,
+    amesim_results: AmesimResults,
+    state_vector: list[float],
+    rhs_diagnostic: TestMqlVariableChamberRhsDiagnostic,
+    chamber_alias: str,
+    piston_alias: str,
+    line_alias: str,
+    time_s: float,
+) -> TestMqlPnch012EnergyEquationDiagnostic:
+    if chamber_alias != "pn_c1_8" or piston_alias != "pn_brp2_8":
+        raise KeyError(f"Unsupported PNCH012 energy diagnostic: {chamber_alias}")
+    if line_alias != "pneumatic_69":
+        raise KeyError(f"Unsupported PNCH012 line diagnostic: {line_alias}")
+
+    snapshot = closure.snapshot_at(time_s, state_vector).pneumatic
+    chamber = closure.pneumatic_closure.components.p4_port3_remote_primary_chamber
+    chamber_properties = snapshot.p4_port3_remote_primary_chamber
+    reference_temperature_k = 298.15
+
+    def amesim_value(data_path: str) -> float:
+        return interpolate_series_value(
+            amesim_results.times,
+            amesim_results.series(data_path),
+            time_s,
+        )
+
+    amesim_dh1 = amesim_value(f"dh1@{line_alias}")
+    amesim_dm1_g_s = amesim_value(f"dm1@{line_alias}")
+    amesim_chamber_temperature = amesim_value(f"temp@{chamber_alias}")
+    amesim_chamber_pressure_abs = (
+        amesim_value(f"press@{chamber_alias}") + AMESIM_REFERENCE_PRESSURE_PA
+    )
+    amesim_mgas_kg = amesim_value(f"mgas1@{chamber_alias}") * 1.0e-3
+    amesim_mass_derivative_fd_g_s = _series_finite_difference_at(
+        times=amesim_results.times,
+        values=amesim_results.series(f"mgas1@{chamber_alias}"),
+        time_s=time_s,
+    )
+    amesim_temperature_derivative_fd = _series_finite_difference_at(
+        times=amesim_results.times,
+        values=amesim_results.series(f"temp@{chamber_alias}"),
+        time_s=time_s,
+    )
+    amesim_volume_rate_fd = (
+        _series_finite_difference_at(
+            times=amesim_results.times,
+            values=amesim_results.series(f"vol@{chamber_alias}"),
+            time_s=time_s,
+        )
+        * 1.0e-6
+    )
+    amesim_piston_volume_rate = amesim_value(f"vvol1@{piston_alias}") * (1.0 / 60000.0)
+    amesim_heat_flow = (
+        chamber.heat_transfer_coefficient
+        * chamber.heat_transfer_area
+        * (chamber.external_temperature - amesim_chamber_temperature)
+    )
+
+    def reference_pn2vol_dtemp(volume_rate_m3_s: float) -> float:
+        mass_flow_kg_s = amesim_dm1_g_s * 1.0e-3
+        reference_offset_flow_w = (
+            chamber.gas.cp * reference_temperature_k
+            - chamber.gas.cv * amesim_chamber_temperature
+        ) * mass_flow_kg_s
+        return (
+            amesim_dh1
+            + reference_offset_flow_w
+            + amesim_heat_flow
+            - amesim_chamber_pressure_abs * volume_rate_m3_s
+        ) / (amesim_mgas_kg * chamber.gas.cv)
+
+    dtemp_from_fd_volume = reference_pn2vol_dtemp(amesim_volume_rate_fd)
+    dtemp_from_piston_volume = reference_pn2vol_dtemp(amesim_piston_volume_rate)
+    python_current_dtemp = (
+        rhs_diagnostic.energy_derivative_w
+        - chamber_properties.u * rhs_diagnostic.mass_derivative_kg_s
+    ) / (chamber.state.m * chamber.gas.cv)
+
+    return TestMqlPnch012EnergyEquationDiagnostic(
+        chamber_alias=chamber_alias,
+        piston_alias=piston_alias,
+        line_alias=line_alias,
+        time_s=time_s,
+        amesim_port_1_enthalpy_flow_w=amesim_dh1,
+        amesim_port_1_mass_flow_g_s=amesim_dm1_g_s,
+        amesim_chamber_mass_derivative_fd_g_s=amesim_mass_derivative_fd_g_s,
+        amesim_chamber_mass_derivative_residual_g_s=(
+            amesim_dm1_g_s - amesim_mass_derivative_fd_g_s
+        ),
+        amesim_chamber_temperature_derivative_fd_k_s=amesim_temperature_derivative_fd,
+        amesim_chamber_volume_rate_fd_m3_s=amesim_volume_rate_fd,
+        amesim_piston_volume_rate_m3_s=amesim_piston_volume_rate,
+        amesim_heat_flow_w=amesim_heat_flow,
+        amesim_boundary_work_fd_volume_w=(
+            -amesim_chamber_pressure_abs * amesim_volume_rate_fd
+        ),
+        amesim_pn2vol_reference_dtemp_fd_volume_k_s=dtemp_from_fd_volume,
+        amesim_pn2vol_reference_dtemp_fd_volume_residual_k_s=(
+            dtemp_from_fd_volume - amesim_temperature_derivative_fd
+        ),
+        amesim_pn2vol_reference_dtemp_piston_volume_k_s=dtemp_from_piston_volume,
+        amesim_pn2vol_reference_dtemp_piston_volume_residual_k_s=(
+            dtemp_from_piston_volume - amesim_temperature_derivative_fd
+        ),
+        python_port_a_mass_flow_kg_s=rhs_diagnostic.port_a_mass_flow_kg_s,
+        python_current_energy_derivative_w=rhs_diagnostic.energy_derivative_w,
+        python_current_chamber_temperature_derivative_k_s=python_current_dtemp,
+        reference_temperature_k=reference_temperature_k,
+    )
+
 def _pnvo_flow_parameter_diagnostic(
     *,
     closure: object,
@@ -930,6 +1395,37 @@ def run_test_mql_pnvo_event_window_diagnostic(
                 time_s=sample_time,
             ),
         )
+        pnl0001_energy_flow_diagnostics = (
+            _pnl0001_energy_flow_diagnostic(
+                closure=closure,
+                amesim_results=amesim_results,
+                state_vector=state_vector_by_sample_time[sample_time],
+                rhs_diagnostic=pnl0001_rhs_diagnostics[0],
+                line_alias="pneumatic_69",
+                chamber_alias="pn_c1_8",
+                node_alias="pnnode4_16",
+                time_s=sample_time,
+            ),
+        )
+        pnch012_rhs_diagnostics = (
+            closure.variable_chamber_rhs_diagnostic(
+                chamber_alias="pn_c1_8",
+                state_vector=state_vector_by_sample_time[sample_time],
+                time_s=sample_time,
+            ),
+        )
+        pnch012_energy_equation_diagnostics = (
+            _pnch012_energy_equation_diagnostic(
+                closure=closure,
+                amesim_results=amesim_results,
+                state_vector=state_vector_by_sample_time[sample_time],
+                rhs_diagnostic=pnch012_rhs_diagnostics[0],
+                chamber_alias="pn_c1_8",
+                piston_alias="pn_brp2_8",
+                line_alias="pneumatic_69",
+                time_s=sample_time,
+            ),
+        )
         sample_diagnostics.append(
             TestMqlPnvoEventWindowSampleDiagnostic(
                 time_s=sample_time,
@@ -939,6 +1435,10 @@ def run_test_mql_pnvo_event_window_diagnostic(
                 pnl0001_rhs_diagnostics=pnl0001_rhs_diagnostics,
                 pnl0001_pressure_loss_diagnostics=pnl0001_pressure_loss_diagnostics,
                 pnvo_flow_parameter_diagnostics=pnvo_flow_parameter_diagnostics,
+                pnl0001_energy_flow_diagnostics=pnl0001_energy_flow_diagnostics,
+                pnch012_energy_equation_diagnostics=(
+                    pnch012_energy_equation_diagnostics
+                ),
             )
         )
     python_values = closure.data_path_values(
@@ -1035,6 +1535,108 @@ def format_test_mql_pnvo_event_window_summary(
                     f"upstream_t={flow_parameter.python_upstream_temperature_k}, "
                     f"amesim_gasvel={flow_parameter.amesim_gas_velocity_m_s}"
                 )
+            for energy_flow in sample.pnl0001_energy_flow_diagnostics:
+                lines.append(
+                    f"    - energy_flow@{energy_flow.line_alias}: "
+                    f"amesim_dh1={energy_flow.amesim_port_1_enthalpy_flow_w}, "
+                    f"amesim_dm1={energy_flow.amesim_port_1_mass_flow_g_s}, "
+                    f"amesim_node_dh2={energy_flow.amesim_node_port_2_enthalpy_flow_w}, "
+                    f"amesim_node_dm2={energy_flow.amesim_node_port_2_mass_flow_g_s}, "
+                    f"amesim_observed_sdh="
+                    f"{energy_flow.amesim_observed_port_enthalpy_sum_w}, "
+                    f"python_port1_u="
+                    f"{energy_flow.python_port_1_internal_energy_flow_w}, "
+                    f"python_port2_u="
+                    f"{energy_flow.python_port_2_internal_energy_flow_w}, "
+                    f"python_dU_dt={energy_flow.python_current_energy_derivative_w}, "
+                    f"direct_h_dU_dt="
+                    f"{energy_flow.python_direct_enthalpy_energy_derivative_w}, "
+                    f"direct_minus_current="
+                    f"{energy_flow.python_direct_minus_current_energy_derivative_w}, "
+                    f"python_p4_node_dh2="
+                    f"{energy_flow.python_p4_node_port_2_enthalpy_flow_w}, "
+                    f"python_p4_node_dm2="
+                    f"{energy_flow.python_p4_node_port_2_mass_flow_g_s}, "
+                    f"python_p4_node_h2="
+                    f"{energy_flow.python_p4_node_port_2_connected_h_j_kg}, "
+                    f"python_p4_node_t2="
+                    f"{energy_flow.python_p4_node_port_2_connected_temperature_k}, "
+                    f"python_p4_node_h2_delta="
+                    f"{energy_flow.python_p4_node_port_2_connected_h_delta_to_line_h_j_kg}, "
+                    f"python_p4_node_counterfactual_dT="
+                    f"{energy_flow.python_p4_node_port_2_counterfactual_dtemp_k_s}, "
+                    f"href_t={energy_flow.reference_temperature_k}, "
+                    f"ref_dh1_est="
+                    f"{energy_flow.amesim_port_1_reference_enthalpy_estimate_w}, "
+                    f"ref_dh1_error="
+                    f"{energy_flow.amesim_port_1_reference_enthalpy_error_w}, "
+                    f"candidate_dm2i="
+                    f"{energy_flow.amesim_candidate_pn2pipefr_dm2i_g_s}, "
+                    f"candidate_dh2i="
+                    f"{energy_flow.amesim_candidate_pn2pipefr_dh2i_w}, "
+                    f"candidate_storage_sdh="
+                    f"{energy_flow.amesim_candidate_storage_enthalpy_sum_w}, "
+                    f"candidate_pn2pipefr_dT="
+                    f"{energy_flow.amesim_candidate_pn2pipefr_dtemp_k_s}, "
+                    f"candidate_pn2pipefr_dT_residual="
+                    f"{energy_flow.amesim_candidate_pn2pipefr_dtemp_residual_k_s}, "
+                    f"python_ref_dh1="
+                    f"{energy_flow.python_reference_port_1_enthalpy_flow_w}, "
+                    f"python_ref_node_dh2="
+                    f"{energy_flow.python_reference_node_port_2_enthalpy_flow_w}, "
+                    f"python_ref_sdh="
+                    f"{energy_flow.python_reference_observed_port_enthalpy_sum_w}, "
+                    f"python_ref_dh1_error="
+                    f"{energy_flow.python_reference_port_1_to_amesim_error_w}, "
+                    f"python_ref_node_dh2_error="
+                    f"{energy_flow.python_reference_node_port_2_to_amesim_error_w}, "
+                    f"python_ref_sdh_error="
+                    f"{energy_flow.python_reference_sum_to_amesim_error_w}, "
+                    f"python_current_dT="
+                    f"{energy_flow.python_current_line_temperature_derivative_k_s}, "
+                    f"amesim_fd_dT="
+                    f"{energy_flow.amesim_line_temperature_derivative_fd_k_s}, "
+                    f"amesim_pn2vol2_dT_node_plus_dh1="
+                    f"{energy_flow.amesim_pn2vol2_dtemp_node_plus_dh1_k_s}, "
+                    f"amesim_pn2vol2_dT_node_minus_dh1="
+                    f"{energy_flow.amesim_pn2vol2_dtemp_node_minus_dh1_k_s}, "
+                    f"python_ref_pn2vol2_dT_node_minus_dh1="
+                    f"{energy_flow.python_reference_pn2vol2_dtemp_node_minus_dh1_k_s}"
+                )
+            for chamber_energy in sample.pnch012_energy_equation_diagnostics:
+                lines.append(
+                    f"    - pnch012_energy@{chamber_energy.chamber_alias}: "
+                    f"amesim_dh1={chamber_energy.amesim_port_1_enthalpy_flow_w}, "
+                    f"amesim_dm1={chamber_energy.amesim_port_1_mass_flow_g_s}, "
+                    f"amesim_fd_dmgas="
+                    f"{chamber_energy.amesim_chamber_mass_derivative_fd_g_s}, "
+                    f"amesim_dmgas_residual="
+                    f"{chamber_energy.amesim_chamber_mass_derivative_residual_g_s}, "
+                    f"amesim_fd_dT="
+                    f"{chamber_energy.amesim_chamber_temperature_derivative_fd_k_s}, "
+                    f"amesim_fd_dvol="
+                    f"{chamber_energy.amesim_chamber_volume_rate_fd_m3_s}, "
+                    f"amesim_piston_dvol="
+                    f"{chamber_energy.amesim_piston_volume_rate_m3_s}, "
+                    f"amesim_dq={chamber_energy.amesim_heat_flow_w}, "
+                    f"amesim_boundary_work_fd="
+                    f"{chamber_energy.amesim_boundary_work_fd_volume_w}, "
+                    f"amesim_pn2vol_ref_dT_fdvol="
+                    f"{chamber_energy.amesim_pn2vol_reference_dtemp_fd_volume_k_s}, "
+                    f"amesim_pn2vol_ref_dT_fdvol_residual="
+                    f"{chamber_energy.amesim_pn2vol_reference_dtemp_fd_volume_residual_k_s}, "
+                    f"amesim_pn2vol_ref_dT_pistonvol="
+                    f"{chamber_energy.amesim_pn2vol_reference_dtemp_piston_volume_k_s}, "
+                    f"amesim_pn2vol_ref_dT_pistonvol_residual="
+                    f"{chamber_energy.amesim_pn2vol_reference_dtemp_piston_volume_residual_k_s}, "
+                    f"python_port_a_dm="
+                    f"{chamber_energy.python_port_a_mass_flow_kg_s}, "
+                    f"python_dU_dt="
+                    f"{chamber_energy.python_current_energy_derivative_w}, "
+                    f"python_current_dT="
+                    f"{chamber_energy.python_current_chamber_temperature_derivative_k_s}, "
+                    f"href_t={chamber_energy.reference_temperature_k}"
+                )
     lines.append("Final comparison:")
     for data_path in diagnostic.data_paths:
         lines.append(

PythonModels/systems/test_mql_closure.py

@@ -1785,6 +1785,17 @@ class TestMqlPn3P4NodeChamberSegmentClosure:
             connected_h=connected_h,
         )
 
+    @staticmethod
+    def _enthalpy_flow_for_node_port(
+        *,
+        flow_kg_s: float,
+        node_h: float,
+        connected_h: float,
+        port_mass_flow_kg_s: float,
+    ) -> float:
+        upstream_h = node_h if flow_kg_s >= 0.0 else connected_h
+        return port_mass_flow_kg_s * upstream_h
+
     def _two_state_orifice_line_flows(
         self,
         *,
@@ -1886,16 +1897,18 @@ class TestMqlPn3P4NodeChamberSegmentClosure:
         return node.balance(
             port_2_temperature_k=primary_properties.T,
             port_2_pressure_pa=primary_properties.p,
-            port_1_enthalpy_flow_w=self._enthalpy_flow_from_node(
+            port_1_enthalpy_flow_w=self._enthalpy_flow_for_node_port(
                 flow_kg_s=port_1.flow_kg_s,
                 node_h=port_1.node_h,
                 connected_h=port_1.connected_h,
+                port_mass_flow_kg_s=port_1.port_mass_flow_kg_s,
             ),
             port_1_mass_flow_g_s=port_1.port_mass_flow_kg_s * 1.0e3,
-            port_3_enthalpy_flow_w=self._enthalpy_flow_from_node(
+            port_3_enthalpy_flow_w=self._enthalpy_flow_for_node_port(
                 flow_kg_s=port_3.flow_kg_s,
                 node_h=port_3.node_h,
                 connected_h=port_3.connected_h,
+                port_mass_flow_kg_s=port_3.port_mass_flow_kg_s,
             ),
             port_3_mass_flow_g_s=port_3.port_mass_flow_kg_s * 1.0e3,
         )
@@ -1912,22 +1925,25 @@ class TestMqlPn3P4NodeChamberSegmentClosure:
         return node.balance(
             port_2_temperature_k=primary_properties.T,
             port_2_pressure_pa=primary_properties.p,
-            port_1_enthalpy_flow_w=self._enthalpy_flow_from_node(
+            port_1_enthalpy_flow_w=self._enthalpy_flow_for_node_port(
                 flow_kg_s=port_1.flow_kg_s,
                 node_h=port_1.node_h,
                 connected_h=port_1.connected_h,
+                port_mass_flow_kg_s=port_1.port_mass_flow_kg_s,
             ),
             port_1_mass_flow_g_s=port_1.port_mass_flow_kg_s * 1.0e3,
-            port_3_enthalpy_flow_w=self._enthalpy_flow_from_node(
+            port_3_enthalpy_flow_w=self._enthalpy_flow_for_node_port(
                 flow_kg_s=port_3.flow_kg_s,
                 node_h=port_3.node_h,
                 connected_h=port_3.connected_h,
+                port_mass_flow_kg_s=port_3.port_mass_flow_kg_s,
             ),
             port_3_mass_flow_g_s=port_3.port_mass_flow_kg_s * 1.0e3,
-            port_4_enthalpy_flow_w=self._enthalpy_flow_from_node(
+            port_4_enthalpy_flow_w=self._enthalpy_flow_for_node_port(
                 flow_kg_s=port_4.flow_kg_s,
                 node_h=port_4.node_h,
                 connected_h=port_4.connected_h,
+                port_mass_flow_kg_s=port_4.port_mass_flow_kg_s,
             ),
             port_4_mass_flow_g_s=port_4.port_mass_flow_kg_s * 1.0e3,
         )

tests/test_run_test_mql_full_state_comparison.py

@@ -11,6 +11,8 @@ from PythonModels.scripts.run_test_mql_full_state_comparison import (
     TestMqlFullStateComparisonExecutionConfig,
     TestMqlFullStateComparisonPathConfig,
     TestMqlFullStateComparisonScriptConfig,
+    TestMqlPnl0001EnergyFlowDiagnostic,
+    TestMqlPnch012EnergyEquationDiagnostic,
     TestMqlPnl0001PressureLossCalibrationDiagnostic,
     TestMqlPnvoEventBoundaryDiagnostic,
     TestMqlPnvoEventWindowDiagnostic,
@@ -284,6 +286,76 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase):
                     python_to_amesim_cm_ratio=0.93,
                 ),
             ),
+            pnl0001_energy_flow_diagnostics=(
+                TestMqlPnl0001EnergyFlowDiagnostic(
+                    line_alias="pneumatic_69",
+                    chamber_alias="pn_c1_8",
+                    node_alias="pnnode4_16",
+                    time_s=0.05,
+                    amesim_port_1_enthalpy_flow_w=-4248.6,
+                    amesim_port_1_mass_flow_g_s=419.9,
+                    amesim_node_port_2_enthalpy_flow_w=-18664.2,
+                    amesim_node_port_2_mass_flow_g_s=438.3,
+                    amesim_observed_port_enthalpy_sum_w=-22912.8,
+                    python_port_1_internal_energy_flow_w=-378000.0,
+                    python_port_2_internal_energy_flow_w=398000.0,
+                    python_current_energy_derivative_w=20000.0,
+                    python_direct_enthalpy_port_1_flow_w=-629000.0,
+                    python_direct_enthalpy_port_2_flow_w=663000.0,
+                    python_direct_enthalpy_energy_derivative_w=34000.0,
+                    python_direct_minus_current_energy_derivative_w=14000.0,
+                    python_p4_node_port_2_enthalpy_flow_w=640000.0,
+                    python_p4_node_port_2_mass_flow_g_s=430.0,
+                    python_p4_node_port_2_connected_h_j_kg=1488372.1,
+                    python_p4_node_port_2_connected_temperature_k=286.6,
+                    python_p4_node_port_2_connected_h_delta_to_line_h_j_kg=-34000.0,
+                    python_p4_node_port_2_counterfactual_dtemp_k_s=-5100.0,
+                    reference_temperature_k=298.15,
+                    amesim_port_1_reference_enthalpy_estimate_w=-4467.0,
+                    amesim_port_1_reference_enthalpy_error_w=-218.4,
+                    amesim_candidate_pn2pipefr_dm2i_g_s=-419.9,
+                    amesim_candidate_pn2pipefr_dh2i_w=4467.0,
+                    amesim_candidate_storage_enthalpy_sum_w=-14197.2,
+                    amesim_candidate_pn2pipefr_dtemp_k_s=-5480.0,
+                    amesim_candidate_pn2pipefr_dtemp_residual_k_s=-5420.3,
+                    python_reference_port_1_enthalpy_flow_w=-10760.0,
+                    python_reference_node_port_2_enthalpy_flow_w=-18500.0,
+                    python_reference_observed_port_enthalpy_sum_w=-29260.0,
+                    python_reference_port_1_to_amesim_error_w=-6511.4,
+                    python_reference_node_port_2_to_amesim_error_w=164.2,
+                    python_reference_sum_to_amesim_error_w=-6347.2,
+                    python_current_line_temperature_derivative_k_s=0.0,
+                    amesim_line_temperature_derivative_fd_k_s=277.3,
+                    amesim_pn2vol2_dtemp_node_plus_dh1_k_s=-9000.0,
+                    amesim_pn2vol2_dtemp_node_minus_dh1_k_s=-5600.0,
+                    python_reference_pn2vol2_dtemp_node_minus_dh1_k_s=-11200.0,
+                ),
+            ),
+            pnch012_energy_equation_diagnostics=(
+                TestMqlPnch012EnergyEquationDiagnostic(
+                    chamber_alias="pn_c1_8",
+                    piston_alias="pn_brp2_8",
+                    line_alias="pneumatic_69",
+                    time_s=0.05,
+                    amesim_port_1_enthalpy_flow_w=-4248.6,
+                    amesim_port_1_mass_flow_g_s=419.9,
+                    amesim_chamber_mass_derivative_fd_g_s=370.8,
+                    amesim_chamber_mass_derivative_residual_g_s=49.1,
+                    amesim_chamber_temperature_derivative_fd_k_s=5134.6,
+                    amesim_chamber_volume_rate_fd_m3_s=0.001004,
+                    amesim_piston_volume_rate_m3_s=2.5e-8,
+                    amesim_heat_flow_w=-0.0,
+                    amesim_boundary_work_fd_volume_w=-289.5,
+                    amesim_pn2vol_reference_dtemp_fd_volume_k_s=8090.0,
+                    amesim_pn2vol_reference_dtemp_fd_volume_residual_k_s=2955.4,
+                    amesim_pn2vol_reference_dtemp_piston_volume_k_s=8107.0,
+                    amesim_pn2vol_reference_dtemp_piston_volume_residual_k_s=2972.4,
+                    python_port_a_mass_flow_kg_s=-0.42,
+                    python_current_energy_derivative_w=-28000.0,
+                    python_current_chamber_temperature_derivative_k_s=-7800.0,
+                    reference_temperature_k=298.15,
+                ),
+            ),
         )
         diagnostic = TestMqlPnvoEventWindowDiagnostic(
             orifice_alias="pn_morifice_1",
@@ -331,6 +403,36 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase):
         self.assertIn("python_cm=0.0147", summary)
         self.assertIn("cm_ratio=0.93", summary)
         self.assertIn("amesim_gasvel=885.9", summary)
+        self.assertIn("energy_flow@pneumatic_69", summary)
+        self.assertIn("amesim_dh1=-4248.6", summary)
+        self.assertIn("amesim_node_dh2=-18664.2", summary)
+        self.assertIn("direct_h_dU_dt=34000.0", summary)
+        self.assertIn("python_p4_node_dh2=640000.0", summary)
+        self.assertIn("python_p4_node_dm2=430.0", summary)
+        self.assertIn("python_p4_node_h2=1488372.1", summary)
+        self.assertIn("python_p4_node_t2=286.6", summary)
+        self.assertIn("python_p4_node_h2_delta=-34000.0", summary)
+        self.assertIn("python_p4_node_counterfactual_dT=-5100.0", summary)
+        self.assertIn("ref_dh1_est=-4467.0", summary)
+        self.assertIn("ref_dh1_error=-218.4", summary)
+        self.assertIn("candidate_dm2i=-419.9", summary)
+        self.assertIn("candidate_dh2i=4467.0", summary)
+        self.assertIn("candidate_storage_sdh=-14197.2", summary)
+        self.assertIn("candidate_pn2pipefr_dT=-5480.0", summary)
+        self.assertIn("python_ref_dh1=-10760.0", summary)
+        self.assertIn("python_ref_node_dh2=-18500.0", summary)
+        self.assertIn("python_ref_sdh=-29260.0", summary)
+        self.assertIn("python_ref_sdh_error=-6347.2", summary)
+        self.assertIn("python_current_dT=0.0", summary)
+        self.assertIn("amesim_fd_dT=277.3", summary)
+        self.assertIn("amesim_pn2vol2_dT_node_plus_dh1=-9000.0", summary)
+        self.assertIn("python_ref_pn2vol2_dT_node_minus_dh1=-11200.0", summary)
+        self.assertIn("pnch012_energy@pn_c1_8", summary)
+        self.assertIn("amesim_fd_dmgas=370.8", summary)
+        self.assertIn("amesim_dmgas_residual=49.1", summary)
+        self.assertIn("amesim_fd_dvol=0.001004", summary)
+        self.assertIn("amesim_pn2vol_ref_dT_fdvol=8090.0", summary)
+        self.assertIn("python_current_dT=-7800.0", summary)
         self.assertIn("Final comparison:", summary)
         self.assertTrue(summary.endswith("\n"))