diff --git a/PythonModels/scripts/run_test_mql_full_state_comparison.py b/PythonModels/scripts/run_test_mql_full_state_comparison.py index 42b86fe..fea6d1b 100644 --- a/PythonModels/scripts/run_test_mql_full_state_comparison.py +++ b/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( diff --git a/PythonModels/systems/test_mql_closure.py b/PythonModels/systems/test_mql_closure.py index a236e40..4aecbae 100644 --- a/PythonModels/systems/test_mql_closure.py +++ b/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, ) diff --git a/tests/test_run_test_mql_full_state_comparison.py b/tests/test_run_test_mql_full_state_comparison.py index ce1cd42..e3428bf 100644 --- a/tests/test_run_test_mql_full_state_comparison.py +++ b/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"))