对齐AMESim管阻孔口与储气耦合

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(

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(

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

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:

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(

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:

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

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():