SystemSimulationApp/app/simulation/solvers/algebraic.py

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,
    AmesimPnvo001FixedOpening,
)
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_EXECUTOR_V2_ENVIRONMENT_VARIABLE = "SIMULATION_CAUSAL_EXECUTOR_V2"
CAUSAL_COORDINATE_KERNEL_ENVIRONMENT_VARIABLE = (
    "SIMULATION_CAUSAL_COORDINATE_KERNEL"
)
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"}


def _causal_executor_v2_environment_enabled() -> bool:
    """Return whether the allocation-light causal executor is enabled."""

    value = os.getenv(CAUSAL_EXECUTOR_V2_ENVIRONMENT_VARIABLE, "1")
    return value.strip().lower() not in {"0", "false", "no", "off"}


def _causal_coordinate_kernel_environment_enabled() -> bool:
    """Return whether the canonical-coordinate causal kernel is enabled."""

    value = os.getenv(CAUSAL_COORDINATE_KERNEL_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 CausalEffortKernelTarget:
    """One canonical effort coordinate extracted from a residual evaluator."""

    coordinate_index: int
    assignment: CausalEffortAssignment
    equation_index: int
    equation_id: str


@dataclass(frozen=True)
class CausalEffortKernelEvaluation:
    """One component call shared by every state anchor that it owns."""

    evaluate: Callable[[], tuple[float, ...]]
    targets: tuple[CausalEffortKernelTarget, ...]


@dataclass(frozen=True)
class CausalEffortKernelStage:
    """Precompiled canonical coordinates and compatibility broadcasts."""

    variable: str
    assignments: tuple[tuple[int, CausalEffortAssignment], ...]
    direct_targets: tuple[tuple[int, CausalEffortAssignment], ...]
    component_evaluations: tuple[CausalEffortKernelEvaluation, ...]


@dataclass(frozen=True)
class CausalFlowKernelStage:
    """Map one existing dependency stage into the canonical workspace."""

    stage: ExplicitFlowStage
    coordinate_indices: tuple[int, ...]


@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 ResistancePnl0001SeriesBinding:
    resistance: AmesimPnor001 | AmesimPnvo001FixedOpening
    resistance_port: str
    resistance_other_port: str
    positive_flow_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._resistance_pnl0001_series_plan = (
            self._build_resistance_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._explicit_flow_unknowns_by_variables = {
            selected: tuple(
                unknown
                for unknown in self.unknowns
                if unknown.variable in selected
            )
            for selected in self._explicit_flow_plans_by_variables
        }
        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_executor_v2_environment_enabled = (
            _causal_executor_v2_environment_enabled()
        )
        self._causal_coordinate_kernel_environment_enabled = (
            _causal_coordinate_kernel_environment_enabled()
        )
        self._causal_compiled_effort_unknown_count = sum(
            len(assignment.members)
            for assignments in self._causal_effort_plan_by_variable.values()
            for assignment in assignments
        )
        self._causal_compiled_flow_assignment_count = sum(
            len(stage.assignments) for stage in self._explicit_flow_plan
        )
        self._causal_external_effort_unknowns = tuple(
            member
            for variable in ("x", "v")
            for assignment in self._causal_effort_plan_by_variable.get(
                variable, ()
            )
            for member in assignment.members
        )
        self._causal_external_x_states = tuple(
            unknown.state
            for unknown in self._causal_external_effort_unknowns
            if unknown.variable == "x"
        )
        self._causal_external_v_states = tuple(
            unknown.state
            for unknown in self._causal_external_effort_unknowns
            if unknown.variable == "v"
        )
        (
            self._causal_effort_kernel_by_variable,
            self._causal_flow_kernel_plan,
            self._causal_coordinate_values,
        ) = self._compile_causal_coordinate_kernel()
        self._causal_logical_effort_coordinate_count = sum(
            len(stage.assignments)
            for stage in self._causal_effort_kernel_by_variable.values()
        )
        self._causal_eliminated_effort_alias_count = max(
            self._causal_compiled_effort_unknown_count
            - self._causal_logical_effort_coordinate_count,
            0,
        )
        self._causal_compatibility_scatter_count = (
            self._causal_compiled_effort_unknown_count
            + self._causal_compiled_flow_assignment_count
        )
        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_v2_fast_solve_count = 0
        self._causal_v2_runtime_validation_failure_count = 0
        self._causal_coordinate_fast_solve_count = 0
        self._causal_last_verified_diagnostics: AlgebraicSolveDiagnostics | None = None
        self._causal_cached_fast_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
        )

    @property
    def causal_executor_v2_enabled(self) -> bool:
        """Whether this run may use the v2 executor (environment opt-out)."""

        return (
            self._causal_executor_v2_environment_enabled
            and self.causal_fast_path_enabled
        )

    @property
    def causal_coordinate_kernel_enabled(self) -> bool:
        """Whether the canonical-coordinate executor may run now."""

        return (
            self._causal_coordinate_kernel_environment_enabled
            and self.causal_executor_v2_enabled
            and bool(self._causal_coordinate_values)
        )

    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
            ),
            "executorV2Configured": self._causal_executor_v2_environment_enabled,
            "executorV2Enabled": self.causal_executor_v2_enabled,
            "executorV2FastSolveCount": self._causal_v2_fast_solve_count,
            "executorV2RuntimeValidationFailureCount": (
                self._causal_v2_runtime_validation_failure_count
            ),
            "coordinateKernelConfigured": (
                self._causal_coordinate_kernel_environment_enabled
            ),
            "coordinateKernelEnabled": self.causal_coordinate_kernel_enabled,
            "coordinateKernelFastSolveCount": (
                self._causal_coordinate_fast_solve_count
            ),
            "compiledEffortUnknownCount": (
                self._causal_compiled_effort_unknown_count
            ),
            "compiledFlowAssignmentCount": (
                self._causal_compiled_flow_assignment_count
            ),
            "compiledAssignmentCount": (
                self._causal_compiled_effort_unknown_count
                + self._causal_compiled_flow_assignment_count
            ),
            "logicalEffortCoordinateCount": (
                self._causal_logical_effort_coordinate_count
            ),
            "eliminatedEffortAliasCount": (
                self._causal_eliminated_effort_alias_count
            ),
            "canonicalCoordinateCount": len(self._causal_coordinate_values),
            "compatibilityScatterCount": (
                self._causal_compatibility_scatter_count
            ),
        }

    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._resistance_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 _compile_causal_coordinate_kernel(
        self,
    ) -> tuple[
        dict[str, CausalEffortKernelStage],
        tuple[CausalFlowKernelStage, ...],
        list[float],
    ]:
        """Compile independent coordinates without changing public port state.

        ``PortState`` remains the compatibility surface consumed by component
        methods.  The workspace stores one value per proven effort equality
        group and one per explicit flow assignment; compatibility aliases are
        populated only after every target in an effort stage has been checked.
        """

        if not self._causal_fast_path_eligible:
            return {}, (), []

        equation_index_by_id = {
            equation.id: index
            for index, equation in enumerate(self._equation_templates)
        }
        effort_stages: dict[str, CausalEffortKernelStage] = {}
        next_coordinate = 0
        for variable in ("p", "x", "v"):
            assignments = self._causal_effort_plan_by_variable.get(variable, ())
            indexed_assignments = tuple(
                (next_coordinate + offset, assignment)
                for offset, assignment in enumerate(assignments)
            )
            next_coordinate += len(indexed_assignments)
            direct_targets: list[tuple[int, CausalEffortAssignment]] = []
            targets_by_component: dict[
                int,
                list[CausalEffortKernelTarget],
            ] = {}
            component_evaluators: dict[int, Callable[[], tuple[float, ...]]] = {}
            for coordinate_index, assignment in indexed_assignments:
                equation_index = equation_index_by_id[
                    assignment.anchor.equation_id
                ]
                kind, evaluation_plan, source = (
                    self._equation_evaluation_locations[equation_index]
                )
                if kind != "component":
                    direct_targets.append((coordinate_index, assignment))
                    continue
                component_plan = evaluation_plan
                key = id(component_plan)
                component_evaluators[key] = component_plan.evaluate
                targets_by_component.setdefault(key, []).append(
                    CausalEffortKernelTarget(
                        coordinate_index=coordinate_index,
                        assignment=assignment,
                        equation_index=source,
                        equation_id=assignment.anchor.equation_id,
                    )
                )
            effort_stages[variable] = CausalEffortKernelStage(
                variable=variable,
                assignments=indexed_assignments,
                direct_targets=tuple(direct_targets),
                component_evaluations=tuple(
                    CausalEffortKernelEvaluation(
                        evaluate=component_evaluators[key],
                        targets=tuple(targets),
                    )
                    for key, targets in targets_by_component.items()
                ),
            )

        flow_stages: list[CausalFlowKernelStage] = []
        for stage in self._explicit_flow_plan:
            coordinate_indices = tuple(
                range(next_coordinate, next_coordinate + len(stage.assignments))
            )
            next_coordinate += len(stage.assignments)
            flow_stages.append(
                CausalFlowKernelStage(
                    stage=stage,
                    coordinate_indices=coordinate_indices,
                )
            )
        return effort_stages, tuple(flow_stages), [0.0] * next_coordinate

    @staticmethod
    def _read_effort_anchor(assignment: CausalEffortAssignment) -> float:
        state = assignment.anchor.unknown.state
        if assignment.variable == "p":
            return state.p
        if assignment.variable == "x":
            return state.x
        return state.v

    @staticmethod
    def _scatter_effort_assignment(
        assignment: CausalEffortAssignment,
        value: float,
    ) -> None:
        if assignment.variable == "p":
            for unknown in assignment.members:
                unknown.state.p = value
            return
        if assignment.variable == "x":
            for unknown in assignment.members:
                unknown.state.x = value
            return
        for unknown in assignment.members:
            unknown.state.v = value

    def _execute_causal_coordinate_effort_plan(
        self,
        variables: tuple[str, ...],
    ) -> bool:
        """Evaluate canonical effort coordinates in component-sized batches."""

        workspace = self._causal_coordinate_values
        for variable in variables:
            stage = self._causal_effort_kernel_by_variable.get(variable)
            if stage is None:
                return False
            for coordinate_index, assignment in stage.direct_targets:
                workspace[coordinate_index] = (
                    self._read_effort_anchor(assignment)
                    - assignment.anchor.evaluate()
                )
            for evaluation in stage.component_evaluations:
                equation_values = evaluation.evaluate()
                for target in evaluation.targets:
                    if target.equation_index >= len(equation_values):
                        raise RuntimeError(
                            "Compiled algebraic equation disappeared at runtime: "
                            f"{target.equation_id}."
                        )
                    workspace[target.coordinate_index] = (
                        self._read_effort_anchor(target.assignment)
                        - float(equation_values[target.equation_index])
                    )
            for coordinate_index, assignment in stage.assignments:
                target = workspace[coordinate_index]
                if not isfinite(target) or (
                    variable == "p" and target <= PRESSURE_LOWER_BOUND_PA
                ):
                    return False
            for coordinate_index, assignment in stage.assignments:
                self._scatter_effort_assignment(
                    assignment,
                    workspace[coordinate_index],
                )
        return True

    def _execute_causal_effort_plan(
        self,
        variables: tuple[str, ...],
    ) -> bool:
        if self.causal_coordinate_kernel_enabled:
            return self._execute_causal_coordinate_effort_plan(variables)
        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
        self._causal_cached_fast_diagnostics = replace(
            diagnostics,
            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 _causal_fast_diagnostics(
        self,
    ) -> AlgebraicSolveDiagnostics:
        cached = self._causal_cached_fast_diagnostics
        if cached is None:
            raise RuntimeError("Causal execution has no verified residual baseline.")
        return cached

    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))
        if self.causal_coordinate_kernel_enabled:
            if self._execute_causal_coordinate_effort_plan(variables):
                return
            self._disable_causal_fast_path("nonFiniteCausalEffortAnchor")
        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 _connection_flow_target_reader(
        self,
        equation,
        unknown: AlgebraicUnknown,
    ) -> Callable[[], float]:
        if equation.relation != "sumToZero" or len(equation.variables) != 2:
            raise ValueError(
                f"Connection flow equation {equation.id} must sum two variables."
            )
        first = self._unknowns_by_id[equation.variables[0]]
        second = self._unknowns_by_id[equation.variables[1]]
        if unknown is first:
            return lambda: 0.0 - (0.0 + second.read())
        if unknown is second:
            return lambda: 0.0 - (first.read() + 0.0)
        raise ValueError(
            f"Connection flow equation {equation.id} does not contain {unknown.id}."
        )

    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._connection_flow_target_reader(equation, unknown),
            )

        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, ...]:
        """Evaluate simultaneous targets while this stage's unknowns are zero."""
        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] = 0.0 - 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)
        reset_unknowns = self._explicit_flow_unknowns_by_variables.get(selected)
        if reset_unknowns is None:
            reset_unknowns = tuple(
                unknown
                for unknown in self.unknowns
                if unknown.variable in selected
            )
        for unknown in reset_unknowns:
            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)
            for assignment, target_value in zip(assignments, values):
                if not isfinite(target_value):
                    continue
                assignment.unknown.write(target_value)
                seeded_ids.add(assignment.unknown.id)
        return seeded_ids

    def _execute_compiled_causal_flow_plan(self) -> str | None:
        """Execute the compile-proven full flow plan without coverage sets."""

        if self.causal_coordinate_kernel_enabled:
            return self._execute_causal_coordinate_flow_plan()

        # Position and velocity are propagated by the mechanical reducer
        # before the pressure-only causal solve.  They are therefore not
        # rewritten below, but remain part of the compiled algebraic contract.
        # Validate that small external boundary explicitly instead of restoring
        # the legacy scan over every pressure/flow/force unknown.
        if any(
            not isfinite(unknown.read())
            for unknown in self._causal_external_effort_unknowns
        ):
            return "nonFiniteCausalExternalEffort"

        reset_unknowns = self._explicit_flow_unknowns_by_variables[
            frozenset(("f", "m_flow"))
        ]
        for unknown in reset_unknowns:
            unknown.write(0.0)

        for stage in self._explicit_flow_plan:
            try:
                values = self._evaluate_explicit_flow_stage(stage)
            except MemoryError:
                raise
            except (ArithmeticError, RuntimeError, ValueError) as exc:
                return f"causalFlowEvaluationFailed:{type(exc).__name__}"
            if len(values) != len(stage.assignments):
                return "causalFlowAssignmentCountMismatch"
            for assignment, target_value in zip(stage.assignments, values):
                if not isfinite(target_value):
                    return "nonFiniteCausalFlowAssignment"
                assignment.unknown.write(target_value)
        return None

    def _execute_causal_coordinate_flow_plan(self) -> str | None:
        """Run flow stages through reusable canonical coordinates."""

        if any(not isfinite(state.x) for state in self._causal_external_x_states):
            return "nonFiniteCausalExternalEffort"
        if any(not isfinite(state.v) for state in self._causal_external_v_states):
            return "nonFiniteCausalExternalEffort"

        return self._execute_causal_coordinate_flow_stages(
            self._causal_flow_kernel_plan,
            self._causal_coordinate_values,
        )

    @staticmethod
    def _execute_causal_coordinate_flow_stages(
        kernel_plan: tuple[CausalFlowKernelStage, ...],
        workspace: list[float],
    ) -> str | None:
        """Execute proven flow stages without per-call result containers."""

        # Residual-based explicit assignments use ``-residual`` and therefore
        # require their target coordinate to be zero.  Keep this compatibility
        # initialization until a component exposes a proven direct target op.
        for kernel_stage in kernel_plan:
            for assignment in kernel_stage.stage.assignments:
                if assignment.unknown.variable == "m_flow":
                    assignment.unknown.state.m_flow = 0.0
                else:
                    assignment.unknown.state.f = 0.0

        for kernel_stage in kernel_plan:
            stage = kernel_stage.stage
            coordinate_indices = kernel_stage.coordinate_indices
            if len(coordinate_indices) != len(stage.assignments):
                return "causalFlowAssignmentCountMismatch"
            try:
                for assignment_index, evaluate in stage.direct_evaluations:
                    workspace[coordinate_indices[assignment_index]] = float(
                        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 "
                                f"runtime: {equation_id}."
                            )
                        workspace[coordinate_indices[assignment_index]] = (
                            0.0 - float(equation_values[equation_index])
                        )
            except MemoryError:
                raise
            except (ArithmeticError, RuntimeError, ValueError) as exc:
                return f"causalFlowEvaluationFailed:{type(exc).__name__}"
            for assignment, coordinate_index in zip(
                stage.assignments,
                coordinate_indices,
            ):
                target_value = workspace[coordinate_index]
                if not isfinite(target_value):
                    return "nonFiniteCausalFlowAssignment"
            for assignment, coordinate_index in zip(
                stage.assignments,
                coordinate_indices,
            ):
                target_value = workspace[coordinate_index]
                if assignment.unknown.variable == "m_flow":
                    assignment.unknown.state.m_flow = target_value
                else:
                    assignment.unknown.state.f = target_value
        return None

    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_resistance_pnl0001_series_plan(
        self,
    ) -> tuple[ResistancePnl0001SeriesBinding, ...]:
        bindings: list[ResistancePnl0001SeriesBinding] = []
        resistance_types = (AmesimPnor001, AmesimPnvo001FixedOpening)
        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, resistance_types) and isinstance(
                second, AmesimPnl0001
            ):
                resistance, resistance_port = first, first_endpoint.port
                pipe, pipe_port = second, second_endpoint.port
            elif isinstance(second, resistance_types) and isinstance(
                first, AmesimPnl0001
            ):
                resistance, resistance_port = second, second_endpoint.port
                pipe, pipe_port = first, first_endpoint.port
            else:
                continue
            if isinstance(pipe, AmesimPnl0002) or pipe_port != "port_1":
                continue
            if isinstance(resistance, AmesimPnor001):
                positive_flow_port, negative_flow_port = "port_1", "port_2"
            else:
                positive_flow_port, negative_flow_port = "port_2", "port_3"
            if resistance_port == positive_flow_port:
                resistance_other_port = negative_flow_port
            elif resistance_port == negative_flow_port:
                resistance_other_port = positive_flow_port
            else:
                continue
            bindings.append(
                ResistancePnl0001SeriesBinding(
                    resistance=resistance,
                    resistance_port=resistance_port,
                    resistance_other_port=resistance_other_port,
                    positive_flow_port=positive_flow_port,
                    pipe=pipe,
                    pipe_port=pipe_port,
                )
            )
        return tuple(bindings)

    def _seed_resistance_pnl0001_series_pressures(self) -> None:
        """Causalize pressure between an orifice/valve and a PNL0001 R port."""

        from scipy.optimize import brentq

        for binding in self._resistance_pnl0001_series_plan:
            resistance = binding.resistance
            resistance_port = binding.resistance_port
            pipe = binding.pipe
            pipe_port = binding.pipe_port

            pressure_a = resistance.get_port(binding.resistance_other_port).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 resistance_port == binding.positive_flow_port:
                    resistance_flow_into_connection = resistance.mass_flow(
                        intermediate_pressure,
                        pressure_a,
                    )
                else:
                    resistance_flow_into_connection = -resistance.mass_flow(
                        pressure_a,
                        intermediate_pressure,
                    )
                pipe_flow_into_connection = pipe.mass_flow(
                    intermediate_pressure,
                    pressure_b,
                    pipe.properties().T,
                )
                return resistance_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,
                    )
                )
            resistance.get_port(resistance_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
        )
        causal_v2_candidate = (
            causal_candidate and self._causal_executor_v2_environment_enabled
        )
        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
                causal_v2_candidate = False
                self._seed_equal_efforts(effort_variables)
        else:
            self._seed_equal_efforts(effort_variables)

        seeded_flow_ids: set[str] | None = None
        contact_bindings: tuple[UnilateralContactBinding, ...] = ()
        if causal_v2_candidate:
            v2_failure_reason = self._execute_compiled_causal_flow_plan()
            if v2_failure_reason is not None:
                self._causal_v2_runtime_validation_failure_count += 1
                self._causal_legacy_fallback_count += 1
                self._disable_causal_fast_path(v2_failure_reason)
                causal_candidate = False
                causal_audit_due = False
                causal_v2_candidate = False
                # Rebuild the ordinary seed from scratch in the same solve.
                # A partial compiled stage must never influence fallback.
                self._seed_equal_efforts(effort_variables)
        if not causal_v2_candidate:
            self._seed_closed_resistance_pressures()
            self._seed_resistance_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:
            if causal_v2_candidate:
                # Compilation proves a disjoint, complete effort/flow
                # partition.  The v2 executors validate each produced value,
                # so no coverage set or full unknown scan is needed here.
                causal_unknowns_are_feasible = True
            else:
                assert seeded_flow_ids is not None
                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
                if causal_v2_candidate:
                    self._causal_v2_fast_solve_count += 1
                    if self.causal_coordinate_kernel_enabled:
                        self._causal_coordinate_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