From cca9d1e883bfd1edc580f58f476678673f230317 Mon Sep 17 00:00:00 2001 From: Codex Date: Sun, 16 Aug 2026 11:34:35 +0000 Subject: [PATCH] fix: seed PNVO pipe series pressure --- app/simulation/solvers/algebraic.py | 85 ++++++++++++------- app/simulation/solvers/algebraic_blocks.py | 8 +- ...est_pressure_flow_solver_initialization.py | 52 ++++++++++++ tests/test_thermofluid_closure_plan.py | 2 +- 4 files changed, 110 insertions(+), 37 deletions(-) diff --git a/app/simulation/solvers/algebraic.py b/app/simulation/solvers/algebraic.py index 40c2840..e3e1200 100644 --- a/app/simulation/solvers/algebraic.py +++ b/app/simulation/solvers/algebraic.py @@ -7,7 +7,10 @@ from math import expm1, isfinite, log, sqrt import os from app.simulation.components.amesim.boundary.sources import AmesimPnpl01 -from app.simulation.components.amesim.flow.orifices import AmesimPnor001 +from app.simulation.components.amesim.flow.orifices import ( + AmesimPnor001, + AmesimPnvo001FixedOpening, +) from app.simulation.components.amesim.flow.pipes import ( AmesimPnl00r, AmesimPnl0001, @@ -179,9 +182,11 @@ class AlgebraicEquationSubset: @dataclass(frozen=True) -class PnorPnl0001SeriesBinding: - orifice: AmesimPnor001 - orifice_port: str +class ResistancePnl0001SeriesBinding: + resistance: AmesimPnor001 | AmesimPnvo001FixedOpening + resistance_port: str + resistance_other_port: str + positive_flow_port: str pipe: AmesimPnl0001 pipe_port: str @@ -315,8 +320,8 @@ class PressureFlowSolver: variable: self._build_effort_equality_groups(variable) for variable in ("p", "x", "v") } - self._pnor_pnl0001_series_plan = ( - self._build_pnor_pnl0001_series_plan() + self._resistance_pnl0001_series_plan = ( + self._build_resistance_pnl0001_series_plan() ) self._closed_resistance_pressure_plan = ( self._build_closed_resistance_pressure_plan() @@ -498,7 +503,7 @@ class PressureFlowSolver: return failed("activeSetCausalizationRequired") if self._closed_resistance_pressure_plan: return failed("specialClosedResistancePressureSeed") - if self._pnor_pnl0001_series_plan: + if self._resistance_pnl0001_series_plan: return failed("specialSeriesPressureSeed") if len(self._unknowns_by_id) != len(self.unknowns): return failed("duplicateAlgebraicUnknown") @@ -1895,68 +1900,84 @@ class PressureFlowSolver: component.get_port(binding.port_name).p = pressure binding.neighbor.get_port(binding.neighbor_port).p = pressure - def _build_pnor_pnl0001_series_plan( + def _build_resistance_pnl0001_series_plan( self, - ) -> tuple[PnorPnl0001SeriesBinding, ...]: - bindings: list[PnorPnl0001SeriesBinding] = [] + ) -> tuple[ResistancePnl0001SeriesBinding, ...]: + bindings: list[ResistancePnl0001SeriesBinding] = [] + resistance_types = (AmesimPnor001, AmesimPnvo001FixedOpening) 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 + if isinstance(first, resistance_types) and isinstance( + second, AmesimPnl0001 + ): + resistance, resistance_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 + elif isinstance(second, resistance_types) and isinstance( + first, AmesimPnl0001 + ): + resistance, resistance_port = second, second_endpoint.port pipe, pipe_port = first, first_endpoint.port else: continue if isinstance(pipe, AmesimPnl0002) or pipe_port != "port_1": continue + if isinstance(resistance, AmesimPnor001): + positive_flow_port, negative_flow_port = "port_1", "port_2" + else: + positive_flow_port, negative_flow_port = "port_2", "port_3" + if resistance_port == positive_flow_port: + resistance_other_port = negative_flow_port + elif resistance_port == negative_flow_port: + resistance_other_port = positive_flow_port + else: + continue bindings.append( - PnorPnl0001SeriesBinding( - orifice=orifice, - orifice_port=orifice_port, + ResistancePnl0001SeriesBinding( + resistance=resistance, + resistance_port=resistance_port, + resistance_other_port=resistance_other_port, + positive_flow_port=positive_flow_port, pipe=pipe, pipe_port=pipe_port, ) ) return tuple(bindings) - def _seed_pnor_pnl0001_series_pressures(self) -> None: - """Causalize the pressure between a PNOR001 and PNL0001 R port.""" + def _seed_resistance_pnl0001_series_pressures(self) -> None: + """Causalize pressure between an orifice/valve and a PNL0001 R port.""" from scipy.optimize import brentq - for binding in self._pnor_pnl0001_series_plan: - orifice = binding.orifice - orifice_port = binding.orifice_port + for binding in self._resistance_pnl0001_series_plan: + resistance = binding.resistance + resistance_port = binding.resistance_port pipe = binding.pipe pipe_port = binding.pipe_port - orifice_other = "port_2" if orifice_port == "port_1" else "port_1" - pressure_a = orifice.get_port(orifice_other).p + pressure_a = resistance.get_port(binding.resistance_other_port).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, + if resistance_port == binding.positive_flow_port: + resistance_flow_into_connection = resistance.mass_flow( intermediate_pressure, + pressure_a, ) else: - orifice_flow_into_connection = orifice.mass_flow( - intermediate_pressure, + resistance_flow_into_connection = -resistance.mass_flow( pressure_a, + intermediate_pressure, ) pipe_flow_into_connection = pipe.mass_flow( intermediate_pressure, pressure_b, pipe.properties().T, ) - return orifice_flow_into_connection + pipe_flow_into_connection + return resistance_flow_into_connection + pipe_flow_into_connection lower_value = mismatch(lower) upper_value = mismatch(upper) @@ -1977,7 +1998,7 @@ class PressureFlowSolver: maxiter=32, ) ) - orifice.get_port(orifice_port).p = pressure + resistance.get_port(resistance_port).p = pressure pipe.get_port(pipe_port).p = pressure def _equation_scale_plan(self, equation) -> EquationScalePlan: @@ -2083,7 +2104,7 @@ class PressureFlowSolver: else: self._seed_equal_efforts(effort_variables) self._seed_closed_resistance_pressures() - self._seed_pnor_pnl0001_series_pressures() + self._seed_resistance_pnl0001_series_pressures() seeded_flow_ids = self._solve_explicit_flow_unknowns() contact_bindings = self._seed_unilateral_contacts() if contact_bindings: diff --git a/app/simulation/solvers/algebraic_blocks.py b/app/simulation/solvers/algebraic_blocks.py index cadbe3a..f0e9077 100644 --- a/app/simulation/solvers/algebraic_blocks.py +++ b/app/simulation/solvers/algebraic_blocks.py @@ -324,7 +324,7 @@ class StreamPressureBlockSolver: return failed(str(parent_reason or "parentCausalPathIneligible")) if solver._closed_resistance_pressure_plan: return failed("specialClosedResistancePressureSeed") - if solver._pnor_pnl0001_series_plan: + if solver._resistance_pnl0001_series_plan: return failed("specialSeriesPressureSeed") selected = self._selected_equation_evaluation if selected is None: @@ -725,10 +725,10 @@ class StreamPressureBlockSolver: f"{binding.neighbor.name}.{binding.neighbor_port}.p", ) ) - for binding in solver._pnor_pnl0001_series_plan: + for binding in solver._resistance_pnl0001_series_plan: special_seed_target_ids.extend( ( - f"{binding.orifice.name}.{binding.orifice_port}.p", + f"{binding.resistance.name}.{binding.resistance_port}.p", f"{binding.pipe.name}.{binding.pipe_port}.p", ) ) @@ -774,7 +774,7 @@ class StreamPressureBlockSolver: # closure deliberately mirrors ``solver.solve(effort_variables=())``: # the primary global solve has already propagated equal pressures. solver._seed_closed_resistance_pressures() - solver._seed_pnor_pnl0001_series_pressures() + solver._seed_resistance_pnl0001_series_pressures() for unknown in self._selected_flow_unknowns: unknown.write(0.0) for stage in self._selected_explicit_flow_plan: diff --git a/tests/test_pressure_flow_solver_initialization.py b/tests/test_pressure_flow_solver_initialization.py index a948892..cfc5071 100644 --- a/tests/test_pressure_flow_solver_initialization.py +++ b/tests/test_pressure_flow_solver_initialization.py @@ -58,6 +58,58 @@ class PressureFlowSolverInitializationTests(unittest.TestCase): "specialSeriesPressureSeed", ) + def test_pnvo_pnl0001_series_pressure_is_seeded_by_flow_balance(self) -> None: + medium = IdealGasMedium() + chamber = AmesimPnch023("source", medium, p0=15.3e6) + source_plug = AmesimPnpl01("source_closed") + valve = AmesimPnvo001SignalOpening( + "valve", + medium, + opening0=1.0, + ) + valve.res.signal = 1.0 + pipe = AmesimPnl0001("pipe", medium, p0=14.0e6) + storage_plug = AmesimPnpl01("storage_closed") + network = SimulationNetwork("pnvo-pnl0001-series") + for component in (chamber, source_plug, valve, pipe, storage_plug): + network.add_component(component) + network.connect("source_closed", "port_1", "source", "port_1") + network.connect("source", "port_2", "valve", "port_3") + network.connect("valve", "port_2", "pipe", "port_1") + network.connect("pipe", "port_2", "storage_closed", "port_1") + chamber.refresh_thermodynamic_ports() + pipe.refresh_thermodynamic_ports() + + solver = PressureFlowSolver(network) + result = solver.solve() + + self.assertTrue(result.success) + self.assertEqual(result.evaluations, 0) + self.assertGreater(valve.port_2.p, pipe.port_2.p) + self.assertLess(valve.port_2.p, chamber.port_2.p) + self.assertAlmostEqual(valve.port_2.p, pipe.port_1.p) + self.assertAlmostEqual(-valve.port_2.m_flow, pipe.port_1.m_flow) + + previous_intermediate_pressure = valve.port_2.p + replay_pressure = 12.0e6 + replay_temperature = 293.15 + replay_mass = ( + medium.density(replay_pressure, replay_temperature) * pipe.volume + ) + pipe.state = VolumeState( + m=replay_mass, + U=replay_mass * medium.specific_internal_energy(replay_temperature), + ) + pipe.refresh_thermodynamic_ports() + + replay = solver.solve() + + self.assertTrue(replay.success) + self.assertEqual(replay.evaluations, 0) + self.assertNotAlmostEqual(valve.port_2.p, previous_intermediate_pressure) + self.assertAlmostEqual(valve.port_2.p, pipe.port_1.p) + self.assertAlmostEqual(-valve.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) diff --git a/tests/test_thermofluid_closure_plan.py b/tests/test_thermofluid_closure_plan.py index 1bcfc27..cdd2074 100644 --- a/tests/test_thermofluid_closure_plan.py +++ b/tests/test_thermofluid_closure_plan.py @@ -360,7 +360,7 @@ class ThermofluidClosurePlanTests(unittest.TestCase): ) ] - self.assertEqual(len(series_solver._pnor_pnl0001_series_plan), 1) + self.assertEqual(len(series_solver._resistance_pnl0001_series_plan), 1) self.assertEqual(len(closed_pipe_solver._closed_resistance_pressure_plan), 2) self.assertEqual(len(resistance_solver._closed_resistance_pressure_plan), 1)