diff --git a/app/main.py b/app/main.py index 96967f0..22dc59b 100644 --- a/app/main.py +++ b/app/main.py @@ -1217,6 +1217,12 @@ def compile_reactflow_network( media_by_component_id[node.id], parameter_values, ) + apply_layout_transform = getattr(component, "apply_layout_transform", None) + if apply_layout_transform is not None: + apply_layout_transform( + rotation=node.data.rotation, + mirrored=node.data.mirrored, + ) validate_component_port_interface(node, component.port_definitions) network.add_component(component) diff --git a/app/simulation/components/amesim/mechanical/pistons.py b/app/simulation/components/amesim/mechanical/pistons.py index 848cf99..0bd19d1 100644 --- a/app/simulation/components/amesim/mechanical/pistons.py +++ b/app/simulation/components/amesim/mechanical/pistons.py @@ -93,8 +93,8 @@ class AmesimPnrp17(AlgebraicComponent): ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_3", "left", order=20), - PortDisplaySpec("port_4", "left", order=30), - PortDisplaySpec("port_2", "right", order=40), + PortDisplaySpec("port_2", "left", order=30), + PortDisplaySpec("port_4", "right", order=40), PortDisplaySpec("port_5", "right", order=50), ), order=60, diff --git a/app/simulation/components/amesim/mechanical/translational.py b/app/simulation/components/amesim/mechanical/translational.py index 1c82b31..3d79b9f 100644 --- a/app/simulation/components/amesim/mechanical/translational.py +++ b/app/simulation/components/amesim/mechanical/translational.py @@ -87,6 +87,7 @@ class AmesimForc(AlgebraicComponent): self.set_parameter_values({}) self.res = self.register_declared_port("res") self.port_2 = self.register_declared_port("port_2") + self._orientation_sign = 1.0 @classmethod def create( @@ -98,6 +99,15 @@ class AmesimForc(AlgebraicComponent): ) -> "AmesimForc": return cls(name=name) + def apply_layout_transform(self, *, rotation: int, mirrored: bool) -> None: + """Apply the AMESim icon direction to the signed force output.""" + + normalized_rotation = int(rotation) % 360 + if normalized_rotation not in {0, 90, 180, 270}: + raise ValueError("FORC rotation must be a multiple of 90 degrees.") + direction = -1.0 if normalized_rotation in {180, 270} else 1.0 + self._orientation_sign = -direction if mirrored else direction + @property def output_force(self) -> float: return float(self.res.signal) @@ -111,7 +121,7 @@ class AmesimForc(AlgebraicComponent): relation="constitutive", variables=(f"{self.name}.port_2.f", f"{self.name}.res.signal"), role="flow", - value=self.port_2.f + self.output_force, + value=self.port_2.f + self._orientation_sign * self.output_force, ), ) @@ -170,8 +180,8 @@ class AmesimMecmas21(DynamicComponent): category_id="mechanical", symbol="amesim_mecmas21", ports=( - PortDisplaySpec("port_1", "left", order=10), - PortDisplaySpec("port_2", "right", order=20), + PortDisplaySpec("port_2", "left", order=10), + PortDisplaySpec("port_1", "right", order=20), ), order=30, ) @@ -426,11 +436,11 @@ class AmesimLstp00a(AlgebraicComponent): assert self._causal_port_2_x is not None penetration = ( self._causal_penetration - + (self.port_2.x - self._causal_port_2_x) - - (self.port_1.x - self._causal_port_1_x) + + (self.port_1.x - self._causal_port_1_x) + - (self.port_2.x - self._causal_port_2_x) ) return -penetration - return self.gap0 - (self.port_2.x - self.port_1.x) + return self.gap0 + (self.port_2.x - self.port_1.x) @property def penetration(self) -> float: @@ -438,7 +448,7 @@ class AmesimLstp00a(AlgebraicComponent): @property def penetration_velocity(self) -> float: - return self.port_2.v - self.port_1.v + return self.port_1.v - self.port_2.v @property def contact_force(self) -> float: @@ -516,7 +526,7 @@ class AmesimLstp00a(AlgebraicComponent): f"{self.name}.port_2.v", ), role="flow", - value=self.port_1.f + force, + value=self.port_1.f - force, ), EquationResidual( id=f"{self.name}:port_2_contact_force", @@ -531,7 +541,7 @@ class AmesimLstp00a(AlgebraicComponent): f"{self.name}.port_2.v", ), role="flow", - value=self.port_2.f - force, + value=self.port_2.f + force, ), ) diff --git a/app/simulation/core/errors.py b/app/simulation/core/errors.py new file mode 100644 index 0000000..bfec131 --- /dev/null +++ b/app/simulation/core/errors.py @@ -0,0 +1,5 @@ +from __future__ import annotations + + +class RecoverableTrialStateError(ValueError): + """A physical-domain failure caused by an integrator trial state.""" diff --git a/app/simulation/core/peng_robinson.py b/app/simulation/core/peng_robinson.py index f822189..51ea44f 100644 --- a/app/simulation/core/peng_robinson.py +++ b/app/simulation/core/peng_robinson.py @@ -1,5 +1,7 @@ from __future__ import annotations +from app.simulation.core.errors import RecoverableTrialStateError + from dataclasses import dataclass from math import acos, cos, isfinite, log, pi, sqrt @@ -69,7 +71,7 @@ class PengRobinsonFluid: def pressure_from_molar_volume(self, temperature: float, molar_volume: float) -> float: self._validate_temperature(temperature) if molar_volume <= self.b_parameter: - raise ValueError("Molar volume must be larger than Peng-Robinson b parameter.") + raise RecoverableTrialStateError("Molar volume must be larger than Peng-Robinson b parameter.") a_alpha = self.attractive_parameter(temperature) b = self.b_parameter repulsive = UNIVERSAL_GAS_CONSTANT * temperature / (molar_volume - b) diff --git a/app/simulation/solvers/algebraic.py b/app/simulation/solvers/algebraic.py index b8fc4c7..643977f 100644 --- a/app/simulation/solvers/algebraic.py +++ b/app/simulation/solvers/algebraic.py @@ -1,5 +1,7 @@ from __future__ import annotations +from collections.abc import Callable + from dataclasses import dataclass from math import expm1, isfinite, log, sqrt @@ -32,11 +34,22 @@ class AlgebraicUnknown: setattr(self.state, self.variable, float(value)) +@dataclass(frozen=True) +class ExplicitFlowAssignment: + equation_id: str + unknown: AlgebraicUnknown + evaluate: Callable[[], float] + +@dataclass(frozen=True) +class EffortAnchor: + unknown: AlgebraicUnknown + evaluate: Callable[[], float] + @dataclass(frozen=True) class EffortEqualityGroup: variable: str members: tuple[AlgebraicUnknown, ...] - anchors: tuple[tuple[AlgebraicUnknown, float], ...] + anchors: tuple[EffortAnchor, ...] @dataclass(frozen=True) @@ -85,6 +98,11 @@ class PressureFlowSolver: self.max_evaluations = max_evaluations self.unknowns = self._build_unknowns() self._unknowns_by_id = {unknown.id: unknown for unknown in self.unknowns} + self._effort_groups = { + variable: self._build_effort_equality_groups(variable) + for variable in ("p", "x", "v") + } + self._explicit_flow_plan = self._build_explicit_flow_plan() self.last_diagnostics: AlgebraicSolveDiagnostics | None = None def _build_unknowns(self) -> tuple[AlgebraicUnknown, ...]: @@ -132,7 +150,7 @@ class PressureFlowSolver: for variable in ("p", "x", "v"): self._seed_equal_effort(variable) - def _effort_equality_groups( + def _build_effort_equality_groups( self, variable: str, ) -> tuple[EffortEqualityGroup, ...]: @@ -195,10 +213,7 @@ class PressureFlowSolver: effort_unknowns[endpoint] ) - anchors_by_root: dict[ - tuple[str, str], - list[tuple[AlgebraicUnknown, float]], - ] = {} + anchors_by_root: dict[tuple[str, str], list[EffortAnchor]] = {} for equations in component_equations.values(): for equation in equations: if equation.relation != "state" or equation.role != "effort": @@ -215,11 +230,11 @@ class PressureFlowSolver: continue endpoint = endpoints[0] unknown = effort_unknowns[endpoint] - target_value = unknown.read() - float(equation.value) - if not isfinite(target_value): - continue anchors_by_root.setdefault(find(endpoint), []).append( - (unknown, target_value) + EffortAnchor( + unknown=unknown, + evaluate=self._equation_value_reader(equation), + ) ) return tuple( @@ -232,9 +247,17 @@ class PressureFlowSolver: ) def _seed_equal_effort(self, variable: str) -> None: - for group in self._effort_equality_groups(variable): + for group in self._effort_groups[variable]: members = group.members - anchors = group.anchors + anchors = tuple( + (anchor.unknown, anchor.unknown.read() - anchor.evaluate()) + for anchor in group.anchors + ) + anchors = tuple( + (unknown, value) + for unknown, value in anchors + if isfinite(value) + ) if anchors: # Keep each state-owned port current even when an invalid model # has conflicting anchors in one equality group. @@ -490,9 +513,9 @@ class PressureFlowSolver: gap0 = float(getattr(component, "gap0", 0.0)) if binding.algebraic_port == 1: - target = component.port_2.x - gap0 - penetration + target = component.port_2.x + gap0 + penetration else: - target = component.port_1.x + gap0 + penetration + target = component.port_1.x - gap0 - penetration for unknown in binding.algebraic_group.members: unknown.write(target) component.set_causal_contact( @@ -515,7 +538,7 @@ class PressureFlowSolver: position_groups = { unknown.id: group - for group in self._effort_equality_groups("x") + for group in self._effort_groups["x"] for unknown in group.members } bindings: list[UnilateralContactBinding] = [] @@ -547,7 +570,7 @@ class PressureFlowSolver: algebraic_group=first_group, neighbor_force=first_neighbor, algebraic_port=1, - force_sign=1.0, + force_sign=-1.0, ) elif not second_group.anchors and second_neighbor is not None: binding = UnilateralContactBinding( @@ -555,7 +578,7 @@ class PressureFlowSolver: algebraic_group=second_group, neighbor_force=second_neighbor, algebraic_port=2, - force_sign=-1.0, + force_sign=1.0, ) else: # With both coordinates state-owned, penetration is a dynamic @@ -574,137 +597,135 @@ class PressureFlowSolver: return tuple(bindings) - def _solve_explicit_flow_unknowns(self) -> set[str]: - """Directly evaluate explicit flow variables before nonlinear closure. + def _flow_unknowns_for_equation(self, equation) -> tuple[AlgebraicUnknown, ...]: + return tuple( + self._unknowns_by_id[variable] + for variable in equation.variables + if variable in self._unknowns_by_id + and self._unknowns_by_id[variable].role == "flow" + ) - Component constitutive equations use the normalized residual form - ``flow_unknown + remainder = 0`` whenever exactly one physical flow - variable is present. Solve those relations by substitution first, - then propagate the known values through component balances and physical - connectors. This covers pneumatic ``m_flow`` variables as well as - mechanical forces ``f`` such as ``FORC`` without asking the nonlinear - optimizer to discover values many orders of magnitude away from zero. + def _equation_value_reader(self, equation) -> Callable[[], float]: + if equation.owner == "connection": + if len(equation.variables) != 2: + raise ValueError( + f"Connection equation {equation.id} must contain two variables." + ) + first = self._unknowns_by_id[equation.variables[0]] + second = self._unknowns_by_id[equation.variables[1]] + if equation.relation == "sumToZero": + return lambda: first.read() + second.read() + if equation.relation == "equal": + return lambda: first.read() - second.read() + raise ValueError( + f"Unsupported connection equation relation: {equation.relation}." + ) - The remaining coupled equations still go through ``least_squares``; - these assignments provide both a consistent initial guess and the - nominal magnitudes used to scale that smaller nonlinear problem. - """ + component = self.network.components[equation.owner_id] + equation_id = equation.id + def read_component_equation() -> float: + for current in component.pressure_flow_equation_residuals(): + if current.id == equation_id: + return float(current.value) + raise RuntimeError( + f"Compiled algebraic equation disappeared at runtime: {equation_id}." + ) + + return read_component_equation + + def _build_explicit_flow_plan(self) -> tuple[ExplicitFlowAssignment, ...]: + """Compile the legacy deterministic flow assignment order once.""" + + assignments: list[ExplicitFlowAssignment] = [] seeded_ids: set[str] = set() - # Mechanical reaction balances can contain null-space forces. Reusing - # an arbitrary least-squares distribution from the preceding RHS call - # makes contact activation history-dependent, so choose deterministic - # zero tear values and rebuild the force chain from current signals, - # states, and pressure loads on every closure. - for unknown in self.unknowns: - if unknown.variable == "f": - unknown.write(0.0) + def append_assignment(equation, unknown: AlgebraicUnknown) -> None: + assignments.append( + ExplicitFlowAssignment( + equation_id=equation.id, + unknown=unknown, + evaluate=self._equation_value_reader(equation), + ) + ) + seeded_ids.add(unknown.id) - # First evaluate constitutive relations that expose one flow unknown - # with unit coefficient. Other variables in the equation (pressure, - # displacement, velocity, or a signal) have already been refreshed for - # the current state and time by the staged system closure. for component in self.network.components.values(): for equation in component.pressure_flow_equation_residuals(): if equation.relation != "constitutive" or equation.role != "flow": continue - flow_unknowns = [ - self._unknowns_by_id[variable] - for variable in equation.variables - if variable in self._unknowns_by_id - and self._unknowns_by_id[variable].role == "flow" - ] + flow_unknowns = self._flow_unknowns_for_equation(equation) if len(flow_unknowns) != 1: continue unknown = flow_unknowns[0] - if unknown.id in seeded_ids: - continue - target_value = unknown.read() - float(equation.value) - if not isfinite(target_value): - continue - unknown.write(target_value) - seeded_ids.add(unknown.id) + if unknown.id not in seeded_ids: + append_assignment(equation, unknown) - # V1/correctness-first implementation: repeatedly solve any balance that - # now has exactly one unknown flow variable left. Rebuilding and - # rescanning the complete residual tuple after every assignment keeps - # propagation deterministic, but costs O(flow unknowns * equations) and - # can dominate long, stiff simulations. A production follow-up should - # compile the assignment/tear order from the static topology once and - # evaluate only each owning component or connection residual here. + equations = self.network.pressure_flow_equation_residuals() while True: propagated = False - for equation in self.network.pressure_flow_equation_residuals(): + for equation in equations: if equation.role != "flow" or equation.relation not in { "constitutive", "sumToZero", }: continue - flow_unknowns = [ - self._unknowns_by_id[variable] - for variable in equation.variables - if variable in self._unknowns_by_id - and self._unknowns_by_id[variable].role == "flow" - ] + flow_unknowns = self._flow_unknowns_for_equation(equation) if not flow_unknowns: continue - variable_names = {unknown.variable for unknown in flow_unknowns} - if len(variable_names) != 1: + if len({unknown.variable for unknown in flow_unknowns}) != 1: continue - unseeded = [ - unknown for unknown in flow_unknowns if unknown.id not in seeded_ids - ] + unseeded = tuple( + unknown + for unknown in flow_unknowns + if unknown.id not in seeded_ids + ) if len(unseeded) != 1: continue - unknown = unseeded[0] - target_value = unknown.read() - float(equation.value) - if not isfinite(target_value): + append_assignment(equation, unseeded[0]) + propagated = True + break + if propagated: + continue + + for equation in equations: + if equation.role != "flow" or equation.relation not in { + "constitutive", + "sumToZero", + }: continue - unknown.write(target_value) - seeded_ids.add(unknown.id) + flow_unknowns = self._flow_unknowns_for_equation(equation) + unseeded = tuple( + unknown + for unknown in flow_unknowns + if unknown.id not in seeded_ids + ) + if len(unseeded) <= 1: + continue + if len({unknown.variable for unknown in flow_unknowns}) != 1: + continue + append_assignment(equation, unseeded[-1]) propagated = True break - if not propagated: - # Causalize one remaining free flow in an otherwise normalized - # linear balance. This is the algebraic equivalent of choosing - # a tear variable: the other free flows retain their current - # guesses and one dependent flow closes the equation exactly. - # It also gives rank-deficient rigid-body reaction balances a - # deterministic starting point before state reduction supplies - # their common acceleration. - for equation in self.network.pressure_flow_equation_residuals(): - if equation.role != "flow" or equation.relation not in { - "constitutive", - "sumToZero", - }: - continue - flow_unknowns = [ - self._unknowns_by_id[variable] - for variable in equation.variables - if variable in self._unknowns_by_id - and self._unknowns_by_id[variable].role == "flow" - ] - unseeded = [ - unknown - for unknown in flow_unknowns - if unknown.id not in seeded_ids - ] - if len(unseeded) <= 1: - continue - if len({unknown.variable for unknown in flow_unknowns}) != 1: - continue - unknown = unseeded[-1] - target_value = unknown.read() - float(equation.value) - if not isfinite(target_value): - continue - unknown.write(target_value) - seeded_ids.add(unknown.id) - propagated = True - break if not propagated: break + return tuple(assignments) + + def _solve_explicit_flow_unknowns(self) -> set[str]: + """Execute the precompiled explicit flow/force causalization plan.""" + + for unknown in self.unknowns: + if unknown.variable == "f": + unknown.write(0.0) + + seeded_ids: set[str] = set() + for assignment in self._explicit_flow_plan: + target_value = assignment.unknown.read() - assignment.evaluate() + if not isfinite(target_value): + continue + assignment.unknown.write(target_value) + seeded_ids.add(assignment.unknown.id) return seeded_ids def _scales(self) -> dict[str, float]: diff --git a/app/simulation/solvers/solver.py b/app/simulation/solvers/solver.py index a5c6029..389649e 100644 --- a/app/simulation/solvers/solver.py +++ b/app/simulation/solvers/solver.py @@ -1,5 +1,7 @@ from __future__ import annotations +from app.simulation.core.errors import RecoverableTrialStateError + import math from dataclasses import dataclass from typing import Callable, Literal, Sequence @@ -535,11 +537,10 @@ def _integrate_scipy_stepwise( ) -> ODESolution: """Initial stepwise integration path for breakpoints and state resets. - Known V1 limitation: an adaptive solver can evaluate a trial state outside - the algebraic or thermodynamic model domain. Such an RHS exception still - aborts the run here; recoverable trial failures are not yet restored to the - last accepted state and retried with a smaller step. This is not specific - to BDF, although implicit Newton/Jacobian probes make it especially visible. + Recoverable physical-domain failures from rejected integrator trial states + restore the last accepted state and rebuild the same solver with a smaller + maximum/first step. Structural, algebraic, and ordinary model errors still + fail immediately. """ import numpy as np from scipy.integrate import BDF, DOP853, LSODA, RK23, RK45, Radau @@ -609,6 +610,9 @@ def _integrate_scipy_stepwise( math.nextafter(segment_end, -math.inf) if is_breakpoint else segment_end ) has_integration_interval = integration_end > last_accepted_time + segment_max_step = float(config.max_step) + recoverable_retry_count = 0 + last_recoverable_error: RecoverableTrialStateError | None = None while has_integration_interval and last_accepted_time < integration_end: if cancel_check(): @@ -619,11 +623,16 @@ def _integrate_scipy_stepwise( solver_options = { "rtol": config.rtol, "atol": config.atol, - "max_step": config.max_step, + "max_step": segment_max_step, } - if config.first_step is not None: + requested_first_step = ( + 0.1 * segment_max_step + if last_recoverable_error is not None + else config.first_step + ) + if requested_first_step is not None: solver_options["first_step"] = min( - config.first_step, + requested_first_step, integration_end - last_accepted_time, ) @@ -639,6 +648,18 @@ def _integrate_scipy_stepwise( status = "cancelled" message = cancellation_message() break + except RecoverableTrialStateError as exc: + recoverable_retry_count += 1 + last_recoverable_error = exc + next_step = 0.5 * segment_max_step + minimum_step = 64.0 * math.ulp(max(abs(last_accepted_time), 1.0)) + if recoverable_retry_count > 16 or next_step <= minimum_step: + status = "failed" + message = str(exc) + error = exc + break + segment_max_step = next_step + continue except Exception as exc: status = "failed" message = str(exc) @@ -646,6 +667,7 @@ def _integrate_scipy_stepwise( break restart_at_transition = False + restart_after_recoverable = False while solver.status == "running": if cancel_check(): status = "cancelled" @@ -664,6 +686,22 @@ def _integrate_scipy_stepwise( "Simulation was stopped before reaching the requested end time." ) break + except RecoverableTrialStateError as exc: + recoverable_retry_count += 1 + last_recoverable_error = exc + attempted_step = segment_max_step + next_step = 0.5 * attempted_step + minimum_step = 64.0 * math.ulp( + max(abs(last_accepted_time), 1.0) + ) + if recoverable_retry_count > 16 or next_step <= minimum_step: + status = "failed" + message = str(exc) + error = exc + break + segment_max_step = next_step + restart_after_recoverable = True + break except Exception as exc: status = "failed" message = str(exc) @@ -672,6 +710,19 @@ def _integrate_scipy_stepwise( integration_progressed = True if solver.status == "failed": + if last_recoverable_error is not None: + recoverable_retry_count += 1 + next_step = 0.5 * segment_max_step + minimum_step = 64.0 * math.ulp( + max(abs(last_accepted_time), 1.0) + ) + if ( + recoverable_retry_count <= 16 + and next_step > minimum_step + ): + segment_max_step = next_step + restart_after_recoverable = True + break status = "failed" message = str(step_message or "Integration step failed.") break @@ -775,6 +826,7 @@ def _integrate_scipy_stepwise( last_accepted_time = step_end_time last_accepted_state = step_end_state + recoverable_retry_count = 0 reported_time = ( float(segment_end) if is_breakpoint and solver.status == "finished" @@ -806,7 +858,11 @@ def _integrate_scipy_stepwise( ) report_step(reported_time) - if status != "completed" or not restart_at_transition: + if status != "completed": + break + if restart_after_recoverable: + continue + if not restart_at_transition: break if status != "completed": diff --git a/frontend/src/App.tsx b/frontend/src/App.tsx index bb25286..e2d8b0b 100644 --- a/frontend/src/App.tsx +++ b/frontend/src/App.tsx @@ -5967,7 +5967,7 @@ function normalizeLoadedPorts( return definition?.ports.map((port) => ({ ...port })) ?? []; } - return rawPorts.flatMap((rawPort, index) => { + const loadedPorts = rawPorts.flatMap((rawPort, index) => { const record = rawPort !== null && typeof rawPort === "object" ? (rawPort as Record) @@ -6004,6 +6004,21 @@ function normalizeLoadedPorts( }, ]; }); + if (!definition) { + return loadedPorts; + } + + const loadedByName = new Map(loadedPorts.map((port) => [port.name, port])); + const registeredNames = new Set(definition.ports.map((port) => port.name)); + return [ + ...definition.ports.map((registered) => ({ + ...(loadedByName.get(registered.name) ?? registered), + // Port placement is versioned catalog metadata, just like the symbol. + // Prefer it so projects saved with an incorrect display order are upgraded. + side: registered.side, + })), + ...loadedPorts.filter((port) => !registeredNames.has(port.name)), + ]; } function isPortNominalRole(value: unknown): value is PortNominalRole { diff --git a/tests/test_amesim_mechanical_public_components.py b/tests/test_amesim_mechanical_public_components.py index 3f74bec..91d5419 100644 --- a/tests/test_amesim_mechanical_public_components.py +++ b/tests/test_amesim_mechanical_public_components.py @@ -17,6 +17,12 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase): def setUp(self) -> None: self.medium = IdealGasMedium() + def test_mecmas21_catalog_ports_follow_amesim_icon_sides(self) -> None: + self.assertEqual( + tuple((port.name, port.side) for port in AmesimMecmas21.DISPLAY.ports), + (("port_2", "left"), ("port_1", "right")), + ) + def test_f000_constrains_mechanical_port_force_to_zero(self) -> None: source = AmesimF000("zero_1") source.port_1.f = 12.5 @@ -38,6 +44,14 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase): self.assertAlmostEqual(residuals[0].value, 0.0) self.assertEqual(converter.component_result_values(), {"force": 20.0}) + converter.apply_layout_transform(rotation=180, mirrored=False) + converter.port_2.f = 20.0 + self.assertAlmostEqual( + converter.pressure_flow_equation_residuals()[0].value, + 0.0, + ) + self.assertEqual(converter.component_result_values(), {"force": 20.0}) + def test_mecmas21_acceleration_uses_connected_port_forces(self) -> None: mass = AmesimMecmas21( "mass_1", @@ -96,12 +110,12 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase): kcont=1000.0, rcont=10.0, ) - contact.port_1.x = 0.0 - contact.port_2.x = 0.002 - contact.port_1.v = 0.0 - contact.port_2.v = 0.1 - contact.port_1.f = -3.0 - contact.port_2.f = 3.0 + contact.port_1.x = 0.002 + contact.port_2.x = 0.0 + contact.port_1.v = 0.1 + contact.port_2.v = 0.0 + contact.port_1.f = 3.0 + contact.port_2.f = -3.0 self.assertAlmostEqual(contact.gap, -0.002) self.assertAlmostEqual(contact.penetration, 0.002) @@ -116,8 +130,8 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase): def test_lstp00a_returns_zero_before_contact(self) -> None: contact = AmesimLstp00a("contact_1", self.medium, gap0=0.001, kcont=1000.0) - contact.port_1.x = 0.0 - contact.port_2.x = 0.0005 + contact.port_1.x = 0.0005 + contact.port_2.x = 0.0 self.assertAlmostEqual(contact.gap, 0.0005) self.assertAlmostEqual(contact.contact_force, 0.0) diff --git a/tests/test_amesim_mechanical_xml.py b/tests/test_amesim_mechanical_xml.py index 57237c1..8cc0ff8 100644 --- a/tests/test_amesim_mechanical_xml.py +++ b/tests/test_amesim_mechanical_xml.py @@ -146,9 +146,9 @@ def signal_force_mass_project() -> ReactFlowProjectPayload: def elastic_contact_project() -> ReactFlowProjectPayload: left_parameters = dict(MECMAS21_DEFAULTS) - left_parameters["x0"] = 0.0 + left_parameters["x0"] = 0.001 right_parameters = dict(MECMAS21_DEFAULTS) - right_parameters["x0"] = 0.001 + right_parameters["x0"] = 0.0 return ReactFlowProjectPayload( name="amesim-mechanical-elastic-contact-smoke", nodes=[ @@ -260,10 +260,10 @@ class AmesimMechanicalXmlTests(unittest.TestCase): self.assertTrue(result["success"], result["message"]) self.assertEqual(result["series"]["time"], [0.0, 0.01, 0.02]) self.assertAlmostEqual(result["series"]["contact_1.force"][0], 1.0) - self.assertAlmostEqual(result["series"]["mass_left.a"][0], 0.5) - self.assertAlmostEqual(result["series"]["mass_right.a"][0], -0.5) - self.assertGreater(result["series"]["mass_left.x"][-1], 0.0) - self.assertLess(result["series"]["mass_right.x"][-1], 0.001) + self.assertAlmostEqual(result["series"]["mass_left.a"][0], -0.5) + self.assertAlmostEqual(result["series"]["mass_right.a"][0], 0.5) + self.assertLess(result["series"]["mass_left.x"][-1], 0.001) + self.assertGreater(result["series"]["mass_right.x"][-1], 0.0) def test_zero_force_mechanical_project_compiles_and_simulates(self) -> None: xml = build_reactflow_system_xml(zero_force_mass_project()) diff --git a/tests/test_amesim_pnrp17_xml.py b/tests/test_amesim_pnrp17_xml.py index d00efb5..3db3504 100644 --- a/tests/test_amesim_pnrp17_xml.py +++ b/tests/test_amesim_pnrp17_xml.py @@ -79,9 +79,9 @@ def pnrp17_coupled_project() -> ReactFlowProjectPayload: "amesim_pnrp17", [ _pneumatic_port("port_1", "left"), - mechanical_port("port_2", "right"), mechanical_port("port_3", "left"), - mechanical_port("port_4", "left"), + mechanical_port("port_2", "left"), + mechanical_port("port_4", "right"), mechanical_port("port_5", "right"), ], {"gi": 0.0, "dp": 0.1, "dr": 0.02, "x0": 0.0}, @@ -131,6 +131,18 @@ def pnrp17_coupled_project() -> ReactFlowProjectPayload: class AmesimPnrp17Tests(unittest.TestCase): + def test_catalog_ports_follow_amesim_pnrp17_icon_sides(self) -> None: + self.assertEqual( + tuple((port.name, port.side) for port in AmesimPnrp17.DISPLAY.ports), + ( + ("port_1", "left"), + ("port_3", "left"), + ("port_2", "left"), + ("port_4", "right"), + ("port_5", "right"), + ), + ) + def test_component_equations_match_pnrp17_geometry_and_signs(self) -> None: medium = IdealGasMedium() piston = AmesimPnrp17( diff --git a/tests/test_component_catalog.py b/tests/test_component_catalog.py index c8e95c6..68b9895 100644 --- a/tests/test_component_catalog.py +++ b/tests/test_component_catalog.py @@ -163,7 +163,7 @@ class ComponentCatalogTests(unittest.TestCase): for parameter in components["amesim_mecmas21"]["parameters"] } self.assertEqual(components["amesim_mecmas21"]["category"]["id"], "mechanical") - self.assertEqual([port["name"] for port in components["amesim_mecmas21"]["ports"]], ["port_1", "port_2"]) + self.assertEqual([port["name"] for port in components["amesim_mecmas21"]["ports"]], ["port_2", "port_1"]) self.assertEqual(mecmas_parameters["mass"]["unit"], "kg") self.assertEqual(mecmas_parameters["Kbmin"]["unit"], "N/m") lstp_parameters = { @@ -185,7 +185,7 @@ class ComponentCatalogTests(unittest.TestCase): for parameter in components["amesim_pnrp17"]["parameters"] } self.assertEqual(components["amesim_pnrp17"]["category"]["id"], "mechanical") - self.assertEqual([port["name"] for port in components["amesim_pnrp17"]["ports"]], ["port_1", "port_3", "port_4", "port_2", "port_5"]) + self.assertEqual([port["name"] for port in components["amesim_pnrp17"]["ports"]], ["port_1", "port_3", "port_2", "port_4", "port_5"]) self.assertEqual(pnrp_parameters["dp"]["unit"], "m") self.assertEqual(lmechn_parameters["v1"]["maximum"], 8.0) self.assertEqual(components["amesim_pnvo001"]["category"]["id"], "flow") diff --git a/tests/test_contact_solver_causalization.py b/tests/test_contact_solver_causalization.py index b999cbb..b44453a 100644 --- a/tests/test_contact_solver_causalization.py +++ b/tests/test_contact_solver_causalization.py @@ -58,7 +58,7 @@ class _PressureCoupledMechanicalLoad(AlgebraicComponent): f"{self.name}.pneumatic.p", ), role="flow", - value=self.mechanical.f - self.pneumatic.p, + value=self.mechanical.f + self.pneumatic.p, ), EquationResidual( id=f"{self.name}:pressure_closure", @@ -163,7 +163,7 @@ class ContactSolverCausalizationTests(unittest.TestCase): self.assertTrue(diagnostics.success, diagnostics.message) self.assertGreater(diagnostics.evaluations, 0) self.assertAlmostEqual(load.pneumatic.p, 41.0, delta=1.0e-3) - self.assertAlmostEqual(load.mechanical.f, 41.0, delta=1.0e-3) + self.assertAlmostEqual(load.mechanical.f, -41.0, delta=1.0e-3) self.assertAlmostEqual(contact.contact_force, 41.0, delta=1.0e-3) self.assertAlmostEqual(contact.penetration, 4.1e-10, delta=1.0e-14) @@ -172,7 +172,7 @@ class ContactSolverCausalizationTests(unittest.TestCase): expected_force = 10.0 - 20.0 * (1.0 - exp(-1.0)) load = _PrescribedMechanicalLoad( "load", - force=expected_force, + force=-expected_force, displacement=0.0, velocity=0.0, ) @@ -191,8 +191,8 @@ class ContactSolverCausalizationTests(unittest.TestCase): mass=1.0, useFriction=0.0, stoptype=4.0, - x0=0.08, - v0=-2.0, + x0=-0.08, + v0=2.0, ) zero = AmesimF000("zero") diff --git a/tests/test_core_solver.py b/tests/test_core_solver.py index 1450f7b..585a4db 100644 --- a/tests/test_core_solver.py +++ b/tests/test_core_solver.py @@ -4,6 +4,7 @@ import types import unittest from unittest.mock import patch +from app.simulation.core.errors import RecoverableTrialStateError from app.simulation.solvers.solver import ( SolveIVPConfig, StateTransition, @@ -71,6 +72,81 @@ class IntegrateOdeTests(unittest.TestCase): self.assertEqual(calls[0]["max_step"], 1.0e-3) self.assertNotIn("first_step", calls[0]) + def test_stepwise_solver_rebuilds_after_recoverable_trial_failure(self) -> None: + import numpy as np + import scipy.integrate + + attempted_max_steps: list[float] = [] + + class RetryBdf: + def __init__(self, fun, t0, y0, t_bound, **kwargs): + self.fun = fun + self.t = float(t0) + self.y = np.asarray(y0, dtype=float) + self.t_bound = float(t_bound) + self.h_abs = float(kwargs["max_step"]) + self.status = "running" + attempted_max_steps.append(self.h_abs) + + def step(self): + if self.h_abs > 0.25: + raise RecoverableTrialStateError("trial state outside domain") + self.t = self.t_bound + self.status = "finished" + return None + + def dense_output(self): + state = self.y.copy() + return lambda _time: state.copy() + + with patch.object(scipy.integrate, "BDF", RetryBdf): + result = integrate_ode( + rhs=lambda _time, _state: [0.0], + initial_state=[1.0], + config=SolveIVPConfig( + t_start=0.0, + t_stop=1.0, + method="BDF", + max_step=1.0, + ), + t_eval=[0.0, 1.0], + cancel_check=lambda: False, + ) + + self.assertTrue(result.success, result.message) + self.assertEqual(attempted_max_steps, [1.0, 0.5, 0.25]) + self.assertEqual(result.t, [0.0, 1.0]) + self.assertEqual(result.y, [[1.0, 1.0]]) + + def test_stepwise_solver_does_not_retry_ordinary_model_errors(self) -> None: + import numpy as np + import scipy.integrate + + attempts = 0 + + class FailingBdf: + def __init__(self, fun, t0, y0, t_bound, **kwargs): + nonlocal attempts + attempts += 1 + self.t = float(t0) + self.y = np.asarray(y0, dtype=float) + self.status = "running" + + def step(self): + raise ValueError("structural model error") + + with patch.object(scipy.integrate, "BDF", FailingBdf): + result = integrate_ode( + rhs=lambda _time, _state: [0.0], + initial_state=[1.0], + config=SolveIVPConfig(t_start=0.0, t_stop=1.0, method="BDF"), + cancel_check=lambda: False, + ) + + self.assertFalse(result.success) + self.assertEqual(result.message, "structural model error") + self.assertEqual(attempts, 1) + def test_scipy_stepwise_solver_can_cancel_before_start(self) -> None: result = integrate_ode( rhs=lambda _time, state: state, diff --git a/tests/test_mechanical_solver_causalization.py b/tests/test_mechanical_solver_causalization.py index 325fb31..0078bd6 100644 --- a/tests/test_mechanical_solver_causalization.py +++ b/tests/test_mechanical_solver_causalization.py @@ -144,10 +144,10 @@ class MechanicalSolverCausalizationTests(unittest.TestCase): Pdis=0.1, discContactOption=1.0, ) - contact.port_1.x = 0.0 - contact.port_2.x = 0.1 - contact.port_1.v = 0.0 - contact.port_2.v = -2.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) @@ -175,7 +175,7 @@ class MechanicalSolverCausalizationTests(unittest.TestCase): def test_causal_contact_survives_unrelated_nonlinear_fallback(self) -> None: medium = IdealGasMedium() source = AmesimForc("contact_force") - source.res.signal = -40.0 + source.res.signal = 40.0 contact = AmesimLstp00a( "contact", medium,