From caca32a513428da36b6f32892c5a967d467d1fb8 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Tue, 11 Aug 2026 12:13:09 +0000 Subject: [PATCH] =?UTF-8?q?=E6=A0=A1=E5=87=86=E7=AC=AC=E4=BA=8C=E6=94=AF?= =?UTF-8?q?=E8=B7=AF=E7=83=AD=E6=B5=81=E4=BD=93=E8=83=BD=E9=87=8F=E4=B8=8E?= =?UTF-8?q?=E7=AE=A1=E8=B7=AF=E6=91=A9=E6=93=A6?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../components/amesim/flow/pipes.py | 116 +++++++++++++++--- .../components/amesim/junctions/nodes.py | 45 +++++-- app/simulation/core/base.py | 13 ++ app/simulation/solvers/mechanical.py | 33 +++-- app/simulation/solvers/stream.py | 20 +++ app/simulation/systems/generic.py | 6 + .../test_amesim_pneumatic_node_components.py | 42 ++++++- .../test_amesim_pnl0002_pnl0003_component.py | 17 +++ tests/test_amesim_pnl00r_component.py | 44 +++++++ tests/test_component_catalog.py | 2 +- tests/test_mechanical_solver_causalization.py | 34 ++++- ...est_pressure_flow_solver_initialization.py | 6 +- 12 files changed, 336 insertions(+), 42 deletions(-) diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index b35b8c2..8ebcc17 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -47,7 +47,7 @@ class AmesimPnl00r(AlgebraicComponent): """ MODEL_TYPE = "amesim_pnl00r" - MODEL_VERSION = "0.1.0" + MODEL_VERSION = "0.2.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -203,13 +203,34 @@ class AmesimPnl00r(AlgebraicComponent): laminar = 64.0 / reynolds_number if reynolds_number <= 2300.0: return laminar - turbulent = 1.0 / ( - -1.8 * log10((self.rr / 3.7) ** 1.11 + 6.9 / reynolds_number) + + # pn2pipefr does not apply the fully rough correction at every + # turbulent Reynolds number. Its saved ff curves first follow the + # hydraulically smooth law and approach the rough asymptote as Re*rr + # grows. Keeping those two limits separate reproduces the AMESim + # curves for both 14 mm and 20 mm test_mql pipes; putting both terms + # directly inside one Haaland logarithm over-predicts PNL0002 friction + # by about 23 percent near Re=57,000. + smooth_turbulent = 1.0 / ( + -1.8 * log10(6.9 / reynolds_number) ) ** 2 + if self.rr <= 0.0: + turbulent = smooth_turbulent + else: + fully_rough = 1.0 / ( + -1.8 * log10((self.rr / 3.7) ** 1.11) + ) ** 2 + roughness_reynolds = reynolds_number * self.rr + roughness_weight = roughness_reynolds * roughness_reynolds / ( + roughness_reynolds * roughness_reynolds + 180.0 * 180.0 + ) + turbulent = smooth_turbulent + roughness_weight * ( + fully_rough - smooth_turbulent + ) if reynolds_number >= 4000.0: return turbulent fraction = (reynolds_number - 2300.0) / 1700.0 - return laminar + fraction * (turbulent - laminar) + return laminar + fraction**0.58 * (turbulent - laminar) def darcy_pressure_drop( self, @@ -279,7 +300,11 @@ class AmesimPnl00r(AlgebraicComponent): density = max(self.medium.density(upstream_pressure, upstream_temperature), 1.0e-12) reynolds = self.reynolds_number(m_flow, upstream_temperature) velocity = m_flow / (density * self.area) - cm = abs(m_flow) / max(self.area * upstream_pressure, 1.0e-18) + cm = ( + abs(m_flow) + * sqrt(upstream_temperature) + / max(self.area * upstream_pressure, 1.0e-18) + ) return { "re": reynolds, "cm": cm, @@ -326,7 +351,7 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): """AMESim PNL0001 C-R pneumatic pipe with compressibility and friction.""" MODEL_TYPE = "amesim_pnl0001" - 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"), @@ -798,7 +823,11 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): "u": props.u, "h": props.h, "re": reynolds, - "cm": abs(flow) / max(self.area * upstream_pressure, 1.0e-18), + "cm": ( + abs(flow) + * sqrt(props.T) + / max(self.area * upstream_pressure, 1.0e-18) + ), "v": flow / (density * self.area), "ff": self.friction_factor(reynolds), } @@ -861,7 +890,7 @@ class AmesimPnl0002(AmesimPnl0001): """AMESim PNL0002 R-C-R pneumatic pipe with one center compliance.""" MODEL_TYPE = "amesim_pnl0002" - MODEL_VERSION = "0.3.0" + MODEL_VERSION = "0.5.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -964,6 +993,20 @@ class AmesimPnl0002(AmesimPnl0001): def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: self._connected_h = dict(connected_h) + def update_flow_temperature_references( + self, + connected_h: Mapping[str, float], + ) -> None: + self._connected_h = dict(connected_h) + + def state_derivative_from_ports( + self, + connected_h: Mapping[str, float], + ) -> list[float]: + # Junctions allocate their energy-balanced outlet enthalpy per port. + # The separate cache is only the temperature input to pn2pipefr. + return super().state_derivative_from_ports(connected_h) + def component_result_values(self) -> Mapping[str, float]: props = self.properties() flow_1 = self.port_mass_flow( @@ -978,10 +1021,45 @@ class AmesimPnl0002(AmesimPnl0001): props.T, port_name="port_2", ) - diagnostic_flow = flow_1 if abs(flow_1) >= abs(flow_2) else flow_2 - upstream_pressure = max(self.port_1.p, self.port_2.p, props.p, 1.0) - density = max(self.medium.density(upstream_pressure, props.T), 1.0e-12) - reynolds = self.reynolds_number(diagnostic_flow, props.T) + resistance_diagnostics: list[tuple[float, float, float, float]] = [] + for port_name, port, flow in ( + ("port_1", self.port_1, flow_1), + ("port_2", self.port_2, flow_2), + ): + if flow >= 0.0: + upstream_pressure = max(port.p, 1.0) + upstream_h = self._connected_h.get(port_name, props.h) + upstream_temperature = max( + self.medium.temperature_from_pressure_enthalpy( + upstream_pressure, + upstream_h, + ), + 1.0, + ) + else: + upstream_pressure = max(props.p, 1.0) + upstream_temperature = props.T + density = max( + self.medium.density(upstream_pressure, upstream_temperature), + 1.0e-12, + ) + reynolds = self.reynolds_number(flow, upstream_temperature) + resistance_diagnostics.append( + ( + reynolds, + ( + abs(flow) + * sqrt(upstream_temperature) + / max(self.area * upstream_pressure, 1.0e-18) + ), + abs(flow) / (density * self.area), + self.friction_factor(reynolds), + ) + ) + reynolds, cm, velocity, friction = ( + sum(values) / len(resistance_diagnostics) + for values in zip(*resistance_diagnostics) + ) return { "m": self.state.m, "U": self.state.U, @@ -991,9 +1069,9 @@ class AmesimPnl0002(AmesimPnl0001): "u": props.u, "h": props.h, "re": reynolds, - "cm": abs(diagnostic_flow) / max(self.area * upstream_pressure, 1.0e-18), - "v": diagnostic_flow / (density * self.area), - "ff": self.friction_factor(reynolds), + "cm": cm, + "v": velocity, + "ff": friction, } def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: @@ -1045,7 +1123,7 @@ class AmesimPnl0003(DynamicComponent): state_size = 4 MODEL_TYPE = "amesim_pnl0003" - 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"), @@ -1318,7 +1396,11 @@ class AmesimPnl0003(DynamicComponent): "h2": port_2.h, "dmctr": center_flow, "re": reynolds, - "cm": abs(center_flow) / max(self.area * max(port_1.p, port_2.p, 1.0), 1.0e-18), + "cm": ( + abs(center_flow) + * sqrt(upstream.T) + / max(self.area * max(port_1.p, port_2.p, 1.0), 1.0e-18) + ), "v": center_flow / (max(upstream.rho, 1.0e-12) * self.area), "ff": self.friction_factor(reynolds), } diff --git a/app/simulation/components/amesim/junctions/nodes.py b/app/simulation/components/amesim/junctions/nodes.py index 56315be..b9414dd 100644 --- a/app/simulation/components/amesim/junctions/nodes.py +++ b/app/simulation/components/amesim/junctions/nodes.py @@ -10,13 +10,20 @@ from app.simulation.core.ports import PortDefinition, PortState class _AmesimPneumaticNode(AlgebraicComponent): - """Shared implementation for AMESim pneumatic junction submodels.""" + """Shared implementation for AMESim pneumatic junction submodels. + + PN3NODE2/P4NODE2 use port 2 as their pressure and temperature reference. + Non-reference outlet ports use that reference temperature. When port 2 is + an outlet, its enthalpy is the residual that closes the junction energy + balance, matching the AMESim dh2 causality. + """ REFERENCE_PORT = "port_2" def __init__(self, name: str) -> None: super().__init__(name=name) self.set_parameter_values({}) + self.temperature_reference_h = 0.0 for definition in self.PORTS: setattr(self, definition.name, self.register_declared_port(definition.name)) @@ -58,6 +65,10 @@ class _AmesimPneumaticNode(AlgebraicComponent): return tuple(residuals) def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: + self.temperature_reference_h = connected_h.get( + self.REFERENCE_PORT, + sum(connected_h.values()) / len(connected_h) if connected_h else 0.0, + ) incoming = [ (port.m_flow, connected_h[name]) for name, port in self.ports.items() @@ -67,19 +78,37 @@ class _AmesimPneumaticNode(AlgebraicComponent): if total_flow > 1e-12: mixed_h = sum(m_flow * h for m_flow, h in incoming) / total_flow else: - mixed_h = connected_h.get( - self.REFERENCE_PORT, - sum(connected_h.values()) / len(connected_h) if connected_h else 0.0, + mixed_h = self.temperature_reference_h + + reference_port = self.get_port(self.REFERENCE_PORT) + for name, port in self.ports.items(): + port.h_outflow = ( + mixed_h + if name == self.REFERENCE_PORT + else self.temperature_reference_h + ) + + if reference_port.m_flow < -1e-12: + energy_without_reference = sum( + port.m_flow + * ( + connected_h[name] + if port.m_flow > 1e-12 + else self.temperature_reference_h + ) + for name, port in self.ports.items() + if name != self.REFERENCE_PORT + ) + reference_port.h_outflow = ( + -energy_without_reference / reference_port.m_flow ) - for port in self.ports.values(): - port.h_outflow = mixed_h class AmesimPn3Node2(_AmesimPneumaticNode): """AMESim PN3NODE2 pneumatic three-port junction.""" MODEL_TYPE = "amesim_pn3node2" - MODEL_VERSION = "0.1.0" + MODEL_VERSION = "0.3.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -115,7 +144,7 @@ class AmesimP4Node2(_AmesimPneumaticNode): """AMESim P4NODE2 pneumatic four-port junction.""" MODEL_TYPE = "amesim_p4node2" - MODEL_VERSION = "0.1.0" + MODEL_VERSION = "0.3.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), diff --git a/app/simulation/core/base.py b/app/simulation/core/base.py index 37dab30..02bfd2e 100644 --- a/app/simulation/core/base.py +++ b/app/simulation/core/base.py @@ -204,6 +204,19 @@ class Component(ABC): return None + def update_flow_temperature_references( + self, + connected_h: Mapping[str, float], + ) -> None: + """Update enthalpy references used only by pressure-flow laws. + + Most components use the normal stream enthalpy for both energy + transport and upstream-property evaluation. AMESim node submodels can + expose a distinct temperature reference, so the default is a no-op. + """ + + return None + def pneumatic_volume_outputs(self) -> Mapping[str, tuple[float, float]]: """Return directed ``volume``/``volume_flow`` values by pneumatic port. diff --git a/app/simulation/solvers/mechanical.py b/app/simulation/solvers/mechanical.py index fb7c1e6..404b50f 100644 --- a/app/simulation/solvers/mechanical.py +++ b/app/simulation/solvers/mechanical.py @@ -60,6 +60,20 @@ class MechanicalConstraintGroup: def _boundary_tolerance(bound: float) -> float: return 1.0e-12 * max(abs(bound), 1.0) + @staticmethod + def _velocity_tolerance(velocity: float) -> float: + """Treat only floating-point-scale motion as stationary at a stop. + + Implicit solvers perturb every state while constructing a numerical + Jacobian. Around an ideal endstop those perturbations must not switch + the unilateral constraint on and off; doing so turns a zero constrained + acceleration into the full outward-force acceleration across a + machine-scale velocity delta. The tolerance is deliberately far below + MECMAS21's physical ``dvel`` threshold so real release motion is kept. + """ + + return 1.0e-12 * max(abs(velocity), 1.0) + def reset_mode(self) -> None: self.mode = "uninitialized" @@ -131,6 +145,7 @@ class MechanicalConstraintGroup: def _static_endstop_side(self, total_force: float) -> str | None: position = self.representative.x velocity = self.representative.v + velocity_tolerance = self._velocity_tolerance(velocity) lower = self.lower_bound upper = self.upper_bound # MECMAS21's dvel is the friction stick threshold. Its discrete @@ -138,14 +153,14 @@ class MechanicalConstraintGroup: if ( lower is not None and position <= lower + self._boundary_tolerance(lower) - and velocity <= 0.0 + and velocity <= velocity_tolerance and total_force <= 0.0 ): return "lower" if ( upper is not None and position >= upper - self._boundary_tolerance(upper) - and velocity >= 0.0 + and velocity >= -velocity_tolerance and total_force >= 0.0 ): return "upper" @@ -489,6 +504,8 @@ class MechanicalStateReducer: position_index = velocity_index + 1 previous_velocity = float(previous_state[velocity_index]) current_velocity = float(current_state[velocity_index]) + previous_velocity_tolerance = group._velocity_tolerance(previous_velocity) + current_velocity_tolerance = group._velocity_tolerance(current_velocity) previous_position = float(previous_state[position_index]) current_position = float(current_state[position_index]) lower = group.lower_bound @@ -496,7 +513,7 @@ class MechanicalStateReducer: if ( lower is not None and previous_position <= lower + group._boundary_tolerance(lower) - and previous_velocity < 0.0 + and previous_velocity < -previous_velocity_tolerance ): candidates.append((previous_time, group, "lower", lower)) elif ( @@ -522,8 +539,8 @@ class MechanicalStateReducer: elif ( lower is not None and previous_position <= lower - and previous_velocity > 0.0 - and current_velocity < 0.0 + and previous_velocity > previous_velocity_tolerance + and current_velocity < -current_velocity_tolerance and current_position <= lower ): turnaround_time = self._locate_turnaround( @@ -551,7 +568,7 @@ class MechanicalStateReducer: if ( upper is not None and previous_position >= upper - group._boundary_tolerance(upper) - and previous_velocity > 0.0 + and previous_velocity > previous_velocity_tolerance ): candidates.append((previous_time, group, "upper", upper)) elif ( @@ -577,8 +594,8 @@ class MechanicalStateReducer: elif ( upper is not None and previous_position >= upper - and previous_velocity < 0.0 - and current_velocity > 0.0 + and previous_velocity < -previous_velocity_tolerance + and current_velocity > current_velocity_tolerance and current_position >= upper ): turnaround_time = self._locate_turnaround( diff --git a/app/simulation/solvers/stream.py b/app/simulation/solvers/stream.py index 573d4d1..c017260 100644 --- a/app/simulation/solvers/stream.py +++ b/app/simulation/solvers/stream.py @@ -63,6 +63,26 @@ class StreamResolver: values[endpoint.component][endpoint.port] = connected_port.h_outflow return values + def connected_temperature_reference_enthalpies( + self, + ) -> dict[str, dict[str, float]]: + """Return connector references used for upstream temperature only.""" + + values: dict[str, dict[str, float]] = { + component.name: {} for component in self.network.components.values() + } + for endpoint, connected in self._connected_endpoint.items(): + connected_component = self.network.components[connected.component] + connected_port = connected_component.get_port(connected.port) + values[endpoint.component][endpoint.port] = float( + getattr( + connected_component, + "temperature_reference_h", + connected_port.h_outflow, + ) + ) + return values + def solve(self) -> tuple[StreamSolveDiagnostics, dict[str, dict[str, float]]]: dynamic_components = [ component diff --git a/app/simulation/systems/generic.py b/app/simulation/systems/generic.py index f50ba17..20fa330 100644 --- a/app/simulation/systems/generic.py +++ b/app/simulation/systems/generic.py @@ -313,8 +313,14 @@ class GenericFluidSystem: for coupling_iteration in range(1, max_coupling_iterations + 1): previous_flows = tuple(port.m_flow for port in physical_ports) stream, connected_h = self.stream_resolver.solve() + temperature_reference_h = ( + self.stream_resolver.connected_temperature_reference_enthalpies() + ) for component in self.dynamic_components: component.update_stream_outflows(connected_h[component.name]) + component.update_flow_temperature_references( + temperature_reference_h[component.name] + ) algebraic = self.pressure_flow_solver.solve() pressure_flow_solve_count += 1 current_flows = tuple(port.m_flow for port in physical_ports) diff --git a/tests/test_amesim_pneumatic_node_components.py b/tests/test_amesim_pneumatic_node_components.py index 8f15ad3..b2d43d7 100644 --- a/tests/test_amesim_pneumatic_node_components.py +++ b/tests/test_amesim_pneumatic_node_components.py @@ -68,7 +68,7 @@ class AmesimPneumaticNodeComponentTests(unittest.TestCase): ), ) - def test_node_stream_outflow_uses_incoming_weighted_mix(self) -> None: + def test_node_uses_port_2_temperature_for_non_reference_outlets(self) -> None: node = AmesimP4Node2("p4_1") node.port_1.m_flow = 0.25 node.port_2.m_flow = 0.75 @@ -84,10 +84,44 @@ class AmesimPneumaticNodeComponentTests(unittest.TestCase): } ) - self.assertEqual(node.port_1.h_outflow, 250.0) + self.assertEqual(node.temperature_reference_h, 300.0) + self.assertEqual(node.port_1.h_outflow, 300.0) self.assertEqual(node.port_2.h_outflow, 250.0) - self.assertEqual(node.port_3.h_outflow, 250.0) - self.assertEqual(node.port_4.h_outflow, 250.0) + self.assertEqual(node.port_3.h_outflow, 300.0) + self.assertEqual(node.port_4.h_outflow, 300.0) + + def test_port_2_outlet_closes_node_energy_balance(self) -> None: + node = AmesimP4Node2("p4_1") + node.port_1.m_flow = -0.01823 + node.port_2.m_flow = -2.56077 + node.port_3.m_flow = 2.579 + node.port_4.m_flow = 0.0 + + node.update_stream_outflows( + { + "port_1": 2_630_000.0, + "port_2": 140_900.0, + "port_3": 63_960.0, + "port_4": 0.0, + } + ) + + self.assertEqual(node.port_1.h_outflow, 140_900.0) + energy_flow = sum( + port.m_flow + * ( + { + "port_1": 2_630_000.0, + "port_2": 140_900.0, + "port_3": 63_960.0, + "port_4": 0.0, + }[name] + if port.m_flow > 0.0 + else port.h_outflow + ) + for name, port in node.ports.items() + ) + self.assertAlmostEqual(energy_flow, 0.0, places=10) if __name__ == "__main__": diff --git a/tests/test_amesim_pnl0002_pnl0003_component.py b/tests/test_amesim_pnl0002_pnl0003_component.py index 98f1a97..3eee9d6 100644 --- a/tests/test_amesim_pnl0002_pnl0003_component.py +++ b/tests/test_amesim_pnl0002_pnl0003_component.py @@ -151,6 +151,23 @@ class AmesimPnl0002ComponentTests(unittest.TestCase): self.assertAlmostEqual(derivative[0], 0.1) + def test_connection_derivative_uses_energy_balanced_node_enthalpy(self) -> None: + pipe = AmesimPnl0002("pnl_2", self.medium) + props = pipe.properties() + pipe.port_1.m_flow = 0.2 + pipe.port_2.m_flow = 0.1 + reference_h = props.h + 10_000.0 + pipe.update_flow_temperature_references( + {"port_1": reference_h, "port_2": reference_h} + ) + + derivative = pipe.state_derivative_from_ports( + {"port_1": props.h, "port_2": props.h} + ) + + self.assertAlmostEqual(derivative[0], 0.3) + self.assertAlmostEqual(derivative[1], 0.3 * props.h) + class AmesimPnl0003ComponentTests(unittest.TestCase): def setUp(self) -> None: diff --git a/tests/test_amesim_pnl00r_component.py b/tests/test_amesim_pnl00r_component.py index a82c83a..a4e33b8 100644 --- a/tests/test_amesim_pnl00r_component.py +++ b/tests/test_amesim_pnl00r_component.py @@ -66,6 +66,50 @@ class AmesimPnl00rComponentTests(unittest.TestCase): self.assertGreater(pipe.friction_factor(100000.0), 0.0) self.assertLess(pipe.friction_factor(100000.0), 0.1) + def test_friction_factor_matches_amesim_smooth_to_rough_transition(self) -> None: + pipe_20mm = AmesimPnl00r( + "pnl_20mm", + self.medium, + diam=0.02, + rr=0.045 / 20.0, + ) + pipe_14mm = AmesimPnl00r( + "pnl_14mm", + self.medium, + diam=0.014, + rr=0.045 / 14.0, + ) + + self.assertAlmostEqual( + pipe_20mm.friction_factor(56_887.5547), + 0.0215740061043, + delta=8.0e-5, + ) + self.assertAlmostEqual( + pipe_20mm.friction_factor(700_686.41), + 0.02400535718, + delta=8.0e-5, + ) + self.assertAlmostEqual( + pipe_14mm.friction_factor(726_799.66), + 0.02658268645, + delta=8.0e-5, + ) + + def test_friction_factor_matches_amesim_transition_regime(self) -> None: + pipe = AmesimPnl00r( + "pnl_20mm", + self.medium, + diam=0.02, + rr=0.045 / 20.0, + ) + + self.assertAlmostEqual( + pipe.friction_factor(3_699.80236), + 0.0391257647759, + delta=4.0e-4, + ) + def test_darcy_pressure_drop_uses_flow_sign(self) -> None: pipe = AmesimPnl00r("pnl_1", self.medium) density = self.medium.density(500000.0, 300.0) diff --git a/tests/test_component_catalog.py b/tests/test_component_catalog.py index 1bce93c..e9ab7b6 100644 --- a/tests/test_component_catalog.py +++ b/tests/test_component_catalog.py @@ -465,7 +465,7 @@ class ComponentCatalogTests(unittest.TestCase): component = self.components[model_type] model_parameters = parameters(model_type) expected_version = ( - "0.3.0" if model_type == "amesim_pnl0002" else "0.2.0" + "0.5.0" if model_type == "amesim_pnl0002" else "0.3.0" ) self.assertEqual(component["modelVersion"], expected_version) self.assertEqual(model_parameters["mode"]["editor"], "choice") diff --git a/tests/test_mechanical_solver_causalization.py b/tests/test_mechanical_solver_causalization.py index 1afca47..fb2e0f0 100644 --- a/tests/test_mechanical_solver_causalization.py +++ b/tests/test_mechanical_solver_causalization.py @@ -396,8 +396,38 @@ class MechanicalSolverCausalizationTests(unittest.TestCase): self.assertAlmostEqual(result.series["mass.v"][0], -0.5 * mass.dvel) self.assertLess(min(result.series["mass.x"]), 0.0) self.assertLessEqual(max(result.series["mass.x"]), 1.0e-15) - self.assertAlmostEqual(result.series["mass.x"][-1], 0.0, places=15) - self.assertAlmostEqual(result.series["mass.v"][-1], 0.0, places=15) + self.assertAlmostEqual(result.series["mass.x"][-1], 0.0, delta=5.0e-15) + self.assertAlmostEqual(result.series["mass.v"][-1], 0.0, delta=1.1e-12) + + def test_ideal_upper_stop_ignores_jacobian_scale_inward_velocity_noise(self) -> None: + system, mass = _single_mass_system( + 100.0, + stoptype=1.0, + x0=0.0, + xmin=-1.0, + xmax=0.0, + ) + initial_state = system.consistent_initial_state_vector() + perturbed_state = [-1.0e-14, initial_state[1]] + + self.assertEqual(system.rhs(0.0, initial_state), [0.0, 0.0]) + self.assertEqual(system.rhs(0.0, perturbed_state), [0.0, 0.0]) + self.assertEqual(mass.v, -1.0e-14) + + def test_ideal_lower_stop_ignores_jacobian_scale_outward_velocity_noise(self) -> None: + system, mass = _single_mass_system( + -100.0, + stoptype=1.0, + x0=0.0, + xmin=0.0, + xmax=1.0, + ) + initial_state = system.consistent_initial_state_vector() + perturbed_state = [1.0e-14, initial_state[1]] + + self.assertEqual(system.rhs(0.0, initial_state), [0.0, 0.0]) + self.assertEqual(system.rhs(0.0, perturbed_state), [0.0, 0.0]) + self.assertEqual(mass.v, 1.0e-14) def test_ideal_stop_rejects_initial_position_outside_limits(self) -> None: system, _mass = _single_mass_system( diff --git a/tests/test_pressure_flow_solver_initialization.py b/tests/test_pressure_flow_solver_initialization.py index 4caff61..adacb56 100644 --- a/tests/test_pressure_flow_solver_initialization.py +++ b/tests/test_pressure_flow_solver_initialization.py @@ -183,9 +183,11 @@ class PressureFlowSolverInitializationTests(unittest.TestCase): first = system.rhs(0.0, state) second = system.rhs(0.0, state) - self.assertGreater(system.max_thermofluid_iterations, 1) + # PN3NODE2 takes its pressure-flow temperature reference from port 2, + # so that closure no longer depends on the flow-weighted energy mix. + self.assertGreaterEqual(system.max_thermofluid_iterations, 1) for first_value, second_value in zip(first, second): - self.assertAlmostEqual(first_value, second_value, places=10) + self.assertAlmostEqual(first_value, second_value, delta=1.0e-8) def test_current_storage_pressure_reseeds_stale_orifice_ports_and_flow(self) -> None: network, medium, high, low, valve = self._near_equal_pressure_network()