from __future__ import annotations from dataclasses import dataclass from app.simulation.examples.test_mql.primitives.mechanical import ( AmesimElasticEndstop, AmesimMassFrictionEndstops, AmesimPistonGeometry, circular_area, mm_to_m, ) from app.simulation.reporting.amesim_results import AmesimResults from app.simulation.reporting.test_mql_variables import ( TestMqlVariableCatalog, build_test_mql_variable_catalog, ) from app.simulation.examples.test_mql.config import TestMqlConfig, TestMqlResolvedComponent MM_TO_M = 1.0e-3 N_PER_MM_TO_N_PER_M = 1.0e3 N_PER_MM_PER_S_TO_N_PER_M_PER_S = 1.0e3 @dataclass(frozen=True) class TestMqlPistonSpec: alias: str piston_diameter_m: float rod_diameter_m: float zero_displacement_m: float piston_area_m2: float rod_area_m2: float annulus_area_m2: float data_paths: tuple[str, ...] def geometry(self) -> AmesimPistonGeometry: return AmesimPistonGeometry( piston_diameter_m=self.piston_diameter_m, rod_diameter_m=self.rod_diameter_m, zero_length_m=self.zero_displacement_m, ) @dataclass(frozen=True) class TestMqlMassEndstopSpec: alias: str mass_kg: float xmin_m: float xmax_m: float min_stiffness_n_per_m: float max_stiffness_n_per_m: float min_damping_n_per_m_per_s: float max_damping_n_per_m_per_s: float min_penetration_m: float max_penetration_m: float stiction_force_n: float coulomb_friction_n: float viscous_friction_n_per_m_per_s: float windage_n_per_m2_per_s2: float stick_velocity_threshold_m_s: float reset_velocity_threshold_m_s: float rest_coeff: float stribeck_constant_m_s: float use_friction: bool stop_type: int initial_velocity_m_s: float initial_displacement_m: float data_paths: tuple[str, ...] def endstop(self) -> AmesimMassFrictionEndstops: return AmesimMassFrictionEndstops( mass_kg=self.mass_kg, lower_limit_m=self.xmin_m, upper_limit_m=self.xmax_m, lower_stiffness_n_per_m=self.min_stiffness_n_per_m, upper_stiffness_n_per_m=self.max_stiffness_n_per_m, lower_damping_n_per_m_per_s=self.min_damping_n_per_m_per_s, upper_damping_n_per_m_per_s=self.max_damping_n_per_m_per_s, viscous_friction_n_per_m_per_s=self.viscous_friction_n_per_m_per_s, coulomb_friction_n=self.coulomb_friction_n, stiction_force_n=self.stiction_force_n, windage_n_per_m2_per_s2=self.windage_n_per_m2_per_s2, ) @dataclass(frozen=True) class TestMqlElasticEndstopSpec: alias: str gap_m: float contact_stiffness_n_per_m: float contact_damping_n_per_m_per_s: float spring_diameter_m: float wire_diameter_m: float data_paths: tuple[str, ...] def endstop(self) -> AmesimElasticEndstop: return AmesimElasticEndstop( contact_stiffness_n_per_m=self.contact_stiffness_n_per_m, contact_damping_n_per_m_per_s=self.contact_damping_n_per_m_per_s, gap0_m=self.gap_m, ) @dataclass(frozen=True) class TestMqlMechanicalNodeSpec: alias: str port_count: int sum_mode: int data_paths: tuple[str, ...] @dataclass(frozen=True) class TestMqlPiecewiseLinearSignalSpec: alias: str t_start_s: float starts: tuple[float, ...] ends: tuple[float, ...] durations_s: tuple[float, ...] stage_count: int is_cyclic: bool data_paths: tuple[str, ...] def output_at(self, time_s: float) -> float: if self.stage_count <= 0: return 0.0 elapsed = max(time_s - self.t_start_s, 0.0) active_durations = self.durations_s[: self.stage_count] total_duration = sum(active_durations) if self.is_cyclic and total_duration > 0.0: elapsed = elapsed % total_duration stage_start_time = 0.0 for index, duration in enumerate(active_durations): stage_end_time = stage_start_time + duration if elapsed < stage_end_time or index == self.stage_count - 1: if duration <= 0.0: return self.ends[index] fraction = (elapsed - stage_start_time) / duration return self.starts[index] + fraction * (self.ends[index] - self.starts[index]) stage_start_time = stage_end_time return self.ends[self.stage_count - 1] @dataclass(frozen=True) class TestMqlForceConnectorSpec: alias: str signal_alias: str target_mass_alias: str data_paths: tuple[str, ...] def force_at( self, time_s: float, signals: dict[str, TestMqlPiecewiseLinearSignalSpec], ) -> float: return signals[self.signal_alias].output_at(time_s) @dataclass(frozen=True) class TestMqlMechanicalAssembly: pistons: dict[str, TestMqlPistonSpec] masses: dict[str, TestMqlMassEndstopSpec] elastic_endstops: dict[str, TestMqlElasticEndstopSpec] mechanical_nodes: dict[str, TestMqlMechanicalNodeSpec] piecewise_signals: dict[str, TestMqlPiecewiseLinearSignalSpec] force_connectors: dict[str, TestMqlForceConnectorSpec] zero_force_sources: tuple[str, ...] @property def component_count(self) -> int: return ( len(self.pistons) + len(self.masses) + len(self.elastic_endstops) + len(self.mechanical_nodes) + len(self.piecewise_signals) + len(self.force_connectors) + len(self.zero_force_sources) ) @property def aliases(self) -> tuple[str, ...]: return tuple( [ *self.pistons, *self.masses, *self.elastic_endstops, *self.mechanical_nodes, *self.piecewise_signals, *self.force_connectors, *self.zero_force_sources, ] ) @dataclass(frozen=True) class TestMqlMechanicalMassState: alias: str velocity_m_s: float displacement_m: float def as_vector(self) -> list[float]: return [self.velocity_m_s, self.displacement_m] @dataclass(frozen=True) class TestMqlMechanicalNodeKinematics: alias: str velocities_m_s: dict[int, float] displacements_m: dict[int, float] @dataclass(frozen=True) class TestMqlPistonKinematics: alias: str port_2_velocity_m_s: float port_2_displacement_m: float port_3_velocity_m_s: float port_3_displacement_m: float @dataclass(frozen=True) class TestMqlMechanicalMassSnapshot: states: tuple[TestMqlMechanicalMassState, ...] node_kinematics_by_alias: dict[str, TestMqlMechanicalNodeKinematics] piston_kinematics_by_alias: dict[str, TestMqlPistonKinematics] @property def state_count(self) -> int: return 2 * len(self.states) class TestMqlMechanicalMassClosure: def __init__(self, assembly: TestMqlMechanicalAssembly) -> None: self.assembly = assembly self.mass_aliases = tuple(assembly.masses) def initial_state_vector(self) -> list[float]: state: list[float] = [] for alias in self.mass_aliases: spec = self.assembly.masses[alias] state.extend([spec.initial_velocity_m_s, spec.initial_displacement_m]) return state def snapshot(self, state_vector: list[float] | None = None) -> TestMqlMechanicalMassSnapshot: values = self.initial_state_vector() if state_vector is None else list(state_vector) if len(values) != 2 * len(self.mass_aliases): raise ValueError("mechanical mass state vector requires two values per mass") states = tuple( TestMqlMechanicalMassState( alias=alias, velocity_m_s=values[2 * index], displacement_m=values[2 * index + 1], ) for index, alias in enumerate(self.mass_aliases) ) node_kinematics = self._node_kinematics_by_alias(states) return TestMqlMechanicalMassSnapshot( states=states, node_kinematics_by_alias=node_kinematics, piston_kinematics_by_alias=self._piston_kinematics_by_alias( states, node_kinematics, ), ) def _node_kinematics_by_alias( self, states: tuple[TestMqlMechanicalMassState, ...], ) -> dict[str, TestMqlMechanicalNodeKinematics]: state_by_alias = {state.alias: state for state in states} front = state_by_alias["mass_friction_endstops_18"] rear = state_by_alias["mass_friction_endstops_19"] return { "dynamic_mechanical_node_alternative_2": TestMqlMechanicalNodeKinematics( alias="dynamic_mechanical_node_alternative_2", velocities_m_s={port: -front.velocity_m_s for port in range(1, 9)}, displacements_m={port: -front.displacement_m for port in range(1, 9)}, ), "dynamic_mechanical_node_alternative_3": TestMqlMechanicalNodeKinematics( alias="dynamic_mechanical_node_alternative_3", velocities_m_s={port: rear.velocity_m_s for port in range(1, 9)}, displacements_m={port: rear.displacement_m for port in range(1, 9)}, ), } def _piston_kinematics_by_alias( self, states: tuple[TestMqlMechanicalMassState, ...], node_kinematics_by_alias: dict[str, TestMqlMechanicalNodeKinematics], ) -> dict[str, TestMqlPistonKinematics]: state_by_alias = {state.alias: state for state in states} rear_node = node_kinematics_by_alias["dynamic_mechanical_node_alternative_3"] piston_bindings = ( ("pn_brp2_8", "mass_friction_endstops_10", 8), ("pn_brp2_9", "mass_friction_endstops_11", 7), ("pn_brp2_10", "mass_friction_endstops_12", 6), ("pn_brp2_11", "mass_friction_endstops_13", 5), ("pn_brp2_12", "mass_friction_endstops_14", 4), ("pn_brp2_13", "mass_friction_endstops_15", 3), ("pn_brp2_14", "mass_friction_endstops_16", 2), ("pn_brp2_15", "mass_friction_endstops_17", 1), ) return { piston_alias: TestMqlPistonKinematics( alias=piston_alias, port_2_velocity_m_s=state_by_alias[mass_alias].velocity_m_s, port_2_displacement_m=state_by_alias[mass_alias].displacement_m, port_3_velocity_m_s=rear_node.velocities_m_s[rear_node_port], port_3_displacement_m=rear_node.displacements_m[rear_node_port], ) for piston_alias, mass_alias, rear_node_port in piston_bindings } def rhs( self, state_vector: list[float], *, force_by_mass_alias: dict[str, float] | None = None, constrained_mass_aliases: set[str] | None = None, ) -> list[float]: snapshot = self.snapshot(state_vector) force_by_mass_alias = force_by_mass_alias or {} constrained_mass_aliases = constrained_mass_aliases or set() derivatives: list[float] = [] for state in snapshot.states: spec = self.assembly.masses[state.alias] mass = spec.endstop() applied_force = force_by_mass_alias.get(state.alias, 0.0) acceleration, velocity = mass.derivatives( velocity_m_s=state.velocity_m_s, displacement_m=state.displacement_m, port_1_force_n=applied_force, ) if state.alias in constrained_mass_aliases and _limit_constraint_holds( spec, state, applied_force, ): acceleration = 0.0 velocity = 0.0 derivatives.extend([acceleration, velocity]) return derivatives def build_test_mql_mechanical_assembly( config: TestMqlConfig | None = None, amesim_results: AmesimResults | None = None, variable_catalog: TestMqlVariableCatalog | None = None, ) -> TestMqlMechanicalAssembly: config = config or TestMqlConfig.from_amesim_specs() if variable_catalog is None and amesim_results is not None: variable_catalog = build_test_mql_variable_catalog(amesim_results) pistons = { component.alias: _build_piston(component, variable_catalog) for component in config.components_by_submodel("PNRP17") } masses = { component.alias: _build_mass(component, variable_catalog, amesim_results) for component in config.components_by_submodel("MECMAS21") } elastic_endstops = { component.alias: _build_elastic_endstop(component, variable_catalog) for component in config.components_by_submodel("LSTP00A") } mechanical_nodes = { component.alias: _build_mechanical_node(component, variable_catalog) for component in config.components_by_submodel("LMECHN1") } piecewise_signals = { component.alias: _build_piecewise_signal(component, variable_catalog) for component in config.components_by_submodel("UD00") } force_connectors = { component.alias: _build_force_connector(component, variable_catalog) for component in config.components_by_submodel("FORC") } zero_force_sources = tuple(component.alias for component in config.components_by_submodel("F000")) return TestMqlMechanicalAssembly( pistons=pistons, masses=masses, elastic_endstops=elastic_endstops, mechanical_nodes=mechanical_nodes, piecewise_signals=piecewise_signals, force_connectors=force_connectors, zero_force_sources=zero_force_sources, ) def _build_piston( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, ) -> TestMqlPistonSpec: geometry = AmesimPistonGeometry( piston_diameter_m=mm_to_m(component.parameter_value("dp")), rod_diameter_m=mm_to_m(component.parameter_value("dr")), zero_length_m=mm_to_m(component.parameter_value("x0")), ) return TestMqlPistonSpec( alias=component.alias, piston_diameter_m=geometry.piston_diameter_m, rod_diameter_m=geometry.rod_diameter_m, zero_displacement_m=geometry.zero_length_m, piston_area_m2=geometry.piston_area_m2, rod_area_m2=geometry.rod_area_m2, annulus_area_m2=geometry.annulus_area_m2, data_paths=_data_paths(variable_catalog, component.alias), ) def _build_mass( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, amesim_results: AmesimResults | None, ) -> TestMqlMassEndstopSpec: return TestMqlMassEndstopSpec( alias=component.alias, mass_kg=component.parameter_value("mass"), xmin_m=component.parameter_value("xmin"), xmax_m=component.parameter_value("xmax"), min_stiffness_n_per_m=n_per_mm_to_n_per_m(component.parameter_value("Kbmin")), max_stiffness_n_per_m=n_per_mm_to_n_per_m(component.parameter_value("Kbmax")), min_damping_n_per_m_per_s=n_per_mm_per_s_to_n_per_m_per_s(component.parameter_value("Dbmin")), max_damping_n_per_m_per_s=n_per_mm_per_s_to_n_per_m_per_s(component.parameter_value("Dbmax")), min_penetration_m=mm_to_m(component.parameter_value("Pdmin")), max_penetration_m=mm_to_m(component.parameter_value("Pdmax")), stiction_force_n=component.parameter_value("fstick"), coulomb_friction_n=component.parameter_value("fcoul"), viscous_friction_n_per_m_per_s=component.parameter_value("rvisc"), windage_n_per_m2_per_s2=component.parameter_value("wind"), stick_velocity_threshold_m_s=component.parameter_value("dvel"), reset_velocity_threshold_m_s=component.parameter_value("restdvel"), rest_coeff=component.parameter_value("restcoeff"), stribeck_constant_m_s=component.parameter_value("astrib"), use_friction=bool(int(component.parameter_value("useFriction"))), stop_type=int(component.parameter_value("stoptype")), initial_velocity_m_s=_initial_value(amesim_results, f"v1@{component.alias}"), initial_displacement_m=_initial_value(amesim_results, f"x1@{component.alias}"), data_paths=_data_paths(variable_catalog, component.alias), ) def _build_elastic_endstop( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, ) -> TestMqlElasticEndstopSpec: return TestMqlElasticEndstopSpec( alias=component.alias, gap_m=mm_to_m(component.parameter_value("gap0")), contact_stiffness_n_per_m=component.parameter_value("kcont"), contact_damping_n_per_m_per_s=component.parameter_value("rcont"), spring_diameter_m=mm_to_m(component.parameter_value("sdiam")), wire_diameter_m=mm_to_m(component.parameter_value("wdiam")), data_paths=_data_paths(variable_catalog, component.alias), ) def _build_mechanical_node( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, ) -> TestMqlMechanicalNodeSpec: return TestMqlMechanicalNodeSpec( alias=component.alias, port_count=int(component.parameter_value("v1")), sum_mode=int(component.parameter_value("sum")), data_paths=_data_paths(variable_catalog, component.alias), ) def _limit_constraint_holds( spec: TestMqlMassEndstopSpec, state: TestMqlMechanicalMassState, applied_force_n: float, ) -> bool: if abs(state.velocity_m_s) > spec.stick_velocity_threshold_m_s: return False at_lower_limit = state.displacement_m <= spec.xmin_m + spec.min_penetration_m at_upper_limit = state.displacement_m >= spec.xmax_m - spec.max_penetration_m return (at_lower_limit and applied_force_n <= 0.0) or ( at_upper_limit and applied_force_n >= 0.0 ) def _build_piecewise_signal( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, ) -> TestMqlPiecewiseLinearSignalSpec: starts = tuple(component.parameter_value(f"start{index}") for index in range(1, 9)) ends = tuple(component.parameter_value(f"end{index}") for index in range(1, 9)) durations = tuple(component.parameter_value(f"t{index}") for index in range(1, 9)) return TestMqlPiecewiseLinearSignalSpec( alias=component.alias, t_start_s=component.parameter_value("tstart"), starts=starts, ends=ends, durations_s=durations, stage_count=int(component.parameter_value("nstages")), is_cyclic=bool(int(component.parameter_value("iscyclic"))), data_paths=_data_paths(variable_catalog, component.alias), ) def _build_force_connector( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, ) -> TestMqlForceConnectorSpec: signal_alias_by_force_connector = { "forcecon_1": "piecewiselinear", "forcecon_2": "piecewiselinear_1", } target_mass_by_force_connector = { "forcecon_1": "mass_friction_endstops_19", "forcecon_2": "mass_friction_endstops_18", } return TestMqlForceConnectorSpec( alias=component.alias, signal_alias=signal_alias_by_force_connector[component.alias], target_mass_alias=target_mass_by_force_connector[component.alias], data_paths=_data_paths(variable_catalog, component.alias), ) def n_per_mm_to_n_per_m(value: float) -> float: return value * N_PER_MM_TO_N_PER_M def n_per_mm_per_s_to_n_per_m_per_s(value: float) -> float: return value * N_PER_MM_PER_S_TO_N_PER_M_PER_S def _initial_value(amesim_results: AmesimResults | None, data_path: str) -> float: if amesim_results is None: return 0.0 return float(amesim_results.series(data_path)[0]) def _data_paths( variable_catalog: TestMqlVariableCatalog | None, alias: str, ) -> tuple[str, ...]: if variable_catalog is None: return () return variable_catalog.data_paths_for_owner(alias)