接入AMESim help文档与PNVO流量诊断

PythonModels/scripts/run_test_mql_full_state_comparison.py

@@ -3,7 +3,7 @@ from __future__ import annotations
 import argparse
 from dataclasses import dataclass, field
 from datetime import UTC, datetime
-from math import nextafter
+from math import nextafter, sqrt
 from pathlib import Path
 
 from PythonModels.core.solver import SolveIVPConfig, integrate_ode
@@ -27,6 +27,7 @@ from PythonModels.systems.test_mql import (
     TestMqlSimulationResult,
     TestMqlSystem,
 )
+from PythonModels.systems.test_mql_pneumatic import AMESIM_REFERENCE_PRESSURE_PA
 
 
 DEFAULT_FULL_STATE_COMPARISON_DATA_PATHS = (
@@ -265,6 +266,45 @@ class TestMqlPnvoEventWindowSegmentDiagnostic:
     message: str
 
 
+@dataclass(frozen=True)
+class TestMqlPnl0001PressureLossCalibrationDiagnostic:
+    line_alias: str
+    chamber_alias: str
+    time_s: float
+    amesim_cm: float
+    amesim_dm1_g_s: float
+    amesim_mass_flow_magnitude_kg_s: float
+    amesim_line_gauge_pressure_pa: float
+    amesim_chamber_gauge_pressure_pa: float
+    amesim_line_temperature_k: float
+    amesim_pressure_drop_pa: float
+    current_darcy_pressure_drop_pa: float
+    pressure_drop_multiplier: float
+    amesim_linear_conductance_kg_s_sqrt_k_per_pa: float
+    python_linear_conductance_kg_s_sqrt_k_per_pa: float
+    python_to_amesim_linear_conductance_ratio: float
+
+
+@dataclass(frozen=True)
+class TestMqlPnvoFlowParameterDiagnostic:
+    orifice_alias: str
+    time_s: float
+    amesim_cm: float
+    amesim_dm2_g_s: float
+    amesim_opening: float
+    amesim_gas_velocity_m_s: float
+    python_opening: float
+    python_flow_coefficient: float
+    python_effective_area_m2: float
+    python_line_pressure_pa: float
+    python_boundary_pressure_pa: float
+    python_upstream_pressure_pa: float
+    python_upstream_temperature_k: float
+    python_mass_flow_kg_s: float
+    python_cm: float
+    python_to_amesim_cm_ratio: float
+
+
 @dataclass(frozen=True)
 class TestMqlPnvoEventWindowSampleDiagnostic:
     time_s: float
@@ -274,6 +314,12 @@ class TestMqlPnvoEventWindowSampleDiagnostic:
     pnl0001_rhs_diagnostics: tuple[TestMqlPnl0001LineRhsDiagnostic, ...] = field(
         default_factory=tuple
     )
+    pnl0001_pressure_loss_diagnostics: tuple[
+        TestMqlPnl0001PressureLossCalibrationDiagnostic, ...
+    ] = field(default_factory=tuple)
+    pnvo_flow_parameter_diagnostics: tuple[
+        TestMqlPnvoFlowParameterDiagnostic, ...
+    ] = field(default_factory=tuple)
 
     def abs_error(self, data_path: str) -> float:
         return abs(
@@ -509,6 +555,166 @@ def format_test_mql_pnvo_event_boundary_summary(
     return "\n".join(lines) + "\n"
 
 
+def _pnl0001_pressure_loss_calibration_diagnostic(
+    *,
+    closure: object,
+    amesim_results: AmesimResults,
+    python_values_by_data_path: dict[str, float],
+    line_alias: str,
+    chamber_alias: str,
+    time_s: float,
+) -> TestMqlPnl0001PressureLossCalibrationDiagnostic:
+    if line_alias != "pneumatic_69" or chamber_alias != "pn_c1_8":
+        raise KeyError(f"Unsupported PNL0001 pressure-loss diagnostic: {line_alias}")
+    line = closure.pneumatic_closure.components.p4_port3_remote_primary_line
+
+    def amesim_value(data_path: str) -> float:
+        return interpolate_series_value(
+            amesim_results.times,
+            amesim_results.series(data_path),
+            time_s,
+        )
+
+    amesim_cm = amesim_value(f"cm@{line_alias}")
+    amesim_dm1_g_s = amesim_value(f"dm1@{line_alias}")
+    amesim_line_gauge_pressure_pa = amesim_value(f"p2@{line_alias}")
+    amesim_chamber_gauge_pressure_pa = amesim_value(f"press@{chamber_alias}")
+    amesim_line_temperature_k = amesim_value(f"t2@{line_alias}")
+    mass_flow_magnitude_kg_s = abs(amesim_dm1_g_s) * 1.0e-3
+    amesim_pressure_drop_pa = abs(
+        amesim_line_gauge_pressure_pa - amesim_chamber_gauge_pressure_pa
+    )
+    current_darcy_pressure_drop_pa = abs(
+        line.darcy_pressure_drop_for_state(
+            mass_flow_kg_s=mass_flow_magnitude_kg_s,
+            pressure_pa=(
+                amesim_line_gauge_pressure_pa + AMESIM_REFERENCE_PRESSURE_PA
+            ),
+            temperature_k=amesim_line_temperature_k,
+        )
+    )
+    if current_darcy_pressure_drop_pa > 0.0:
+        pressure_drop_multiplier = (
+            amesim_pressure_drop_pa / current_darcy_pressure_drop_pa
+        )
+    else:
+        pressure_drop_multiplier = (
+            float("inf") if amesim_pressure_drop_pa > 0.0 else 1.0
+        )
+    amesim_linear_conductance = _pnl0001_linear_conductance(
+        dm1_g_s=amesim_dm1_g_s,
+        temperature_k=amesim_line_temperature_k,
+        pressure_drop_pa=amesim_pressure_drop_pa,
+    )
+    python_pressure_drop_pa = abs(
+        python_values_by_data_path[f"p2@{line_alias}"]
+        - python_values_by_data_path[f"press@{chamber_alias}"]
+    )
+    python_linear_conductance = _pnl0001_linear_conductance(
+        dm1_g_s=python_values_by_data_path[f"dm1@{line_alias}"],
+        temperature_k=python_values_by_data_path[f"t2@{line_alias}"],
+        pressure_drop_pa=python_pressure_drop_pa,
+    )
+    if amesim_linear_conductance > 0.0:
+        conductance_ratio = python_linear_conductance / amesim_linear_conductance
+    else:
+        conductance_ratio = (
+            float("inf") if python_linear_conductance > 0.0 else 1.0
+        )
+    return TestMqlPnl0001PressureLossCalibrationDiagnostic(
+        line_alias=line_alias,
+        chamber_alias=chamber_alias,
+        time_s=time_s,
+        amesim_cm=amesim_cm,
+        amesim_dm1_g_s=amesim_dm1_g_s,
+        amesim_mass_flow_magnitude_kg_s=mass_flow_magnitude_kg_s,
+        amesim_line_gauge_pressure_pa=amesim_line_gauge_pressure_pa,
+        amesim_chamber_gauge_pressure_pa=amesim_chamber_gauge_pressure_pa,
+        amesim_line_temperature_k=amesim_line_temperature_k,
+        amesim_pressure_drop_pa=amesim_pressure_drop_pa,
+        current_darcy_pressure_drop_pa=current_darcy_pressure_drop_pa,
+        pressure_drop_multiplier=pressure_drop_multiplier,
+        amesim_linear_conductance_kg_s_sqrt_k_per_pa=(
+            amesim_linear_conductance
+        ),
+        python_linear_conductance_kg_s_sqrt_k_per_pa=(
+            python_linear_conductance
+        ),
+        python_to_amesim_linear_conductance_ratio=conductance_ratio,
+    )
+
+
+def _pnl0001_linear_conductance(
+    *,
+    dm1_g_s: float,
+    temperature_k: float,
+    pressure_drop_pa: float,
+) -> float:
+    if temperature_k <= 0.0:
+        raise ValueError("temperature_k must be positive")
+    if pressure_drop_pa <= 0.0:
+        return 0.0
+    return abs(dm1_g_s) * 1.0e-3 * sqrt(temperature_k) / pressure_drop_pa
+
+
+def _pnvo_flow_parameter_diagnostic(
+    *,
+    closure: object,
+    amesim_results: AmesimResults,
+    state_vector: list[float],
+    orifice_alias: str,
+    time_s: float,
+) -> TestMqlPnvoFlowParameterDiagnostic:
+    if orifice_alias != "pn_morifice_1":
+        raise KeyError(f"Unsupported PNVO flow parameter diagnostic: {orifice_alias}")
+    snapshot = closure.snapshot_at(time_s, state_vector).pneumatic
+    orifice = closure.pneumatic_closure.components.p4_port3_remote_orifice
+    if orifice.name != orifice_alias:
+        raise KeyError(f"Unexpected PNVO diagnostic orifice: {orifice.name}")
+    line_properties = snapshot.p4_port3_remote_orifice_line_port_2
+    boundary_properties = snapshot.p4_port3_remote_primary_line
+    upstream_properties = (
+        line_properties if line_properties.p >= boundary_properties.p else boundary_properties
+    )
+    python_mass_flow_kg_s = snapshot.p4_port3_remote_orifice_to_node_flow
+    denominator = (
+        orifice.flow_coefficient * orifice.effective_area * upstream_properties.p
+    )
+    python_cm = (
+        abs(python_mass_flow_kg_s) * sqrt(upstream_properties.T) / denominator
+        if denominator > 0.0 and upstream_properties.T > 0.0
+        else 0.0
+    )
+
+    def amesim_value(data_path: str) -> float:
+        return interpolate_series_value(
+            amesim_results.times,
+            amesim_results.series(data_path),
+            time_s,
+        )
+
+    amesim_cm = amesim_value(f"cm@{orifice_alias}")
+    cm_ratio = python_cm / amesim_cm if amesim_cm != 0.0 else 0.0
+    return TestMqlPnvoFlowParameterDiagnostic(
+        orifice_alias=orifice_alias,
+        time_s=time_s,
+        amesim_cm=amesim_cm,
+        amesim_dm2_g_s=amesim_value(f"dm2@{orifice_alias}"),
+        amesim_opening=amesim_value(f"xv@{orifice_alias}"),
+        amesim_gas_velocity_m_s=amesim_value(f"gasvel@{orifice_alias}"),
+        python_opening=orifice.opening,
+        python_flow_coefficient=orifice.flow_coefficient,
+        python_effective_area_m2=orifice.effective_area,
+        python_line_pressure_pa=line_properties.p,
+        python_boundary_pressure_pa=boundary_properties.p,
+        python_upstream_pressure_pa=upstream_properties.p,
+        python_upstream_temperature_k=upstream_properties.T,
+        python_mass_flow_kg_s=python_mass_flow_kg_s,
+        python_cm=python_cm,
+        python_to_amesim_cm_ratio=cm_ratio,
+    )
+
+
 def run_test_mql_pnvo_event_window_diagnostic(
     config: TestMqlFullStateComparisonScriptConfig | None = None,
     *,
@@ -679,6 +885,25 @@ def run_test_mql_pnvo_event_window_diagnostic(
                 time_s=sample_time,
             ),
         )
+        pnl0001_pressure_loss_diagnostics = (
+            _pnl0001_pressure_loss_calibration_diagnostic(
+                closure=closure,
+                amesim_results=amesim_results,
+                python_values_by_data_path=python_sample_values,
+                line_alias="pneumatic_69",
+                chamber_alias="pn_c1_8",
+                time_s=sample_time,
+            ),
+        )
+        pnvo_flow_parameter_diagnostics = (
+            _pnvo_flow_parameter_diagnostic(
+                closure=closure,
+                amesim_results=amesim_results,
+                state_vector=state_vector_by_sample_time[sample_time],
+                orifice_alias=orifice_alias,
+                time_s=sample_time,
+            ),
+        )
         sample_diagnostics.append(
             TestMqlPnvoEventWindowSampleDiagnostic(
                 time_s=sample_time,
@@ -686,6 +911,8 @@ def run_test_mql_pnvo_event_window_diagnostic(
                 python_values_by_data_path=python_sample_values,
                 amesim_values_by_data_path=amesim_sample_values,
                 pnl0001_rhs_diagnostics=pnl0001_rhs_diagnostics,
+                pnl0001_pressure_loss_diagnostics=pnl0001_pressure_loss_diagnostics,
+                pnvo_flow_parameter_diagnostics=pnvo_flow_parameter_diagnostics,
             )
         )
     python_values = closure.data_path_values(
@@ -749,6 +976,35 @@ def format_test_mql_pnvo_event_window_summary(
                     f"dm_dt={rhs.mass_derivative_kg_s}, "
                     f"dU_dt={rhs.energy_derivative_w}"
                 )
+            for pressure_loss in sample.pnl0001_pressure_loss_diagnostics:
+                lines.append(
+                    f"    - pressure_loss@{pressure_loss.line_alias}: "
+                    f"cm={pressure_loss.amesim_cm}, "
+                    f"amesim_dm1={pressure_loss.amesim_dm1_g_s}, "
+                    f"amesim_dp={pressure_loss.amesim_pressure_drop_pa}, "
+                    f"darcy_dp={pressure_loss.current_darcy_pressure_drop_pa}, "
+                    f"dp_multiplier={pressure_loss.pressure_drop_multiplier}, "
+                    f"amesim_linear_k="
+                    f"{pressure_loss.amesim_linear_conductance_kg_s_sqrt_k_per_pa}, "
+                    f"python_linear_k="
+                    f"{pressure_loss.python_linear_conductance_kg_s_sqrt_k_per_pa}, "
+                    f"linear_k_ratio="
+                    f"{pressure_loss.python_to_amesim_linear_conductance_ratio}"
+                )
+            for flow_parameter in sample.pnvo_flow_parameter_diagnostics:
+                lines.append(
+                    f"    - flow_parameter@{flow_parameter.orifice_alias}: "
+                    f"amesim_cm={flow_parameter.amesim_cm}, "
+                    f"python_cm={flow_parameter.python_cm}, "
+                    f"cm_ratio={flow_parameter.python_to_amesim_cm_ratio}, "
+                    f"amesim_dm2={flow_parameter.amesim_dm2_g_s}, "
+                    f"python_m={flow_parameter.python_mass_flow_kg_s}, "
+                    f"opening={flow_parameter.python_opening}, "
+                    f"effective_area={flow_parameter.python_effective_area_m2}, "
+                    f"upstream_p={flow_parameter.python_upstream_pressure_pa}, "
+                    f"upstream_t={flow_parameter.python_upstream_temperature_k}, "
+                    f"amesim_gasvel={flow_parameter.amesim_gas_velocity_m_s}"
+                )
     lines.append("Final comparison:")
     for data_path in diagnostic.data_paths:
         lines.append(