From 1193b2dceec5d5d2f6d9a7df9027a8f1d5df46a9 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Fri, 17 Jul 2026 03:28:36 +0000 Subject: [PATCH] =?UTF-8?q?=E5=AE=9E=E7=8E=B0test=5Fmql=20PNL0001=E7=AE=A1?= =?UTF-8?q?=E8=B7=AF=E5=8A=A8=E6=80=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../components/amesim_pneumatic_line.py | 284 ++++++++++++++++++ PythonModels/systems/test_mql.py | 90 ++++++ PythonModels/systems/test_mql_closure.py | 176 +++++++++++ .../systems/test_mql_line_parameters.py | 162 ++++++++++ .../systems/test_mql_pneumatic_lines.py | 50 +++ tests/test_amesim_pnl0001_pipe.py | 71 +++++ tests/test_test_mql_pnl0001.py | 82 +++++ tests/test_test_mql_pnl0001_segment.py | 80 +++++ 8 files changed, 995 insertions(+) create mode 100644 PythonModels/components/amesim_pneumatic_line.py create mode 100644 PythonModels/systems/test_mql_line_parameters.py create mode 100644 PythonModels/systems/test_mql_pneumatic_lines.py create mode 100644 tests/test_amesim_pnl0001_pipe.py create mode 100644 tests/test_test_mql_pnl0001.py create mode 100644 tests/test_test_mql_pnl0001_segment.py diff --git a/PythonModels/components/amesim_pneumatic_line.py b/PythonModels/components/amesim_pneumatic_line.py new file mode 100644 index 0000000..08f6a35 --- /dev/null +++ b/PythonModels/components/amesim_pneumatic_line.py @@ -0,0 +1,284 @@ +from __future__ import annotations + +from dataclasses import dataclass +from math import log10, pi + +from PythonModels.components.amesim_pneumatic import ( + HELIUM_PNEUMATIC_GAS, + AmesimPneumaticGas, + diameter_mm_to_area_m2, +) +from PythonModels.core.base import DynamicComponent +from PythonModels.core.medium import ThermodynamicProperties +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +@dataclass(frozen=True) +class AmesimPnl0001Diagnostics: + mass_flow_kg_s: float + reynolds_number: float + gas_velocity_m_s: float + friction_factor: float + pressure_drop_pa: float + + +class AmesimPnl0001Pipe(DynamicComponent): + """Physical first-pass implementation of AMESim ``PNL0001`` (C-R). + + Port 2 owns the lumped gas storage. Port 1 is connected through a Darcy + resistance. Both connection mass flows use the PythonModels convention: + positive values enter the pipe storage. + + AMESim's proprietary pressure-loss calibration is not available in the + archive. This implementation therefore uses an explicit Darcy-Weisbach + law while preserving the real geometry, state count, mass/energy balance, + heat-transfer parameter, and observable diagnostics. + """ + + def __init__( + self, + name: str, + *, + diameter_mm: float, + length_m: float, + relative_roughness: float, + polytropic_constant: float = 1.35, + heat_transfer_coefficient: float = 0.0, + external_temperature_k: float = 293.15, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + p0: float = 101_325.0, + T0: float = 293.15, + ) -> None: + if diameter_mm <= 0.0: + raise ValueError("diameter_mm must be positive") + if length_m <= 0.0: + raise ValueError("length_m must be positive") + if relative_roughness < 0.0: + raise ValueError("relative_roughness must be non-negative") + if polytropic_constant <= 0.0: + raise ValueError("polytropic_constant must be positive") + if heat_transfer_coefficient < 0.0: + raise ValueError("heat_transfer_coefficient must be non-negative") + if external_temperature_k <= 0.0: + raise ValueError("external_temperature_k must be positive") + + super().__init__(name=name) + self.diameter = diameter_mm * 1.0e-3 + self.length = length_m + self.relative_roughness = relative_roughness + self.polytropic_constant = polytropic_constant + self.heat_transfer_coefficient = heat_transfer_coefficient + self.external_temperature = external_temperature_k + self.gas = gas + self.area = diameter_mm_to_area_m2(diameter_mm) + self.volume = self.area * self.length + self.heat_transfer_area = pi * self.diameter * self.length + + rho0 = gas.density(p0, T0) + mass0 = rho0 * self.volume + self.state = VolumeState( + m=mass0, + U=mass0 * gas.specific_internal_energy(T0), + ) + self.port_1 = PortState() + self.port_2 = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + if self.state.m <= 0.0: + raise ValueError("pipe mass must stay positive") + temperature = self.gas.temperature_from_internal_energy( + self.state.U / self.state.m + ) + density = self.state.m / self.volume + pressure = self.gas.pressure(density, temperature) + properties = ThermodynamicProperties( + p=pressure, + T=temperature, + rho=density, + u=self.state.U / self.state.m, + h=self.gas.specific_enthalpy(temperature), + ) + self.port_2.p = pressure + self.port_2.h_outflow = properties.h + return properties + + def gas_mass_g(self) -> float: + return self.state.m * 1.0e3 + + def resistance_mass_flow( + self, + *, + port_1_pressure_pa: float, + port_1_temperature_k: float, + ) -> float: + """Return mass flow from port 1 into the port-2 storage in kg/s.""" + if port_1_pressure_pa <= 0.0: + raise ValueError("port_1_pressure_pa must be positive") + if port_1_temperature_k <= 0.0: + raise ValueError("port_1_temperature_k must be positive") + + internal = self.properties() + pressure_difference = port_1_pressure_pa - internal.p + if pressure_difference == 0.0: + return 0.0 + upstream_pressure = max(port_1_pressure_pa, internal.p) + upstream_temperature = ( + port_1_temperature_k if pressure_difference > 0.0 else internal.T + ) + density = self.gas.density(upstream_pressure, upstream_temperature) + magnitude = self._mass_flow_for_pressure_drop( + abs(pressure_difference), + density=density, + temperature=upstream_temperature, + ) + return magnitude if pressure_difference > 0.0 else -magnitude + + def diagnostics( + self, + *, + mass_flow_kg_s: float, + temperature_k: float | None = None, + ) -> AmesimPnl0001Diagnostics: + properties = self.properties() + temperature = temperature_k or properties.T + reynolds = self._reynolds_number(mass_flow_kg_s, temperature) + friction_factor = self._friction_factor(reynolds) + velocity = mass_flow_kg_s / (properties.rho * self.area) + pressure_drop = self._darcy_pressure_drop( + mass_flow_kg_s, + density=properties.rho, + temperature=temperature, + ) + return AmesimPnl0001Diagnostics( + mass_flow_kg_s=mass_flow_kg_s, + reynolds_number=reynolds, + gas_velocity_m_s=velocity, + friction_factor=friction_factor, + pressure_drop_pa=pressure_drop, + ) + + def derivatives_from_connections( + self, + *, + port_1_m_flow: float, + connected_h_1: float, + port_2_m_flow: float, + connected_h_2: float, + ) -> VolumeState: + internal = self.properties() + inlet_h_1 = self.connection_inlet_enthalpy( + port_m_flow=port_1_m_flow, + connected_h=connected_h_1, + internal_h=internal.h, + ) + inlet_h_2 = self.connection_inlet_enthalpy( + port_m_flow=port_2_m_flow, + connected_h=connected_h_2, + internal_h=internal.h, + ) + heat_flow = ( + self.heat_transfer_coefficient + * self.heat_transfer_area + * (self.external_temperature - internal.T) + ) + return VolumeState( + m=port_1_m_flow + port_2_m_flow, + U=port_1_m_flow * inlet_h_1 + port_2_m_flow * inlet_h_2 + heat_flow, + ) + + def _mass_flow_for_pressure_drop( + self, + pressure_drop_pa: float, + *, + density: float, + temperature: float, + ) -> float: + if pressure_drop_pa <= 0.0: + return 0.0 + upper = 1.0e-9 + while self._darcy_pressure_drop( + upper, + density=density, + temperature=temperature, + ) < pressure_drop_pa: + upper *= 10.0 + if upper > 1.0e3: + raise ValueError("unable to bracket PNL0001 resistance flow") + lower = 0.0 + for _ in range(80): + middle = 0.5 * (lower + upper) + if self._darcy_pressure_drop( + middle, + density=density, + temperature=temperature, + ) < pressure_drop_pa: + lower = middle + else: + upper = middle + return 0.5 * (lower + upper) + + def _darcy_pressure_drop( + self, + mass_flow_kg_s: float, + *, + density: float, + temperature: float, + ) -> float: + if mass_flow_kg_s == 0.0: + return 0.0 + reynolds = self._reynolds_number(mass_flow_kg_s, temperature) + friction_factor = self._friction_factor(reynolds) + velocity = mass_flow_kg_s / (density * self.area) + magnitude = ( + friction_factor + * (self.length / self.diameter) + * density + * velocity + * velocity + / 2.0 + ) + return magnitude if mass_flow_kg_s > 0.0 else -magnitude + + def _reynolds_number(self, mass_flow_kg_s: float, temperature: float) -> float: + viscosity = helium_dynamic_viscosity(temperature) + return 4.0 * abs(mass_flow_kg_s) / (pi * self.diameter * viscosity) + + def _friction_factor(self, reynolds_number: float) -> float: + if reynolds_number <= 0.0: + return 64_000_000.0 + laminar = 64.0 / reynolds_number + if reynolds_number <= 2_300.0: + return laminar + turbulent = 1.0 / ( + -1.8 + * log10( + (self.relative_roughness / 3.7) ** 1.11 + + 6.9 / reynolds_number + ) + ) ** 2 + if reynolds_number >= 4_000.0: + return turbulent + fraction = (reynolds_number - 2_300.0) / 1_700.0 + return laminar + fraction * (turbulent - laminar) + + +def helium_dynamic_viscosity(temperature_k: float) -> float: + """Sutherland approximation centered on the test_mql initial condition.""" + if temperature_k <= 0.0: + raise ValueError("temperature_k must be positive") + reference_temperature = 293.15 + reference_viscosity = 2.0e-5 + sutherland_constant = 79.4 + return ( + reference_viscosity + * (temperature_k / reference_temperature) ** 1.5 + * (reference_temperature + sutherland_constant) + / (temperature_k + sutherland_constant) + ) diff --git a/PythonModels/systems/test_mql.py b/PythonModels/systems/test_mql.py index 55dd178..42eb77c 100644 --- a/PythonModels/systems/test_mql.py +++ b/PythonModels/systems/test_mql.py @@ -3859,6 +3859,7 @@ class TestMqlSystem: def __init__(self, archive_path: Path | None = None) -> None: self.archive_path = archive_path or Path(__file__).resolve().parents[2] / AMESIM_ARCHIVE_RELATIVE_PATH self.network = SimulationNetwork(name=MODEL_NAME) + self.pnl0001_assembly = self._build_pnl0001_assembly() self.pneumatic_assembly = self._build_pneumatic_assembly() pneumatic_components = self._pneumatic_components_by_alias() for spec in COMPONENT_SPECS: @@ -3881,6 +3882,13 @@ class TestMqlSystem: return build_test_mql_pneumatic_assembly() + def _build_pnl0001_assembly(self): + from PythonModels.systems.test_mql_pneumatic_lines import ( + build_test_mql_pnl0001_assembly, + ) + + return build_test_mql_pnl0001_assembly(self.archive_path) + def _pneumatic_components_by_alias(self) -> dict[str, Component]: return { **self.pneumatic_assembly.fixed_chambers, @@ -3893,6 +3901,10 @@ class TestMqlSystem: def typed_pneumatic_component_count(self) -> int: return len(self._pneumatic_components_by_alias()) + @property + def typed_pnl0001_line_count(self) -> int: + return len(self.pnl0001_assembly.lines) + def pneumatic_state_vector(self) -> list[float]: return self.network.initial_state_vector() @@ -4071,6 +4083,84 @@ class TestMqlSystem: t_eval=t_eval, ) + def pnl0001_chamber_segment_closure_from_spec( + self, + spec, + *, + inlet_node_pressure_pa: float, + outlet_pressure_pa: float, + inlet_node_temperature_k: float = 293.15, + outlet_temperature_k: float = 293.15, + ): + """Insert the topology-derived inlet PNL0001 into a chamber segment.""" + from PythonModels.components.amesim_pneumatic import ( + AmesimPneumaticOrifice, + AmesimPneumaticVolume, + ) + from PythonModels.systems.test_mql_closure import ( + TestMqlPneumaticBoundaryCondition, + TestMqlPnl0001ChamberSegmentClosure, + TestMqlPnl0001ChamberSegmentComponents, + ) + + inlet_line = self.pnl0001_assembly.lines[spec.inlet_line_alias] + volume = self.network.components[spec.volume_alias] + inlet_orifice = self.network.components[spec.inlet_orifice_alias] + outlet_orifice = self.network.components[spec.outlet_orifice_alias] + if not isinstance(volume, AmesimPneumaticVolume): + raise TypeError(f"{spec.volume_alias} is not an AMESim pneumatic volume") + if not isinstance(inlet_orifice, AmesimPneumaticOrifice): + raise TypeError( + f"{spec.inlet_orifice_alias} is not an AMESim pneumatic orifice" + ) + if not isinstance(outlet_orifice, AmesimPneumaticOrifice): + raise TypeError( + f"{spec.outlet_orifice_alias} is not an AMESim pneumatic orifice" + ) + return TestMqlPnl0001ChamberSegmentClosure( + components=TestMqlPnl0001ChamberSegmentComponents( + inlet_line=inlet_line, + volume=volume, + inlet_orifice=inlet_orifice, + outlet_orifice=outlet_orifice, + spec=spec, + ), + inlet_node=TestMqlPneumaticBoundaryCondition( + pressure_pa=inlet_node_pressure_pa, + temperature_k=inlet_node_temperature_k, + ), + outlet_boundary=TestMqlPneumaticBoundaryCondition( + pressure_pa=outlet_pressure_pa, + temperature_k=outlet_temperature_k, + ), + ) + + def simulate_pnl0001_chamber_segment_from_spec( + self, + spec, + *, + inlet_node_pressure_pa: float, + outlet_pressure_pa: float, + inlet_node_temperature_k: float = 293.15, + outlet_temperature_k: float = 293.15, + config: SolveIVPConfig | None = None, + t_eval: list[float] | None = None, + ): + closure = self.pnl0001_chamber_segment_closure_from_spec( + spec, + inlet_node_pressure_pa=inlet_node_pressure_pa, + outlet_pressure_pa=outlet_pressure_pa, + inlet_node_temperature_k=inlet_node_temperature_k, + outlet_temperature_k=outlet_temperature_k, + ) + run_config = config or SolveIVPConfig(t_stop=1.0e-4, max_step=1.0e-5) + return integrate_ode( + rhs=lambda t, state: closure.rhs(state), + initial_state=closure.initial_state_vector(), + config=run_config, + t_eval=t_eval, + ) + def pneumatic_branch_closure_from_spec(self, spec): return self.pneumatic_branch_closure( name=spec.name, diff --git a/PythonModels/systems/test_mql_closure.py b/PythonModels/systems/test_mql_closure.py index 9a5651b..af948ef 100644 --- a/PythonModels/systems/test_mql_closure.py +++ b/PythonModels/systems/test_mql_closure.py @@ -8,6 +8,7 @@ from PythonModels.components.amesim_pneumatic import ( AmesimPneumaticOrifice, AmesimPneumaticVolume, ) +from PythonModels.components.amesim_pneumatic_line import AmesimPnl0001Pipe from PythonModels.core.medium import ThermodynamicProperties from PythonModels.core.ports import PortState from PythonModels.core.state import VolumeState @@ -98,6 +99,26 @@ class TestMqlPneumaticChamberSegmentSnapshot: outlet_flow: float +@dataclass(frozen=True) +class TestMqlPnl0001ChamberSegmentComponents: + inlet_line: AmesimPnl0001Pipe + volume: AmesimPneumaticVolume + inlet_orifice: AmesimPneumaticOrifice + outlet_orifice: AmesimPneumaticOrifice + spec: TestMqlPneumaticChamberSegmentSpec + + +@dataclass(frozen=True) +class TestMqlPnl0001ChamberSegmentSnapshot: + inlet_line: ThermodynamicProperties + chamber: ThermodynamicProperties + inlet_node: ThermodynamicProperties + outlet_boundary: ThermodynamicProperties + node_to_line_flow: float + line_to_chamber_flow: float + outlet_flow: float + + @dataclass(frozen=True) class TestMqlPneumaticBranchComponents: name: str @@ -379,3 +400,158 @@ class TestMqlPneumaticChamberSegmentClosure: internal_h=snapshot.chamber.h, ) return derivative.as_vector() + + +class TestMqlPnl0001ChamberSegmentClosure: + """Fixed chamber segment with the topology-derived inlet PNL0001 state. + + The inlet line has a closed causal boundary here: port-1 pressure and + temperature come from the PN3 node boundary, while port-2 flow comes from + the fixed orifice. The outlet line remains a boundary until its three-line + PN3 node balance is assembled. + """ + + def __init__( + self, + *, + components: TestMqlPnl0001ChamberSegmentComponents, + inlet_node: TestMqlPneumaticBoundaryCondition, + outlet_boundary: TestMqlPneumaticBoundaryCondition, + ) -> None: + self.components = components + self.inlet_node = inlet_node + self.outlet_boundary = outlet_boundary + + def initial_state_vector(self) -> list[float]: + return [ + *self.components.inlet_line.get_state_vector(), + *self.components.volume.get_state_vector(), + ] + + def apply_state_vector(self, values: list[float]) -> None: + if len(values) != 4: + raise ValueError("PNL0001/chamber segment state vector requires four values") + self.components.inlet_line.set_state_vector(values[:2]) + self.components.volume.set_state_vector(values[2:]) + + def snapshot( + self, + state_vector: list[float] | None = None, + ) -> TestMqlPnl0001ChamberSegmentSnapshot: + if state_vector is not None: + self.apply_state_vector(state_vector) + line = self.components.inlet_line.properties() + chamber = self.components.volume.properties() + inlet_node = self.inlet_node.properties(self.components.inlet_line.gas) + outlet_boundary = self.outlet_boundary.properties(self.components.volume.gas) + node_to_line_flow = self.components.inlet_line.resistance_mass_flow( + port_1_pressure_pa=inlet_node.p, + port_1_temperature_k=inlet_node.T, + ) + inlet_temperature = line.T if line.p >= chamber.p else chamber.T + line_to_chamber_flow = self.components.inlet_orifice.mass_flow( + line.p, + chamber.p, + inlet_temperature, + ) + outlet_temperature = ( + chamber.T if chamber.p >= outlet_boundary.p else outlet_boundary.T + ) + outlet_flow = self.components.outlet_orifice.mass_flow( + chamber.p, + outlet_boundary.p, + outlet_temperature, + ) + snapshot = TestMqlPnl0001ChamberSegmentSnapshot( + inlet_line=line, + chamber=chamber, + inlet_node=inlet_node, + outlet_boundary=outlet_boundary, + node_to_line_flow=node_to_line_flow, + line_to_chamber_flow=line_to_chamber_flow, + outlet_flow=outlet_flow, + ) + self._write_port_states(snapshot) + return snapshot + + @staticmethod + def _port( + component: AmesimPneumaticVolume | AmesimPneumaticOrifice, + port_name: str, + ) -> PortState: + return TestMqlPneumaticChamberSegmentClosure._port(component, port_name) + + def _write_port_states( + self, + snapshot: TestMqlPnl0001ChamberSegmentSnapshot, + ) -> None: + spec = self.components.spec + line = self.components.inlet_line + chamber = self.components.volume + inlet_orifice = self.components.inlet_orifice + outlet_orifice = self.components.outlet_orifice + + line.port_1.p = snapshot.inlet_node.p + line.port_1.m_flow = snapshot.node_to_line_flow + line.port_1.h_outflow = snapshot.inlet_line.h + line.port_2.p = snapshot.inlet_line.p + line.port_2.m_flow = -snapshot.line_to_chamber_flow + + inlet_boundary_port = self._port( + inlet_orifice, + spec.inlet_orifice_boundary_port, + ) + inlet_volume_port = self._port(inlet_orifice, spec.inlet_orifice_volume_port) + inlet_boundary_port.p = snapshot.inlet_line.p + inlet_boundary_port.m_flow = snapshot.line_to_chamber_flow + inlet_boundary_port.h_outflow = snapshot.inlet_line.h + inlet_volume_port.p = snapshot.chamber.p + inlet_volume_port.m_flow = -snapshot.line_to_chamber_flow + inlet_volume_port.h_outflow = snapshot.chamber.h + + chamber_inlet_port = self._port(chamber, spec.volume_inlet_port) + chamber_outlet_port = self._port(chamber, spec.volume_outlet_port) + chamber_inlet_port.m_flow = snapshot.line_to_chamber_flow + chamber_outlet_port.m_flow = -snapshot.outlet_flow + + outlet_volume_port = self._port( + outlet_orifice, + spec.outlet_orifice_volume_port, + ) + outlet_boundary_port = self._port( + outlet_orifice, + spec.outlet_orifice_boundary_port, + ) + outlet_volume_port.p = snapshot.chamber.p + outlet_volume_port.m_flow = snapshot.outlet_flow + outlet_volume_port.h_outflow = snapshot.chamber.h + outlet_boundary_port.p = snapshot.outlet_boundary.p + outlet_boundary_port.m_flow = -snapshot.outlet_flow + outlet_boundary_port.h_outflow = snapshot.outlet_boundary.h + + def rhs(self, state_vector: list[float]) -> list[float]: + snapshot = self.snapshot(state_vector) + line_derivative = self.components.inlet_line.derivatives_from_connections( + port_1_m_flow=snapshot.node_to_line_flow, + connected_h_1=snapshot.inlet_node.h, + port_2_m_flow=-snapshot.line_to_chamber_flow, + connected_h_2=snapshot.chamber.h, + ) + chamber = self.components.volume + spec = self.components.spec + chamber_derivative = chamber.derivatives_from_two_connections( + port_a_m_flow=self._port(chamber, "port_1").m_flow, + connected_h_a=( + snapshot.inlet_line.h + if spec.volume_inlet_port == "port_1" + else snapshot.outlet_boundary.h + ), + port_b_m_flow=self._port(chamber, "port_2").m_flow, + connected_h_b=( + snapshot.inlet_line.h + if spec.volume_inlet_port == "port_2" + else snapshot.outlet_boundary.h + ), + internal_h=snapshot.chamber.h, + ) + return [*line_derivative.as_vector(), *chamber_derivative.as_vector()] diff --git a/PythonModels/systems/test_mql_line_parameters.py b/PythonModels/systems/test_mql_line_parameters.py new file mode 100644 index 0000000..fef40c7 --- /dev/null +++ b/PythonModels/systems/test_mql_line_parameters.py @@ -0,0 +1,162 @@ +from __future__ import annotations + +import re +import tarfile +from dataclasses import dataclass +from pathlib import Path + +from PythonModels.systems.test_mql import CONNECTION_SPECS, GLOBAL_PARAMETERS +from PythonModels.systems.test_mql_config import resolve_numeric_expression + + +AMESIM_REFERENCE_PRESSURE_PA = 101_300.0 + + +@dataclass(frozen=True) +class TestMqlPnl0001Spec: + alias: str + source_component: str + source_port: str + target_component: str + target_port: str + diameter_mm: float + length_m: float + relative_roughness: float + polytropic_constant: float + heat_transfer_coefficient: float + external_temperature_k: float + gas_type_index: int + mode: int + initial_temperature_k: float + initial_gauge_pressure_pa: float + + @property + def initial_absolute_pressure_pa(self) -> float: + return self.initial_gauge_pressure_pa + AMESIM_REFERENCE_PRESSURE_PA + + +def load_test_mql_pnl0001_specs( + archive_path: str | Path, + *, + cir_member: str = "test_mql_.cir", +) -> tuple[TestMqlPnl0001Spec, ...]: + """Load resolved PNL0001 geometry and initial states from the AMESim source.""" + with tarfile.open(archive_path) as archive: + cir_file = archive.extractfile(cir_member) + if cir_file is None: + raise ValueError(f"Missing AMESim circuit member: {cir_member}") + cir_text = cir_file.read().decode("latin1") + + numeric_globals = { + name: value + for name, expression in GLOBAL_PARAMETERS.items() + if (value := resolve_numeric_expression(expression, {})) is not None + } + connections = { + str(connection["alias"]): connection + for connection in CONNECTION_SPECS + if connection["submodel"] == "PNL0001" + } + specs = [] + for block in re.findall(r".*?", cir_text, flags=re.DOTALL): + if _optional_text(block, "SUB_NAME") != "PNL0001": + continue + alias = _required_text(block, "ALIAS") + connection = connections.get(alias) + if connection is None: + raise ValueError(f"PNL0001 line {alias!r} is absent from CONNECTION_SPECS") + real_parameters = _parameter_expressions(block, "RPARAM") + integer_parameters = _parameter_expressions(block, "IPARAM") + state_values = _evar_values(block) + specs.append( + TestMqlPnl0001Spec( + alias=alias, + source_component=str(connection["source_component"]), + source_port=str(connection["source_port"]), + target_component=str(connection["target_component"]), + target_port=str(connection["target_port"]), + diameter_mm=_required_numeric( + alias, "diam", real_parameters, numeric_globals + ), + length_m=_required_numeric(alias, "le", real_parameters, numeric_globals), + relative_roughness=_required_numeric( + alias, "rr", real_parameters, numeric_globals + ), + polytropic_constant=_required_numeric( + alias, "k", real_parameters, numeric_globals + ), + heat_transfer_coefficient=_required_numeric( + alias, "kth", real_parameters, numeric_globals + ), + external_temperature_k=_required_numeric( + alias, "extemp", real_parameters, numeric_globals + ), + gas_type_index=int( + _required_numeric(alias, "gi", integer_parameters, numeric_globals) + ), + mode=int( + _required_numeric(alias, "mode", integer_parameters, numeric_globals) + ), + initial_temperature_k=_required_numeric( + alias, "t2", state_values, numeric_globals + ), + initial_gauge_pressure_pa=_required_numeric( + alias, "p2", state_values, numeric_globals + ), + ) + ) + if set(connections) != {spec.alias for spec in specs}: + missing = sorted(set(connections) - {spec.alias for spec in specs}) + raise ValueError(f"Missing PNL0001 parameter blocks: {missing}") + return tuple(specs) + + +def _parameter_expressions(block: str, tag_name: str) -> dict[str, str]: + parameters = {} + for parameter_block in re.findall( + rf"<{tag_name}>.*?", + block, + flags=re.DOTALL, + ): + parameters[_required_text(parameter_block, "VARNAME")] = _required_text( + parameter_block, + "VALUE", + ) + return parameters + + +def _evar_values(block: str) -> dict[str, str]: + values = {} + for variable_block in re.findall(r".*?", block, flags=re.DOTALL): + value = _optional_text(variable_block, "VALUE") + if value: + values[_required_text(variable_block, "VARNAME")] = value + return values + + +def _required_numeric( + alias: str, + name: str, + expressions: dict[str, str], + variables: dict[str, float], +) -> float: + if name not in expressions: + raise ValueError(f"Missing {name!r} on PNL0001 line {alias!r}") + value = resolve_numeric_expression(expressions[name], variables) + if value is None: + raise ValueError( + f"Cannot resolve {name!r}={expressions[name]!r} on PNL0001 line {alias!r}" + ) + return value + + +def _required_text(block: str, tag_name: str) -> str: + value = _optional_text(block, tag_name) + if value is None: + raise ValueError(f"Missing AMESim circuit element: {tag_name}") + return value + + +def _optional_text(block: str, tag_name: str) -> str | None: + match = re.search(rf"<{tag_name}>(.*?)", block, flags=re.DOTALL) + return match.group(1).strip() if match is not None else None diff --git a/PythonModels/systems/test_mql_pneumatic_lines.py b/PythonModels/systems/test_mql_pneumatic_lines.py new file mode 100644 index 0000000..9d4db9a --- /dev/null +++ b/PythonModels/systems/test_mql_pneumatic_lines.py @@ -0,0 +1,50 @@ +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +from PythonModels.components.amesim_pneumatic import ( + HELIUM_PNEUMATIC_GAS, + AmesimPneumaticGas, +) +from PythonModels.components.amesim_pneumatic_line import AmesimPnl0001Pipe +from PythonModels.systems.test_mql_line_parameters import ( + TestMqlPnl0001Spec, + load_test_mql_pnl0001_specs, +) + + +@dataclass(frozen=True) +class TestMqlPnl0001Assembly: + specs: tuple[TestMqlPnl0001Spec, ...] + lines: dict[str, AmesimPnl0001Pipe] + + def spec(self, alias: str) -> TestMqlPnl0001Spec: + for spec in self.specs: + if spec.alias == alias: + return spec + raise KeyError(alias) + + +def build_test_mql_pnl0001_assembly( + archive_path: str | Path, + *, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, +) -> TestMqlPnl0001Assembly: + specs = load_test_mql_pnl0001_specs(archive_path) + lines = { + spec.alias: AmesimPnl0001Pipe( + name=spec.alias, + diameter_mm=spec.diameter_mm, + length_m=spec.length_m, + relative_roughness=spec.relative_roughness, + polytropic_constant=spec.polytropic_constant, + heat_transfer_coefficient=spec.heat_transfer_coefficient, + external_temperature_k=spec.external_temperature_k, + gas=gas, + p0=spec.initial_absolute_pressure_pa, + T0=spec.initial_temperature_k, + ) + for spec in specs + } + return TestMqlPnl0001Assembly(specs=specs, lines=lines) diff --git a/tests/test_amesim_pnl0001_pipe.py b/tests/test_amesim_pnl0001_pipe.py new file mode 100644 index 0000000..2c6e78f --- /dev/null +++ b/tests/test_amesim_pnl0001_pipe.py @@ -0,0 +1,71 @@ +from __future__ import annotations + +import unittest + +from PythonModels.components.amesim_pneumatic_line import AmesimPnl0001Pipe + + +class AmesimPnl0001PipeTests(unittest.TestCase): + def setUp(self) -> None: + self.pipe = AmesimPnl0001Pipe( + name="pneumatic_96", + diameter_mm=14.0, + length_m=1.0, + relative_roughness=0.045 / 14.0, + p0=15.3e6, + T0=293.15, + ) + + def test_initial_state_uses_real_pipe_volume_and_two_states(self) -> None: + properties = self.pipe.properties() + + self.assertAlmostEqual(self.pipe.volume, 1.539380400258999e-4) + self.assertEqual(len(self.pipe.get_state_vector()), 2) + self.assertAlmostEqual(properties.p, 15.3e6, delta=1.0e-5) + self.assertAlmostEqual(properties.T, 293.15) + self.assertAlmostEqual(self.pipe.gas_mass_g(), 3.716965219188, places=10) + + def test_resistance_flow_follows_pressure_gradient(self) -> None: + forward = self.pipe.resistance_mass_flow( + port_1_pressure_pa=15.31e6, + port_1_temperature_k=293.15, + ) + reverse = self.pipe.resistance_mass_flow( + port_1_pressure_pa=15.29e6, + port_1_temperature_k=293.15, + ) + + self.assertGreater(forward, 0.0) + self.assertLess(reverse, 0.0) + self.assertAlmostEqual(abs(forward), abs(reverse), delta=abs(forward) * 0.01) + + def test_diagnostics_reproduce_laminar_friction_contract(self) -> None: + diagnostics = self.pipe.diagnostics(mass_flow_kg_s=8.90603914774626e-6) + + self.assertGreater(diagnostics.reynolds_number, 40.0) + self.assertLess(diagnostics.reynolds_number, 45.0) + self.assertAlmostEqual( + diagnostics.friction_factor, + 64.0 / diagnostics.reynolds_number, + ) + self.assertGreater(diagnostics.gas_velocity_m_s, 0.0) + self.assertGreater(diagnostics.pressure_drop_pa, 0.0) + + def test_connection_derivative_preserves_mass_and_stream_direction(self) -> None: + internal = self.pipe.properties() + derivative = self.pipe.derivatives_from_connections( + port_1_m_flow=0.2, + connected_h_1=internal.h + 1000.0, + port_2_m_flow=-0.1, + connected_h_2=internal.h - 1000.0, + ) + + self.assertAlmostEqual(derivative.m, 0.1) + self.assertAlmostEqual( + derivative.U, + 0.2 * (internal.h + 1000.0) - 0.1 * internal.h, + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_test_mql_pnl0001.py b/tests/test_test_mql_pnl0001.py new file mode 100644 index 0000000..0d2730c --- /dev/null +++ b/tests/test_test_mql_pnl0001.py @@ -0,0 +1,82 @@ +from __future__ import annotations + +import unittest +from pathlib import Path + +from PythonModels.reporting.amesim_results import load_test_mql_amesim_results +from PythonModels.systems.test_mql_line_parameters import load_test_mql_pnl0001_specs +from PythonModels.systems.test_mql_pneumatic_lines import ( + build_test_mql_pnl0001_assembly, +) + + +REPO_ROOT = Path(__file__).resolve().parents[1] +TEST_MQL_AME = REPO_ROOT / "AmesimModels" / "test_mql.ame" + + +class TestMqlPnl0001Tests(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + cls.specs = load_test_mql_pnl0001_specs(TEST_MQL_AME) + cls.assembly = build_test_mql_pnl0001_assembly(TEST_MQL_AME) + cls.results = load_test_mql_amesim_results(TEST_MQL_AME) + + def test_loads_all_real_pnl0001_parameters_from_cir(self) -> None: + self.assertEqual(len(self.specs), 20) + spec = self.assembly.spec("pneumatic_96") + + self.assertEqual(spec.source_component, "pn_node3_8") + self.assertEqual(spec.target_component, "pn_orifice_18") + self.assertEqual(spec.diameter_mm, 14.0) + self.assertEqual(spec.length_m, 1.0) + self.assertAlmostEqual(spec.relative_roughness, 0.045 / 14.0) + self.assertEqual(spec.polytropic_constant, 1.35) + self.assertEqual(spec.heat_transfer_coefficient, 0.0) + self.assertEqual(spec.external_temperature_k, 293.15) + self.assertEqual(spec.gas_type_index, 1) + self.assertEqual(spec.mode, 2) + self.assertAlmostEqual(spec.initial_gauge_pressure_pa, 15_198_700.0) + self.assertAlmostEqual(spec.initial_absolute_pressure_pa, 15_300_000.0) + + def test_builds_twenty_two_state_physical_line_components(self) -> None: + self.assertEqual(len(self.assembly.lines), 20) + self.assertTrue( + all(len(line.get_state_vector()) == 2 for line in self.assembly.lines.values()) + ) + + def test_pneumatic_96_initial_observables_match_amesim_baseline(self) -> None: + pipe = self.assembly.lines["pneumatic_96"] + + self.assertAlmostEqual( + pipe.properties().p - 101_300.0, + self.results.series("p2@pneumatic_96")[0], + delta=1.0e-5, + ) + self.assertAlmostEqual( + pipe.properties().T, + self.results.series("t2@pneumatic_96")[0], + ) + self.assertAlmostEqual( + pipe.gas_mass_g(), + self.results.series("mgas@pneumatic_96")[0], + delta=0.005, + ) + + def test_baseline_mass_balance_identifies_port_two_input_sign(self) -> None: + mass = self.results.series("mgas@pneumatic_96") + dm1 = self.results.series("dm1@pneumatic_96") + dm2_peer = self.results.series("dm2@pn_orifice_18") + index = 500 + finite_difference = (mass[index] - mass[index - 1]) / ( + self.results.times[index] - self.results.times[index - 1] + ) + + self.assertAlmostEqual( + finite_difference, + dm2_peer[index] - dm1[index], + delta=3.0e-5, + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_test_mql_pnl0001_segment.py b/tests/test_test_mql_pnl0001_segment.py new file mode 100644 index 0000000..3e251ee --- /dev/null +++ b/tests/test_test_mql_pnl0001_segment.py @@ -0,0 +1,80 @@ +from __future__ import annotations + +import unittest + +from PythonModels.core.solver import SolveIVPConfig +from PythonModels.systems.test_mql import TestMqlSystem + + +class TestMqlPnl0001ChamberSegmentTests(unittest.TestCase): + def setUp(self) -> None: + self.system = TestMqlSystem() + self.spec = self.system.discover_pneumatic_branch_topology().chamber_segment_specs[0] + + def test_system_assembles_all_twenty_physical_pnl0001_lines(self) -> None: + self.assertEqual(self.system.typed_pnl0001_line_count, 20) + self.assertIn("pneumatic_96", self.system.pnl0001_assembly.lines) + + def test_closure_inserts_topology_derived_inlet_line(self) -> None: + closure = self.system.pnl0001_chamber_segment_closure_from_spec( + self.spec, + inlet_node_pressure_pa=16.0e6, + outlet_pressure_pa=15.3e6, + ) + + snapshot = closure.snapshot() + rhs = closure.rhs(closure.initial_state_vector()) + + self.assertEqual(closure.components.inlet_line.name, "pneumatic_96") + self.assertEqual(len(closure.initial_state_vector()), 4) + self.assertEqual(len(rhs), 4) + self.assertGreater(snapshot.node_to_line_flow, 0.0) + self.assertAlmostEqual(snapshot.line_to_chamber_flow, 0.0, delta=1.0e-7) + self.assertAlmostEqual(snapshot.outlet_flow, 0.0, delta=1.0e-7) + self.assertAlmostEqual(rhs[0] + rhs[2], snapshot.node_to_line_flow) + + def test_internal_line_orifice_flow_cancels_from_total_mass_balance(self) -> None: + closure = self.system.pnl0001_chamber_segment_closure_from_spec( + self.spec, + inlet_node_pressure_pa=15.3e6, + outlet_pressure_pa=14.0e6, + ) + state = closure.initial_state_vector() + state[0] *= 1.01 + state[1] *= 1.01 + + snapshot = closure.snapshot(state) + rhs = closure.rhs(state) + + self.assertGreater(snapshot.line_to_chamber_flow, 0.0) + self.assertAlmostEqual( + rhs[0] + rhs[2], + snapshot.node_to_line_flow - snapshot.outlet_flow, + delta=1.0e-12, + ) + self.assertAlmostEqual( + closure.components.inlet_line.port_2.m_flow, + -snapshot.line_to_chamber_flow, + ) + + def test_simulates_four_state_line_chamber_segment(self) -> None: + solution = self.system.simulate_pnl0001_chamber_segment_from_spec( + self.spec, + inlet_node_pressure_pa=15.31e6, + outlet_pressure_pa=15.29e6, + config=SolveIVPConfig( + t_start=0.0, + t_stop=1.0e-7, + max_step=1.0e-8, + ), + t_eval=[0.0, 5.0e-8, 1.0e-7], + ) + + self.assertTrue(solution.success) + self.assertEqual(len(solution.t), 3) + self.assertEqual(len(solution.y), 4) + self.assertEqual([len(row) for row in solution.y], [3, 3, 3, 3]) + + +if __name__ == "__main__": + unittest.main()