From 568be2e632ff0c9bf3384a85f9c4b4de576bdba3 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Fri, 17 Jul 2026 08:54:38 +0000 Subject: [PATCH] =?UTF-8?q?=E5=AE=9E=E7=8E=B0test=5Fmql=20PNL0003=E4=B8=8E?= =?UTF-8?q?PNL00R=E7=AE=A1=E8=B7=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../components/amesim_pneumatic_line.py | 421 +++++++++++++++--- PythonModels/systems/test_mql.py | 24 + .../systems/test_mql_line_parameters.py | 185 +++++++- .../systems/test_mql_pneumatic_lines.py | 79 +++- .../test_amesim_pneumatic_line_components.py | 109 +++++ tests/test_test_mql_pnl0003_pnl00r.py | 94 ++++ 6 files changed, 844 insertions(+), 68 deletions(-) create mode 100644 tests/test_amesim_pneumatic_line_components.py create mode 100644 tests/test_test_mql_pnl0003_pnl00r.py diff --git a/PythonModels/components/amesim_pneumatic_line.py b/PythonModels/components/amesim_pneumatic_line.py index 08f6a35..83dcb5a 100644 --- a/PythonModels/components/amesim_pneumatic_line.py +++ b/PythonModels/components/amesim_pneumatic_line.py @@ -8,7 +8,7 @@ from PythonModels.components.amesim_pneumatic import ( AmesimPneumaticGas, diameter_mm_to_area_m2, ) -from PythonModels.core.base import DynamicComponent +from PythonModels.core.base import AlgebraicComponent, DynamicComponent from PythonModels.core.medium import ThermodynamicProperties from PythonModels.core.ports import PortState from PythonModels.core.state import VolumeState @@ -23,7 +23,89 @@ class AmesimPnl0001Diagnostics: pressure_drop_pa: float -class AmesimPnl0001Pipe(DynamicComponent): +class _DarcyPipeResistanceMixin: + diameter: float + length: float + relative_roughness: float + area: float + + 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 pneumatic pipe 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) + + +class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent): """Physical first-pass implementation of AMESim ``PNL0001`` (C-R). Port 2 owns the lumped gas storage. Port 1 is connected through a Darcy @@ -193,80 +275,289 @@ class AmesimPnl0001Pipe(DynamicComponent): 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( +class AmesimPnl0003Pipe(_DarcyPipeResistanceMixin, DynamicComponent): + """First-pass AMESim ``PNL0003`` (C-R-C) pipe. + + The two pipe-end compliances are represented as equal half-volume gas + stores connected by the same auditable Darcy resistance used for PNL0001. + Center flow is positive from port 1 storage to port 2 storage. + """ + + state_size = 4 + + def __init__( self, - mass_flow_kg_s: float, + name: str, *, - density: float, - temperature: float, - ) -> float: - if mass_flow_kg_s == 0.0: + 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, + p1_0: float = 101_325.0, + T1_0: float = 293.15, + p2_0: float = 101_325.0, + T2_0: 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.compliance_volume = self.volume / 2.0 + self.heat_transfer_area = pi * self.diameter * self.length + + self.state_1 = self._initial_state(p1_0, T1_0) + self.state_2 = self._initial_state(p2_0, T2_0) + self.port_1 = PortState() + self.port_2 = PortState() + + def _initial_state(self, pressure: float, temperature: float) -> VolumeState: + rho = self.gas.density(pressure, temperature) + mass = rho * self.compliance_volume + return VolumeState( + m=mass, + U=mass * self.gas.specific_internal_energy(temperature), + ) + + def get_state_vector(self) -> list[float]: + return [*self.state_1.as_vector(), *self.state_2.as_vector()] + + def set_state_vector(self, values: list[float]) -> None: + if len(values) != 4: + raise ValueError("PNL0003 state vector requires four values") + self.state_1 = VolumeState.from_vector(values[:2]) + self.state_2 = VolumeState.from_vector(values[2:]) + + def properties_1(self) -> ThermodynamicProperties: + properties = self._properties(self.state_1) + self.port_1.p = properties.p + self.port_1.h_outflow = properties.h + return properties + + def properties_2(self) -> ThermodynamicProperties: + properties = self._properties(self.state_2) + self.port_2.p = properties.p + self.port_2.h_outflow = properties.h + return properties + + def _properties(self, state: VolumeState) -> ThermodynamicProperties: + if state.m <= 0.0: + raise ValueError("pipe mass must stay positive") + temperature = self.gas.temperature_from_internal_energy(state.U / state.m) + density = state.m / self.compliance_volume + pressure = self.gas.pressure(density, temperature) + return ThermodynamicProperties( + p=pressure, + T=temperature, + rho=density, + u=state.U / state.m, + h=self.gas.specific_enthalpy(temperature), + ) + + def gas_mass_g(self) -> float: + return (self.state_1.m + self.state_2.m) * 1.0e3 + + def resistance_mass_flow(self) -> float: + """Return center mass flow from port 1 storage to port 2 storage.""" + port_1 = self.properties_1() + port_2 = self.properties_2() + pressure_difference = port_1.p - port_2.p + if pressure_difference == 0.0: return 0.0 + upstream = port_1 if pressure_difference > 0.0 else port_2 + magnitude = self._mass_flow_for_pressure_drop( + abs(pressure_difference), + density=upstream.rho, + temperature=upstream.T, + ) + return magnitude if pressure_difference > 0.0 else -magnitude + + def diagnostics( + self, + *, + mass_flow_kg_s: float, + temperature_k: float | None = None, + ) -> AmesimPnl0001Diagnostics: + port_1 = self.properties_1() + port_2 = self.properties_2() + temperature = temperature_k or (port_1.T if mass_flow_kg_s >= 0.0 else port_2.T) + density = port_1.rho if mass_flow_kg_s >= 0.0 else port_2.rho 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 + pressure_drop = self._darcy_pressure_drop( + mass_flow_kg_s, + density=density, + 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, + ) -> tuple[VolumeState, VolumeState]: + port_1 = self.properties_1() + port_2 = self.properties_2() + center_flow = self.resistance_mass_flow() + heat_flow_each = ( + self.heat_transfer_coefficient + * self.heat_transfer_area + * (self.external_temperature - 0.5 * (port_1.T + port_2.T)) / 2.0 ) - return magnitude if mass_flow_kg_s > 0.0 else -magnitude + port_1_external_h = self.connection_inlet_enthalpy( + port_m_flow=port_1_m_flow, + connected_h=connected_h_1, + internal_h=port_1.h, + ) + port_2_external_h = self.connection_inlet_enthalpy( + port_m_flow=port_2_m_flow, + connected_h=connected_h_2, + internal_h=port_2.h, + ) + port_1_center_h = self.connection_inlet_enthalpy( + port_m_flow=-center_flow, + connected_h=port_2.h, + internal_h=port_1.h, + ) + port_2_center_h = self.connection_inlet_enthalpy( + port_m_flow=center_flow, + connected_h=port_1.h, + internal_h=port_2.h, + ) + return ( + VolumeState( + m=port_1_m_flow - center_flow, + U=( + port_1_m_flow * port_1_external_h + - center_flow * port_1_center_h + + heat_flow_each + ), + ), + VolumeState( + m=port_2_m_flow + center_flow, + U=( + port_2_m_flow * port_2_external_h + + center_flow * port_2_center_h + + heat_flow_each + ), + ), + ) - 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) +class AmesimPnl00rPipe(_DarcyPipeResistanceMixin, AlgebraicComponent): + """First-pass AMESim ``PNL00R`` (R) pipe resistance.""" + + def __init__( + self, + name: str, + *, + diameter_mm: float, + length_m: float, + relative_roughness: float, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + ) -> 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") + + super().__init__(name=name) + self.diameter = diameter_mm * 1.0e-3 + self.length = length_m + self.relative_roughness = relative_roughness + self.gas = gas + self.area = diameter_mm_to_area_m2(diameter_mm) + self.port_1 = PortState() + self.port_2 = PortState() + + def mass_flow( + self, + *, + port_1_pressure_pa: float, + port_1_temperature_k: float, + port_2_pressure_pa: float, + port_2_temperature_k: float, + ) -> float: + """Return mass flow from port 1 to port 2 in kg/s.""" + if port_1_pressure_pa <= 0.0 or port_2_pressure_pa <= 0.0: + raise ValueError("port pressures must be positive") + if port_1_temperature_k <= 0.0 or port_2_temperature_k <= 0.0: + raise ValueError("port temperatures must be positive") + pressure_difference = port_1_pressure_pa - port_2_pressure_pa + if pressure_difference == 0.0: + return 0.0 + upstream_pressure = max(port_1_pressure_pa, port_2_pressure_pa) + upstream_temperature = ( + port_1_temperature_k + if pressure_difference > 0.0 + else port_2_temperature_k + ) + 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, + pressure_pa: float, + temperature_k: float, + ) -> AmesimPnl0001Diagnostics: + density = self.gas.density(pressure_pa, temperature_k) + reynolds = self._reynolds_number(mass_flow_kg_s, temperature_k) + friction_factor = self._friction_factor(reynolds) + velocity = mass_flow_kg_s / (density * self.area) + pressure_drop = self._darcy_pressure_drop( + mass_flow_kg_s, + density=density, + temperature=temperature_k, + ) + 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 helium_dynamic_viscosity(temperature_k: float) -> float: diff --git a/PythonModels/systems/test_mql.py b/PythonModels/systems/test_mql.py index 3b26199..fcb0328 100644 --- a/PythonModels/systems/test_mql.py +++ b/PythonModels/systems/test_mql.py @@ -3860,6 +3860,8 @@ class TestMqlSystem: 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.pnl0003_assembly = self._build_pnl0003_assembly() + self.pnl00r_assembly = self._build_pnl00r_assembly() self.node3_assembly = self._build_node3_assembly() self.pneumatic_assembly = self._build_pneumatic_assembly() pneumatic_components = self._pneumatic_components_by_alias() @@ -3890,6 +3892,20 @@ class TestMqlSystem: return build_test_mql_pnl0001_assembly(self.archive_path) + def _build_pnl0003_assembly(self): + from PythonModels.systems.test_mql_pneumatic_lines import ( + build_test_mql_pnl0003_assembly, + ) + + return build_test_mql_pnl0003_assembly(self.archive_path) + + def _build_pnl00r_assembly(self): + from PythonModels.systems.test_mql_pneumatic_lines import ( + build_test_mql_pnl00r_assembly, + ) + + return build_test_mql_pnl00r_assembly(self.archive_path) + @staticmethod def _build_node3_assembly(): from PythonModels.systems.test_mql_nodes import build_test_mql_node3_assembly @@ -3912,6 +3928,14 @@ class TestMqlSystem: def typed_pnl0001_line_count(self) -> int: return len(self.pnl0001_assembly.lines) + @property + def typed_pnl0003_line_count(self) -> int: + return len(self.pnl0003_assembly.lines) + + @property + def typed_pnl00r_line_count(self) -> int: + return len(self.pnl00r_assembly.lines) + @property def typed_node3_count(self) -> int: return len(self.node3_assembly) diff --git a/PythonModels/systems/test_mql_line_parameters.py b/PythonModels/systems/test_mql_line_parameters.py index fef40c7..ff4c87d 100644 --- a/PythonModels/systems/test_mql_line_parameters.py +++ b/PythonModels/systems/test_mql_line_parameters.py @@ -35,6 +35,48 @@ class TestMqlPnl0001Spec: return self.initial_gauge_pressure_pa + AMESIM_REFERENCE_PRESSURE_PA +@dataclass(frozen=True) +class TestMqlPnl0003Spec: + 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_1_k: float + initial_gauge_pressure_1_pa: float + initial_temperature_2_k: float + initial_gauge_pressure_2_pa: float + + @property + def initial_absolute_pressure_1_pa(self) -> float: + return self.initial_gauge_pressure_1_pa + AMESIM_REFERENCE_PRESSURE_PA + + @property + def initial_absolute_pressure_2_pa(self) -> float: + return self.initial_gauge_pressure_2_pa + AMESIM_REFERENCE_PRESSURE_PA + + +@dataclass(frozen=True) +class TestMqlPnl00rSpec: + alias: str + source_component: str + source_port: str + target_component: str + target_port: str + diameter_mm: float + length_m: float + relative_roughness: float + gas_type_index: int + + def load_test_mql_pnl0001_specs( archive_path: str | Path, *, @@ -111,6 +153,145 @@ def load_test_mql_pnl0001_specs( return tuple(specs) +def load_test_mql_pnl0003_specs( + archive_path: str | Path, + *, + cir_member: str = "test_mql_.cir", +) -> tuple[TestMqlPnl0003Spec, ...]: + """Load resolved PNL0003 geometry and both compliance initial states.""" + 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"] == "PNL0003" + } + specs = [] + for block in re.findall(r".*?", cir_text, flags=re.DOTALL): + if _optional_text(block, "SUB_NAME") != "PNL0003": + continue + alias = _required_text(block, "ALIAS") + connection = connections.get(alias) + if connection is None: + raise ValueError(f"PNL0003 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( + TestMqlPnl0003Spec( + 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_1_k=_required_numeric( + alias, "t1", state_values, numeric_globals + ), + initial_gauge_pressure_1_pa=_required_numeric( + alias, "p1", state_values, numeric_globals + ), + initial_temperature_2_k=_required_numeric( + alias, "t2", state_values, numeric_globals + ), + initial_gauge_pressure_2_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 PNL0003 parameter blocks: {missing}") + return tuple(specs) + + +def load_test_mql_pnl00r_specs( + archive_path: str | Path, + *, + cir_member: str = "test_mql_.cir", +) -> tuple[TestMqlPnl00rSpec, ...]: + """Load resolved PNL00R geometry 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"] == "PNL00R" + } + specs = [] + for block in re.findall(r".*?", cir_text, flags=re.DOTALL): + if _optional_text(block, "SUB_NAME") != "PNL00R": + continue + alias = _required_text(block, "ALIAS") + connection = connections.get(alias) + if connection is None: + raise ValueError(f"PNL00R line {alias!r} is absent from CONNECTION_SPECS") + real_parameters = _parameter_expressions(block, "RPARAM") + integer_parameters = _parameter_expressions(block, "IPARAM") + specs.append( + TestMqlPnl00rSpec( + 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 + ), + gas_type_index=int( + _required_numeric(alias, "gi", integer_parameters, 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 PNL00R parameter blocks: {missing}") + return tuple(specs) + + def _parameter_expressions(block: str, tag_name: str) -> dict[str, str]: parameters = {} for parameter_block in re.findall( @@ -141,11 +322,11 @@ def _required_numeric( variables: dict[str, float], ) -> float: if name not in expressions: - raise ValueError(f"Missing {name!r} on PNL0001 line {alias!r}") + raise ValueError(f"Missing {name!r} on 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}" + f"Cannot resolve {name!r}={expressions[name]!r} on line {alias!r}" ) return value diff --git a/PythonModels/systems/test_mql_pneumatic_lines.py b/PythonModels/systems/test_mql_pneumatic_lines.py index 9d4db9a..50ebc54 100644 --- a/PythonModels/systems/test_mql_pneumatic_lines.py +++ b/PythonModels/systems/test_mql_pneumatic_lines.py @@ -7,10 +7,18 @@ from PythonModels.components.amesim_pneumatic import ( HELIUM_PNEUMATIC_GAS, AmesimPneumaticGas, ) -from PythonModels.components.amesim_pneumatic_line import AmesimPnl0001Pipe +from PythonModels.components.amesim_pneumatic_line import ( + AmesimPnl0001Pipe, + AmesimPnl0003Pipe, + AmesimPnl00rPipe, +) from PythonModels.systems.test_mql_line_parameters import ( TestMqlPnl0001Spec, + TestMqlPnl0003Spec, + TestMqlPnl00rSpec, load_test_mql_pnl0001_specs, + load_test_mql_pnl0003_specs, + load_test_mql_pnl00r_specs, ) @@ -26,6 +34,30 @@ class TestMqlPnl0001Assembly: raise KeyError(alias) +@dataclass(frozen=True) +class TestMqlPnl0003Assembly: + specs: tuple[TestMqlPnl0003Spec, ...] + lines: dict[str, AmesimPnl0003Pipe] + + def spec(self, alias: str) -> TestMqlPnl0003Spec: + for spec in self.specs: + if spec.alias == alias: + return spec + raise KeyError(alias) + + +@dataclass(frozen=True) +class TestMqlPnl00rAssembly: + specs: tuple[TestMqlPnl00rSpec, ...] + lines: dict[str, AmesimPnl00rPipe] + + def spec(self, alias: str) -> TestMqlPnl00rSpec: + for spec in self.specs: + if spec.alias == alias: + return spec + raise KeyError(alias) + + def build_test_mql_pnl0001_assembly( archive_path: str | Path, *, @@ -48,3 +80,48 @@ def build_test_mql_pnl0001_assembly( for spec in specs } return TestMqlPnl0001Assembly(specs=specs, lines=lines) + + +def build_test_mql_pnl0003_assembly( + archive_path: str | Path, + *, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, +) -> TestMqlPnl0003Assembly: + specs = load_test_mql_pnl0003_specs(archive_path) + lines = { + spec.alias: AmesimPnl0003Pipe( + 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, + p1_0=spec.initial_absolute_pressure_1_pa, + T1_0=spec.initial_temperature_1_k, + p2_0=spec.initial_absolute_pressure_2_pa, + T2_0=spec.initial_temperature_2_k, + ) + for spec in specs + } + return TestMqlPnl0003Assembly(specs=specs, lines=lines) + + +def build_test_mql_pnl00r_assembly( + archive_path: str | Path, + *, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, +) -> TestMqlPnl00rAssembly: + specs = load_test_mql_pnl00r_specs(archive_path) + lines = { + spec.alias: AmesimPnl00rPipe( + name=spec.alias, + diameter_mm=spec.diameter_mm, + length_m=spec.length_m, + relative_roughness=spec.relative_roughness, + gas=gas, + ) + for spec in specs + } + return TestMqlPnl00rAssembly(specs=specs, lines=lines) diff --git a/tests/test_amesim_pneumatic_line_components.py b/tests/test_amesim_pneumatic_line_components.py new file mode 100644 index 0000000..66a101b --- /dev/null +++ b/tests/test_amesim_pneumatic_line_components.py @@ -0,0 +1,109 @@ +from __future__ import annotations + +import unittest + +from PythonModels.components.amesim_pneumatic_line import ( + AmesimPnl0003Pipe, + AmesimPnl00rPipe, +) + + +class AmesimPnl0003PipeTests(unittest.TestCase): + def setUp(self) -> None: + self.pipe = AmesimPnl0003Pipe( + name="pneumatic_88", + diameter_mm=20.0, + length_m=0.3, + relative_roughness=0.045 / 20.0, + p1_0=15.3e6, + p2_0=15.3e6, + T1_0=293.15, + T2_0=293.15, + ) + + def test_initial_state_uses_two_half_volume_compliances(self) -> None: + port_1 = self.pipe.properties_1() + port_2 = self.pipe.properties_2() + + self.assertAlmostEqual(self.pipe.volume, 9.424777960769381e-5) + self.assertAlmostEqual(self.pipe.compliance_volume, self.pipe.volume / 2.0) + self.assertEqual(len(self.pipe.get_state_vector()), 4) + self.assertAlmostEqual(port_1.p, 15.3e6, delta=1.0e-5) + self.assertAlmostEqual(port_2.p, 15.3e6, delta=1.0e-5) + self.assertAlmostEqual(port_1.T, 293.15) + self.assertAlmostEqual(port_2.T, 293.15) + + def test_center_resistance_flow_follows_end_pressure_gradient(self) -> None: + state = self.pipe.get_state_vector() + state[0] *= 1.01 + state[1] *= 1.01 + self.pipe.set_state_vector(state) + + forward = self.pipe.resistance_mass_flow() + state[0] /= 1.01 * 1.01 + state[1] /= 1.01 * 1.01 + state[2] *= 1.01 + state[3] *= 1.01 + self.pipe.set_state_vector(state) + reverse = self.pipe.resistance_mass_flow() + + self.assertGreater(forward, 0.0) + self.assertLess(reverse, 0.0) + + def test_connection_derivatives_conserve_internal_center_flow_mass(self) -> None: + state = self.pipe.get_state_vector() + state[0] *= 1.01 + state[1] *= 1.01 + self.pipe.set_state_vector(state) + + d1, d2 = self.pipe.derivatives_from_connections( + port_1_m_flow=0.2, + connected_h_1=self.pipe.properties_1().h + 1000.0, + port_2_m_flow=-0.1, + connected_h_2=self.pipe.properties_2().h - 1000.0, + ) + + self.assertAlmostEqual(d1.m + d2.m, 0.1) + + +class AmesimPnl00rPipeTests(unittest.TestCase): + def setUp(self) -> None: + self.pipe = AmesimPnl00rPipe( + name="pneumatic_100", + diameter_mm=14.0, + length_m=1.0, + relative_roughness=0.045 / 14.0, + ) + + def test_stateless_resistance_flow_follows_pressure_gradient(self) -> None: + forward = self.pipe.mass_flow( + port_1_pressure_pa=15.31e6, + port_1_temperature_k=293.15, + port_2_pressure_pa=15.29e6, + port_2_temperature_k=293.15, + ) + reverse = self.pipe.mass_flow( + port_1_pressure_pa=15.29e6, + port_1_temperature_k=293.15, + port_2_pressure_pa=15.31e6, + port_2_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_expose_darcy_terms(self) -> None: + diagnostics = self.pipe.diagnostics( + mass_flow_kg_s=1.0e-4, + pressure_pa=15.3e6, + temperature_k=293.15, + ) + + self.assertGreater(diagnostics.reynolds_number, 0.0) + self.assertGreater(diagnostics.gas_velocity_m_s, 0.0) + self.assertGreater(diagnostics.pressure_drop_pa, 0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_test_mql_pnl0003_pnl00r.py b/tests/test_test_mql_pnl0003_pnl00r.py new file mode 100644 index 0000000..a47bb19 --- /dev/null +++ b/tests/test_test_mql_pnl0003_pnl00r.py @@ -0,0 +1,94 @@ +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 import TestMqlSystem +from PythonModels.systems.test_mql_line_parameters import ( + load_test_mql_pnl0003_specs, + load_test_mql_pnl00r_specs, +) +from PythonModels.systems.test_mql_pneumatic_lines import ( + build_test_mql_pnl0003_assembly, + build_test_mql_pnl00r_assembly, +) + + +REPO_ROOT = Path(__file__).resolve().parents[1] +TEST_MQL_AME = REPO_ROOT / "AmesimModels" / "test_mql.ame" + + +class TestMqlPnl0003AndPnl00rTests(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + cls.pnl0003_specs = load_test_mql_pnl0003_specs(TEST_MQL_AME) + cls.pnl00r_specs = load_test_mql_pnl00r_specs(TEST_MQL_AME) + cls.pnl0003_assembly = build_test_mql_pnl0003_assembly(TEST_MQL_AME) + cls.pnl00r_assembly = build_test_mql_pnl00r_assembly(TEST_MQL_AME) + cls.results = load_test_mql_amesim_results(TEST_MQL_AME) + + def test_loads_all_real_pnl0003_parameters_from_cir(self) -> None: + self.assertEqual(len(self.pnl0003_specs), 8) + spec = self.pnl0003_assembly.spec("pneumatic_88") + + self.assertEqual(spec.source_component, "pn_node3_9") + self.assertEqual(spec.target_component, "pn_morifice_9") + self.assertEqual(spec.diameter_mm, 20.0) + self.assertEqual(spec.length_m, 0.3) + self.assertAlmostEqual(spec.relative_roughness, 0.045 / 20.0) + self.assertEqual(spec.gas_type_index, 1) + self.assertEqual(spec.mode, 2) + self.assertAlmostEqual(spec.initial_gauge_pressure_1_pa, 15_198_700.0) + self.assertAlmostEqual(spec.initial_gauge_pressure_2_pa, 15_198_700.0) + self.assertAlmostEqual(spec.initial_absolute_pressure_1_pa, 15_300_000.0) + self.assertAlmostEqual(spec.initial_absolute_pressure_2_pa, 15_300_000.0) + + def test_loads_all_real_pnl00r_parameters_from_cir(self) -> None: + self.assertEqual(len(self.pnl00r_specs), 4) + spec = self.pnl00r_assembly.spec("pneumatic_100") + + self.assertEqual(spec.source_component, "pn_node3_9") + self.assertEqual(spec.target_component, "pn_node3_10") + 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.gas_type_index, 1) + + def test_builds_pnl0003_and_pnl00r_physical_line_components(self) -> None: + self.assertEqual(len(self.pnl0003_assembly.lines), 8) + self.assertEqual(len(self.pnl00r_assembly.lines), 4) + self.assertTrue( + all(len(line.get_state_vector()) == 4 for line in self.pnl0003_assembly.lines.values()) + ) + + def test_pneumatic_88_initial_observables_match_amesim_baseline(self) -> None: + pipe = self.pnl0003_assembly.lines["pneumatic_88"] + + self.assertAlmostEqual( + pipe.properties_1().p - 101_300.0, + self.results.series("p1@pneumatic_88")[0], + delta=1.0e-5, + ) + self.assertAlmostEqual( + pipe.properties_2().p - 101_300.0, + self.results.series("p2@pneumatic_88")[0], + delta=1.0e-5, + ) + self.assertAlmostEqual( + pipe.gas_mass_g(), + self.results.series("mgas@pneumatic_88")[0], + delta=0.005, + ) + + def test_system_exposes_real_pnl0003_and_pnl00r_assemblies(self) -> None: + system = TestMqlSystem() + + self.assertEqual(system.typed_pnl0003_line_count, 8) + self.assertEqual(system.typed_pnl00r_line_count, 4) + self.assertIn("pneumatic_88", system.pnl0003_assembly.lines) + self.assertIn("pneumatic_100", system.pnl00r_assembly.lines) + + +if __name__ == "__main__": + unittest.main()