from __future__ import annotations from collections.abc import Callable, Mapping from dataclasses import dataclass, replace from math import expm1, isfinite, log, sqrt import os 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.performance import profile_phase from app.simulation.systems.network import SimulationNetwork # The pneumatic constitutive laws require strictly positive absolute pressure. # Keep the optimizer's open lower bound at zero: adaptive explicit integrators # can legitimately probe positive sub-pascal trial states before rejecting a # high-stiffness step. A 1 Pa bound made those otherwise valid seeds fail in # SciPy before the residuals were evaluated. PRESSURE_LOWER_BOUND_PA = 0.0 CAUSAL_FAST_PATH_ENVIRONMENT_VARIABLE = "SIMULATION_CAUSAL_FAST_PATH" CAUSAL_FAST_PATH_AUDIT_INTERVAL = 64 def _causal_fast_path_environment_enabled() -> bool: value = os.getenv(CAUSAL_FAST_PATH_ENVIRONMENT_VARIABLE, "1") return value.strip().lower() not in {"0", "false", "no", "off"} class AlgebraicSolveError(RuntimeError): def __init__( self, message: str, diagnostics: "AlgebraicSolveDiagnostics", *, scope_kind: str = "network", scope_components: tuple[str, ...] = (), ) -> None: super().__init__(message) self.diagnostics = diagnostics self.scope_kind = scope_kind self.scope_components = scope_components @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] equation_id: str @dataclass(frozen=True) class CausalEffortAssignment: """One uniquely state-anchored effort equality group.""" variable: str members: tuple[AlgebraicUnknown, ...] anchor: EffortAnchor @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 AlgebraicEquationBlock: """One connected component of the equation/unknown incidence graph.""" unknown_indices: tuple[int, ...] equation_indices: tuple[int, ...] @dataclass(frozen=True) class ScopedComponentEquationEvaluation: evaluate: Callable[[], tuple[float, ...]] targets: tuple[tuple[int, int, str], ...] @dataclass(frozen=True) class ScopedConnectionEquationEvaluation: target: int evaluate: Callable[[], float] @dataclass(frozen=True) class AlgebraicEquationSubset: """Compiled residual evaluation for a union of independent blocks.""" unknown_indices: tuple[int, ...] equation_indices: tuple[int, ...] component_evaluations: tuple[ScopedComponentEquationEvaluation, ...] connection_evaluations: tuple[ScopedConnectionEquationEvaluation, ...] jacobian_sparsity: object def equation_values(self) -> tuple[float, ...]: values: list[float | None] = [None] * len(self.equation_indices) for evaluation in self.component_evaluations: component_values = evaluation.evaluate() for target, source, equation_id in evaluation.targets: if source >= len(component_values): raise RuntimeError( "Compiled algebraic equation disappeared at runtime: " f"{equation_id}." ) values[target] = float(component_values[source]) for evaluation in self.connection_evaluations: values[evaluation.target] = float(evaluation.evaluate()) if any(value is None for value in values): raise RuntimeError("Scoped algebraic evaluation returned no value.") return tuple(float(value) for value in values) @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 residual_evaluations: int = 0 jacobian_mode: str = "seeded" dense_fallback_used: bool = False nonlinear_block_count: int = 0 nonlinear_block_unknown_count: int = 0 block_fallback_used: bool = False block_fallback_reason: str | None = None residual_verified_this_solve: bool = True causal_fast_path_used: bool = False 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, "residualEvaluations": self.residual_evaluations, "jacobianMode": self.jacobian_mode, "denseFallbackUsed": self.dense_fallback_used, "nonlinearBlockCount": self.nonlinear_block_count, "nonlinearBlockUnknownCount": self.nonlinear_block_unknown_count, "blockFallbackUsed": self.block_fallback_used, "blockFallbackReason": self.block_fallback_reason, "residualVerifiedThisSolve": self.residual_verified_this_solve, "causalFastPathUsed": self.causal_fast_path_used, } 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, scope_kind: str = "network", ) -> None: self.network = network self.residual_tolerance = residual_tolerance self.max_evaluations = max_evaluations self.scope_kind = scope_kind 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._jacobian_sparsity, self._jacobian_sparsity_is_trusted, self._jacobian_sparsity_fallback_reason, ) = self._build_jacobian_sparsity() self._equation_evaluation_locations = ( self._build_equation_evaluation_locations() ) ( self._equation_blocks, self._equation_blocks_are_trusted, self._equation_blocks_fallback_reason, ) = self._build_equation_blocks() ( self._causal_effort_plan_by_variable, self._causal_flow_unknown_ids, self._causal_flow_equation_ids, self._causal_fast_path_eligible, self._causal_fast_path_fallback_reason, ) = self._build_causal_execution_plan() self._causal_fast_path_environment_enabled = ( _causal_fast_path_environment_enabled() ) self._causal_runtime_disabled_reason: str | None = None self._causal_audit_interval = CAUSAL_FAST_PATH_AUDIT_INTERVAL self._causal_audit_required = True self._causal_solves_since_audit = 0 self._causal_fast_solve_count = 0 self._causal_full_residual_audit_count = 0 self._causal_audit_failure_count = 0 self._causal_legacy_fallback_count = 0 self._causal_last_verified_diagnostics: AlgebraicSolveDiagnostics | None = None self.last_diagnostics: AlgebraicSolveDiagnostics | None = None @property def equation_templates(self) -> tuple[EquationResidual, ...]: """Immutable equation metadata compiled for this solver.""" return self._equation_templates @property def jacobian_sparsity(self): """Compiled equation/unknown dependency pattern for finite differences.""" return self._jacobian_sparsity @property def jacobian_sparsity_is_trusted(self) -> bool: """Whether every algebraic row and column has declared structure.""" return self._jacobian_sparsity_is_trusted @property def jacobian_sparsity_fallback_reason(self) -> str | None: """Reason sparse finite differences are disabled for this solver.""" return self._jacobian_sparsity_fallback_reason @property def equation_blocks(self) -> tuple[AlgebraicEquationBlock, ...]: """Trusted connected blocks in the algebraic incidence graph.""" return self._equation_blocks @property def equation_blocks_are_trusted(self) -> bool: """Whether scoped nonlinear fallback is safe for this network.""" return self._equation_blocks_are_trusted @property def equation_blocks_fallback_reason(self) -> str | None: """Reason nonlinear block pruning is disabled for this network.""" return self._equation_blocks_fallback_reason @property def causal_fast_path_eligible(self) -> bool: """Whether the compiled equations form one strictly causal program.""" return self._causal_fast_path_eligible @property def causal_fast_path_enabled(self) -> bool: """Whether new solves may currently use the causal fast path.""" return ( self._causal_fast_path_environment_enabled and self._causal_fast_path_eligible and self._causal_runtime_disabled_reason is None ) def causal_execution_diagnostics(self) -> dict[str, object]: disabled_reason = self._causal_runtime_disabled_reason if not self._causal_fast_path_environment_enabled: disabled_reason = "disabledByEnvironment" elif not self._causal_fast_path_eligible: disabled_reason = self._causal_fast_path_fallback_reason last_verified = self._causal_last_verified_diagnostics return { "eligible": self._causal_fast_path_eligible, "enabled": self.causal_fast_path_enabled, "fallbackReason": self._causal_fast_path_fallback_reason, "disabledReason": disabled_reason, "fastSolveCount": self._causal_fast_solve_count, "fullResidualAuditCount": self._causal_full_residual_audit_count, "auditFailureCount": self._causal_audit_failure_count, "legacyFallbackCount": self._causal_legacy_fallback_count, "auditInterval": self._causal_audit_interval, "solvesSinceAudit": self._causal_solves_since_audit, "lastVerifiedMaxScaledResidual": ( last_verified.max_scaled_residual if last_verified is not None else None ), } def request_causal_audit(self) -> None: """Require the next eligible solve to verify every residual.""" self._causal_audit_required = True def _disable_causal_fast_path(self, reason: str) -> None: if self._causal_runtime_disabled_reason is None: self._causal_runtime_disabled_reason = reason self._causal_audit_required = True def _build_causal_execution_plan( self, ) -> tuple[ dict[str, tuple[CausalEffortAssignment, ...]], frozenset[str], frozenset[str], bool, str | None, ]: """Prove that effort propagation plus staged flow assignment is complete. The fast path is deliberately narrower than ordinary sparse fallback. It accepts only audited built-in residual contracts, one state anchor per effort equality tree, and a flow plan whose dependency stages never use the arbitrary cycle-breaking seed retained by the compatibility solver. """ empty = ({}, frozenset(), frozenset()) def failed(reason: str): return (*empty, False, reason) if self.scope_kind != "network": return failed("nonGlobalAlgebraicScope") if not self._jacobian_sparsity_is_trusted: return failed( self._jacobian_sparsity_fallback_reason or "untrustedAlgebraicStructure" ) if self._unilateral_contact_plan: return failed("activeSetCausalizationRequired") if self._closed_resistance_pressure_plan: return failed("specialClosedResistancePressureSeed") if self._pnor_pnl0001_series_plan: return failed("specialSeriesPressureSeed") if len(self._unknowns_by_id) != len(self.unknowns): return failed("duplicateAlgebraicUnknown") if len({equation.id for equation in self._equation_templates}) != len( self._equation_templates ): return failed("duplicateAlgebraicEquation") effort_unknown_ids = frozenset( unknown.id for unknown in self.unknowns if unknown.role == "effort" ) effort_groups = tuple( group for variable_groups in self._effort_groups.values() for group in variable_groups ) effort_member_ids = frozenset( unknown.id for group in effort_groups for unknown in group.members ) if effort_member_ids != effort_unknown_ids: return failed("incompleteEffortGroupCoverage") if any(len(group.anchors) != 1 for group in effort_groups): return failed("effortGroupDoesNotHaveOneAnchor") group_by_unknown_id = { unknown.id: group for group in effort_groups for unknown in group.members } state_equations_by_group = {id(group): 0 for group in effort_groups} equal_equations_by_group = {id(group): 0 for group in effort_groups} effort_equations = tuple( equation for equation in self._equation_templates if equation.role == "effort" ) for equation in effort_equations: referenced_ids = tuple( dict.fromkeys( variable for variable in equation.variables if variable in effort_unknown_ids ) ) if equation.relation == "state": if len(referenced_ids) != 1: return failed("invalidEffortStateEquation") group = group_by_unknown_id[referenced_ids[0]] if equation.id != group.anchors[0].equation_id: return failed("effortAnchorEquationMismatch") state_equations_by_group[id(group)] += 1 continue if equation.relation == "equal": if len(referenced_ids) != 2: return failed("invalidEffortEqualityEquation") first_group = group_by_unknown_id[referenced_ids[0]] second_group = group_by_unknown_id[referenced_ids[1]] if first_group is not second_group: return failed("crossGroupEffortEquality") equal_equations_by_group[id(first_group)] += 1 continue return failed("unsupportedEffortEquationRelation") for group in effort_groups: if state_equations_by_group[id(group)] != 1: return failed("invalidEffortStateEquationCount") if equal_equations_by_group[id(group)] != len(group.members) - 1: return failed("effortEqualityGroupIsNotATree") effort_plan_by_variable = { variable: tuple( CausalEffortAssignment( variable=variable, members=group.members, anchor=group.anchors[0], ) for group in self._effort_groups[variable] ) for variable in ("p", "x", "v") } flow_unknown_ids = frozenset( unknown.id for unknown in self.unknowns if unknown.role == "flow" ) flow_equations = tuple( equation for equation in self._equation_templates if equation.role == "flow" ) if any( equation.relation not in {"constitutive", "sumToZero"} for equation in flow_equations ): return failed("unsupportedFlowEquationRelation") if len(effort_equations) + len(flow_equations) != len( self._equation_templates ): return failed("unsupportedAlgebraicEquationRole") equation_by_id = { equation.id: equation for equation in self._equation_templates } assigned_unknown_ids: list[str] = [] assigned_equation_ids: list[str] = [] seeded_unknown_ids: set[str] = set() for stage in self._explicit_flow_plan: stage_unknown_ids: set[str] = set() for assignment in stage.assignments: equation = equation_by_id.get(assignment.equation_id) if equation is None or equation.role != "flow": return failed("unknownExplicitFlowEquation") dependency_ids = { unknown.id for unknown in self._flow_unknowns_for_equation(equation) } if dependency_ids - seeded_unknown_ids != {assignment.unknown.id}: return failed("cyclicExplicitFlowPlan") if ( assignment.unknown.id in seeded_unknown_ids or assignment.unknown.id in stage_unknown_ids ): return failed("duplicateExplicitFlowAssignment") stage_unknown_ids.add(assignment.unknown.id) assigned_unknown_ids.append(assignment.unknown.id) assigned_equation_ids.append(assignment.equation_id) seeded_unknown_ids.update(stage_unknown_ids) flow_equation_ids = frozenset( equation.id for equation in flow_equations ) if ( frozenset(assigned_unknown_ids) != flow_unknown_ids or len(assigned_unknown_ids) != len(flow_unknown_ids) ): return failed("incompleteExplicitFlowCoverage") if ( frozenset(assigned_equation_ids) != flow_equation_ids or len(assigned_equation_ids) != len(flow_equation_ids) ): return failed("incompleteExplicitFlowEquationCoverage") return ( effort_plan_by_variable, flow_unknown_ids, flow_equation_ids, True, None, ) def _execute_causal_effort_plan( self, variables: tuple[str, ...], ) -> bool: for variable in variables: assignments = self._causal_effort_plan_by_variable.get(variable) if assignments is None: return False for assignment in assignments: anchor = assignment.anchor target = anchor.unknown.read() - anchor.evaluate() if not isfinite(target) or ( variable == "p" and target <= PRESSURE_LOWER_BOUND_PA ): return False for unknown in assignment.members: unknown.write(target) return True def _causal_audit_is_due(self) -> bool: return ( self._causal_audit_required or self._causal_last_verified_diagnostics is None or self._causal_solves_since_audit >= self._causal_audit_interval ) def _record_causal_audit( self, diagnostics: AlgebraicSolveDiagnostics, ) -> None: self._causal_full_residual_audit_count += 1 self._causal_solves_since_audit = 0 self._causal_audit_required = False self._causal_last_verified_diagnostics = diagnostics def _causal_fast_diagnostics( self, ) -> AlgebraicSolveDiagnostics: verified = self._causal_last_verified_diagnostics if verified is None: raise RuntimeError("Causal execution has no verified residual baseline.") return replace( verified, message=( "Compiled causal pressure-flow program completed; residuals " "reuse the latest full audit." ), evaluations=0, residual_evaluations=0, dense_fallback_used=False, nonlinear_block_count=0, nonlinear_block_unknown_count=0, block_fallback_used=False, block_fallback_reason=None, residual_verified_this_solve=False, causal_fast_path_used=True, ) def _build_jacobian_sparsity(self): """Compile the residual dependency contract into one CSR pattern. Component authoring requires ``EquationResidual.variables`` to list every algebraic port value read by the residual. A missing row/column or a reference to a physical algebraic variable that is not present in this solver makes that contract structurally incomplete, so nonlinear fallback stays on the legacy dense finite-difference path. """ try: from scipy.sparse import csr_matrix except ImportError as exc: raise RuntimeError( "Topology-driven simulation requires SciPy; install requirements.txt." ) from exc unknown_indices = { unknown.id: index for index, unknown in enumerate(self.unknowns) } row_indices: list[int] = [] column_indices: list[int] = [] declared_algebraic_names = {"p", "m_flow", "x", "v", "f"} unresolved_algebraic_reference = False for row_index, equation in enumerate(self._equation_templates): for variable in equation.variables: unknown = self._unknowns_by_id.get(variable) if unknown is not None: row_indices.append(row_index) column_indices.append(unknown_indices[unknown.id]) continue parts = variable.rsplit(".", 2) if len(parts) == 3 and parts[-1] in declared_algebraic_names: unresolved_algebraic_reference = True pattern = csr_matrix( ( [True] * len(row_indices), (row_indices, column_indices), ), shape=(len(self._equation_templates), len(self.unknowns)), dtype=bool, ) if unresolved_algebraic_reference: return pattern, False, "unresolvedAlgebraicVariable" if pattern.shape[0] and any(pattern.getnnz(axis=1) == 0): return pattern, False, "equationWithoutDeclaredUnknown" if pattern.shape[1] and any(pattern.getnnz(axis=0) == 0): return pattern, False, "unknownWithoutDeclaredEquation" if any( not type(item.component).__module__.startswith( "app.simulation.components." ) for item in self._component_equation_plan ): # Built-in component declarations are covered by the repository's # structural and numerical dependency tests. An external component # may omit a read while still leaving every row/column non-empty, # which cannot be detected from shape checks alone. Keep such # networks on the original dense finite-difference contract. return pattern, False, "untrustedCustomComponent" return pattern, True, None def _build_equation_evaluation_locations( self, ) -> tuple[tuple[str, object, int], ...]: locations: list[tuple[str, object, int]] = [] for component_plan in self._component_equation_plan: locations.extend( ("component", component_plan, local_index) for local_index, _template in enumerate(component_plan.templates) ) locations.extend( ("connection", connection_plan, 0) for connection_plan in self._connection_equation_plan ) if len(locations) != len(self._equation_templates): raise RuntimeError("Compiled algebraic evaluation plan is inconsistent.") return tuple(locations) def _build_equation_blocks( self, ) -> tuple[tuple[AlgebraicEquationBlock, ...], bool, str | None]: """Partition a trusted square incidence graph into exact blocks. A scoped optimizer calls component equation evaluators directly. Keep that optimization limited to the audited component package: an external/custom component may have undeclared reads or evaluation side effects even when its structural metadata happens to look complete. """ if not self._jacobian_sparsity_is_trusted: return ( (), False, self._jacobian_sparsity_fallback_reason or "untrustedAlgebraicStructure", ) if any( not type(item.component).__module__.startswith( "app.simulation.components." ) for item in self._component_equation_plan ): return (), False, "untrustedCustomComponent" unknown_index_by_id = { unknown.id: index for index, unknown in enumerate(self.unknowns) } if len(unknown_index_by_id) != len(self.unknowns): return (), False, "duplicateAlgebraicUnknown" equation_unknowns: list[tuple[int, ...]] = [] equations_by_unknown: list[list[int]] = [ [] for _unknown in self.unknowns ] for equation_index, equation in enumerate(self._equation_templates): dependencies = tuple( dict.fromkeys( unknown_index_by_id[variable] for variable in equation.variables if variable in unknown_index_by_id ) ) if not dependencies: return (), False, "equationWithoutDeclaredUnknown" equation_unknowns.append(dependencies) for unknown_index in dependencies: equations_by_unknown[unknown_index].append(equation_index) if any(not attached for attached in equations_by_unknown): return (), False, "unknownWithoutDeclaredEquation" blocks: list[AlgebraicEquationBlock] = [] visited_unknowns: set[int] = set() visited_equations: set[int] = set() for root_unknown in range(len(self.unknowns)): if root_unknown in visited_unknowns: continue block_unknowns: set[int] = set() block_equations: set[int] = set() pending_unknowns = [root_unknown] while pending_unknowns: unknown_index = pending_unknowns.pop() if unknown_index in visited_unknowns: continue visited_unknowns.add(unknown_index) block_unknowns.add(unknown_index) for equation_index in equations_by_unknown[unknown_index]: if equation_index not in visited_equations: visited_equations.add(equation_index) block_equations.add(equation_index) for dependency in equation_unknowns[equation_index]: if dependency not in visited_unknowns: pending_unknowns.append(dependency) unknown_indices = tuple(sorted(block_unknowns)) equation_indices = tuple(sorted(block_equations)) if len(unknown_indices) != len(equation_indices): return (), False, "nonSquareEquationBlock" blocks.append( AlgebraicEquationBlock( unknown_indices=unknown_indices, equation_indices=equation_indices, ) ) if len(visited_equations) != len(self._equation_templates): return (), False, "unreachableEquationBlock" return tuple(blocks), True, None def _compile_equation_subset( self, blocks: tuple[AlgebraicEquationBlock, ...], ) -> AlgebraicEquationSubset: """Compile one batched residual for the union of independent blocks.""" equation_indices = tuple( sorted( equation_index for block in blocks for equation_index in block.equation_indices ) ) unknown_indices = tuple( sorted( unknown_index for block in blocks for unknown_index in block.unknown_indices ) ) equation_target = { equation_index: target for target, equation_index in enumerate(equation_indices) } component_targets: dict[int, list[tuple[int, int, str]]] = {} component_plans: dict[int, ComponentEquationEvaluation] = {} connection_evaluations: list[ScopedConnectionEquationEvaluation] = [] for equation_index in equation_indices: kind, evaluation_plan, source = self._equation_evaluation_locations[ equation_index ] target = equation_target[equation_index] if kind == "component": key = id(evaluation_plan) component_plan = evaluation_plan component_plans[key] = component_plan component_targets.setdefault(key, []).append( ( target, source, self._equation_templates[equation_index].id, ) ) else: connection_plan = evaluation_plan connection_evaluations.append( ScopedConnectionEquationEvaluation( target=target, evaluate=connection_plan.evaluate, ) ) component_evaluations = tuple( ScopedComponentEquationEvaluation( evaluate=component_plans[key].evaluate, targets=tuple(targets), ) for key, targets in component_targets.items() ) jacobian_sparsity = self._jacobian_sparsity[ list(equation_indices), : ][:, list(unknown_indices)].tocsr() if ( any(jacobian_sparsity.getnnz(axis=1) == 0) or any(jacobian_sparsity.getnnz(axis=0) == 0) ): raise RuntimeError("Scoped algebraic Jacobian is structurally empty.") return AlgebraicEquationSubset( unknown_indices=unknown_indices, equation_indices=equation_indices, component_evaluations=component_evaluations, connection_evaluations=tuple(connection_evaluations), jacobian_sparsity=jacobian_sparsity, ) def scale_context(self) -> dict[str, float]: """Capture global normalization data for a compatible block solve.""" scales = self._scales() positive_pressures = [ unknown.read() for unknown in self._unknowns_by_variable["p"] if unknown.read() > 0.0 ] scales["fallback_pressure"] = ( sum(positive_pressures) / len(positive_pressures) if positive_pressures else scales["p"] ) return scales 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), equation_id=equation.id, ) ) 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] ), } @profile_phase("simulation.pressure_flow", minimum_mode="audit") def solve( self, *, effort_variables: tuple[str, ...] = ("p", "x", "v"), scale_context: Mapping[str, float] | None = None, ) -> 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 causal_candidate = ( self.causal_fast_path_enabled and effort_variables == ("p",) and scale_context is None ) causal_audit_due = ( self._causal_audit_is_due() if causal_candidate else False ) for component in self._causal_contact_components: component.clear_causal_contact() if causal_candidate: causal_efforts_are_valid = self._execute_causal_effort_plan( effort_variables ) if not causal_efforts_are_valid: self._causal_legacy_fallback_count += 1 self._disable_causal_fast_path("nonFiniteCausalEffortAnchor") causal_candidate = False causal_audit_due = False self._seed_equal_efforts(effort_variables) else: self._seed_equal_efforts(effort_variables) self._seed_closed_resistance_pressures() self._seed_pnor_pnl0001_series_pressures() seeded_flow_ids = self._solve_explicit_flow_unknowns() contact_bindings = self._seed_unilateral_contacts() if contact_bindings: seeded_flow_ids.update(self._solve_explicit_flow_unknowns(("f",))) self._refresh_unilateral_contacts(contact_bindings) if causal_candidate: causal_unknowns_are_feasible = ( not contact_bindings and seeded_flow_ids == self._causal_flow_unknown_ids and all( isfinite(unknown.read()) and ( unknown.variable != "p" or unknown.read() > PRESSURE_LOWER_BOUND_PA ) for unknown in self.unknowns ) ) if not causal_unknowns_are_feasible: self._causal_legacy_fallback_count += 1 self._disable_causal_fast_path("causalRuntimeGateFailed") causal_candidate = False causal_audit_due = False elif not causal_audit_due: diagnostics = self._causal_fast_diagnostics() self._causal_fast_solve_count += 1 self._causal_solves_since_audit += 1 self.last_diagnostics = diagnostics return diagnostics scales = ( { name: float(scale_context[name]) for name in ("p", "m_flow", "x", "v", "f") } if scale_context is not None else 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 > PRESSURE_LOWER_BOUND_PA ) for unknown, value in seeded_unknown_values ) seeded_state_is_valid = ( seeded_unknowns_are_feasible and all(isfinite(value) for value in seeded_scaled) and seeded_max_scaled_residual <= self.residual_tolerance ) if causal_candidate and causal_audit_due and not seeded_state_is_valid: self._causal_audit_failure_count += 1 self._causal_legacy_fallback_count += 1 self._disable_causal_fast_path("causalResidualAuditFailed") causal_candidate = False causal_audit_due = False if seeded_state_is_valid: 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, ), ) if causal_candidate and causal_audit_due: self._record_causal_audit(diagnostics) 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 = ( float(scale_context["fallback_pressure"]) if scale_context is not None else ( 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( [ ( PRESSURE_LOWER_BOUND_PA / 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)) residual_evaluations = 0 def scaled_residuals(values): nonlocal residual_evaluations residual_evaluations += 1 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, ) optimizer_arguments = { "bounds": (lower, upper), "x_scale": "jac", "ftol": 1e-10, "xtol": 1e-10, "gtol": 1e-10, "max_nfev": self.max_evaluations, } def evaluate_result(result): assign(result.x) self._refresh_unilateral_contacts(contact_bindings) current_equation_values = self._pressure_flow_equation_values() current_scaled = [ abs(value / scale) for value, scale in zip( current_equation_values, equation_scales, ) ] current_max_scaled_residual = max(current_scaled, default=0.0) residuals_converged = ( all(isfinite(value) for value in current_scaled) and current_max_scaled_residual <= self.residual_tolerance ) optimizer_status_is_acceptable = ( bool(result.success) or int(result.status) == 0 ) return ( residuals_converged and optimizer_status_is_acceptable, current_equation_values, current_max_scaled_residual, ) nonlinear_blocks: tuple[AlgebraicEquationBlock, ...] = () nonlinear_block_unknown_count = 0 block_fallback_used = False block_fallback_reason: str | None = None total_optimizer_evaluations = 0 if ( self._equation_blocks_are_trusted and self._jacobian_sparsity_is_trusted ): infeasible_unknown_indices = { index for index, (unknown, value) in enumerate(seeded_unknown_values) if not isfinite(value) or ( unknown.variable == "p" and value <= PRESSURE_LOWER_BOUND_PA ) } nonlinear_blocks = tuple( block for block in self._equation_blocks if any( not isfinite(seeded_scaled[equation_index]) or seeded_scaled[equation_index] > self.residual_tolerance for equation_index in block.equation_indices ) or any( unknown_index in infeasible_unknown_indices for unknown_index in block.unknown_indices ) ) nonlinear_block_unknown_count = sum( len(block.unknown_indices) for block in nonlinear_blocks ) else: block_fallback_reason = ( self._equation_blocks_fallback_reason or self._jacobian_sparsity_fallback_reason or "untrustedAlgebraicStructure" ) if contact_bindings: # Contact activity is the stronger runtime reason even when its # causal projection also makes the static graph appear rectangular. block_fallback_reason = "activeCausalContact" block_subset: AlgebraicEquationSubset | None = None if nonlinear_blocks: if contact_bindings: # Contact causalization mutates coordinates and component-local # active-set caches during a residual call. The first version # deliberately retains the proven global dense path whenever a # contact is active; independent-contact pruning can be added # only with an explicit side-effect dependency contract. block_fallback_reason = "activeCausalContact" elif nonlinear_block_unknown_count >= len(self.unknowns): block_fallback_reason = "fullScopeNonlinearBlock" else: try: block_subset = self._compile_equation_subset( nonlinear_blocks ) except MemoryError: raise except Exception as exc: block_fallback_reason = ( f"blockCompilationFailed:{type(exc).__name__}" ) if block_subset is not None: subset_unknowns = tuple( self.unknowns[index] for index in block_subset.unknown_indices ) subset_x0 = x0[list(block_subset.unknown_indices)] subset_lower = lower[list(block_subset.unknown_indices)] subset_upper = upper[list(block_subset.unknown_indices)] subset_equation_scales = tuple( equation_scales[index] for index in block_subset.equation_indices ) def assign_subset(values) -> None: for unknown, value in zip(subset_unknowns, values): unknown.write(float(value) * variable_scale(unknown)) def scaled_subset_residuals(values): nonlocal residual_evaluations residual_evaluations += 1 assign_subset(values) return np.asarray( [ value / scale for value, scale in zip( block_subset.equation_values(), subset_equation_scales, ) ], dtype=float, ) block_result = None try: block_result = least_squares( scaled_subset_residuals, subset_x0, bounds=(subset_lower, subset_upper), jac_sparsity=block_subset.jacobian_sparsity, x_scale="jac", ftol=1e-10, xtol=1e-10, gtol=1e-10, max_nfev=self.max_evaluations, ) total_optimizer_evaluations += int(block_result.nfev) assign_subset(block_result.x) block_equation_values = self._pressure_flow_equation_values() block_scaled = [ abs(value / scale) for value, scale in zip( block_equation_values, equation_scales, ) ] block_unknowns_are_feasible = all( isfinite(unknown.read()) and ( unknown.variable != "p" or unknown.read() > PRESSURE_LOWER_BOUND_PA ) for unknown in self.unknowns ) block_success = ( block_unknowns_are_feasible and all(isfinite(value) for value in block_scaled) and max(block_scaled, default=0.0) <= self.residual_tolerance and ( bool(block_result.success) or int(block_result.status) == 0 ) ) if block_success: diagnostics = AlgebraicSolveDiagnostics( success=True, message=str(block_result.message), evaluations=total_optimizer_evaluations, pressure_scale=pressure_scale, flow_scale=flow_scale, max_scaled_residual=max(block_scaled, default=0.0), max_raw_residual=max( (abs(value) for value in block_equation_values), default=0.0, ), residual_evaluations=residual_evaluations, jacobian_mode="blockSparse", dense_fallback_used=False, nonlinear_block_count=len(nonlinear_blocks), nonlinear_block_unknown_count=( nonlinear_block_unknown_count ), block_fallback_used=False, block_fallback_reason=None, ) self.last_diagnostics = diagnostics return diagnostics block_fallback_reason = ( "blockResidualNotConverged" if block_result is not None else "blockOptimizerFailed" ) except MemoryError: assign(x0) self._refresh_unilateral_contacts(contact_bindings) raise except Exception as exc: block_fallback_reason = ( f"blockSolveFailed:{type(exc).__name__}" ) except BaseException: assign(x0) self._refresh_unilateral_contacts(contact_bindings) raise # A scoped solve is strictly an optimization. Its candidate must # never seed the compatibility fallback, including candidates from # blocks that converged before another block failed. block_fallback_used = True assign(x0) self._refresh_unilateral_contacts(contact_bindings) elif seeded_max_scaled_residual > self.residual_tolerance or not ( seeded_unknowns_are_feasible and all(isfinite(value) for value in seeded_scaled) ): block_fallback_used = True sparse_is_safe = ( self._jacobian_sparsity_is_trusted and not contact_bindings and self._jacobian_sparsity.nnz > 0 ) dense_fallback_used = False result = None success = False equation_values = seeded_values max_scaled_residual = seeded_max_scaled_residual if sparse_is_safe: try: sparse_result = least_squares( scaled_residuals, x0, jac_sparsity=self._jacobian_sparsity, **optimizer_arguments, ) except (ArithmeticError, RuntimeError, ValueError): sparse_result = None if sparse_result is not None: total_optimizer_evaluations += int(sparse_result.nfev) try: ( success, equation_values, max_scaled_residual, ) = evaluate_result(sparse_result) except (ArithmeticError, RuntimeError, ValueError): success = False result = sparse_result if not success: # The sparse LSMR path is an optimization, not a new numerical # contract. Restart the legacy dense solve from the exact # original seed; a failed sparse candidate must not influence # the fallback result through shared PortState objects. dense_fallback_used = True assign(x0) self._refresh_unilateral_contacts(contact_bindings) result = least_squares( scaled_residuals, x0, **optimizer_arguments, ) total_optimizer_evaluations += int(result.nfev) ( success, equation_values, max_scaled_residual, ) = evaluate_result(result) jacobian_mode = ( "sparseThenDense" if dense_fallback_used else "sparse" ) else: result = least_squares( scaled_residuals, x0, **optimizer_arguments, ) total_optimizer_evaluations = int(result.nfev) ( success, equation_values, max_scaled_residual, ) = evaluate_result(result) jacobian_mode = "dense" assert result is not None diagnostics = AlgebraicSolveDiagnostics( success=success, message=str(result.message), evaluations=total_optimizer_evaluations, 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, ), residual_evaluations=residual_evaluations, jacobian_mode=jacobian_mode, dense_fallback_used=dense_fallback_used, nonlinear_block_count=len(nonlinear_blocks), nonlinear_block_unknown_count=nonlinear_block_unknown_count, block_fallback_used=block_fallback_used, block_fallback_reason=block_fallback_reason, ) self.last_diagnostics = diagnostics if not success: raise AlgebraicSolveError( "Pressure-flow equations did not converge to the requested tolerance.", diagnostics, scope_kind=self.scope_kind, scope_components=tuple(self.network.components), ) return diagnostics