lujingze/SystemSimulationApp
1
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

接入PNL0001事件窗口RHS诊断

Browse Source
This commit is contained in:
huojiarong committed 2026-07-22 04:21:38 +00:00
1 parent 2120901909
commit f0f8f40c7a
4 files changed
+597 -11
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PythonModels/scripts/run_test_mql_full_state_comparison.py
+314 -2
View File
@@ -6,7 +6,7 @@ from datetime import UTC, datetime
from math import nextafter
from pathlib import Path
from PythonModels.core.solver import SolveIVPConfig
from PythonModels.core.solver import SolveIVPConfig, integrate_ode
from PythonModels.reporting.amesim_results import AmesimResults, load_test_mql_amesim_results
from PythonModels.reporting.test_mql_comparison import (
TestMqlComparisonResult,
@@ -22,7 +22,11 @@ from PythonModels.reporting.test_mql_output_validation import (
compare_validated_test_mql_output,
validate_test_mql_output,
)
from PythonModels.systems.test_mql import TestMqlSimulationResult, TestMqlSystem
from PythonModels.systems.test_mql import (
TestMqlPnl0001LineRhsDiagnostic,
TestMqlSimulationResult,
TestMqlSystem,
)
DEFAULT_FULL_STATE_COMPARISON_DATA_PATHS = (
@@ -248,6 +252,54 @@ class TestMqlPnvoEventBoundaryDiagnostic:
)
@dataclass(frozen=True)
class TestMqlPnvoEventWindowSegmentDiagnostic:
t_start: float
t_stop: float
method: str
rtol: float
atol: float
max_step: float
rhs_evaluations: int
success: bool
message: str
@dataclass(frozen=True)
class TestMqlPnvoEventWindowSampleDiagnostic:
time_s: float
data_paths: tuple[str, ...]
python_values_by_data_path: dict[str, float]
amesim_values_by_data_path: dict[str, float]
pnl0001_rhs_diagnostics: tuple[TestMqlPnl0001LineRhsDiagnostic, ...] = field(
default_factory=tuple
)
def abs_error(self, data_path: str) -> float:
return abs(
self.python_values_by_data_path[data_path]
- self.amesim_values_by_data_path[data_path]
)
@dataclass(frozen=True)
class TestMqlPnvoEventWindowDiagnostic:
orifice_alias: str
event_time_s: float
final_time_s: float
data_paths: tuple[str, ...]
segment_diagnostics: tuple[TestMqlPnvoEventWindowSegmentDiagnostic, ...]
sample_diagnostics: tuple[TestMqlPnvoEventWindowSampleDiagnostic, ...]
python_values_by_data_path: dict[str, float]
amesim_values_by_data_path: dict[str, float]
def abs_error(self, data_path: str) -> float:
return abs(
self.python_values_by_data_path[data_path]
- self.amesim_values_by_data_path[data_path]
)
@dataclass(frozen=True)
class TestMqlFullStateComparisonPathConfig:
archive_path: Path = field(
@@ -457,6 +509,257 @@ def format_test_mql_pnvo_event_boundary_summary(
return "\n".join(lines) + "\n"
def run_test_mql_pnvo_event_window_diagnostic(
config: TestMqlFullStateComparisonScriptConfig | None = None,
*,
orifice_alias: str = "pn_morifice_1",
) -> TestMqlPnvoEventWindowDiagnostic:
config = config or TestMqlFullStateComparisonScriptConfig()
system = TestMqlSystem(archive_path=config.paths.archive_path)
control = system.pneumatic_assembly.variable_orifice_controls[orifice_alias]
event_time_s = control.step.step_time_s
final_time_s = 0.05
spec = system.discover_pneumatic_branch_topology().chamber_segment_specs[0]
closure_kwargs = {
"inlet_node_pressure_pa": config.execution.inlet_node_pressure_pa,
"resistance_boundary_pressure_pa": (
config.execution.resistance_boundary_pressure_pa
),
"inlet_node_temperature_k": config.execution.inlet_node_temperature_k,
"resistance_boundary_temperature_k": (
config.execution.resistance_boundary_temperature_k
),
}
closure = system.full_state_closure_from_spec(spec, **closure_kwargs)
state_vector = closure.initial_state_vector()
solver_template = config.execution.solver
segments = (
SolveIVPConfig(
t_start=solver_template.t_start,
t_stop=nextafter(event_time_s, 0.0),
method=solver_template.method,
rtol=solver_template.rtol,
atol=solver_template.atol,
max_step=solver_template.max_step,
),
SolveIVPConfig(
t_start=event_time_s,
t_stop=0.040001,
method="Radau",
rtol=solver_template.rtol,
atol=solver_template.atol,
max_step=1.0e-7,
first_step=1.0e-10,
),
SolveIVPConfig(
t_start=0.040001,
t_stop=0.04001,
method="Radau",
rtol=solver_template.rtol,
atol=solver_template.atol,
max_step=1.0e-6,
),
SolveIVPConfig(
t_start=0.04001,
t_stop=0.0401,
method="Radau",
rtol=solver_template.rtol,
atol=solver_template.atol,
max_step=1.0e-5,
),
SolveIVPConfig(
t_start=0.0401,
t_stop=0.041,
method="Radau",
rtol=solver_template.rtol,
atol=solver_template.atol,
max_step=1.0e-5,
),
SolveIVPConfig(
t_start=0.041,
t_stop=0.048,
method="BDF",
rtol=solver_template.rtol,
atol=solver_template.atol,
max_step=1.0e-5,
),
SolveIVPConfig(
t_start=0.048,
t_stop=final_time_s,
method="BDF",
rtol=1.0e-5,
atol=1.0e-8,
max_step=1.0e-5,
),
)
data_paths = (
"press@pn_c1_8",
"temp@pn_c1_8",
"vol@pn_c1_8",
"mgas1@pn_c1_8",
"p2@pneumatic_69",
"t2@pneumatic_69",
"mgas@pneumatic_69",
"re@pneumatic_69",
"v@pneumatic_69",
"ff@pneumatic_69",
"dm1@pneumatic_69",
f"xv@{orifice_alias}",
f"dm2@{orifice_alias}",
)
sample_times = (event_time_s, 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:
rhs_evaluations = 0
segment_sample_times = tuple(
sample_time
for sample_time in sample_times
if segment.t_start <= sample_time <= segment.t_stop
)
t_eval = tuple(
sorted({segment.t_start, segment.t_stop, *segment_sample_times})
)
def rhs(time_s, values):
nonlocal rhs_evaluations
rhs_evaluations += 1
return closure.rhs_at(time_s, values)
solution = integrate_ode(
rhs=rhs,
initial_state=state_vector,
config=segment,
t_eval=list(t_eval),
)
segment_diagnostics.append(
TestMqlPnvoEventWindowSegmentDiagnostic(
t_start=segment.t_start,
t_stop=segment.t_stop,
method=segment.method,
rtol=segment.rtol,
atol=segment.atol,
max_step=segment.max_step,
rhs_evaluations=rhs_evaluations,
success=bool(solution.success),
message=str(solution.message),
)
)
for sample_time in segment_sample_times:
if sample_time in solution.t:
sample_index = list(solution.t).index(sample_time)
state_vector_by_sample_time[sample_time] = [
row[sample_index] for row in solution.y
]
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)
sample_diagnostics = []
for sample_time in sample_times:
if sample_time not in state_vector_by_sample_time:
continue
python_sample_values = closure.data_path_values(
time_s=sample_time,
state_vector=state_vector_by_sample_time[sample_time],
data_paths=data_paths,
)
amesim_sample_values = {
data_path: interpolate_series_value(
amesim_results.times,
amesim_results.series(data_path),
sample_time,
)
for data_path in data_paths
}
pnl0001_rhs_diagnostics = (
closure.pnl0001_line_rhs_diagnostic(
line_alias="pneumatic_69",
state_vector=state_vector_by_sample_time[sample_time],
time_s=sample_time,
),
)
sample_diagnostics.append(
TestMqlPnvoEventWindowSampleDiagnostic(
time_s=sample_time,
data_paths=data_paths,
python_values_by_data_path=python_sample_values,
amesim_values_by_data_path=amesim_sample_values,
pnl0001_rhs_diagnostics=pnl0001_rhs_diagnostics,
)
)
python_values = closure.data_path_values(
time_s=final_time_s,
state_vector=state_vector,
data_paths=data_paths,
)
amesim_values = {
data_path: interpolate_series_value(
amesim_results.times,
amesim_results.series(data_path),
final_time_s,
)
for data_path in data_paths
}
return TestMqlPnvoEventWindowDiagnostic(
orifice_alias=orifice_alias,
event_time_s=event_time_s,
final_time_s=final_time_s,
data_paths=data_paths,
segment_diagnostics=tuple(segment_diagnostics),
sample_diagnostics=tuple(sample_diagnostics),
python_values_by_data_path=python_values,
amesim_values_by_data_path=amesim_values,
)
def format_test_mql_pnvo_event_window_summary(
diagnostic: TestMqlPnvoEventWindowDiagnostic,
) -> str:
lines = [
"Model: test_mql",
f"Mode: PNVO event window diagnostic ({diagnostic.orifice_alias})",
f"Event time: {diagnostic.event_time_s}",
f"Final time: {diagnostic.final_time_s}",
"Segments:",
]
for segment in diagnostic.segment_diagnostics:
lines.append(
f" - {segment.t_start} -> {segment.t_stop}: "
f"method={segment.method}, rtol={segment.rtol}, "
f"atol={segment.atol}, max_step={segment.max_step}, "
f"rhs={segment.rhs_evaluations}, success={segment.success}"
)
if diagnostic.sample_diagnostics:
lines.append("Sample comparisons:")
for sample in diagnostic.sample_diagnostics:
lines.append(f" t={sample.time_s}")
for data_path in sample.data_paths:
lines.append(
f" - {data_path}: "
f"python={sample.python_values_by_data_path[data_path]}, "
f"amesim={sample.amesim_values_by_data_path[data_path]}, "
f"abs_error={sample.abs_error(data_path)}"
)
for rhs in sample.pnl0001_rhs_diagnostics:
lines.append(
f" - rhs@{rhs.line_alias}: "
f"chamber_to_line={rhs.chamber_to_line_flow_kg_s}, "
f"node_to_line={rhs.node_to_line_flow_kg_s}, "
f"dm_dt={rhs.mass_derivative_kg_s}, "
f"dU_dt={rhs.energy_derivative_w}"
)
lines.append("Final comparison:")
for data_path in diagnostic.data_paths:
lines.append(
f" - {data_path}: "
f"python={diagnostic.python_values_by_data_path[data_path]}, "
f"amesim={diagnostic.amesim_values_by_data_path[data_path]}, "
f"abs_error={diagnostic.abs_error(data_path)}"
)
return "\n".join(lines) + "\n"
def format_test_mql_full_state_comparison_summary(
run: TestMqlFullStateComparisonRun,
) -> str:
@@ -609,11 +912,20 @@ def main() -> None:
action="store_true",
help="compare the STEP0 left-state/right-opening boundary at t=0.04 s",
)
parser.add_argument(
"--pnvo-event-window",
action="store_true",
help="run the segmented PNVO opening window through the t=0.05 s save point",
)
args = parser.parse_args()
if args.pnvo_event_boundary:
diagnostic = run_test_mql_pnvo_event_boundary_diagnostic()
print(format_test_mql_pnvo_event_boundary_summary(diagnostic), end="")
return
if args.pnvo_event_window:
diagnostic = run_test_mql_pnvo_event_window_diagnostic()
print(format_test_mql_pnvo_event_window_summary(diagnostic), end="")
return
run, output_dir = run_test_mql_full_state_comparison()
print(format_test_mql_full_state_comparison_summary(run), end="")
PythonModels/systems/test_mql.py
+121 -9
View File
@@ -3874,6 +3874,23 @@ class TestMqlVariableChamberRhsDiagnostic:
energy_derivative_w: float
@dataclass(frozen=True)
class TestMqlPnl0001LineRhsDiagnostic:
line_alias: str
chamber_alias: str
line_pressure_pa: float
line_temperature_k: float
chamber_pressure_pa: float
chamber_temperature_k: float
chamber_to_line_flow_kg_s: float
node_to_line_flow_kg_s: float
mass_derivative_kg_s: float
port_1_energy_flow_w: float
port_2_energy_flow_w: float
thermal_energy_flow_w: float
energy_derivative_w: float
_PISTON_FORCE_BINDINGS = (
("mass_friction_endstops_10", "pn_brp2_8", "p4_port3_remote_primary_chamber"),
("mass_friction_endstops_11", "pn_brp2_9", "p4_primary_chamber"),
@@ -3977,8 +3994,13 @@ _PRIMARY_CHAMBER_DIAGNOSTIC_FIELDS = {
}
_PNL0001_MASS_FLOW_DIAGNOSTIC_FIELDS = {
"pneumatic_69": "p4_port3_remote_chamber_to_line_flow",
_PNL0001_LINE_DIAGNOSTIC_FIELDS = {
"pneumatic_69": (
"p4_port3_remote_primary_line",
"p4_port3_remote_primary_chamber",
"p4_port3_remote_chamber_to_line_flow",
"p4_port3_remote_node_to_primary_line_flow",
),
}
_VARIABLE_ORIFICE_MASS_FLOW_DIAGNOSTIC_FIELDS = {
@@ -4146,6 +4168,65 @@ class TestMqlFullStateClosure:
),
)
def pnl0001_line_rhs_diagnostic(
self,
*,
line_alias: str,
state_vector: list[float],
time_s: float = 0.0,
) -> TestMqlPnl0001LineRhsDiagnostic:
if line_alias not in _PNL0001_LINE_DIAGNOSTIC_FIELDS:
raise KeyError(line_alias)
(
line_field,
chamber_field,
chamber_flow_field,
node_flow_field,
) = _PNL0001_LINE_DIAGNOSTIC_FIELDS[line_alias]
chamber_alias = _PNCH012_ALIAS_BY_SNAPSHOT_FIELD[chamber_field]
snapshot = self.snapshot_at(time_s, state_vector)
line = getattr(self.pneumatic_closure.components, line_field)
line_properties = getattr(snapshot.pneumatic, line_field)
chamber_properties = getattr(snapshot.pneumatic, chamber_field)
chamber_to_line_flow = getattr(snapshot.pneumatic, chamber_flow_field)
node_to_line_flow = getattr(snapshot.pneumatic, node_flow_field)
line_derivative = line.derivatives_from_connections(
port_1_m_flow=chamber_to_line_flow,
connected_h_1=chamber_properties.h,
port_2_m_flow=node_to_line_flow,
connected_h_2=line_properties.h,
)
port_1_inlet_u = (
chamber_properties.h / line.gas.gamma
if chamber_to_line_flow > 0.0
else line_properties.u
)
port_2_inlet_u = (
line_properties.h / line.gas.gamma
if node_to_line_flow > 0.0
else line_properties.u
)
thermal_energy_flow = (
line.heat_transfer_coefficient
* line.heat_transfer_area
* (line.external_temperature - line_properties.T)
)
return TestMqlPnl0001LineRhsDiagnostic(
line_alias=line_alias,
chamber_alias=chamber_alias,
line_pressure_pa=line_properties.p,
line_temperature_k=line_properties.T,
chamber_pressure_pa=chamber_properties.p,
chamber_temperature_k=chamber_properties.T,
chamber_to_line_flow_kg_s=chamber_to_line_flow,
node_to_line_flow_kg_s=node_to_line_flow,
mass_derivative_kg_s=line_derivative.m,
port_1_energy_flow_w=chamber_to_line_flow * port_1_inlet_u,
port_2_energy_flow_w=node_to_line_flow * port_2_inlet_u,
thermal_energy_flow_w=thermal_energy_flow,
energy_derivative_w=line_derivative.U,
)
def data_path_values(
self,
*,
@@ -4173,14 +4254,43 @@ class TestMqlFullStateClosure:
values_by_data_path = {}
for data_path in selected_paths:
signal, alias = self._split_data_path(data_path)
if alias in _PNL0001_MASS_FLOW_DIAGNOSTIC_FIELDS:
if signal != "dm1":
if alias in _PNL0001_LINE_DIAGNOSTIC_FIELDS:
(
line_field,
_chamber_field,
chamber_flow_field,
_node_flow_field,
) = _PNL0001_LINE_DIAGNOSTIC_FIELDS[alias]
if signal == "dm1":
canonical_flow_kg_s = getattr(snapshot.pneumatic, chamber_flow_field)
values_by_data_path[data_path] = -canonical_flow_kg_s * 1.0e3
elif signal == "p2":
line_properties = getattr(snapshot.pneumatic, line_field)
values_by_data_path[data_path] = pressure_to_amesim_gauge_pa(
line_properties.p
)
elif signal == "t2":
line_properties = getattr(snapshot.pneumatic, line_field)
values_by_data_path[data_path] = line_properties.T
elif signal == "mgas":
line = getattr(self.pneumatic_closure.components, line_field)
values_by_data_path[data_path] = line.gas_mass_g()
elif signal in {"re", "v", "ff"}:
line = getattr(self.pneumatic_closure.components, line_field)
line_properties = getattr(snapshot.pneumatic, line_field)
canonical_flow_kg_s = getattr(snapshot.pneumatic, chamber_flow_field)
diagnostics = line.diagnostics(
mass_flow_kg_s=-canonical_flow_kg_s,
temperature_k=line_properties.T,
)
if signal == "re":
values_by_data_path[data_path] = diagnostics.reynolds_number
elif signal == "v":
values_by_data_path[data_path] = diagnostics.gas_velocity_m_s
else:
values_by_data_path[data_path] = diagnostics.friction_factor
else:
raise KeyError(data_path)
canonical_flow_kg_s = getattr(
snapshot.pneumatic,
_PNL0001_MASS_FLOW_DIAGNOSTIC_FIELDS[alias],
)
values_by_data_path[data_path] = -canonical_flow_kg_s * 1.0e3
elif alias in self.pneumatic_assembly.variable_orifices:
orifice = self.pneumatic_assembly.variable_orifices[alias]
if signal == "xv":
@@ -4207,6 +4317,8 @@ class TestMqlFullStateClosure:
values_by_data_path[data_path] = chamber_properties.T
elif signal == "vol":
values_by_data_path[data_path] = chamber.volume_cm3()
elif signal in {"mgas", "mgas1"}:
values_by_data_path[data_path] = chamber.gas_mass_g()
else:
raise KeyError(data_path)
elif alias in self.mechanical_closure.assembly.pistons: