from __future__ import annotations from collections.abc import Callable from dataclasses import dataclass from math import expm1, isfinite, log, sqrt from app.simulation.components.amesim.boundary.sources import AmesimPnpl01 from app.simulation.components.amesim.flow.orifices import AmesimPnor001 from app.simulation.components.amesim.flow.pipes import ( AmesimPnl00r, AmesimPnl0001, AmesimPnl0002, ) from app.simulation.core.equations import EquationResidual from app.simulation.core.ports import PortState, VariableRole from app.simulation.systems.network import SimulationNetwork class AlgebraicSolveError(RuntimeError): def __init__(self, message: str, diagnostics: "AlgebraicSolveDiagnostics") -> None: super().__init__(message) self.diagnostics = diagnostics @dataclass(frozen=True) class AlgebraicUnknown: component: str port: str variable: str role: VariableRole state: PortState @property def id(self) -> str: return f"{self.component}.{self.port}.{self.variable}" def read(self) -> float: return float(getattr(self.state, self.variable)) def write(self, value: float) -> None: setattr(self.state, self.variable, float(value)) @dataclass(frozen=True) class ExplicitFlowAssignment: equation_id: str unknown: AlgebraicUnknown evaluate: Callable[[], float] | None component: object | None = None equation_index: int | None = None @dataclass(frozen=True) class ExplicitFlowStage: assignments: tuple[ExplicitFlowAssignment, ...] direct_evaluations: tuple[tuple[int, Callable[[], float]], ...] component_evaluations: tuple["ExplicitFlowComponentEvaluation", ...] @dataclass(frozen=True) class ExplicitFlowComponentEvaluation: component: object evaluate: Callable[[], tuple[float, ...]] assignment_indices: tuple[int, ...] equation_indices: tuple[int, ...] equation_ids: tuple[str, ...] @dataclass(frozen=True) class EquationScalePlan: variable_names: tuple[str, ...] force_unknowns: tuple[AlgebraicUnknown, ...] @dataclass(frozen=True) class EffortAnchor: unknown: AlgebraicUnknown evaluate: Callable[[], float] @dataclass(frozen=True) class ConnectionEquationEvaluation: template: EquationResidual evaluate: Callable[[], float] @dataclass(frozen=True) class ComponentEquationEvaluation: component: object evaluate: Callable[[], tuple[float, ...]] templates: tuple[EquationResidual, ...] @dataclass(frozen=True) class PnorPnl0001SeriesBinding: orifice: AmesimPnor001 orifice_port: str pipe: AmesimPnl0001 pipe_port: str @dataclass(frozen=True) class ClosedResistancePressureBinding: component: object port_name: str neighbor: object neighbor_port: str pressure_source_port: str | None @dataclass(frozen=True) class EffortEqualityGroup: variable: str members: tuple[AlgebraicUnknown, ...] anchors: tuple[EffortAnchor, ...] @dataclass(frozen=True) class UnilateralContactBinding: component: object algebraic_group: EffortEqualityGroup neighbor_force: AlgebraicUnknown algebraic_port: int force_sign: float @dataclass(frozen=True) class AlgebraicSolveDiagnostics: success: bool message: str evaluations: int pressure_scale: float flow_scale: float max_scaled_residual: float max_raw_residual: float def as_dict(self) -> dict[str, object]: return { "success": self.success, "message": self.message, "evaluations": self.evaluations, "pressureScale": self.pressure_scale, "flowScale": self.flow_scale, "maxScaledResidual": self.max_scaled_residual, "maxRawResidual": self.max_raw_residual, } class PressureFlowSolver: """Solve the acausal pressure-flow subsystem for a compiled network.""" def __init__( self, network: SimulationNetwork, *, residual_tolerance: float = 1e-7, max_evaluations: int = 500, ) -> None: self.network = network self.residual_tolerance = residual_tolerance self.max_evaluations = max_evaluations self.unknowns = self._build_unknowns() self._unknowns_by_id = {unknown.id: unknown for unknown in self.unknowns} self._unknowns_by_variable = { variable: tuple( unknown for unknown in self.unknowns if unknown.variable == variable ) for variable in ("p", "m_flow", "x", "v", "f") } self._component_equation_owners = tuple(network.components.values()) self._component_equation_plan = tuple( ComponentEquationEvaluation( component=component, evaluate=component.pressure_flow_equation_values, templates=component.pressure_flow_equation_residuals(), ) for component in self._component_equation_owners ) self._component_equation_plans_by_id = { item.component.name: item for item in self._component_equation_plan } self._estimated_flow_components = tuple( component for component in self._component_equation_owners if hasattr(component, "K_eff") ) self._causal_contact_components = tuple( component for component in self._component_equation_owners if getattr(component, "clear_causal_contact", None) is not None ) self._connection_equation_plan = tuple( ConnectionEquationEvaluation( template=equation, evaluate=self._equation_value_reader(equation), ) for equation in network.connection_equation_residuals() ) self._equation_templates = tuple( equation for item in self._component_equation_plan for equation in item.templates ) + tuple(item.template for item in self._connection_equation_plan) self._effort_groups = { variable: self._build_effort_equality_groups(variable) for variable in ("p", "x", "v") } self._pnor_pnl0001_series_plan = ( self._build_pnor_pnl0001_series_plan() ) self._closed_resistance_pressure_plan = ( self._build_closed_resistance_pressure_plan() ) self._unilateral_contact_plan = self._build_unilateral_contact_plan() self._explicit_flow_plan = self._build_explicit_flow_plan() self._explicit_flow_plans_by_variables = { selected: self._filter_explicit_flow_plan(selected) for selected in (frozenset(("f", "m_flow")), frozenset(("f",))) } self._equation_scale_plans = self._build_equation_scale_plans( self._equation_templates ) self.last_diagnostics: AlgebraicSolveDiagnostics | None = None def _build_unknowns(self) -> tuple[AlgebraicUnknown, ...]: unknowns: list[AlgebraicUnknown] = [] for component in self.network.components.values(): for definition in component.active_port_definitions: if definition.kind != "physical": continue state = component.get_port(definition.name) for variable in definition.variables: if variable.role not in {"effort", "flow"}: continue unknowns.append( AlgebraicUnknown( component=component.name, port=definition.name, variable=variable.name, role=variable.role, state=state, ) ) return tuple(unknowns) @staticmethod def _port_key(variable: str, expected_variable: str) -> tuple[str, str] | None: try: component_name, port_name, variable_name = variable.rsplit(".", 2) except ValueError: return None if variable_name != expected_variable: return None return component_name, port_name def _seed_equal_efforts( self, variables: tuple[str, ...] = ("p", "x", "v"), ) -> None: """Lift state-owned efforts across their complete equality groups. Dynamic components refresh their own ports before each closure, while connected algebraic ports retain values from the preceding RHS evaluation. State equations expose the current effort as ``port.variable - target``; use that target as the authoritative anchor for every connected/equal pressure, displacement, and velocity port before evaluating explicit flow laws. """ self.propagate_equal_efforts(variables) def propagate_equal_efforts(self, variables: tuple[str, ...]) -> None: """Propagate selected state-owned efforts without solving flows. Piston geometry needs current mechanical ``x``/``v`` before swept volume propagation, but pressure and flow equations can wait until the connected chamber has refreshed that volume. """ unknown = sorted(set(variables) - set(self._effort_groups)) if unknown: raise ValueError("Unsupported effort variables: " + ", ".join(unknown)) for variable in variables: self._seed_equal_effort(variable) def _build_effort_equality_groups( self, variable: str, ) -> tuple[EffortEqualityGroup, ...]: effort_unknowns = { (unknown.component, unknown.port): unknown for unknown in self.unknowns if unknown.variable == variable } if not effort_unknowns: return () parent = {key: key for key in effort_unknowns} def find(key: tuple[str, str]) -> tuple[str, str]: root = key while parent[root] != root: root = parent[root] while parent[key] != key: next_key = parent[key] parent[key] = root key = next_key return root def union(first: tuple[str, str], second: tuple[str, str]) -> None: first_root = find(first) second_root = find(second) if first_root != second_root: parent[second_root] = first_root for connection in self.network.connections: if connection.kind != "physical": continue first = connection.endpoint_a.key second = connection.endpoint_b.key if first in effort_unknowns and second in effort_unknowns: union(first, second) component_equations = { item.component.name: item.templates for item in self._component_equation_plan } for equations in component_equations.values(): for equation in equations: if equation.relation != "equal" or equation.role != "effort": continue endpoints = [ endpoint for equation_variable in equation.variables if ( (endpoint := self._port_key(equation_variable, variable)) in effort_unknowns ) ] for endpoint in endpoints[1:]: union(endpoints[0], endpoint) members_by_root: dict[tuple[str, str], list[AlgebraicUnknown]] = {} for endpoint in effort_unknowns: members_by_root.setdefault(find(endpoint), []).append( effort_unknowns[endpoint] ) 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": continue endpoints = [ endpoint for equation_variable in equation.variables if ( (endpoint := self._port_key(equation_variable, variable)) in effort_unknowns ) ] if len(endpoints) != 1: continue endpoint = endpoints[0] unknown = effort_unknowns[endpoint] anchors_by_root.setdefault(find(endpoint), []).append( EffortAnchor( unknown=unknown, evaluate=self._equation_value_reader(equation), ) ) return tuple( EffortEqualityGroup( variable=variable, members=tuple(members), anchors=tuple(anchors_by_root.get(root, ())), ) for root, members in members_by_root.items() ) def _seed_equal_effort(self, variable: str) -> None: for group in self._effort_groups[variable]: members = group.members 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. for unknown, target_value in anchors: unknown.write(target_value) anchor_values = [value for _unknown, value in anchors] effort_scale = max([abs(value) for value in anchor_values] + [1.0]) if max(anchor_values) - min(anchor_values) > 1.0e-9 * effort_scale: # A conflicting multi-storage group is structurally invalid; # leave it for the residual solver/preparation diagnostics. continue target_value = sum(anchor_values) / len(anchor_values) for unknown in members: unknown.write(target_value) continue if variable == "p": seed = next( ( unknown.read() for unknown in members if unknown.read() > 0.0 ), None, ) if seed is None: continue else: seed = members[0].read() for unknown in members: if variable != "p" or unknown.read() <= 0.0: unknown.write(seed) def _connected_flow_unknown( self, component_name: str, port_name: str, variable: str, ) -> AlgebraicUnknown | None: endpoint_key = (component_name, port_name) for connection in self.network.connections: if connection.kind != "physical": continue if connection.endpoint_a.key == endpoint_key: other = connection.endpoint_b elif connection.endpoint_b.key == endpoint_key: other = connection.endpoint_a else: continue return self._unknowns_by_id.get( f"{other.component}.{other.port}.{variable}" ) return None @staticmethod def _bisect_contact_root( value_at, lower: float, upper: float, target: float, ) -> float | None: lower_value = float(value_at(lower)) - target upper_value = float(value_at(upper)) - target tolerance = 1.0e-13 * max(abs(target), 1.0) if abs(lower_value) <= tolerance: return lower if abs(upper_value) <= tolerance: return upper if not isfinite(lower_value) or not isfinite(upper_value): return None if (lower_value < 0.0) == (upper_value < 0.0): return None for _iteration in range(100): middle = 0.5 * (lower + upper) middle_value = float(value_at(middle)) - target if abs(middle_value) <= tolerance: return middle if (lower_value < 0.0) == (middle_value < 0.0): lower = middle lower_value = middle_value else: upper = middle upper_value = middle_value return 0.5 * (lower + upper) def _contact_penetration_for_force( self, component, requested_force: float, current_penetration: float, ) -> float | None: """Invert one LSTP force law and select the root nearest its current state.""" if not isfinite(requested_force): return None option = int(getattr(component, "discContactOption", 2.0)) if option != 1: requested_force = max(requested_force, 0.0) stiffness = max(float(getattr(component, "kcont", 0.0)), 0.0) damping = max(float(getattr(component, "rcont", 0.0)), 0.0) damping_length = float(getattr(component, "Pdis", 0.0)) relative_velocity = float(getattr(component, "penetration_velocity")) damping_term = damping * relative_velocity current_penetration = ( max(float(current_penetration), 0.0) if isfinite(current_penetration) else 0.0 ) force_tolerance = 1.0e-12 * max(abs(requested_force), 1.0) def raw_force(penetration: float) -> float: if penetration <= 0.0: return 0.0 damping_fraction = ( -expm1(-penetration / damping_length) if damping_length > 0.0 else 1.0 ) return stiffness * penetration + damping_term * damping_fraction def contact_force(penetration: float) -> float: force = raw_force(penetration) return force if option == 1 else max(force, 0.0) candidates: list[float] = [] def add_candidate(penetration: float | None) -> None: if penetration is None or not isfinite(penetration) or penetration < 0.0: return if abs(contact_force(penetration) - requested_force) > force_tolerance: return if not any( abs(penetration - candidate) <= 1.0e-12 * max(abs(penetration), abs(candidate), 1.0e-18) for candidate in candidates ): candidates.append(penetration) add_candidate(current_penetration) add_candidate(0.0) if option != 1 and requested_force == 0.0: return min( candidates or [0.0], key=lambda penetration: abs(penetration - current_penetration), ) if damping_length <= 0.0: if stiffness > 0.0: penetration = (requested_force - damping_term) / stiffness if penetration > 0.0: add_candidate(penetration) elif abs(requested_force - damping_term) <= force_tolerance: add_candidate(max(current_penetration, 1.0e-18)) elif stiffness > 0.0: critical_penetration: float | None = None if damping_term < -stiffness * damping_length: critical_penetration = damping_length * log( -damping_term / (stiffness * damping_length) ) add_candidate(critical_penetration) upper = max( current_penetration, damping_length, abs(requested_force) / stiffness, critical_penetration or 0.0, 1.0e-18, ) for _iteration in range(100): upper_value = raw_force(upper) if isfinite(upper_value) and upper_value >= requested_force: break upper *= 2.0 else: upper = float("nan") if isfinite(upper): if critical_penetration is not None: add_candidate( self._bisect_contact_root( raw_force, 0.0, critical_penetration, requested_force, ) ) add_candidate( self._bisect_contact_root( raw_force, critical_penetration, upper, requested_force, ) ) else: add_candidate( self._bisect_contact_root( raw_force, 0.0, upper, requested_force, ) ) elif damping_term != 0.0: upper = max(current_penetration, damping_length, 1.0e-18) for _iteration in range(100): upper_value = raw_force(upper) crossed = ( upper_value >= requested_force if damping_term > 0.0 else upper_value <= requested_force ) if isfinite(upper_value) and crossed: add_candidate( self._bisect_contact_root( raw_force, 0.0, upper, requested_force, ) ) break upper *= 2.0 if not candidates: return None return min( candidates, key=lambda penetration: abs(penetration - current_penetration), ) def _apply_unilateral_contact_binding( self, binding: UnilateralContactBinding, ) -> bool: component = binding.component requested_force = binding.force_sign * binding.neighbor_force.read() if int(getattr(component, "discContactOption", 2.0)) != 1: requested_force = max(requested_force, 0.0) cached_penetration = getattr(component, "_causal_penetration", None) penetration = self._contact_penetration_for_force( component, requested_force, ( float(cached_penetration) if cached_penetration is not None else float(getattr(component, "penetration")) ), ) if penetration is None: component.clear_causal_contact() return False gap0 = float(getattr(component, "gap0", 0.0)) if binding.algebraic_port == 1: target = component.port_2.x + gap0 + penetration else: target = component.port_1.x - gap0 - penetration for unknown in binding.algebraic_group.members: unknown.write(target) component.set_causal_contact( penetration=penetration, force=requested_force, ) return True def _refresh_unilateral_contacts( self, bindings: tuple[UnilateralContactBinding, ...], ) -> None: for binding in bindings: self._apply_unilateral_contact_binding(binding) def _build_unilateral_contact_plan( self, ) -> tuple[UnilateralContactBinding, ...]: """Compile contacts that can eliminate one algebraic coordinate.""" position_groups = { unknown.id: group for group in self._effort_groups["x"] for unknown in group.members } bindings: list[UnilateralContactBinding] = [] for component in self.network.components.values(): if component.model_type != "amesim_lstp00a": continue first_neighbor = self._connected_flow_unknown( component.name, "port_1", "f", ) second_neighbor = self._connected_flow_unknown( component.name, "port_2", "f", ) first_group = position_groups.get(f"{component.name}.port_1.x") second_group = position_groups.get(f"{component.name}.port_2.x") if ( first_group is None or second_group is None or first_group is second_group ): continue if not first_group.anchors and first_neighbor is not None: binding = UnilateralContactBinding( component=component, algebraic_group=first_group, neighbor_force=first_neighbor, algebraic_port=1, force_sign=-1.0, ) elif not second_group.anchors and second_neighbor is not None: binding = UnilateralContactBinding( component=component, algebraic_group=second_group, neighbor_force=second_neighbor, algebraic_port=2, force_sign=1.0, ) else: # With both coordinates state-owned, penetration is a dynamic # result rather than an algebraic active-set choice. continue bindings.append(binding) return tuple(bindings) def _seed_unilateral_contacts( self, ) -> tuple[UnilateralContactBinding, ...]: """Apply compiled local contact eliminations for the current state.""" bindings: list[UnilateralContactBinding] = [] bound_group_ids: set[int] = set() for binding in self._unilateral_contact_plan: group_id = id(binding.algebraic_group) if group_id in bound_group_ids: # One relative contact law may eliminate a free coordinate. # Any other contact sharing that coordinate must remain in the # nonlinear system or the projections would overwrite each # other and make root selection order-dependent. continue if self._apply_unilateral_contact_binding(binding): bindings.append(binding) bound_group_ids.add(group_id) return tuple(bindings) 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" ) 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}." ) evaluation = self._component_equation_plans_by_id[equation.owner_id] component = evaluation.component equation_id = equation.id equation_ids = tuple(current.id for current in evaluation.templates) try: equation_index = equation_ids.index(equation_id) except ValueError as exc: raise RuntimeError( f"Compiled algebraic equation disappeared at runtime: {equation_id}." ) from exc def read_component_equation() -> float: current_values = evaluation.evaluate() if equation_index >= len(current_values): raise RuntimeError( "Compiled algebraic equation disappeared at runtime: " f"{equation_id}." ) return float(current_values[equation_index]) return read_component_equation def _pressure_flow_equation_values(self) -> tuple[float, ...]: """Evaluate live equation values through the compiled topology.""" values: list[float] = [] for item in self._component_equation_plan: current_values = item.evaluate() if len(current_values) != len(item.templates): raise RuntimeError( "Compiled algebraic equation count changed at runtime for " f"{item.component.name}." ) values.extend(current_values) values.extend(item.evaluate() for item in self._connection_equation_plan) return tuple(values) def _pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: """Materialize public residual objects only when explicitly requested.""" return tuple( EquationResidual( id=template.id, owner=template.owner, owner_id=template.owner_id, relation=template.relation, variables=template.variables, value=value, role=template.role, ) for template, value in zip( self._equation_templates, self._pressure_flow_equation_values(), ) ) def _explicit_flow_assignment( self, equation, unknown: AlgebraicUnknown, ) -> ExplicitFlowAssignment: if equation.owner == "connection": return ExplicitFlowAssignment( equation_id=equation.id, unknown=unknown, evaluate=self._equation_value_reader(equation), ) evaluation = self._component_equation_plans_by_id[equation.owner_id] component = evaluation.component equations = evaluation.templates equation_ids = tuple(current.id for current in equations) try: equation_index = equation_ids.index(equation.id) except ValueError as exc: raise RuntimeError( f"Compiled algebraic equation disappeared at runtime: {equation.id}." ) from exc return ExplicitFlowAssignment( equation_id=equation.id, unknown=unknown, evaluate=None, component=component, equation_index=equation_index, ) def _build_explicit_flow_plan(self) -> tuple[ExplicitFlowStage, ...]: """Compile flow causalization into independent dependency stages.""" stages: list[ExplicitFlowStage] = [] seeded_ids: set[str] = set() initial_assignments: list[ExplicitFlowAssignment] = [] for evaluation in self._component_equation_plan: for equation in evaluation.templates: if equation.relation != "constitutive" or equation.role != "flow": continue 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 initial_assignments.append( self._explicit_flow_assignment(equation, unknown) ) seeded_ids.add(unknown.id) if initial_assignments: stages.append(self._compile_explicit_flow_stage(tuple(initial_assignments))) equations = self._equation_templates while True: stage_assignments: list[ExplicitFlowAssignment] = [] stage_unknown_ids: set[str] = set() for equation in equations: if equation.role != "flow" or equation.relation not in { "constitutive", "sumToZero", }: continue flow_unknowns = self._flow_unknowns_for_equation(equation) if not flow_unknowns: continue if len({unknown.variable for unknown in flow_unknowns}) != 1: continue unseeded = tuple( unknown for unknown in flow_unknowns if unknown.id not in seeded_ids ) if len(unseeded) != 1: continue unknown = unseeded[0] if unknown.id in stage_unknown_ids: continue stage_assignments.append( self._explicit_flow_assignment(equation, unknown) ) stage_unknown_ids.add(unknown.id) if stage_assignments: stages.append( self._compile_explicit_flow_stage(tuple(stage_assignments)) ) seeded_ids.update(stage_unknown_ids) continue fallback_assignment: ExplicitFlowAssignment | None = None for equation in equations: if equation.role != "flow" or equation.relation not in { "constitutive", "sumToZero", }: continue 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 unknown = unseeded[-1] fallback_assignment = self._explicit_flow_assignment( equation, unknown, ) seeded_ids.add(unknown.id) break if fallback_assignment is None: break stages.append(self._compile_explicit_flow_stage((fallback_assignment,))) return tuple(stages) @staticmethod def _compile_explicit_flow_stage( assignments: tuple[ExplicitFlowAssignment, ...], ) -> ExplicitFlowStage: direct_evaluations: list[tuple[int, Callable[[], float]]] = [] assignments_by_component: dict[ object, list[tuple[int, ExplicitFlowAssignment]], ] = {} for assignment_index, assignment in enumerate(assignments): if assignment.component is None: if assignment.evaluate is None: raise RuntimeError( "Explicit flow assignment has no compiled evaluator: " f"{assignment.equation_id}." ) direct_evaluations.append((assignment_index, assignment.evaluate)) continue assignments_by_component.setdefault(assignment.component, []).append( (assignment_index, assignment) ) component_evaluations: list[ExplicitFlowComponentEvaluation] = [] for component, component_assignments in assignments_by_component.items(): if any( assignment.equation_index is None for _assignment_index, assignment in component_assignments ): raise RuntimeError( "Explicit component flow assignment has no equation index." ) component_evaluations.append( ExplicitFlowComponentEvaluation( component=component, evaluate=component.pressure_flow_equation_values, assignment_indices=tuple( assignment_index for assignment_index, _assignment in component_assignments ), equation_indices=tuple( int(assignment.equation_index) for _assignment_index, assignment in component_assignments ), equation_ids=tuple( assignment.equation_id for _assignment_index, assignment in component_assignments ), ) ) return ExplicitFlowStage( assignments=assignments, direct_evaluations=tuple(direct_evaluations), component_evaluations=tuple(component_evaluations), ) def _filter_explicit_flow_plan( self, selected: frozenset[str], ) -> tuple[ExplicitFlowStage, ...]: return tuple( self._compile_explicit_flow_stage( tuple( assignment for assignment in stage.assignments if assignment.unknown.variable in selected ) ) for stage in self._explicit_flow_plan ) @staticmethod def _evaluate_explicit_flow_stage( stage: ExplicitFlowStage, ) -> tuple[float, ...]: assignments = stage.assignments values: list[float | None] = [None] * len(assignments) for assignment_index, evaluate in stage.direct_evaluations: values[assignment_index] = evaluate() for evaluation in stage.component_evaluations: equation_values = evaluation.evaluate() for assignment_index, equation_index, equation_id in zip( evaluation.assignment_indices, evaluation.equation_indices, evaluation.equation_ids, ): if equation_index >= len(equation_values): raise RuntimeError( "Compiled algebraic equation disappeared at runtime: " f"{equation_id}." ) values[assignment_index] = float(equation_values[equation_index]) if any(value is None for value in values): raise RuntimeError("Explicit flow evaluation plan returned no value.") return tuple(float(value) for value in values) def _solve_explicit_flow_unknowns( self, variables: tuple[str, ...] = ("f", "m_flow"), ) -> set[str]: """Execute staged flow/force assignments without repeated equations.""" selected = frozenset(variables) for unknown in self.unknowns: if unknown.variable in selected: unknown.write(0.0) seeded_ids: set[str] = set() plan = self._explicit_flow_plans_by_variables.get(selected) if plan is None: plan = self._filter_explicit_flow_plan(selected) for stage in plan: assignments = stage.assignments values = self._evaluate_explicit_flow_stage(stage) targets = tuple( ( assignment, assignment.unknown.read() - value, ) for assignment, value in zip(assignments, values) ) for assignment, target_value in targets: if not isfinite(target_value): continue assignment.unknown.write(target_value) seeded_ids.add(assignment.unknown.id) return seeded_ids def _build_closed_resistance_pressure_plan( self, ) -> tuple[ClosedResistancePressureBinding, ...]: """Compile sealed resistance ends whose zero-flow pressure is known.""" bindings: list[ClosedResistancePressureBinding] = [] connected: dict[tuple[str, str], tuple[str, str]] = {} for connection in self.network.connections: if connection.kind != "physical" or connection.domain != "pneumatic": continue first, second = connection.endpoints connected[first.key] = second.key connected[second.key] = first.key for component in self.network.components.values(): if not isinstance(component, (AmesimPnl00r, AmesimPnl0001)): continue for port_name in component.ports: neighbor_key = connected.get((component.name, port_name)) if neighbor_key is None: continue neighbor = self.network.components[neighbor_key[0]] if not isinstance(neighbor, AmesimPnpl01): continue if isinstance(component, AmesimPnl0002): pressure_source_port = None elif isinstance(component, AmesimPnl0001): if port_name != "port_1": continue pressure_source_port = None else: pressure_source_port = ( "port_2" if port_name == "port_1" else "port_1" ) bindings.append( ClosedResistancePressureBinding( component=component, port_name=port_name, neighbor=neighbor, neighbor_port=neighbor_key[1], pressure_source_port=pressure_source_port, ) ) return tuple(bindings) def _seed_closed_resistance_pressures(self) -> None: """Seed a sealed resistance end at its zero-flow pressure. A PNPL01 fixes flow, not pressure. Starting a dead-ended Darcy branch with the plug-side pressure at the medium reference can otherwise put the nonlinear solver on the singular square-root part of the inverse flow law. At zero flow, these AMESim pipe resistances have exactly zero pressure drop, which gives a deterministic and physically exact seed. """ for binding in self._closed_resistance_pressure_plan: component = binding.component pressure = ( component.properties().p if binding.pressure_source_port is None else component.get_port(binding.pressure_source_port).p ) component.get_port(binding.port_name).p = pressure binding.neighbor.get_port(binding.neighbor_port).p = pressure def _build_pnor_pnl0001_series_plan( self, ) -> tuple[PnorPnl0001SeriesBinding, ...]: bindings: list[PnorPnl0001SeriesBinding] = [] for connection in self.network.connections: first_endpoint, second_endpoint = connection.endpoints first = self.network.components[first_endpoint.component] second = self.network.components[second_endpoint.component] if isinstance(first, AmesimPnor001) and isinstance(second, AmesimPnl0001): orifice, orifice_port = first, first_endpoint.port pipe, pipe_port = second, second_endpoint.port elif isinstance(second, AmesimPnor001) and isinstance(first, AmesimPnl0001): orifice, orifice_port = second, second_endpoint.port pipe, pipe_port = first, first_endpoint.port else: continue if isinstance(pipe, AmesimPnl0002) or pipe_port != "port_1": continue bindings.append( PnorPnl0001SeriesBinding( orifice=orifice, orifice_port=orifice_port, pipe=pipe, pipe_port=pipe_port, ) ) return tuple(bindings) def _seed_pnor_pnl0001_series_pressures(self) -> None: """Causalize the pressure between a PNOR001 and PNL0001 R port.""" from scipy.optimize import brentq for binding in self._pnor_pnl0001_series_plan: orifice = binding.orifice orifice_port = binding.orifice_port pipe = binding.pipe pipe_port = binding.pipe_port orifice_other = "port_2" if orifice_port == "port_1" else "port_1" pressure_a = orifice.get_port(orifice_other).p pressure_b = pipe.properties().p lower = min(pressure_a, pressure_b) upper = max(pressure_a, pressure_b) def mismatch(intermediate_pressure: float) -> float: if orifice_port == "port_2": orifice_flow_into_connection = -orifice.mass_flow( pressure_a, intermediate_pressure, ) else: orifice_flow_into_connection = orifice.mass_flow( intermediate_pressure, pressure_a, ) pipe_flow_into_connection = pipe.mass_flow( intermediate_pressure, pressure_b, pipe.properties().T, ) return orifice_flow_into_connection + pipe_flow_into_connection lower_value = mismatch(lower) upper_value = mismatch(upper) if lower_value == 0.0: pressure = lower elif upper_value == 0.0: pressure = upper elif (lower_value < 0.0) == (upper_value < 0.0): continue else: pressure = float( brentq( mismatch, lower, upper, xtol=1.0e-6, rtol=1.0e-12, maxiter=32, ) ) orifice.get_port(orifice_port).p = pressure pipe.get_port(pipe_port).p = pressure def _equation_scale_plan(self, equation) -> EquationScalePlan: return EquationScalePlan( variable_names=tuple( variable.rsplit(".", 1)[-1] for variable in equation.variables ), force_unknowns=tuple( unknown for variable in equation.variables if (unknown := self._unknowns_by_id.get(variable)) is not None and unknown.variable == "f" ), ) def _build_equation_scale_plans( self, equations: tuple[EquationResidual, ...], ) -> tuple[EquationScalePlan, ...]: return tuple( self._equation_scale_plan(equation) for equation in equations ) def _scales(self) -> dict[str, float]: pressure_values = [ unknown.read() for unknown in self._unknowns_by_variable["p"] ] pressure_scale = max( [ abs(value) for value in pressure_values if value > 0.0 ] + [1e5] ) estimated_flows = [ abs(float(getattr(component, "K_eff"))) * sqrt(pressure_scale) for component in self._estimated_flow_components ] mass_flow_scale = max( estimated_flows + [ abs(unknown.read()) for unknown in self._unknowns_by_variable["m_flow"] ] + [1e-3] ) return { "p": pressure_scale, "m_flow": mass_flow_scale, "x": max( [abs(unknown.read()) for unknown in self._unknowns_by_variable["x"]] + [1.0] ), "v": max( [abs(unknown.read()) for unknown in self._unknowns_by_variable["v"]] + [1.0] ), "f": max( [abs(unknown.read()) for unknown in self._unknowns_by_variable["f"]] + [1.0] ), } def solve( self, *, effort_variables: tuple[str, ...] = ("p", "x", "v"), ) -> AlgebraicSolveDiagnostics: try: import numpy as np from scipy.optimize import least_squares except ImportError as exc: raise RuntimeError( "Topology-driven simulation requires SciPy; install requirements.txt." ) from exc for component in self._causal_contact_components: component.clear_causal_contact() self._seed_equal_efforts(effort_variables) self._seed_closed_resistance_pressures() self._seed_pnor_pnl0001_series_pressures() self._solve_explicit_flow_unknowns() contact_bindings = self._seed_unilateral_contacts() if contact_bindings: self._solve_explicit_flow_unknowns(("f",)) self._refresh_unilateral_contacts(contact_bindings) scales = self._scales() pressure_scale = scales["p"] flow_scale = scales["m_flow"] seeded_values = self._pressure_flow_equation_values() def initial_equation_scale( equation: EquationResidual, value: float, scale_plan: EquationScalePlan, ) -> float: variable_names = scale_plan.variable_names if equation.role == "flow": force_scales = [ max(abs(unknown.read()), 1.0) for unknown in scale_plan.force_unknowns ] if force_scales: # Freeze force scaling per equation. A 1e17 N source must # not hide an unrelated 40 N piston/contact imbalance in a # different mechanical branch. return max(force_scales + [abs(value), 1.0]) return flow_scale if equation.role == "effort": if "x" in variable_names: return scales["x"] if "v" in variable_names: return scales["v"] return pressure_scale return max([scales.get(name, 1.0) for name in variable_names] + [1.0]) equation_scales = tuple( initial_equation_scale(equation, value, scale_plan) for equation, value, scale_plan in zip( self._equation_templates, seeded_values, self._equation_scale_plans, ) ) seeded_scaled = [ abs(value / scale) for value, scale in zip(seeded_values, equation_scales) ] seeded_max_scaled_residual = max(seeded_scaled, default=0.0) seeded_unknown_values = [ (unknown, unknown.read()) for unknown in self.unknowns ] seeded_unknowns_are_feasible = all( isfinite(value) and (unknown.variable != "p" or value >= 1.0) for unknown, value in seeded_unknown_values ) if ( seeded_unknowns_are_feasible and all(isfinite(value) for value in seeded_scaled) and seeded_max_scaled_residual <= self.residual_tolerance ): diagnostics = AlgebraicSolveDiagnostics( success=True, message="Seeded pressure-flow state satisfies the residual tolerance.", evaluations=0, pressure_scale=pressure_scale, flow_scale=flow_scale, max_scaled_residual=seeded_max_scaled_residual, max_raw_residual=max( (abs(value) for value in seeded_values), default=0.0, ), ) self.last_diagnostics = diagnostics return diagnostics unknown_scales = { unknown.id: ( max(abs(unknown.read()), 1.0) if unknown.variable == "f" else scales[unknown.variable] ) for unknown in self.unknowns } def variable_scale(unknown: AlgebraicUnknown) -> float: return unknown_scales[unknown.id] positive_pressures = [ unknown.read() for unknown in self._unknowns_by_variable["p"] if unknown.read() > 0.0 ] fallback_pressure = ( sum(positive_pressures) / len(positive_pressures) if positive_pressures else pressure_scale ) # A causal contact retains its small relative penetration around the # current absolute port coordinates. Keep that local coordinate during # nonlinear fallback: the contact law remains responsive to optimizer # increments, while a sub-ULP penetration is not lost by subtracting two # large absolute displacements. x0 = np.asarray( [ ( unknown.read() if unknown.variable != "p" or unknown.read() > 0.0 else fallback_pressure ) / variable_scale(unknown) for unknown in self.unknowns ], dtype=float, ) lower = np.asarray( [ 1.0 / pressure_scale if unknown.variable == "p" else -np.inf for unknown in self.unknowns ] ) upper = np.full(len(self.unknowns), np.inf) def assign(values) -> None: for unknown, value in zip(self.unknowns, values): unknown.write(float(value) * variable_scale(unknown)) def scaled_residuals(values): assign(values) self._refresh_unilateral_contacts(contact_bindings) equation_values = self._pressure_flow_equation_values() return np.asarray( [ value / scale for value, scale in zip(equation_values, equation_scales) ], dtype=float, ) result = least_squares( scaled_residuals, x0, bounds=(lower, upper), x_scale="jac", ftol=1e-10, xtol=1e-10, gtol=1e-10, max_nfev=self.max_evaluations, ) assign(result.x) self._refresh_unilateral_contacts(contact_bindings) equation_values = self._pressure_flow_equation_values() scaled = [ abs(value / scale) for value, scale in zip(equation_values, equation_scales) ] max_scaled_residual = max(scaled, default=0.0) residuals_converged = ( all(isfinite(value) for value in scaled) and max_scaled_residual <= self.residual_tolerance ) optimizer_status_is_acceptable = bool(result.success) or int(result.status) == 0 success = residuals_converged and optimizer_status_is_acceptable diagnostics = AlgebraicSolveDiagnostics( success=success, message=str(result.message), evaluations=int(result.nfev), pressure_scale=pressure_scale, flow_scale=flow_scale, max_scaled_residual=max_scaled_residual, max_raw_residual=max( (abs(value) for value in equation_values), default=0.0, ), ) self.last_diagnostics = diagnostics if not success: raise AlgebraicSolveError( "Pressure-flow equations did not converge to the requested tolerance.", diagnostics, ) return diagnostics