diff --git a/app/simulation/components/amesim/flow/orifices.py b/app/simulation/components/amesim/flow/orifices.py index 3b5b58e..2a9c56d 100644 --- a/app/simulation/components/amesim/flow/orifices.py +++ b/app/simulation/components/amesim/flow/orifices.py @@ -51,13 +51,12 @@ _PNVO001_FLOW_COEFFICIENT_GROUP = ParameterGroupDisplaySpec( class AmesimPnor001(AlgebraicComponent): """AMESim PNOR001 constant-flow-coefficient pneumatic orifice. - This public component preserves the PNOR001 catalog/XML contract and uses a - finite bidirectional compressible-orifice approximation. The Siemens - `pn2rcqfix_` details remain a later calibration target. + This public component preserves the PNOR001 catalog/XML contract and uses + real-gas pressure-ratio flow with AMESim-style near-equal-pressure smoothing. """ MODEL_TYPE = "amesim_pnor001" - MODEL_VERSION = "0.2.0" + MODEL_VERSION = "0.3.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -188,6 +187,7 @@ class AmesimPnor001(AlgebraicComponent): self.port_1.h_outflow = initial_h self.port_2 = self.register_declared_port("port_2") self.port_2.h_outflow = initial_h + self._connected_h: dict[str, float] = {} @staticmethod def _integer_parameter(name: str, value: float) -> int: @@ -243,16 +243,113 @@ class AmesimPnor001(AlgebraicComponent): def _upstream_temperature(self, port_name: str) -> float: port = self.get_port(port_name) + inlet_h = self._connected_h.get(port_name, port.h_outflow) return max( self.medium.temperature_from_pressure_enthalpy( max(port.p, 1.0), - port.h_outflow, + inlet_h, ), 1.0, ) + @staticmethod + def _subsonic_mass_flow_parameter( + *, + pressure_ratio: float, + gamma_s: float, + density: float, + upstream_temperature: float, + upstream_pressure: float, + ) -> float: + expansion = ( + pressure_ratio ** (2.0 * gamma_s) + - pressure_ratio ** (1.0 + gamma_s) + ) + return sqrt( + max( + 2.0 + / (1.0 - gamma_s) + * density + * upstream_temperature + / upstream_pressure + * expansion, + 0.0, + ) + ) + + def _one_way_flow_characteristics( + self, + *, + upstream_pressure: float, + downstream_pressure: float, + upstream_temperature: float, + ) -> tuple[float, float]: + p_up = max(upstream_pressure, 1.0) + p_down = max(min(downstream_pressure, p_up), 0.0) + T_up = max(upstream_temperature, 1.0) + gamma_s = self.medium.isentropic_density_pressure_factor( + p_up, + T_up, + p_down, + ) + gamma_s = min(max(gamma_s, 1.0e-9), 1.0 - 1.0e-9) + density = max(self.medium.density(p_up, T_up), 1.0e-12) + pressure_ratio = max(p_down / p_up, 0.0) + critical_ratio = (2.0 * gamma_s / (gamma_s + 1.0)) ** ( + 1.0 / (1.0 - gamma_s) + ) + if pressure_ratio <= critical_ratio: + effective_pressure_ratio = critical_ratio + mass_flow_parameter = ( + sqrt(2.0 / (1.0 + gamma_s) * density * T_up / p_up) + * (2.0 * gamma_s / (gamma_s + 1.0)) + ** (gamma_s / (1.0 - gamma_s)) + ) + gas_velocity = sqrt(2.0 / (1.0 + gamma_s) * p_up / density) + else: + effective_pressure_ratio = pressure_ratio + mass_flow_parameter = self._subsonic_mass_flow_parameter( + pressure_ratio=pressure_ratio, + gamma_s=gamma_s, + density=density, + upstream_temperature=T_up, + upstream_pressure=p_up, + ) + gas_velocity = sqrt( + max( + 2.0 + / (1.0 - gamma_s) + * p_up + / density + * (1.0 - pressure_ratio ** (1.0 - gamma_s)), + 0.0, + ) + ) + + reference = self._subsonic_mass_flow_parameter( + pressure_ratio=_PN_PRESSURE_RATIO_ACCURACY, + gamma_s=gamma_s, + density=density, + upstream_temperature=T_up, + upstream_pressure=p_up, + ) + if mass_flow_parameter > 0.0 and reference > 0.0: + argument = ( + _PN_LAMINAR_SMOOTHING_GAIN + * abs(mass_flow_parameter / reference) + * log(effective_pressure_ratio) + / log(_PN_PRESSURE_RATIO_ACCURACY) + ) + smoothing_factor = tanh(max(argument, 0.0)) + mass_flow_parameter *= smoothing_factor + gas_velocity *= smoothing_factor + return mass_flow_parameter, gas_velocity + def mass_flow(self, p_1: float, p_2: float) -> float: - if p_1 == p_2 or self.effective_area == 0.0: + if ( + isclose(p_1, p_2, rel_tol=0.0, abs_tol=1.0e-8) + or self.effective_area == 0.0 + ): return 0.0 if p_1 > p_2: return self._one_way_mass_flow( @@ -274,43 +371,41 @@ class AmesimPnor001(AlgebraicComponent): upstream_temperature: float, ) -> float: p_up = max(upstream_pressure, 1.0) - p_down = max(min(downstream_pressure, p_up), 0.0) T_up = max(upstream_temperature, 1.0) - gamma = max(self.medium.gamma, 1.000001) - pressure_ratio = max(p_down / p_up, 0.0) - critical_ratio = (2.0 / (gamma + 1.0)) ** (gamma / (gamma - 1.0)) - if pressure_ratio <= critical_ratio: - flow_factor = sqrt(gamma / (self.medium.R_gas * T_up)) * ( - 2.0 / (gamma + 1.0) - ) ** ((gamma + 1.0) / (2.0 * (gamma - 1.0))) - else: - expansion = pressure_ratio ** (2.0 / gamma) - pressure_ratio ** ( - (gamma + 1.0) / gamma - ) - flow_factor = sqrt( - max( - 2.0 - * gamma - * expansion - / (self.medium.R_gas * T_up * (gamma - 1.0)), - 0.0, - ) - ) - return self.effective_cq * self.effective_area * p_up * flow_factor + mass_flow_parameter, _ = self._one_way_flow_characteristics( + upstream_pressure=p_up, + downstream_pressure=downstream_pressure, + upstream_temperature=T_up, + ) + return ( + self.effective_cq + * self.effective_area + * p_up + * mass_flow_parameter + / sqrt(T_up) + ) def component_result_values(self) -> Mapping[str, float]: p_1 = max(self.port_1.p, 1.0) p_2 = max(self.port_2.p, 1.0) - m_flow = abs(self.mass_flow(self.port_1.p, self.port_2.p)) - upstream_pressure = max(p_1, p_2) - upstream_temperature = self._upstream_temperature( - "port_1" if p_1 >= p_2 else "port_2" + if p_1 >= p_2: + upstream_port_name = "port_1" + upstream_pressure = p_1 + downstream_pressure = p_2 + flow_direction = 1.0 + else: + upstream_port_name = "port_2" + upstream_pressure = p_2 + downstream_pressure = p_1 + flow_direction = -1.0 + mass_flow_parameter, gas_velocity = self._one_way_flow_characteristics( + upstream_pressure=upstream_pressure, + downstream_pressure=downstream_pressure, + upstream_temperature=self._upstream_temperature(upstream_port_name), ) - density = max(self.medium.density(upstream_pressure, upstream_temperature), 1.0e-12) - area = max(self.effective_area, 1.0e-18) return { - "cm": m_flow / (self.effective_cq * area * upstream_pressure), - "gasvel": m_flow / (density * area), + "cm": mass_flow_parameter, + "gasvel": flow_direction * gas_velocity, } def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: @@ -344,6 +439,7 @@ class AmesimPnor001(AlgebraicComponent): ) def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: + self._connected_h = dict(connected_h) self.port_1.h_outflow = connected_h["port_2"] self.port_2.h_outflow = connected_h["port_1"] @@ -593,7 +689,10 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): ) def mass_flow(self, p_2: float, p_3: float) -> float: - if p_2 == p_3 or self.effective_area == 0.0: + if ( + isclose(p_2, p_3, rel_tol=1.0e-7, abs_tol=1.0e-9) + or self.effective_area == 0.0 + ): return 0.0 if p_2 > p_3: return self._one_way_mass_flow( diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index d635614..75553f5 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -1,7 +1,7 @@ from __future__ import annotations from collections.abc import Mapping -from math import isclose, log10, pi, sqrt +from math import isclose, log, log10, pi, sqrt, tanh from app.simulation.components.amesim.gases import ( AMESIM_GAS_INDEX_PARAMETER, @@ -246,7 +246,7 @@ class AmesimPnl00r(AlgebraicComponent): while self.darcy_pressure_drop(upper, density=density, temperature=temperature) < pressure_drop: upper *= 10.0 if upper > 1.0e3: - raise ValueError("unable to bracket PNL00R resistance flow") + return 1.0e3 lower = 0.0 for _ in range(48): middle = 0.5 * (lower + upper) @@ -257,7 +257,7 @@ class AmesimPnl00r(AlgebraicComponent): return 0.5 * (lower + upper) def mass_flow(self, p_1: float, p_2: float) -> float: - if p_1 == p_2: + if isclose(p_1, p_2, rel_tol=1.0e-7, abs_tol=1.0e-9): return 0.0 pressure_difference = p_1 - p_2 upstream_pressure = max(p_1, p_2, 1.0) @@ -645,7 +645,7 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): ) < pressure_drop: upper *= 10.0 if upper > 1.0e3: - raise ValueError("unable to bracket PNL0001 resistance flow") + return 1.0e3 lower = 0.0 for _ in range(48): middle = 0.5 * (lower + upper) @@ -659,16 +659,126 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): upper = middle return 0.5 * (lower + upper) + def _one_way_pn2pipefr_mass_flow( + self, + *, + upstream_pressure: float, + downstream_pressure: float, + upstream_temperature: float, + resistance_length: float, + ) -> float: + """AMESim pn2pipefr-style compressible friction flow.""" + + p_up = max(float(upstream_pressure), 1.0) + p_down = max(min(float(downstream_pressure), p_up), 0.0) + T_up = max(float(upstream_temperature), 1.0) + if resistance_length <= 0.0: + raise ValueError("Pipe resistance length must be positive.") + + gamma_s = self.medium.isentropic_density_pressure_factor( + p_up, + T_up, + p_down, + ) + gamma_s = min(max(gamma_s, 1.0e-9), 1.0 - 1.0e-9) + density = max(self.medium.density(p_up, T_up), 1.0e-12) + pressure_ratio = max(p_down / p_up, 0.0) + critical_ratio = (2.0 * gamma_s / (gamma_s + 1.0)) ** ( + 1.0 / (1.0 - gamma_s) + ) + + def mass_flow_parameter(ratio: float) -> float: + if ratio <= critical_ratio: + value = ( + sqrt(2.0 / (1.0 + gamma_s) * density * T_up / p_up) + * (2.0 * gamma_s / (gamma_s + 1.0)) + ** (gamma_s / (1.0 - gamma_s)) + ) + effective_ratio = critical_ratio + else: + expansion = ratio ** (2.0 * gamma_s) - ratio ** (1.0 + gamma_s) + value = sqrt( + max( + 2.0 + / (1.0 - gamma_s) + * density + * T_up + / p_up + * expansion, + 0.0, + ) + ) + effective_ratio = ratio + + accuracy = 0.9999 + reference_expansion = ( + accuracy ** (2.0 * gamma_s) + - accuracy ** (1.0 + gamma_s) + ) + reference = sqrt( + max( + 2.0 + / (1.0 - gamma_s) + * density + * T_up + / p_up + * reference_expansion, + 0.0, + ) + ) + if value > 0.0 and reference > 0.0: + smoothing_argument = ( + 12.0 + * abs(value / reference) + * log(effective_ratio) + / log(accuracy) + ) + value *= tanh(max(smoothing_argument, 0.0)) + return value + + def target_flow(mass_flow: float) -> float: + reynolds = self.reynolds_number(mass_flow, T_up) + friction = self.friction_factor(reynolds) + flow_coefficient = sqrt( + self.diam / (resistance_length * friction) + ) + return ( + flow_coefficient + * self.area + * p_up + * mass_flow_parameter(pressure_ratio) + / sqrt(T_up) + ) + + flow_coefficient = sqrt(self.diam / (resistance_length * 0.02)) + magnitude = ( + flow_coefficient + * self.area + * p_up + * mass_flow_parameter(pressure_ratio) + / sqrt(T_up) + ) + for _iteration in range(16): + next_magnitude = target_flow(magnitude) + if abs(next_magnitude - magnitude) <= max( + 1.0e-12, + abs(magnitude) * 1.0e-9, + ): + return next_magnitude + magnitude = 0.5 * (magnitude + next_magnitude) + return magnitude + def mass_flow(self, p_1: float, p_2: float, temperature: float) -> float: - if p_1 == p_2: + if isclose(p_1, p_2, rel_tol=0.0, abs_tol=1.0e-8): return 0.0 pressure_difference = p_1 - p_2 - upstream_pressure = max(p_1, p_2, 1.0) - density = max(self.medium.density(upstream_pressure, temperature), 1.0e-12) - magnitude = self._mass_flow_for_pressure_drop( - abs(pressure_difference), - density=density, - temperature=temperature, + upstream_temperature = max(float(temperature), 1.0) + resistance_length = getattr(self, "resistance_length", self.le) + magnitude = self._one_way_pn2pipefr_mass_flow( + upstream_pressure=max(p_1, p_2), + downstream_pressure=min(p_1, p_2), + upstream_temperature=upstream_temperature, + resistance_length=resistance_length, ) return magnitude if pressure_difference > 0.0 else -magnitude @@ -1115,7 +1225,7 @@ class AmesimPnl0003(DynamicComponent): while self.darcy_pressure_drop(upper, density=density, temperature=temperature) < pressure_drop: upper *= 10.0 if upper > 1.0e3: - raise ValueError("unable to bracket PNL0003 resistance flow") + return 1.0e3 lower = 0.0 for _ in range(48): middle = 0.5 * (lower + upper) @@ -1129,7 +1239,7 @@ class AmesimPnl0003(DynamicComponent): port_1 = self._properties(self.state_1) port_2 = self._properties(self.state_2) pressure_difference = port_1.p - port_2.p - if pressure_difference == 0.0: + if isclose(port_1.p, port_2.p, rel_tol=1.0e-7, abs_tol=1.0e-9): return 0.0 upstream = port_1 if pressure_difference > 0.0 else port_2 magnitude = self._mass_flow_for_pressure_drop( diff --git a/app/simulation/components/amesim/media/mediums.py b/app/simulation/components/amesim/media/mediums.py index 7ef9290..d9c54df 100644 --- a/app/simulation/components/amesim/media/mediums.py +++ b/app/simulation/components/amesim/media/mediums.py @@ -4,6 +4,7 @@ from collections.abc import Callable from dataclasses import dataclass from typing import ClassVar +from app.simulation.core.errors import RecoverableTrialStateError from app.simulation.core.medium import ( GasMedium, IdealGasMedium, @@ -214,7 +215,9 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: if m <= 0.0: - raise ValueError("Mass must stay positive when recovering temperature.") + raise RecoverableTrialStateError( + "Mass must stay positive when recovering temperature." + ) return self.temperature_from_internal_energy(U / m) def properties_from_mU( @@ -224,7 +227,9 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): V: float, ) -> ThermodynamicProperties: if m <= 0.0: - raise ValueError("Mass must stay positive when recovering temperature.") + raise RecoverableTrialStateError( + "Mass must stay positive when recovering temperature." + ) if V <= 0.0: raise ValueError("Volume must stay positive.") density = m / V diff --git a/app/simulation/core/medium.py b/app/simulation/core/medium.py index 9108ec4..614d418 100644 --- a/app/simulation/core/medium.py +++ b/app/simulation/core/medium.py @@ -3,6 +3,8 @@ from __future__ import annotations from dataclasses import dataclass from typing import Protocol +from app.simulation.core.errors import RecoverableTrialStateError + @dataclass(frozen=True) class ThermodynamicProperties: @@ -195,7 +197,9 @@ class IdealGasMedium: def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: if m <= 0.0: - raise ValueError("Mass must stay positive when recovering temperature.") + raise RecoverableTrialStateError( + "Mass must stay positive when recovering temperature." + ) return self.temperature_from_internal_energy(U / m) def pressure(self, m: float, T: float, V: float) -> float: diff --git a/app/simulation/core/peng_robinson.py b/app/simulation/core/peng_robinson.py index 2365a0f..a711e3d 100644 --- a/app/simulation/core/peng_robinson.py +++ b/app/simulation/core/peng_robinson.py @@ -299,12 +299,12 @@ class PengRobinsonFluid: @staticmethod def _validate_temperature(temperature: float) -> None: if temperature <= 0.0: - raise ValueError("Temperature must be positive.") + raise RecoverableTrialStateError("Temperature must be positive.") @classmethod def _validate_pressure_temperature(cls, pressure: float, temperature: float) -> None: if pressure <= 0.0: - raise ValueError("Pressure must be positive.") + raise RecoverableTrialStateError("Pressure must be positive.") cls._validate_temperature(temperature) HELIUM_PR = PengRobinsonFluid( diff --git a/app/simulation/solvers/algebraic.py b/app/simulation/solvers/algebraic.py index 643977f..6130a86 100644 --- a/app/simulation/solvers/algebraic.py +++ b/app/simulation/solvers/algebraic.py @@ -5,6 +5,13 @@ from collections.abc import Callable from dataclasses import dataclass from math import expm1, isfinite, log, sqrt +from app.simulation.components.amesim.boundary.sources import AmesimPnpl01 +from app.simulation.components.amesim.flow.orifices import AmesimPnor001 +from app.simulation.components.amesim.flow.pipes import ( + AmesimPnl00r, + AmesimPnl0001, + AmesimPnl0002, +) from app.simulation.core.ports import PortState, VariableRole from app.simulation.systems.network import SimulationNetwork @@ -716,7 +723,7 @@ class PressureFlowSolver: """Execute the precompiled explicit flow/force causalization plan.""" for unknown in self.unknowns: - if unknown.variable == "f": + if unknown.variable in {"f", "m_flow"}: unknown.write(0.0) seeded_ids: set[str] = set() @@ -728,6 +735,109 @@ class PressureFlowSolver: seeded_ids.add(assignment.unknown.id) return seeded_ids + def _seed_closed_resistance_pressures(self) -> None: + """Seed a sealed resistance end at its zero-flow pressure. + + A PNPL01 fixes flow, not pressure. Starting a dead-ended Darcy branch + with the plug-side pressure at the medium reference can otherwise put + the nonlinear solver on the singular square-root part of the inverse + flow law. At zero flow, these AMESim pipe resistances have exactly zero + pressure drop, which gives a deterministic and physically exact seed. + """ + + connected: dict[tuple[str, str], tuple[str, str]] = {} + for connection in self.network.connections: + if connection.kind != "physical" or connection.domain != "pneumatic": + continue + first, second = connection.endpoints + connected[first.key] = second.key + connected[second.key] = first.key + + for component in self.network.components.values(): + if not isinstance(component, (AmesimPnl00r, AmesimPnl0001)): + continue + for port_name in component.ports: + neighbor_key = connected.get((component.name, port_name)) + if neighbor_key is None: + continue + neighbor = self.network.components[neighbor_key[0]] + if not isinstance(neighbor, AmesimPnpl01): + continue + if isinstance(component, AmesimPnl0002): + pressure = component.properties().p + elif isinstance(component, AmesimPnl0001): + if port_name != "port_1": + continue + pressure = component.properties().p + else: + other_port_name = "port_2" if port_name == "port_1" else "port_1" + pressure = component.get_port(other_port_name).p + component.get_port(port_name).p = pressure + neighbor.get_port(neighbor_key[1]).p = pressure + + def _seed_pnor_pnl0001_series_pressures(self) -> None: + """Causalize the pressure between a PNOR001 and PNL0001 R port.""" + + for connection in self.network.connections: + first_endpoint, second_endpoint = connection.endpoints + first = self.network.components[first_endpoint.component] + second = self.network.components[second_endpoint.component] + if isinstance(first, AmesimPnor001) and isinstance(second, AmesimPnl0001): + orifice, orifice_port = first, first_endpoint.port + pipe, pipe_port = second, second_endpoint.port + elif isinstance(second, AmesimPnor001) and isinstance(first, AmesimPnl0001): + orifice, orifice_port = second, second_endpoint.port + pipe, pipe_port = first, first_endpoint.port + else: + continue + if isinstance(pipe, AmesimPnl0002) or pipe_port != "port_1": + continue + + orifice_other = "port_2" if orifice_port == "port_1" else "port_1" + pressure_a = orifice.get_port(orifice_other).p + pressure_b = pipe.properties().p + lower = min(pressure_a, pressure_b) + upper = max(pressure_a, pressure_b) + + def mismatch(intermediate_pressure: float) -> float: + if orifice_port == "port_2": + orifice_flow_into_connection = -orifice.mass_flow( + pressure_a, + intermediate_pressure, + ) + else: + orifice_flow_into_connection = orifice.mass_flow( + intermediate_pressure, + pressure_a, + ) + pipe_flow_into_connection = pipe.mass_flow( + intermediate_pressure, + pressure_b, + pipe.properties().T, + ) + return orifice_flow_into_connection + pipe_flow_into_connection + + lower_value = mismatch(lower) + upper_value = mismatch(upper) + if lower_value == 0.0: + pressure = lower + elif upper_value == 0.0: + pressure = upper + elif (lower_value < 0.0) == (upper_value < 0.0): + continue + else: + for _iteration in range(64): + middle = 0.5 * (lower + upper) + middle_value = mismatch(middle) + if (middle_value < 0.0) == (lower_value < 0.0): + lower = middle + lower_value = middle_value + else: + upper = middle + pressure = 0.5 * (lower + upper) + orifice.get_port(orifice_port).p = pressure + pipe.get_port(pipe_port).p = pressure + def _scales(self) -> dict[str, float]: pressure_scale = max( [ @@ -783,6 +893,8 @@ class PressureFlowSolver: clear_causal_contact() self._seed_equal_efforts() + self._seed_closed_resistance_pressures() + self._seed_pnor_pnl0001_series_pressures() self._solve_explicit_flow_unknowns() contact_bindings = self._seed_unilateral_contacts() if contact_bindings: diff --git a/app/simulation/solvers/pneumatic_storage.py b/app/simulation/solvers/pneumatic_storage.py new file mode 100644 index 0000000..33621a9 --- /dev/null +++ b/app/simulation/solvers/pneumatic_storage.py @@ -0,0 +1,256 @@ +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Iterable, Sequence + +from app.simulation.components.amesim.flow.pipes import ( + AmesimPnl0001, + AmesimPnl0002, + AmesimPnl0003, +) +from app.simulation.solvers.mechanical import MechanicalConstraintGroup +from app.simulation.systems.network import Endpoint, SimulationNetwork + +if TYPE_CHECKING: + from app.simulation.core.base import DynamicComponent + from app.simulation.solvers.mechanical import MechanicalStateReducer + + +@dataclass(frozen=True) +class PneumaticStoragePartition: + """One fixed-volume ``[mass, internal energy]`` pressure-state partition.""" + + component: DynamicComponent + state_offset: int + volume: float + + @property + def key(self) -> tuple[str, int]: + return self.component.name, self.state_offset + + +@dataclass(frozen=True) +class IdealPneumaticStorageGroup: + partitions: tuple[PneumaticStoragePartition, ...] + + @property + def names(self) -> tuple[str, ...]: + return tuple(partition.component.name for partition in self.partitions) + + +def pneumatic_storage_partition( + network: SimulationNetwork, + endpoint: Endpoint, +) -> PneumaticStoragePartition | None: + """Map an AMESim pressure-state port to its fixed gas-volume state slice. + + PNL0001 exposes its C side at ``port_2``. PNL0003 exposes one compliance + at each end. PNL0002 has resistances at both external ports, so its center + compliance is intentionally not returned here. + """ + + component = network.components[endpoint.component] + if isinstance(component, AmesimPnl0003): + if endpoint.port == "port_1": + return PneumaticStoragePartition( + component=component, + state_offset=0, + volume=component.compliance_volume, + ) + if endpoint.port == "port_2": + return PneumaticStoragePartition( + component=component, + state_offset=2, + volume=component.compliance_volume, + ) + return None + if isinstance(component, AmesimPnl0001) and not isinstance( + component, AmesimPnl0002 + ): + if endpoint.port == "port_2": + return PneumaticStoragePartition( + component=component, + state_offset=0, + volume=component.volume, + ) + return None + + +def ideal_storage_group_is_reducible( + network: SimulationNetwork, + storage_endpoints: Iterable[Endpoint], +) -> bool: + partitions = [ + pneumatic_storage_partition(network, endpoint) + for endpoint in storage_endpoints + ] + if not partitions or any(partition is None for partition in partitions): + return False + unique = {partition.key: partition for partition in partitions if partition} + if len(unique) < 2: + return False + media = {id(partition.component.medium) for partition in unique.values()} + return len(media) == 1 and all(partition.volume > 0.0 for partition in unique.values()) + + +def _pressure_storage_endpoint_groups( + network: SimulationNetwork, +) -> tuple[tuple[Endpoint, ...], ...]: + pneumatic_endpoints = { + Endpoint(component.name, definition.name) + for component in network.components.values() + for definition in component.port_definitions + if definition.kind == "physical" and definition.domain == "pneumatic" + } + parent = {endpoint: endpoint for endpoint in pneumatic_endpoints} + + def find(endpoint: Endpoint) -> Endpoint: + root = endpoint + while parent[root] != root: + root = parent[root] + while parent[endpoint] != endpoint: + next_endpoint = parent[endpoint] + parent[endpoint] = root + endpoint = next_endpoint + return root + + def union(first: Endpoint, second: Endpoint) -> None: + first_root = find(first) + second_root = find(second) + if first_root != second_root: + parent[second_root] = first_root + + for connection in network.connections: + first, second = connection.endpoints + if first in pneumatic_endpoints and second in pneumatic_endpoints: + union(first, second) + + storage_endpoints: set[Endpoint] = set() + for component in network.components.values(): + for equation in component.pressure_flow_equation_residuals(): + pressure_endpoints = [ + Endpoint(component.name, variable.rsplit(".", 2)[1]) + for variable in equation.variables + if variable.startswith(f"{component.name}.") and variable.endswith(".p") + ] + if equation.relation == "equal": + for endpoint in pressure_endpoints[1:]: + union(pressure_endpoints[0], endpoint) + elif equation.relation == "state": + storage_endpoints.update(pressure_endpoints) + + by_root: dict[Endpoint, list[Endpoint]] = {} + for endpoint in storage_endpoints: + by_root.setdefault(find(endpoint), []).append(endpoint) + return tuple(tuple(endpoints) for endpoints in by_root.values()) + + +class IdealPneumaticStorageReducer: + """Project supported ideal C-C connections onto one thermodynamic state. + + AMESim permits compatible pipe compliances to share an ideal pneumatic + junction. The public ODE solver keeps the original result states but + projects their mass and energy densities together and distributes the + group's total derivative by physical volume. This removes the redundant + pressure constraint without adding a fictitious resistance. + """ + + def __init__( + self, + network: SimulationNetwork, + mechanical_state_reducer: MechanicalStateReducer, + ) -> None: + self.network = network + self.mechanical_state_reducer = mechanical_state_reducer + self.groups = self._build_groups() + self._component_offsets = self._build_component_offsets() + + def _build_groups(self) -> tuple[IdealPneumaticStorageGroup, ...]: + groups: list[IdealPneumaticStorageGroup] = [] + for endpoints in _pressure_storage_endpoint_groups(self.network): + partitions = [ + pneumatic_storage_partition(self.network, endpoint) + for endpoint in endpoints + ] + unique = { + partition.key: partition + for partition in partitions + if partition is not None + } + if len(unique) > 1 and len(unique) == len( + {endpoint.component for endpoint in endpoints} + ): + group = IdealPneumaticStorageGroup(tuple(unique.values())) + if ideal_storage_group_is_reducible(self.network, endpoints): + groups.append(group) + return tuple(groups) + + def _build_component_offsets(self) -> dict[str, int]: + offsets: dict[str, int] = {} + cursor = 0 + for entry in self.mechanical_state_reducer.state_entries: + if isinstance(entry, MechanicalConstraintGroup): + cursor += 2 + else: + offsets[entry.name] = cursor + cursor += entry.state_size + return offsets + + def _global_offset(self, partition: PneumaticStoragePartition) -> int: + return self._component_offsets[partition.component.name] + partition.state_offset + + def synchronize_state_vector( + self, + values: Sequence[float], + *, + validate: bool = False, + ) -> list[float]: + projected = [float(value) for value in values] + for group in self.groups: + volumes = [partition.volume for partition in group.partitions] + offsets = [self._global_offset(partition) for partition in group.partitions] + mass_densities = [ + projected[offset] / volume + for offset, volume in zip(offsets, volumes) + ] + energy_densities = [ + projected[offset + 1] / volume + for offset, volume in zip(offsets, volumes) + ] + if validate: + mass_scale = max([abs(value) for value in mass_densities] + [1.0]) + energy_scale = max([abs(value) for value in energy_densities] + [1.0]) + if ( + max(mass_densities) - min(mass_densities) > 1.0e-9 * mass_scale + or max(energy_densities) - min(energy_densities) + > 1.0e-9 * energy_scale + ): + raise ValueError( + "Ideally coupled AMESim pipe compliances require consistent " + "initial pressure and temperature: " + ", ".join(group.names) + ) + total_volume = sum(volumes) + mass_density = sum(projected[offset] for offset in offsets) / total_volume + energy_density = ( + sum(projected[offset + 1] for offset in offsets) / total_volume + ) + for offset, volume in zip(offsets, volumes): + projected[offset] = mass_density * volume + projected[offset + 1] = energy_density * volume + return projected + + def coupled_derivatives(self, values: Sequence[float]) -> list[float]: + derivatives = [float(value) for value in values] + for group in self.groups: + volumes = [partition.volume for partition in group.partitions] + offsets = [self._global_offset(partition) for partition in group.partitions] + total_volume = sum(volumes) + total_mass_derivative = sum(derivatives[offset] for offset in offsets) + total_energy_derivative = sum( + derivatives[offset + 1] for offset in offsets + ) + for offset, volume in zip(offsets, volumes): + fraction = volume / total_volume + derivatives[offset] = total_mass_derivative * fraction + derivatives[offset + 1] = total_energy_derivative * fraction + return derivatives diff --git a/app/simulation/systems/generic.py b/app/simulation/systems/generic.py index 0871ca7..dfffb0f 100644 --- a/app/simulation/systems/generic.py +++ b/app/simulation/systems/generic.py @@ -9,6 +9,10 @@ from app.simulation.core.base import DynamicComponent from app.simulation.core.metadata import ResultVariableMetadata from app.simulation.solvers.algebraic import PressureFlowSolver from app.simulation.solvers.mechanical import MechanicalStateReducer +from app.simulation.solvers.pneumatic_storage import ( + IdealPneumaticStorageReducer, + ideal_storage_group_is_reducible, +) from app.simulation.solvers.pneumatic_volume import PneumaticVolumeResolver from app.simulation.solvers.solver import ODESolution, SolveIVPConfig, integrate_ode from app.simulation.solvers.signal import SignalResolver @@ -182,13 +186,19 @@ def simulation_preparation_issues( for endpoint in pressure_ports: storage_ports[endpoint] = component.name - storages_by_group: dict[Endpoint, set[str]] = {} + storages_by_group: dict[Endpoint, dict[Endpoint, str]] = {} for endpoint, component_name in storage_ports.items(): - storages_by_group.setdefault(effort_groups.find(endpoint), set()).add( - component_name - ) - for storage_names in storages_by_group.values(): + storages_by_group.setdefault(effort_groups.find(endpoint), {})[ + endpoint + ] = component_name + for storage_endpoints in storages_by_group.values(): + storage_names = set(storage_endpoints.values()) if len(storage_names) > 1: + if ideal_storage_group_is_reducible( + network, + storage_endpoints, + ): + continue issues.append( SimulationPreparationIssue( "IDEAL_STORAGE_COUPLING_UNSUPPORTED", @@ -238,6 +248,10 @@ class GenericFluidSystem: network, self.dynamic_components, ) + self.pneumatic_storage_reducer = IdealPneumaticStorageReducer( + network, + self.mechanical_state_reducer, + ) self.pressure_flow_solver = PressureFlowSolver(network) self.pneumatic_volume_resolver = PneumaticVolumeResolver(network) self.signal_resolver = SignalResolver(network) @@ -250,10 +264,15 @@ class GenericFluidSystem: self.pneumatic_volume_propagation_count = 0 def initial_state_vector(self) -> list[float]: - return self.mechanical_state_reducer.initial_state_vector() + return self.pneumatic_storage_reducer.synchronize_state_vector( + self.mechanical_state_reducer.initial_state_vector(), + validate=True, + ) def apply_state_vector(self, values: list[float]) -> None: - self.mechanical_state_reducer.apply_state_vector(values) + self.mechanical_state_reducer.apply_state_vector( + self.pneumatic_storage_reducer.synchronize_state_vector(values) + ) def _close_current_state(self, time: float) -> dict[str, dict[str, float]]: signal = self.signal_resolver.solve(time) @@ -298,7 +317,9 @@ class GenericFluidSystem: def rhs(self, _time: float, state_vector: list[float]) -> list[float]: self.apply_state_vector(state_vector) connected_h = self._close_current_state(_time) - return self.mechanical_state_reducer.state_derivatives(connected_h) + return self.pneumatic_storage_reducer.coupled_derivatives( + self.mechanical_state_reducer.state_derivatives(connected_h) + ) def _append_current_state(self, series: dict[str, list[float]]) -> None: for component in self.network.components.values(): diff --git a/tests/test_amesim_helium_medium.py b/tests/test_amesim_helium_medium.py index f909368..803eb9c 100644 --- a/tests/test_amesim_helium_medium.py +++ b/tests/test_amesim_helium_medium.py @@ -10,6 +10,7 @@ from app.simulation.components.amesim.media.mediums import ( AMESIM_HELIUM_PROPERTY_MODELS, AmesimHeliumPengRobinsonMedium, ) +from app.simulation.core.errors import RecoverableTrialStateError from app.simulation.core.medium import IdealGasMedium from app.simulation.core.peng_robinson import HELIUM_PR from app.simulation.registry import COMPONENT_MODEL_REGISTRY @@ -109,6 +110,16 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): delta=2.0e-8, ) + def test_negative_mass_is_a_recoverable_integrator_trial_state(self) -> None: + media = (IdealGasMedium(), AmesimHeliumPengRobinsonMedium()) + + for medium in media: + with self.subTest(medium=medium.name): + with self.assertRaises(RecoverableTrialStateError): + medium.temperature_from_mass_internal_energy(-1.0, 1.0) + with self.assertRaises(RecoverableTrialStateError): + medium.properties_from_mU(-1.0, 1.0, 1.0) + def test_amesim_2404_real_gas_isentropic_factor_reference(self) -> None: medium = AmesimHeliumPengRobinsonMedium() diff --git a/tests/test_amesim_pnl0001_component.py b/tests/test_amesim_pnl0001_component.py index 11d8f22..e6bdd43 100644 --- a/tests/test_amesim_pnl0001_component.py +++ b/tests/test_amesim_pnl0001_component.py @@ -88,6 +88,20 @@ class AmesimPnl0001ComponentTests(unittest.TestCase): self.assertLess(reverse, 0.0) self.assertAlmostEqual(abs(forward), abs(reverse), delta=abs(forward) * 0.02) + def test_high_pressure_near_equal_flow_is_continuous_outside_roundoff(self) -> None: + pipe = AmesimPnl0001("pnl_1", self.medium, p0=15.3e6) + + small = pipe.mass_flow(15.3e6 + 0.1, 15.3e6, 293.15) + larger = pipe.mass_flow(15.3e6 + 1.0, 15.3e6, 293.15) + + self.assertGreater(small, 0.0) + self.assertGreater(larger, small) + + def test_pressure_roundoff_does_not_create_a_fictitious_flow(self) -> None: + pipe = AmesimPnl0001("pnl_1", self.medium, p0=15.3e6) + + self.assertEqual(pipe.mass_flow(15.3e6, 15.3e6 + 2.0e-9, 293.15), 0.0) + def test_pressure_flow_residuals_do_not_force_storage_mass_balance(self) -> None: pipe = AmesimPnl0001("pnl_1", self.medium) pipe.port_1.p = 101000.0 diff --git a/tests/test_amesim_pnor001_component.py b/tests/test_amesim_pnor001_component.py index 3bd48be..8766b96 100644 --- a/tests/test_amesim_pnor001_component.py +++ b/tests/test_amesim_pnor001_component.py @@ -62,6 +62,17 @@ class AmesimPnor001ComponentTests(unittest.TestCase): self.assertLess(reverse, 0.0) self.assertAlmostEqual(forward, -reverse) + def test_high_pressure_near_equal_flow_is_continuous(self) -> None: + orifice = AmesimPnor001("pnor_1", self.medium) + orifice.port_1.p = 15.3e6 + orifice.port_2.p = 15.3e6 + + small = orifice.mass_flow(15.3e6 + 0.1, 15.3e6) + larger = orifice.mass_flow(15.3e6 + 1.0, 15.3e6) + + self.assertGreater(small, 0.0) + self.assertGreater(larger, small) + def test_cv_and_kv_modes_produce_positive_equivalent_area(self) -> None: cv_orifice = AmesimPnor001( "pnor_cv", diff --git a/tests/test_component_catalog.py b/tests/test_component_catalog.py index 10afd1d..f96dbb0 100644 --- a/tests/test_component_catalog.py +++ b/tests/test_component_catalog.py @@ -421,7 +421,10 @@ class ComponentCatalogTests(unittest.TestCase): with self.subTest(model_type=model_type): component = self.components[model_type] model_parameters = parameters(model_type) - self.assertEqual(component["modelVersion"], "0.2.0") + self.assertEqual( + component["modelVersion"], + "0.3.0" if model_type == "amesim_pnor001" else "0.2.0", + ) self.assertEqual(model_parameters["flowset"]["editor"], "choice") self.assertEqual(model_parameters["flowset"]["options"], flow_options) self.assertEqual( diff --git a/tests/test_ideal_pneumatic_storage_coupling.py b/tests/test_ideal_pneumatic_storage_coupling.py new file mode 100644 index 0000000..4b79844 --- /dev/null +++ b/tests/test_ideal_pneumatic_storage_coupling.py @@ -0,0 +1,58 @@ +from __future__ import annotations + +import unittest + +from app.simulation.components.amesim.boundary.sources import AmesimPnpl01 +from app.simulation.components.amesim.flow.pipes import AmesimPnl0001, AmesimPnl0003 +from app.simulation.core.medium import IdealGasMedium +from app.simulation.solvers.solver import SolveIVPConfig +from app.simulation.systems.generic import GenericFluidSystem, simulation_preparation_issues +from app.simulation.systems.network import SimulationNetwork + + +class IdealPneumaticStorageCouplingTests(unittest.TestCase): + def build_network(self) -> SimulationNetwork: + medium = IdealGasMedium() + network = SimulationNetwork("ideal-pipe-compliance-coupling") + for component in ( + AmesimPnpl01("closed_1"), + AmesimPnl0001("pnl0001", medium, p0=15.3e6, T0=293.15), + AmesimPnl0003( + "pnl0003", + medium, + p1_0=15.3e6, + T1_0=293.15, + p2_0=15.3e6, + T2_0=293.15, + ), + AmesimPnpl01("closed_2"), + ): + network.add_component(component) + network.connect("closed_1", "port_1", "pnl0001", "port_1") + network.connect("pnl0001", "port_2", "pnl0003", "port_1") + network.connect("pnl0003", "port_2", "closed_2", "port_1") + return network + + def test_supported_amesim_pipe_compliances_are_not_rejected(self) -> None: + issues = simulation_preparation_issues(self.build_network()) + + self.assertNotIn( + "IDEAL_STORAGE_COUPLING_UNSUPPORTED", + {issue.code for issue in issues}, + ) + + def test_coupled_pipe_compliances_run_without_fictitious_resistance(self) -> None: + system = GenericFluidSystem(self.build_network()) + + result = system.simulate( + SolveIVPConfig(t_stop=0.001, method="BDF", max_step=0.0002), + sample_step=0.001, + ) + + self.assertTrue(result.success) + self.assertAlmostEqual(result.final["pnl0001.p"], 15.3e6, delta=1.0e-4) + self.assertAlmostEqual(result.final["pnl0003.p1"], 15.3e6, delta=1.0e-4) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_pressure_flow_solver_initialization.py b/tests/test_pressure_flow_solver_initialization.py index bc844d6..74a72d7 100644 --- a/tests/test_pressure_flow_solver_initialization.py +++ b/tests/test_pressure_flow_solver_initialization.py @@ -6,11 +6,15 @@ from unittest.mock import patch from app.simulation.components.amesim.boundary.sources import AmesimPnpl01 from app.simulation.components.amesim.flow.orifices import ( + AmesimPnor001, AmesimPnvo001SignalOpening, ) +from app.simulation.components.amesim.flow.pipes import AmesimPnl00r +from app.simulation.components.amesim.flow.pipes import AmesimPnl0001 from app.simulation.components.amesim.media.mediums import ( AmesimHeliumPengRobinsonMedium, ) +from app.simulation.components.amesim.storage.chambers import AmesimPnch023 from app.simulation.components.experimental.storage.cylinder import Cylinder from app.simulation.components.experimental.storage.tank import Tank from app.simulation.core.medium import IdealGasMedium @@ -21,6 +25,50 @@ from app.simulation.systems.network import SimulationNetwork class PressureFlowSolverInitializationTests(unittest.TestCase): + def test_pnor_pnl0001_series_pressure_is_seeded_by_flow_balance(self) -> None: + medium = IdealGasMedium() + chamber = AmesimPnch023("source", medium, p0=15.3e6) + source_plug = AmesimPnpl01("source_closed") + orifice = AmesimPnor001("orifice", medium) + pipe = AmesimPnl0001("pipe", medium, p0=14.0e6) + storage_plug = AmesimPnpl01("storage_closed") + network = SimulationNetwork("pnor-pnl0001-series") + for component in (chamber, source_plug, orifice, pipe, storage_plug): + network.add_component(component) + network.connect("source_closed", "port_1", "source", "port_1") + network.connect("source", "port_2", "orifice", "port_1") + network.connect("orifice", "port_2", "pipe", "port_1") + network.connect("pipe", "port_2", "storage_closed", "port_1") + chamber.refresh_thermodynamic_ports() + pipe.refresh_thermodynamic_ports() + + result = PressureFlowSolver(network).solve() + + self.assertTrue(result.success) + self.assertGreater(orifice.port_2.p, pipe.port_2.p) + self.assertLess(orifice.port_2.p, chamber.port_2.p) + self.assertAlmostEqual(orifice.port_2.p, pipe.port_1.p) + self.assertAlmostEqual(-orifice.port_2.m_flow, pipe.port_1.m_flow) + + def test_dead_ended_pnl00r_is_seeded_at_zero_flow_pressure(self) -> None: + medium = IdealGasMedium() + pipe = AmesimPnl00r("resistance", medium) + plug = AmesimPnpl01("closed") + pipe.port_1.p = 15.3e6 + pipe.port_2.p = 100_000.0 + network = SimulationNetwork("dead-ended-resistance") + network.add_component(pipe) + network.add_component(plug) + network.connect("resistance", "port_2", "closed", "port_1") + + result = PressureFlowSolver(network).solve() + + self.assertTrue(result.success) + self.assertEqual(result.evaluations, 0) + self.assertAlmostEqual(pipe.port_2.p, pipe.port_1.p) + self.assertEqual(pipe.port_1.m_flow, 0.0) + self.assertEqual(pipe.port_2.m_flow, 0.0) + @staticmethod def _near_equal_pressure_network() -> tuple[ SimulationNetwork,