Fix pneumatic node zero-flow reversal

app/simulation/components/amesim/junctions/nodes.py

@@ -9,6 +9,24 @@ from app.simulation.core.medium import IdealGasMedium
 from app.simulation.core.ports import PortDefinition, PortState
 
 
+_REFERENCE_OUTFLOW_REGULARIZATION_RATIO = 0.05
+
+
+def _regularized_inverse_outflow(flow: float, transition_flow: float) -> float:
+    """Return a C1 inverse that tends to zero as a negative flow vanishes."""
+
+    if flow >= 0.0:
+        return 0.0
+    transition_flow = max(float(transition_flow), 1.0e-12)
+    if -flow >= transition_flow:
+        return 1.0 / flow
+    return (
+        flow
+        * (2.0 * transition_flow * transition_flow - flow * flow)
+        / transition_flow**4
+    )
+
+
 class _AmesimPneumaticNode(AlgebraicComponent):
     """Shared implementation for AMESim pneumatic junction submodels.
 
@@ -102,7 +120,7 @@ class _AmesimPneumaticNode(AlgebraicComponent):
                 else self.temperature_reference_h
             )
 
-        if reference_port.m_flow < -1e-12:
+        if reference_port.m_flow < 0.0:
             energy_without_reference = sum(
                 port.m_flow
                 * (
@@ -113,8 +131,29 @@ class _AmesimPneumaticNode(AlgebraicComponent):
                 for name, port in self.ports.items()
                 if name != self.REFERENCE_PORT
             )
+            non_reference_flow_scale = sum(
+                abs(port.m_flow)
+                for name, port in self.ports.items()
+                if name != self.REFERENCE_PORT
+            )
+            transition_flow = (
+                _REFERENCE_OUTFLOW_REGULARIZATION_RATIO
+                * non_reference_flow_scale
+            )
+            # Port 2 carries AMESim's residual-energy causality.  Exact
+            # division is singular when its outflow reverses through zero, so
+            # use a C1 band that matches the exact balance at its boundary and
+            # tends to the mixed enthalpy at zero flow.
+            inverse_flow = _regularized_inverse_outflow(
+                reference_port.m_flow,
+                transition_flow,
+            )
+            energy_residual_at_mixed_h = (
+                energy_without_reference
+                + reference_port.m_flow * mixed_h
+            )
             reference_port.h_outflow = (
-                -energy_without_reference / reference_port.m_flow
+                mixed_h - energy_residual_at_mixed_h * inverse_flow
             )
 
 

tests/test_amesim_pneumatic_node_components.py

@@ -123,6 +123,37 @@ class AmesimPneumaticNodeComponentTests(unittest.TestCase):
         )
         self.assertAlmostEqual(energy_flow, 0.0, places=10)
 
+    def test_port_2_outflow_enthalpy_remains_finite_through_flow_reversal(self) -> None:
+        node = AmesimP4Node2("p4_1")
+        connected_h = {
+            "port_1": 300.0,
+            "port_2": 300.0,
+            "port_3": 200.0,
+            "port_4": 100.0,
+        }
+
+        outflow_enthalpies: list[float] = []
+        for reference_flow in (-1.0e-6, 0.0, 1.0e-6):
+            node.port_1.m_flow = -1.0
+            node.port_2.m_flow = reference_flow
+            node.port_3.m_flow = -reference_flow
+            node.port_4.m_flow = 1.0
+
+            node.update_stream_outflows(connected_h)
+
+            self.assertAlmostEqual(
+                sum(port.m_flow for port in node.ports.values()),
+                0.0,
+                places=15,
+            )
+            outflow_enthalpies.append(node.port_2.h_outflow)
+
+        self.assertTrue(all(abs(value) < 1.0e4 for value in outflow_enthalpies))
+        self.assertLess(
+            max(outflow_enthalpies) - min(outflow_enthalpies),
+            1.0,
+        )
+
 
 if __name__ == "__main__":
     unittest.main()

tests/test_mechanical_solver_causalization.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from math import exp
+from math import exp, isfinite
 import unittest
 
 from app.simulation.components.amesim.mechanical.translational import (
@@ -130,7 +130,93 @@ def _single_mass_system(
     return GenericFluidSystem(network), mass
 
 
+def _high_stiffness_contact_impact_system() -> GenericFluidSystem:
+    """Return the reduced mechanical core of the four-branch stop impact.
+
+    The production model has four 50 kg piston groups coupled through LSTP00A
+    contacts to one 170000 kg support group.  When the support reaches its
+    ideal upper stop, its velocity is reset while each contact still carries
+    the finite, extremely fast LSTP damping transient.  One branch is enough
+    to retain that stiffness and state-transition interaction in a fast test.
+    """
+
+    medium = IdealGasMedium()
+    zero_left = AmesimF000("zero_left")
+    moving_mass = AmesimMecmas21(
+        "moving_mass",
+        medium,
+        mass=50.0,
+        useFriction=1.0,
+        stoptype=4.0,
+        x0=0.0090025,
+        v0=1.0,
+    )
+    contact = AmesimLstp00a(
+        "contact",
+        medium,
+        gap0=0.0,
+        kcont=1.0e11,
+        rcont=1.0e11,
+        Pdis=1.0e-7,
+        discContactOption=1.0,
+    )
+    support_mass = AmesimMecmas21(
+        "support_mass",
+        medium,
+        mass=170000.0,
+        useFriction=1.0,
+        stoptype=1.0,
+        xmin=0.0,
+        xmax=0.01,
+        x0=0.009,
+        v0=1.0,
+    )
+    zero_right = AmesimF000("zero_right")
+
+    network = SimulationNetwork("high-stiffness-contact-impact")
+    for component in (
+        zero_left,
+        moving_mass,
+        contact,
+        support_mass,
+        zero_right,
+    ):
+        network.add_component(component)
+    network.connect("zero_left", "port_1", "moving_mass", "port_1")
+    network.connect("moving_mass", "port_2", "contact", "port_1")
+    network.connect("contact", "port_2", "support_mass", "port_1")
+    network.connect("support_mass", "port_2", "zero_right", "port_1")
+    return GenericFluidSystem(network)
+
+
 class MechanicalSolverCausalizationTests(unittest.TestCase):
+    def test_high_stiffness_lstp_survives_ideal_stop_transition(self) -> None:
+        system = _high_stiffness_contact_impact_system()
+
+        result = system.simulate(
+            SolveIVPConfig(
+                t_start=0.0,
+                t_stop=0.002,
+                max_step=0.002,
+                method="BDF",
+                rtol=1.0e-6,
+            ),
+            sample_step=0.00025,
+        )
+
+        self.assertTrue(result.success, result.message)
+        totals = result.diagnostics["integration"]["totals"]
+        self.assertEqual(totals["stateTransitionCount"], 1)
+        self.assertEqual(totals["solverStartCount"], 2)
+        self.assertEqual(totals["recoverableRetryCount"], 0)
+
+        self.assertAlmostEqual(result.series["support_mass.x"][-1], 0.01)
+        self.assertAlmostEqual(result.series["support_mass.v"][-1], 0.0, delta=1.0e-9)
+        self.assertAlmostEqual(result.series["moving_mass.v"][-1], 0.0, delta=1.0e-4)
+        self.assertGreater(max(result.series["contact.force"]), 1.0e10)
+        for values in result.series.values():
+            self.assertTrue(all(isfinite(value) for value in values))
+
     def test_lstp_contact_uses_exponential_damping_ramp_and_negative_force_option(
         self,
     ) -> None: