from __future__ import annotations from math import exp import unittest from app.simulation.components.amesim.mechanical.translational import ( AmesimF000, AmesimForc, AmesimLstp00a, AmesimMecmas21, ) from app.simulation.core.base import AlgebraicComponent from app.simulation.core.equations import EquationResidual from app.simulation.core.medium import IdealGasMedium from app.simulation.core.ports import PortDefinition from app.simulation.solvers.algebraic import PressureFlowSolver from app.simulation.solvers.mechanical import MechanicalConstraintGroup from app.simulation.solvers.solver import SolveIVPConfig from app.simulation.systems.generic import GenericFluidSystem from app.simulation.systems.network import SimulationNetwork class _AnchoredMechanicalForce(AlgebraicComponent): """Test boundary whose modest force must survive unrelated large scales.""" PORTS = (PortDefinition.mechanical_translational("port_1"),) def __init__(self, name: str, force: float) -> None: super().__init__(name=name) self.force = float(force) self.port_1 = self.register_declared_port("port_1") def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: return ( EquationResidual( id=f"{self.name}:x_state", owner="component", owner_id=self.name, relation="state", variables=(f"{self.name}.port_1.x",), role="effort", value=self.port_1.x, ), EquationResidual( id=f"{self.name}:v_state", owner="component", owner_id=self.name, relation="state", variables=(f"{self.name}.port_1.v",), role="effort", value=self.port_1.v, ), EquationResidual( id=f"{self.name}:force_state", owner="component", owner_id=self.name, relation="state", variables=(f"{self.name}.port_1.f",), role="flow", value=self.port_1.f - self.force, ), ) class _RigidMechanicalLink(AlgebraicComponent): """Massless link whose position constraint supplies the rigid coordinate.""" PORTS = ( PortDefinition.mechanical_translational("port_1"), PortDefinition.mechanical_translational("port_2"), ) def __init__(self, name: str) -> None: super().__init__(name=name) self.port_1 = self.register_declared_port("port_1") self.port_2 = self.register_declared_port("port_2") def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: return ( EquationResidual( id=f"{self.name}:x_equal", owner="component", owner_id=self.name, relation="equal", variables=(f"{self.name}.port_1.x", f"{self.name}.port_2.x"), role="effort", value=self.port_1.x - self.port_2.x, ), EquationResidual( id=f"{self.name}:v_equal", owner="component", owner_id=self.name, relation="equal", variables=(f"{self.name}.port_1.v", f"{self.name}.port_2.v"), role="effort", value=self.port_1.v - self.port_2.v, ), ) def _single_mass_system( applied_force: float, *, stoptype: float = 4.0, x0: float = 0.0, xmin: float = -1.0, xmax: float = 1.0, ) -> tuple[GenericFluidSystem, AmesimMecmas21]: medium = IdealGasMedium() source = AmesimForc("force") source.res.signal = applied_force mass = AmesimMecmas21( "mass", medium, mass=2.0, useFriction=0.0, stoptype=stoptype, x0=x0, v0=0.0, xmin=xmin, xmax=xmax, ) zero = AmesimF000("zero") network = SimulationNetwork("single-mass-causalization") for component in (source, mass, zero): network.add_component(component) network.connect("force", "port_2", "mass", "port_1") network.connect("mass", "port_2", "zero", "port_1") return GenericFluidSystem(network), mass class MechanicalSolverCausalizationTests(unittest.TestCase): def test_lstp_contact_uses_exponential_damping_ramp_and_negative_force_option( self, ) -> None: medium = IdealGasMedium() contact = AmesimLstp00a( "contact", medium, gap0=0.0, kcont=100.0, rcont=10.0, Pdis=0.1, discContactOption=1.0, ) contact.port_1.x = 0.1 contact.port_2.x = 0.0 contact.port_1.v = -2.0 contact.port_2.v = 0.0 expected = 10.0 - 20.0 * (1.0 - exp(-1.0)) self.assertAlmostEqual(contact.contact_force, expected, places=12) clipped = AmesimLstp00a( "clipped_contact", medium, gap0=0.0, kcont=100.0, rcont=10.0, Pdis=0.1, discContactOption=2.0, ) clipped.port_1.x = contact.port_1.x clipped.port_2.x = contact.port_2.x clipped.port_1.v = contact.port_1.v clipped.port_2.v = contact.port_2.v self.assertEqual(clipped.contact_force, 0.0) contact.set_causal_contact(penetration=0.1, force=expected) self.assertAlmostEqual(contact.contact_force, expected, places=12) clipped.set_causal_contact(penetration=0.1, force=expected) self.assertEqual(clipped.contact_force, 0.0) def test_causal_contact_survives_unrelated_nonlinear_fallback(self) -> None: medium = IdealGasMedium() source = AmesimForc("contact_force") source.res.signal = 40.0 contact = AmesimLstp00a( "contact", medium, gap0=0.0, kcont=1.0e11, rcont=0.0, Pdis=1.0e-7, discContactOption=1.0, ) mass = AmesimMecmas21( "mass", medium, mass=2.0, useFriction=0.0, x0=1.0e9, ) zero = AmesimF000("zero") unrelated = _AnchoredMechanicalForce("unrelated", 7.0) network = SimulationNetwork("causal-contact-with-nonlinear-fallback") for component in (source, contact, mass, zero, unrelated): network.add_component(component) network.connect("contact_force", "port_2", "contact", "port_1") network.connect("contact", "port_2", "mass", "port_1") network.connect("mass", "port_2", "zero", "port_1") diagnostics = PressureFlowSolver(network).solve() self.assertTrue(diagnostics.success, diagnostics.message) self.assertGreater(diagnostics.evaluations, 0) self.assertAlmostEqual(contact.penetration, 4.0e-10, places=20) self.assertAlmostEqual(contact.contact_force, 40.0, places=8) self.assertAlmostEqual(unrelated.port_1.f, 7.0, places=8) def test_elastic_mass_endstop_applies_contact_force_option(self) -> None: medium = IdealGasMedium() parameters = { "mass": 2.0, "useFriction": 0.0, "stoptype": 2.0, "x0": 0.1, "xmax": 0.0, "Kbmax": 100.0, "Dbmax": 10.0, "Pdmax": 0.01, "v0": -2.0, } negative_allowed = AmesimMecmas21( "negative_allowed", medium, discContactOption=1.0, **parameters, ) clipped = AmesimMecmas21( "clipped", medium, discContactOption=2.0, **parameters, ) self.assertAlmostEqual(negative_allowed._upper_limit_force(), -10.0) self.assertAlmostEqual(negative_allowed.acceleration(), 5.0) self.assertEqual(clipped._upper_limit_force(), 0.0) self.assertEqual(clipped.acceleration(), 0.0) def test_large_explicit_force_does_not_mask_small_local_force_residual(self) -> None: medium = IdealGasMedium() source = AmesimForc("large_force") source.res.signal = 1.0e17 mass = AmesimMecmas21( "large_mass", medium, mass=90_000.0, useFriction=0.0, ) zero = AmesimF000("large_zero") local_force = _AnchoredMechanicalForce("local_force", 40.0) network = SimulationNetwork("large-and-local-force-scales") for component in (source, mass, zero, local_force): network.add_component(component) network.connect("large_force", "port_2", "large_mass", "port_1") network.connect("large_mass", "port_2", "large_zero", "port_1") diagnostics = PressureFlowSolver(network).solve() self.assertTrue(diagnostics.success) self.assertEqual(source.port_2.f, -1.0e17) self.assertEqual(mass.port_1.f, 1.0e17) self.assertAlmostEqual(local_force.port_1.f, 40.0, places=9) residuals = { equation.id: equation.value for equation in network.pressure_flow_equation_residuals() } self.assertLess(abs(residuals["local_force:force_state"]), 1.0e-9) def test_rigidly_connected_masses_share_state_and_acceleration(self) -> None: medium = IdealGasMedium() source = AmesimForc("force") source.res.signal = 100.0 first_mass = AmesimMecmas21( "first_mass", medium, mass=2.0, useFriction=0.0, x0=0.25, v0=0.5, ) second_mass = AmesimMecmas21( "second_mass", medium, mass=3.0, useFriction=0.0, x0=0.25, v0=0.5, ) link = _RigidMechanicalLink("rigid_link") zero = AmesimF000("zero") network = SimulationNetwork("rigid-mass-group") for component in (source, first_mass, link, second_mass, zero): network.add_component(component) network.connect("force", "port_2", "first_mass", "port_1") network.connect("first_mass", "port_2", "rigid_link", "port_1") network.connect("rigid_link", "port_2", "second_mass", "port_1") network.connect("second_mass", "port_2", "zero", "port_1") system = GenericFluidSystem(network) initial_state = system.consistent_initial_state_vector() derivatives = system.rhs(0.0, initial_state) self.assertEqual(len(initial_state), 2) self.assertEqual(initial_state, [0.5, 0.25]) self.assertAlmostEqual(derivatives[0], 20.0, places=12) self.assertAlmostEqual(first_mass.acceleration(), 20.0, places=12) self.assertAlmostEqual(second_mass.acceleration(), 20.0, places=12) result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.01, max_step=0.001), sample_step=0.005, ) self.assertTrue(result.success, result.message) self.assertEqual(result.series["first_mass.x"], result.series["second_mass.x"]) self.assertEqual(result.series["first_mass.v"], result.series["second_mass.v"]) self.assertEqual(result.series["first_mass.a"], result.series["second_mass.a"]) for acceleration in result.series["first_mass.a"]: self.assertAlmostEqual(acceleration, 20.0, places=9) def test_ideal_upper_stop_locks_mass_under_outward_force(self) -> None: system, _mass = _single_mass_system( 100.0, stoptype=1.0, x0=0.0, xmin=-1.0, xmax=0.0, ) result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.01, max_step=0.001), sample_step=0.005, ) self.assertTrue(result.success, result.message) for value in result.series["mass.x"]: self.assertAlmostEqual(value, 0.0, places=12) for value in result.series["mass.v"]: self.assertAlmostEqual(value, 0.0, places=12) for value in result.series["mass.a"]: self.assertAlmostEqual(value, 0.0, places=12) # Ideal-contact reaction is an internal constraint force. AMESim's # Fmax output is reserved for the elastic (stoptype=2) endstop. self.assertEqual(result.series["mass.Fmax"], [0.0, 0.0, 0.0]) def test_ideal_upper_stop_releases_mass_under_inward_force(self) -> None: system, _mass = _single_mass_system( -100.0, stoptype=1.0, x0=0.0, xmin=-1.0, xmax=0.0, ) result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.01, max_step=0.001), sample_step=0.005, ) self.assertTrue(result.success, result.message) self.assertAlmostEqual(result.series["mass.a"][0], -50.0, places=12) self.assertEqual(result.series["mass.Fmax"], [0.0, 0.0, 0.0]) self.assertLess(result.series["mass.v"][-1], 0.0) self.assertLess(result.series["mass.x"][-1], 0.0) def test_ideal_upper_stop_releases_subthreshold_inward_velocity(self) -> None: system, mass = _single_mass_system( 100.0, stoptype=1.0, x0=0.0, xmin=-1.0, xmax=0.0, ) mass.v = -0.5 * mass.dvel mass.refresh_thermodynamic_ports() initial_state = system.consistent_initial_state_vector() derivatives = system.rhs(0.0, initial_state) self.assertAlmostEqual(derivatives[0], 50.0, places=12) self.assertAlmostEqual(derivatives[1], -0.5 * mass.dvel, places=18) result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=1.0e-6, max_step=1.0e-6), sample_step=5.0e-9, ) self.assertTrue(result.success, result.message) 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) def test_ideal_stop_rejects_initial_position_outside_limits(self) -> None: system, _mass = _single_mass_system( 100.0, stoptype=1.0, x0=0.01, xmin=-1.0, xmax=0.0, ) with self.assertRaisesRegex(ValueError, "outside the discrete endstop limits"): system.consistent_initial_state_vector() def test_rhs_trial_state_does_not_commit_ideal_stop_mode(self) -> None: system, _mass = _single_mass_system( 100.0, stoptype=1.0, x0=0.0, xmin=-1.0, xmax=0.01, ) initial_state = system.consistent_initial_state_vector() group = system.mechanical_state_reducer.groups[0] committed_mode = group.mode trial_derivatives = system.rhs(0.0, [0.0, 0.02]) self.assertEqual(trial_derivatives, [0.0, 0.0]) self.assertEqual(group.mode, committed_mode) accepted_derivatives = system.rhs(0.0, initial_state) self.assertEqual(group.mode, committed_mode) self.assertAlmostEqual(accepted_derivatives[0], 50.0, places=12) self.assertAlmostEqual(accepted_derivatives[1], 0.0, places=12) def test_ideal_upper_stop_projects_outward_velocity_at_step_start(self) -> None: system, mass = _single_mass_system( 100.0, stoptype=1.0, x0=0.01, xmin=-1.0, xmax=0.01, ) mass.v = 1.0 mass.refresh_thermodynamic_ports() result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.005, max_step=0.001), sample_step=0.001, ) self.assertTrue(result.success, result.message) for value in result.series["mass.x"]: self.assertAlmostEqual(value, 0.01, places=12) for value in result.series["mass.v"]: self.assertAlmostEqual(value, 0.0, places=12) for value in result.series["mass.a"]: self.assertAlmostEqual(value, 0.0, places=12) def test_ideal_upper_stop_projects_a_high_speed_impact(self) -> None: system, _mass = _single_mass_system( 100.0, stoptype=1.0, x0=0.0, xmin=-1.0, xmax=0.01, ) result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.04, max_step=0.01), sample_step=0.005, ) self.assertTrue(result.success, result.message) self.assertLessEqual(max(result.series["mass.x"]), 0.01 + 1.0e-12) after_impact = [ index for index, time in enumerate(result.series["time"]) if time >= 0.02 - 1.0e-10 ] self.assertTrue(after_impact) for index in after_impact: self.assertAlmostEqual(result.series["mass.x"][index], 0.01, places=12) self.assertAlmostEqual(result.series["mass.v"][index], 0.0, places=12) self.assertAlmostEqual(result.series["mass.a"][index], 0.0, places=12) def test_restitution_upper_stop_rebounds_a_high_speed_impact(self) -> None: system, mass = _single_mass_system( 0.0, stoptype=3.0, x0=-0.01, xmin=-1.0, xmax=0.0, ) mass.v = 2.0 mass.restdvel = 0.1 mass.restcoeff = 0.25 mass.refresh_thermodynamic_ports() result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.012, max_step=0.01), sample_step=0.002, ) self.assertTrue(result.success, result.message) self.assertLessEqual(max(result.series["mass.x"]), 1.0e-12) rebound_indices = [ index for index, velocity in enumerate(result.series["mass.v"]) if velocity < 0.0 ] self.assertTrue(rebound_indices) for index in rebound_indices: self.assertAlmostEqual(result.series["mass.v"][index], -0.5, places=12) self.assertAlmostEqual(result.series["mass.x"][-1], -0.0035, places=10) def test_restitution_upper_stop_locks_at_velocity_threshold(self) -> None: system, mass = _single_mass_system( 100.0, stoptype=3.0, x0=0.0, xmin=-1.0, xmax=0.0, ) mass.restdvel = 0.1 mass.restcoeff = 0.8 mass.v = mass.restdvel mass.refresh_thermodynamic_ports() result = system.simulate( SolveIVPConfig(t_start=0.0, t_stop=0.01, max_step=0.001), sample_step=0.005, ) self.assertTrue(result.success, result.message) for value in result.series["mass.x"]: self.assertAlmostEqual(value, 0.0, places=12) for value in result.series["mass.v"]: self.assertAlmostEqual(value, 0.0, places=12) for value in result.series["mass.a"]: self.assertAlmostEqual(value, 0.0, places=12) def test_plastic_stop_wins_over_restitution_at_shared_boundary(self) -> None: medium = IdealGasMedium() plastic = AmesimMecmas21( "plastic", medium, stoptype=1.0, xmin=-1.0, xmax=0.0, useFriction=0.0, ) restitution = AmesimMecmas21( "restitution", medium, stoptype=3.0, xmin=-1.0, xmax=0.0, restdvel=0.1, restcoeff=0.8, useFriction=0.0, ) group = MechanicalConstraintGroup((plastic, restitution)) self.assertEqual(group.impact_velocity("upper", 2.0), 0.0) if __name__ == "__main__": unittest.main()