SystemSimulationApp/app/simulation/ir/compiler.py

"""Compile :class:`GenericFluidSystem` into callback-free System IR v2.

The current component implementations are Python reference kernels.  This
compiler therefore records their complete, stable numeric layout and execution
plans, but deliberately marks every component kernel ``reference_only`` until
the C-02 pure numeric kernel contract exists.  Object references and callbacks
are used only while inspecting the already-constructed system; none are stored
in the returned program.
"""

from __future__ import annotations

from collections.abc import Iterable, Mapping, Sequence
from dataclasses import fields, is_dataclass, replace
from math import isfinite
from typing import TYPE_CHECKING, Any

from app.simulation.core.base import Component, DynamicComponent
from app.simulation.solvers.causal_ir import compile_causal_numeric_ir
from app.simulation.solvers.mechanical import MechanicalConstraintGroup

from .schema import (
    CURRENT_SYSTEM_IR_VERSION,
    SYSTEM_NUMERIC_IR_COMPILER_ID,
    SYSTEM_NUMERIC_IR_COMPILER_VERSION,
    IRAlgebraicUnknown,
    IRBlockKind,
    IRBufferKind,
    IRBufferSpec,
    IRCSRMatrix,
    IRCSRPattern,
    IRCacheKind,
    IRCapabilityIssue,
    IRCapabilityLevel,
    IRCapabilityReport,
    IRCausalEffortStageRef,
    IRCausalPlan,
    IRCheckFiniteOperation,
    IRComponentCapability,
    IRComponentInstance,
    IRConnectionRule,
    IRConnectionSpec,
    IRConnectionVariable,
    IRConvergenceSpec,
    IRCopyOperation,
    IRDType,
    IRDiagnosticSeverity,
    IREffortBroadcastOperation,
    IREntryPoint,
    IREntryPointKind,
    IREquationOwner,
    IREquationRelation,
    IREventDirection,
    IREventSpec,
    IRExecutionBlock,
    IRFailurePolicy,
    IRFillOperation,
    IRFiniteDifferenceColumn,
    IRFlowAssignmentOperation,
    IRJacobianPlan,
    IRKernelAvailability,
    IRKernelCallOperation,
    IRKernelPhase,
    IRKernelPhaseSpec,
    IRKernelSpec,
    IRLinearCombinationOperation,
    IRMediumSpec,
    IRModeSpec,
    IRModeValueSpec,
    IROutputSpec,
    IRPortKind,
    IRPortNominalRole,
    IRPortSpec,
    IRPortVariable,
    IRPositiveFlowDirection,
    IRPressureFlowBlock,
    IRPressureFlowEquation,
    IRPressureFlowPlan,
    IRPressureFlowScope,
    IRPressureFlowScopeKind,
    IRSlotRef,
    IRStage,
    IRStageKind,
    IRStateMapKind,
    IRStateMapOperation,
    IRStateReducer,
    IRStepKind,
    IRStepRef,
    IRStreamPlan,
    IRStreamSCC,
    IRThermofluidPlan,
    IRTransactionPlan,
    IRValueSpec,
    IRVariableRole,
    SystemIR,
    operation_read_slots,
    operation_write_slots,
)

if TYPE_CHECKING:
    from app.simulation.systems.generic import GenericFluidSystem


_FLOAT_BUFFERS = frozenset(
    kind
    for kind in IRBufferKind
    if kind not in {IRBufferKind.MODE, IRBufferKind.WORK_INT}
)


def _unique_slots(slots: Iterable[IRSlotRef]) -> tuple[IRSlotRef, ...]:
    return tuple(dict.fromkeys(slots))


class _Slots:
    """Deterministic buffer allocator plus value dictionary builder."""

    def __init__(self) -> None:
        self.initial: dict[IRBufferKind, list[float | int]] = {
            kind: [] for kind in IRBufferKind
        }
        self.values: list[IRValueSpec] = []

    def add(
        self,
        buffer: IRBufferKind,
        value_id: str,
        *,
        initial: float | int = 0.0,
        semantic: str,
        role: str,
        quantity: str = "dimensionless",
        unit: str = "",
        scale: float = 1.0,
        lower_bound: float | None = None,
        upper_bound: float | None = None,
        owner_component_index: int | None = None,
    ) -> IRSlotRef:
        values = self.initial[buffer]
        slot = IRSlotRef(buffer, len(values))
        if buffer in _FLOAT_BUFFERS:
            value = float(initial)
            if not isfinite(value):
                raise ValueError(f"Initial IR value {value_id!r} must be finite.")
            values.append(value)
        else:
            values.append(int(initial))
        numeric_scale = float(scale)
        if not isfinite(numeric_scale) or numeric_scale <= 0.0:
            numeric_scale = 1.0
        self.values.append(
            IRValueSpec(
                value_id=value_id,
                slot=slot,
                semantic=semantic,
                role=role,
                quantity=quantity or "dimensionless",
                unit=unit or "",
                scale=numeric_scale,
                lower_bound=lower_bound,
                upper_bound=upper_bound,
                owner_component_index=owner_component_index,
            )
        )
        return slot

    def buffers(self) -> tuple[IRBufferSpec, ...]:
        result: list[IRBufferSpec] = []
        for kind in IRBufferKind:
            values = self.initial[kind]
            if kind in _FLOAT_BUFFERS:
                result.append(
                    IRBufferSpec(
                        kind=kind,
                        dtype=IRDType.FLOAT64,
                        size=len(values),
                        initial_float_values=tuple(float(item) for item in values),
                    )
                )
            else:
                result.append(
                    IRBufferSpec(
                        kind=kind,
                        dtype=IRDType.INT32,
                        size=len(values),
                        initial_int_values=tuple(int(item) for item in values),
                    )
                )
        return tuple(result)


def _csr_matrix_from_rows(
    row_count: int,
    column_count: int,
    rows: Sequence[Mapping[int, float]],
) -> IRCSRMatrix:
    if len(rows) != row_count:
        raise ValueError("CSR row data does not match the declared row count.")
    row_pointers = [0]
    column_indices: list[int] = []
    values: list[float] = []
    for row in rows:
        for column, value in sorted(row.items()):
            if not 0 <= int(column) < column_count:
                raise ValueError("CSR column is outside the declared matrix shape.")
            numeric = float(value)
            if numeric == 0.0:
                continue
            column_indices.append(int(column))
            values.append(numeric)
        row_pointers.append(len(column_indices))
    return IRCSRMatrix(
        IRCSRPattern(
            row_count=row_count,
            column_count=column_count,
            row_pointers=tuple(row_pointers),
            column_indices=tuple(column_indices),
        ),
        tuple(values),
    )


def _csr_pattern(matrix: object) -> IRCSRPattern:
    csr = matrix.tocsr().astype(bool)
    csr.sort_indices()
    return IRCSRPattern(
        row_count=int(csr.shape[0]),
        column_count=int(csr.shape[1]),
        row_pointers=tuple(int(item) for item in csr.indptr),
        column_indices=tuple(int(item) for item in csr.indices),
    )


def _state_names(component: DynamicComponent) -> tuple[str, ...]:
    size = int(component.state_size)
    model_type = str(component.model_type)
    if model_type == "amesim_mecmas21" and size == 2:
        return ("v", "x")
    if hasattr(component, "medium") and size >= 2 and size % 2 == 0:
        if size == 2:
            return ("m", "U")
        names: list[str] = []
        for partition in range(size // 2):
            names.extend((f"m_{partition + 1}", f"U_{partition + 1}"))
        return tuple(names)
    return tuple(f"state_{index}" for index in range(size))


def _state_metadata(name: str) -> tuple[str, str]:
    if name == "v":
        return "velocity", "m/s"
    if name == "x":
        return "length", "m"
    if name == "m" or name.startswith("m_"):
        return "mass", "kg"
    if name == "U" or name.startswith("U_"):
        return "internal_energy", "J"
    return "dimensionless", ""


def _mechanical_group_key(
    group: MechanicalConstraintGroup,
) -> tuple[str, ...]:
    """Return the topology-stable identity used only during compilation."""

    return tuple(sorted(component.name for component in group.components))


def _component_phases(component: Component) -> tuple[IRKernelPhase, ...]:
    phases = [IRKernelPhase.PRIMAL]
    if component.pressure_flow_equation_residuals():
        phases.append(IRKernelPhase.RESIDUAL)
    if isinstance(component, DynamicComponent):
        phases.extend((IRKernelPhase.PROPERTY, IRKernelPhase.DERIVATIVE))
    if getattr(component, "signal_event_times", None) is not None:
        phases.extend((IRKernelPhase.EVENT, IRKernelPhase.RESET))
    if component.RESULT_VARIABLES or any(
        variable.result_visible
        for port in component.active_port_definitions
        for variable in port.variables
    ):
        phases.append(IRKernelPhase.OUTPUT)
    return tuple(dict.fromkeys(phases))


def _numeric_constants(value: object, prefix: str = "") -> tuple[tuple[str, float], ...]:
    """Flatten stable finite numeric dataclass constants without object identity."""

    result: list[tuple[str, float]] = []
    if is_dataclass(value) and not isinstance(value, type):
        for item in fields(value):
            current = getattr(value, item.name)
            name = f"{prefix}.{item.name}" if prefix else item.name
            if isinstance(current, bool):
                continue
            if isinstance(current, (int, float)) and isfinite(float(current)):
                result.append((name, float(current)))
            elif is_dataclass(current) and not isinstance(current, type):
                result.extend(_numeric_constants(current, name))
    fluid = getattr(value, "fluid", None)
    if fluid is not None and fluid is not value and is_dataclass(fluid):
        result.extend(_numeric_constants(fluid, "fluid"))
    return tuple(dict.fromkeys(result))


def _medium_identity(medium: object) -> tuple[str, str, str]:
    medium_type = type(medium)
    substance = str(getattr(medium_type, "SUBSTANCE_ID", medium_type.__name__))
    method = str(getattr(medium_type, "PROPERTY_METHOD_ID", "reference"))
    implementation = f"{medium_type.__module__}.{medium_type.__qualname__}"
    return f"{substance}:{method}", implementation, str(
        getattr(medium_type, "IMPLEMENTATION_VERSION", "python-reference-1")
    )


def _stage(
    stage_id: str,
    kind: IRStageKind,
    operations: Sequence[Any],
) -> IRStage:
    operation_tuple = tuple(operations)
    return IRStage(
        stage_id=stage_id,
        kind=kind,
        operations=operation_tuple,
        declared_read_slots=_unique_slots(
            slot
            for operation in operation_tuple
            for slot in operation_read_slots(operation)
        ),
        declared_write_slots=_unique_slots(
            slot
            for operation in operation_tuple
            for slot in operation_write_slots(operation)
        ),
    )


def _step(index: int, kind: IRStepKind = IRStepKind.STAGE) -> IRStepRef:
    return IRStepRef(kind, int(index))


def compile_system_ir(
    system: "GenericFluidSystem",
    *,
    model_version: str = "unversioned",
) -> SystemIR:
    """Compile an already-constructed generic system into validated IR v2."""

    # Local import keeps the public IR package independent of GenericFluidSystem
    # construction and prevents a systems.generic -> ir -> systems.generic cycle.
    from app.simulation.systems.generic import GenericFluidSystem

    if not isinstance(system, GenericFluidSystem):
        raise TypeError("compile_system_ir() requires a GenericFluidSystem instance.")

    network = system.network
    components = tuple(network.components.values())
    component_index = {component.name: index for index, component in enumerate(components)}
    connections = tuple(network.connections)
    connection_index = {connection.id: index for index, connection in enumerate(connections)}
    slots = _Slots()
    capability_issues: list[IRCapabilityIssue] = []

    time_slot = slots.add(
        IRBufferKind.TIME,
        "time",
        semantic="independent_variable",
        role="input",
        quantity="time",
        unit="s",
    )

    # The reducer's initialization is the authoritative state order and also
    # performs the same supported initial consistency projection as simulate().
    initial_state = tuple(float(item) for item in system.initial_state_vector())
    state_input_slots: list[IRSlotRef] = []
    derivative_output_slots: list[IRSlotRef] = []
    for index, value in enumerate(initial_state):
        state_input_slots.append(
            slots.add(
                IRBufferKind.STATE_INPUT,
                f"solver.state.{index}",
                initial=value,
                semantic="solver_state",
                role="state",
            )
        )
        derivative_output_slots.append(
            slots.add(
                IRBufferKind.DERIVATIVE_OUTPUT,
                f"solver.derivative.{index}",
                semantic="solver_derivative",
                role="derivative",
            )
        )

    dynamic_components = tuple(system.dynamic_components)
    local_state_by_component: dict[str, tuple[IRSlotRef, ...]] = {}
    local_derivative_by_component: dict[str, tuple[IRSlotRef, ...]] = {}
    local_state_row_by_component: dict[str, tuple[int, ...]] = {}
    local_derivative_column_by_component: dict[str, tuple[int, ...]] = {}
    local_state_slots: list[IRSlotRef] = []
    local_derivative_slots: list[IRSlotRef] = []
    for component in dynamic_components:
        state_values = tuple(float(item) for item in component.get_state_vector())
        names = _state_names(component)
        component_states: list[IRSlotRef] = []
        component_derivatives: list[IRSlotRef] = []
        state_rows: list[int] = []
        derivative_columns: list[int] = []
        for local_index, (name, value) in enumerate(zip(names, state_values)):
            quantity, unit = _state_metadata(name)
            state_rows.append(len(local_state_slots))
            derivative_columns.append(len(local_derivative_slots))
            state_slot = slots.add(
                IRBufferKind.LOCAL_STATE,
                f"{component.name}.state.{name}",
                initial=value,
                semantic="component_state",
                role="state",
                quantity=quantity,
                unit=unit,
                scale=max(abs(value), 1.0),
                owner_component_index=component_index[component.name],
            )
            derivative_slot = slots.add(
                IRBufferKind.LOCAL_DERIVATIVE,
                f"{component.name}.derivative.{name}",
                semantic="component_derivative",
                role="derivative",
                quantity=quantity,
                unit=unit,
                owner_component_index=component_index[component.name],
            )
            component_states.append(state_slot)
            component_derivatives.append(derivative_slot)
            local_state_slots.append(state_slot)
            local_derivative_slots.append(derivative_slot)
        local_state_by_component[component.name] = tuple(component_states)
        local_derivative_by_component[component.name] = tuple(component_derivatives)
        local_state_row_by_component[component.name] = tuple(state_rows)
        local_derivative_column_by_component[component.name] = tuple(derivative_columns)

    # Build exact scatter/gather matrices from the real reduced state entries.
    scatter_rows: list[dict[int, float]] = [
        {} for _ in range(len(local_state_slots))
    ]
    gather_rows: list[dict[int, float]] = [
        {} for _ in range(len(initial_state))
    ]
    solver_cursor = 0
    group_solver_offset: dict[tuple[str, ...], int] = {}
    for entry in system.mechanical_state_reducer.state_entries:
        if isinstance(entry, MechanicalConstraintGroup):
            entry_size = 2
            group_solver_offset[_mechanical_group_key(entry)] = solver_cursor
            representative = entry.representative
            for member in entry.components:
                for local_row, solver_column in zip(
                    local_state_row_by_component[member.name],
                    range(solver_cursor, solver_cursor + entry_size),
                ):
                    scatter_rows[local_row][solver_column] = 1.0
            for output_row, local_column in zip(
                range(solver_cursor, solver_cursor + entry_size),
                local_derivative_column_by_component[representative.name],
            ):
                gather_rows[output_row][local_column] = 1.0
        else:
            entry_size = int(entry.state_size)
            for local_row, solver_column in zip(
                local_state_row_by_component[entry.name],
                range(solver_cursor, solver_cursor + entry_size),
            ):
                scatter_rows[local_row][solver_column] = 1.0
            for output_row, local_column in zip(
                range(solver_cursor, solver_cursor + entry_size),
                local_derivative_column_by_component[entry.name],
            ):
                gather_rows[output_row][local_column] = 1.0
        solver_cursor += entry_size
    if solver_cursor != len(initial_state):
        raise RuntimeError("Reduced state entries do not cover the solver vector.")

    # Ideal C-C storage coupling keeps public state coordinates but projects
    # derivatives by physical volume.  Fold that exact projection into gather.
    for group in system.pneumatic_storage_reducer.groups:
        offsets = [
            int(system.pneumatic_storage_reducer._global_offset(partition))
            for partition in group.partitions
        ]
        volumes = [float(partition.volume) for partition in group.partitions]
        total_volume = sum(volumes)
        # ``GenericFluidSystem.apply_state_vector`` first projects every
        # ideally connected storage partition onto one common mass/energy
        # density, then scatters that projected vector into components.  The
        # IR matrix must include the same projection (not merely the matching
        # derivative projection below), otherwise arbitrary solver trial
        # vectors would reach different component states in Python and IR.
        for partition, volume in zip(group.partitions, volumes):
            target_rows = local_state_row_by_component[partition.component.name]
            fraction = volume / total_volume
            for field_offset in (0, 1):
                target_row = target_rows[partition.state_offset + field_offset]
                scatter_rows[target_row] = {
                    source_offset + field_offset: fraction
                    for source_offset in offsets
                }
        for offset, volume in zip(offsets, volumes):
            fraction = volume / total_volume
            for field_offset in (0, 1):
                combined: dict[int, float] = {}
                for source_offset in offsets:
                    for column, weight in gather_rows[source_offset + field_offset].items():
                        combined[column] = combined.get(column, 0.0) + fraction * weight
                gather_rows[offset + field_offset] = combined

    state_scatter = _csr_matrix_from_rows(
        len(local_state_slots), len(initial_state), scatter_rows
    )
    derivative_gather = _csr_matrix_from_rows(
        len(initial_state), len(local_derivative_slots), gather_rows
    )
    absolute_tolerances = tuple(
        float(item)
        for item in system.mechanical_state_reducer.absolute_tolerances(
            1.0e-8,
            mechanical=1.0e-12,
            mode="legacy",
        )
    )
    state_reducer = IRStateReducer(
        solver_state_count=len(initial_state),
        local_state_slots=tuple(local_state_slots),
        raw_derivative_slots=tuple(local_derivative_slots),
        state_scatter=state_scatter,
        derivative_gather=derivative_gather,
        initial_state=initial_state,
        absolute_tolerances=absolute_tolerances,
    )

    # Normalized component parameter slots.
    parameter_slots_by_component: dict[str, tuple[IRSlotRef, ...]] = {}
    for component in components:
        values = component.parameter_values
        component_slots: list[IRSlotRef] = []
        for definition in component.PARAMETERS:
            value = float(values[definition.name])
            component_slots.append(
                slots.add(
                    IRBufferKind.PARAMETER,
                    f"{component.name}.parameter.{definition.name}",
                    initial=value,
                    semantic="component_parameter",
                    role="parameter",
                    quantity=definition.quantity,
                    unit=definition.unit,
                    scale=max(abs(value), abs(float(definition.default)), 1.0),
                    lower_bound=definition.minimum,
                    upper_bound=definition.maximum,
                    owner_component_index=component_index[component.name],
                )
            )
        parameter_slots_by_component[component.name] = tuple(component_slots)

    # Port and connector tables use actual active-port order from compilation.
    port_specs: list[IRPortSpec] = []
    port_index_by_id: dict[str, int] = {}
    port_slot_by_variable_id: dict[str, IRSlotRef] = {}
    port_slots_by_component: dict[str, list[IRSlotRef]] = {
        component.name: [] for component in components
    }
    port_indices_by_component: dict[str, list[int]] = {
        component.name: [] for component in components
    }
    for component in components:
        owner = component_index[component.name]
        for definition in component.active_port_definitions:
            port = component.get_port(definition.name)
            variables: list[IRPortVariable] = []
            port_id = f"{component.name}.{definition.name}"
            for variable in definition.variables:
                variable_id = f"{port_id}.{variable.name}"
                buffer = (
                    IRBufferKind.SIGNAL
                    if definition.kind == "signal"
                    else IRBufferKind.ALGEBRAIC
                )
                value = float(getattr(port, variable.name))
                lower_bound = 0.0 if variable.name == "p" else None
                slot = slots.add(
                    buffer,
                    variable_id,
                    initial=value,
                    semantic="port_variable",
                    role=variable.role,
                    quantity=variable.quantity,
                    unit=variable.unit,
                    scale=max(abs(value), 1.0),
                    lower_bound=lower_bound,
                    owner_component_index=owner,
                )
                port_slot_by_variable_id[variable_id] = slot
                port_slots_by_component[component.name].append(slot)
                variables.append(
                    IRPortVariable(
                        variable_id=variable_id,
                        name=variable.name,
                        role=IRVariableRole(variable.role),
                        connection_rule=IRConnectionRule(variable.connection_rule),
                        quantity=variable.quantity or "dimensionless",
                        unit=variable.unit or "",
                        result_visible=bool(variable.result_visible),
                        slot=slot,
                    )
                )
            port_index = len(port_specs)
            port_index_by_id[port_id] = port_index
            port_indices_by_component[component.name].append(port_index)
            port_specs.append(
                IRPortSpec(
                    port_id=port_id,
                    component_index=owner,
                    name=definition.name,
                    kind=IRPortKind(definition.kind),
                    domain=definition.domain,
                    nominal_role=IRPortNominalRole(definition.nominal_role),
                    positive_flow_direction=(
                        IRPositiveFlowDirection(definition.positive_flow_direction)
                        if definition.positive_flow_direction is not None
                        else None
                    ),
                    variables=tuple(variables),
                )
            )

    connection_specs: list[IRConnectionSpec] = []
    for connection in connections:
        endpoint_a_id = str(connection.endpoint_a)
        endpoint_b_id = str(connection.endpoint_b)
        endpoint_a_port = port_specs[port_index_by_id[endpoint_a_id]]
        endpoint_b_port = port_specs[port_index_by_id[endpoint_b_id]]
        variables = tuple(
            IRConnectionVariable(
                name=first.name,
                rule=first.connection_rule,
                endpoint_a_slot=first.slot,
                endpoint_b_slot=second.slot,
            )
            for first, second in zip(
                endpoint_a_port.variables,
                endpoint_b_port.variables,
                strict=True,
            )
        )
        connection_specs.append(
            IRConnectionSpec(
                connection_id=connection.id,
                kind=IRPortKind(connection.kind),
                domain=connection.domain,
                endpoint_a_port_index=port_index_by_id[endpoint_a_id],
                endpoint_b_port_index=port_index_by_id[endpoint_b_id],
                variables=variables,
            )
        )

    # Project-scoped media are grouped by semantic implementation and numeric
    # constants, never by Python object address.
    medium_groups: dict[
        tuple[str, str, str, tuple[tuple[str, float], ...]], list[int]
    ] = {}
    medium_objects: dict[
        tuple[str, str, str, tuple[tuple[str, float], ...]], object
    ] = {}
    for index, component in enumerate(components):
        medium = getattr(component, "medium", None)
        if medium is None:
            continue
        medium_id, implementation_id, implementation_version = _medium_identity(medium)
        constants = _numeric_constants(medium)
        key = (medium_id, implementation_id, implementation_version, constants)
        medium_groups.setdefault(key, []).append(index)
        medium_objects.setdefault(key, medium)
    medium_specs: list[IRMediumSpec] = []
    medium_variant_counts: dict[str, int] = {}
    for medium_id, _implementation_id, _implementation_version, _constants in medium_groups:
        medium_variant_counts[medium_id] = medium_variant_counts.get(medium_id, 0) + 1
    medium_variant_offsets: dict[str, int] = {}
    for key, member_indices in medium_groups.items():
        medium_id, implementation_id, implementation_version, constants = key
        variant_offset = medium_variant_offsets.get(medium_id, 0)
        medium_variant_offsets[medium_id] = variant_offset + 1
        compiled_medium_id = (
            medium_id
            if medium_variant_counts[medium_id] == 1
            else f"{medium_id}:variant-{variant_offset + 1}"
        )
        parameter_slots = tuple(
            slots.add(
                IRBufferKind.CONSTANT,
                f"medium.{compiled_medium_id}.{name}",
                initial=value,
                semantic="medium_constant",
                role="constant",
                scale=max(abs(value), 1.0),
            )
            for name, value in constants
        )
        medium = medium_objects[key]
        medium_specs.append(
            IRMediumSpec(
                medium_id=compiled_medium_id,
                name=str(getattr(medium, "name", medium_id)),
                implementation_id=implementation_id,
                implementation_version=implementation_version,
                parameter_slots=parameter_slots,
                component_indices=tuple(member_indices),
            )
        )

    # Discrete mechanical modes are shared by every inertia in one rigid group.
    modes: list[IRModeSpec] = []
    mode_slots_by_component: dict[str, list[IRSlotRef]] = {
        component.name: [] for component in components
    }
    mode_slot_by_group: dict[tuple[str, ...], IRSlotRef] = {}
    event_groups: list[MechanicalConstraintGroup] = []
    for group in system.mechanical_state_reducer.groups:
        if not group.discrete_endstop_components:
            continue
        member_names = tuple(sorted(component.name for component in group.components))
        mode_id = "mechanical_group:" + ",".join(member_names) + ".mode"
        mode_slot = slots.add(
            IRBufferKind.MODE,
            mode_id,
            initial=-1,
            semantic="mechanical_constraint_mode",
            role="mode",
        )
        mode_slot_by_group[_mechanical_group_key(group)] = mode_slot
        event_groups.append(group)
        owners = tuple(component_index[name] for name in member_names)
        for name in member_names:
            mode_slots_by_component[name].append(mode_slot)
        modes.append(
            IRModeSpec(
                mode_id=mode_id,
                slot=mode_slot,
                owner_component_indices=owners,
                values=(
                    IRModeValueSpec(-1, "uninitialized"),
                    IRModeValueSpec(0, "free"),
                    IRModeValueSpec(1, "lower"),
                    IRModeValueSpec(2, "upper"),
                ),
                initial_value=-1,
            )
        )

    # Allocate result projection sources before component records so each
    # instance can point at its immutable output indices.
    result_metadata = network.result_variable_metadata()
    output_indices_by_component: dict[str, list[int]] = {
        component.name: [] for component in components
    }
    output_source_slots: list[IRSlotRef] = []
    output_slots: list[IRSlotRef] = []
    for output_index, metadata in enumerate(result_metadata):
        output_indices_by_component[metadata.component_id].append(output_index)
        if metadata.scope == "port":
            source_id = f"{metadata.component_id}.{metadata.port_name}.{metadata.name}"
            source_slot = port_slot_by_variable_id[source_id]
        else:
            source_slot = slots.add(
                IRBufferKind.ALGEBRAIC,
                f"{metadata.key}.output_source",
                semantic="component_result_source",
                role="output_source",
                quantity=metadata.quantity,
                unit=metadata.unit,
                owner_component_index=component_index[metadata.component_id],
            )
        output_source_slots.append(source_slot)
        output_slots.append(
            slots.add(
                IRBufferKind.RESULT_OUTPUT,
                f"result.{metadata.key}",
                semantic="result_output",
                role="output",
                quantity=metadata.quantity,
                unit=metadata.unit,
                owner_component_index=component_index[metadata.component_id],
            )
        )

    # One reference kernel per component implementation plus an explicit
    # reference-only system executor for topology-wide operations.
    kernels: list[IRKernelSpec] = []
    component_kernel_index: dict[str, int] = {}
    components_by_kernel: dict[
        tuple[str, str, str, int, int, int], list[Component]
    ] = {}
    for component in components:
        component_type = type(component)
        model_type = str(component.model_type)
        model_model_version = str(component_type.MODEL_VERSION or "unversioned")
        implementation = f"{component_type.__module__}.{component_type.__qualname__}"
        key = (
            model_type,
            model_model_version,
            implementation,
            len(component.PARAMETERS),
            int(component.state_size) if isinstance(component, DynamicComponent) else 0,
            len(mode_slots_by_component[component.name]),
        )
        components_by_kernel.setdefault(key, []).append(component)
    for key, instances in components_by_kernel.items():
        (
            model_type,
            model_model_version,
            implementation,
            parameter_count,
            state_count,
            mode_count,
        ) = key
        phases = tuple(
            dict.fromkeys(
                phase for component in instances for phase in _component_phases(component)
            )
        )
        kernel_index = len(kernels)
        kernels.append(
            IRKernelSpec(
                kernel_id=(
                    f"reference:{model_type}@{model_model_version}"
                    f":{implementation}:p{parameter_count}:s{state_count}:m{mode_count}"
                ),
                model_type=model_type,
                model_version=model_model_version,
                implementation_version=implementation,
                availability=IRKernelAvailability.REFERENCE_ONLY,
                unavailable_reason="C-02 pure numeric kernel contract is not declared.",
                phases=tuple(IRKernelPhaseSpec(phase) for phase in phases),
                parameter_count=parameter_count,
                state_count=state_count,
                mode_count=mode_count,
                workspace_float_count=0,
                workspace_int_count=0,
            )
        )
        for component in instances:
            component_kernel_index[component.name] = kernel_index
    system_kernel_index = len(kernels)
    kernels.append(
        IRKernelSpec(
            kernel_id="reference:generic-fluid-system@2",
            model_type="generic_fluid_system",
            model_version="2.0.0",
            implementation_version=SYSTEM_NUMERIC_IR_COMPILER_VERSION,
            availability=IRKernelAvailability.REFERENCE_ONLY,
            unavailable_reason="Whole-system native execution is not implemented.",
            phases=tuple(
                IRKernelPhaseSpec(phase)
                for phase in (
                    IRKernelPhase.PRIMAL,
                    IRKernelPhase.RESIDUAL,
                    IRKernelPhase.DERIVATIVE,
                    IRKernelPhase.EVENT,
                    IRKernelPhase.RESET,
                    IRKernelPhase.JACOBIAN,
                    IRKernelPhase.OUTPUT,
                )
            ),
            parameter_count=0,
            state_count=len(initial_state),
            mode_count=len(modes),
            workspace_float_count=0,
            workspace_int_count=0,
        )
    )

    component_instances = tuple(
        IRComponentInstance(
            instance_id=component.name,
            kernel_index=component_kernel_index[component.name],
            parameter_slots=parameter_slots_by_component[component.name],
            state_slots=local_state_by_component.get(component.name, ()),
            derivative_slots=local_derivative_by_component.get(component.name, ()),
            mode_slots=tuple(mode_slots_by_component[component.name]),
            port_indices=tuple(port_indices_by_component[component.name]),
            port_slots=tuple(port_slots_by_component[component.name]),
            output_indices=tuple(output_indices_by_component[component.name]),
            workspace_float_slots=(),
            workspace_int_slots=(),
        )
        for component in components
    )

    # Pressure-flow structure and work slots.
    pressure_solver = system.pressure_flow_solver
    algebraic_scales = pressure_solver._scales()
    unknown_index_by_id = {
        unknown.id: index for index, unknown in enumerate(pressure_solver.unknowns)
    }
    algebraic_unknowns: list[IRAlgebraicUnknown] = []
    for unknown in pressure_solver.unknowns:
        scale = float(algebraic_scales.get(unknown.variable, 1.0))
        algebraic_unknowns.append(
            IRAlgebraicUnknown(
                unknown_id=unknown.id,
                component_index=component_index[unknown.component],
                port_index=port_index_by_id[f"{unknown.component}.{unknown.port}"],
                variable=unknown.variable,
                role=IRVariableRole(unknown.role),
                slot=port_slot_by_variable_id[unknown.id],
                scale=max(scale, 1.0e-300),
                lower_bound=0.0 if unknown.variable == "p" else None,
                upper_bound=None,
            )
        )

    equation_index_by_id = {
        equation.id: index
        for index, equation in enumerate(pressure_solver.equation_templates)
    }
    residual_slots: list[IRSlotRef] = []
    pressure_equations: list[IRPressureFlowEquation] = []
    for equation_index_value, equation in enumerate(pressure_solver.equation_templates):
        residual_slot = slots.add(
            IRBufferKind.WORK_FLOAT,
            f"pressure_flow.residual.{equation.id}",
            semantic="pressure_flow_residual",
            role="residual",
        )
        residual_slots.append(residual_slot)
        variable_slots: list[IRSlotRef] = []
        for variable_id in equation.variables:
            slot = port_slot_by_variable_id.get(variable_id)
            if slot is not None:
                variable_slots.append(slot)
                continue
            if variable_id.endswith(".state"):
                owner_name = variable_id[: -len(".state")]
                variable_slots.extend(local_state_by_component.get(owner_name, ()))
                continue
            capability_issues.append(
                IRCapabilityIssue(
                    code="IR_EQUATION_VARIABLE_UNRESOLVED",
                    severity=IRDiagnosticSeverity.ERROR,
                    scope_id=equation.id,
                    message=f"Equation variable {variable_id!r} has no numeric slot.",
                )
            )
        if equation.owner == "component":
            owner = component_index[equation.owner_id]
        else:
            owner = connection_index[equation.owner_id]
        if equation.role == "flow":
            scale = algebraic_scales.get("m_flow", 1.0)
            if any(variable.endswith(".f") for variable in equation.variables):
                scale = algebraic_scales.get("f", 1.0)
        elif equation.role == "effort":
            suffixes = {item.rsplit(".", 1)[-1] for item in equation.variables}
            scale = (
                algebraic_scales.get("x", 1.0)
                if "x" in suffixes
                else algebraic_scales.get("v", 1.0)
                if "v" in suffixes
                else algebraic_scales.get("p", 1.0)
            )
        else:
            scale = 1.0
        pressure_equations.append(
            IRPressureFlowEquation(
                equation_id=equation.id,
                owner=IREquationOwner(equation.owner),
                owner_index=owner,
                relation=IREquationRelation(equation.relation),
                role=IRVariableRole(equation.role) if equation.role is not None else None,
                variable_slots=tuple(variable_slots),
                residual_slot=residual_slot,
                scale=max(float(scale), 1.0e-300),
            )
        )

    pressure_blocks: list[IRPressureFlowBlock] = []
    block_key_to_index: dict[tuple[tuple[int, ...], tuple[int, ...]], int] = {}

    def add_pressure_block(
        unknown_indices: Sequence[int],
        equation_indices: Sequence[int],
        pattern: object,
    ) -> int:
        unknown_tuple = tuple(int(item) for item in unknown_indices)
        equation_tuple = tuple(int(item) for item in equation_indices)
        key = (unknown_tuple, equation_tuple)
        existing = block_key_to_index.get(key)
        if existing is not None:
            return existing
        index = len(pressure_blocks)
        block_key_to_index[key] = index
        pressure_blocks.append(
            IRPressureFlowBlock(
                block_id=f"pressure_flow.block.{index}",
                unknown_indices=unknown_tuple,
                equation_indices=equation_tuple,
                jacobian_pattern=_csr_pattern(pattern),
            )
        )
        return index

    for block in pressure_solver.equation_blocks:
        pattern = pressure_solver.jacobian_sparsity[
            list(block.equation_indices), :
        ][:, list(block.unknown_indices)]
        add_pressure_block(block.unknown_indices, block.equation_indices, pattern)

    # Preserve the callback-free structural half of causal IR v1.
    causal_compilation = compile_causal_numeric_ir(pressure_solver)
    causal_plans: list[IRCausalPlan] = []
    causal_program = causal_compilation.ir.program if causal_compilation.supported else None
    causal_canonical_slots: list[IRSlotRef] = []
    if causal_program is not None:
        for canonical in causal_program.canonical_slots:
            causal_canonical_slots.append(
                slots.add(
                    IRBufferKind.WORK_FLOAT,
                    f"causal_v1.canonical.{canonical.id}",
                    semantic="causal_coordinate",
                    role=canonical.kind,
                )
            )

    # Main and causal stages are assembled together so stage references in the
    # causal plan are ordinary v2 indices.
    stages: list[IRStage] = []

    def add_stage(stage_id: str, kind: IRStageKind, operations: Sequence[Any]) -> int:
        index = len(stages)
        stages.append(_stage(stage_id, kind, operations))
        return index

    state_stage = add_stage(
        "rhs.state_scatter",
        IRStageKind.STATE_REDUCE,
        (
            IRStateMapOperation(
                IRStateMapKind.SCATTER,
                tuple(state_input_slots),
                tuple(local_state_slots),
            ),
        ),
    )

    # Signal output functions, followed by directed connection propagation.
    signal_operations: list[Any] = []
    for binding in system.signal_resolver._output_bindings:
        component = binding.component
        writes = tuple(
            variable.slot
            for port_index_value in port_indices_by_component[component.name]
            for variable in port_specs[port_index_value].variables
            if port_specs[port_index_value].kind is IRPortKind.SIGNAL
            and port_specs[port_index_value].nominal_role is IRPortNominalRole.OUTPUT
        )
        signal_operations.append(
            IRKernelCallOperation(
                kernel_index=component_kernel_index[component.name],
                component_index=component_index[component.name],
                phase=IRKernelPhase.PRIMAL,
                read_slots=(time_slot, *parameter_slots_by_component[component.name]),
                write_slots=writes,
            )
        )
    for connection in connections:
        if connection.kind != "signal":
            continue
        endpoint_a = port_specs[port_index_by_id[str(connection.endpoint_a)]]
        endpoint_b = port_specs[port_index_by_id[str(connection.endpoint_b)]]
        source, target = (
            (endpoint_a, endpoint_b)
            if endpoint_a.nominal_role is IRPortNominalRole.OUTPUT
            else (endpoint_b, endpoint_a)
        )
        signal_operations.append(IRCopyOperation(source.variables[0].slot, target.variables[0].slot))
    signal_stage = add_stage("rhs.signal", IRStageKind.SIGNAL, signal_operations)

    effort_slots = tuple(
        unknown.slot
        for unknown in algebraic_unknowns
        if unknown.variable in {"x", "v"}
    )
    mechanical_equivalence_stage = add_stage(
        "rhs.mechanical_equivalence",
        IRStageKind.MECHANICAL_EQUIVALENCE,
        (
            IRKernelCallOperation(
                system_kernel_index,
                None,
                IRKernelPhase.PRIMAL,
                tuple(local_state_slots),
                effort_slots,
            ),
        ),
    )

    pneumatic_volume_slots = tuple(
        variable.slot
        for port in port_specs
        if port.kind is IRPortKind.PHYSICAL and port.domain == "pneumatic"
        for variable in port.variables
        if variable.name in {"volume", "volume_flow"}
    )
    volume_operations: list[Any] = [IRFillOperation(pneumatic_volume_slots, 0.0)]
    volume_output_names: dict[str, tuple[str, ...]] = {}
    for component in system.pneumatic_volume_resolver._output_components:
        output_names = tuple(component.pneumatic_volume_outputs())
        volume_output_names[component.name] = output_names
        writes = tuple(
            port_slot_by_variable_id[f"{component.name}.{port_name}.{variable}"]
            for port_name in output_names
            for variable in ("volume", "volume_flow")
        )
        volume_operations.append(
            IRKernelCallOperation(
                component_kernel_index[component.name],
                component_index[component.name],
                IRKernelPhase.PRIMAL,
                (
                    *local_state_by_component.get(component.name, ()),
                    *parameter_slots_by_component[component.name],
                    *port_slots_by_component[component.name],
                ),
                writes,
            )
        )
    connected_endpoint: dict[str, str] = {}
    for connection in connections:
        if connection.kind == "physical" and connection.domain == "pneumatic":
            connected_endpoint[str(connection.endpoint_a)] = str(connection.endpoint_b)
            connected_endpoint[str(connection.endpoint_b)] = str(connection.endpoint_a)
    for component_name, port_names in volume_output_names.items():
        for port_name in port_names:
            source_endpoint = f"{component_name}.{port_name}"
            target_endpoint = connected_endpoint.get(source_endpoint)
            if target_endpoint is None:
                continue
            for variable in ("volume", "volume_flow"):
                volume_operations.append(
                    IRCopyOperation(
                        port_slot_by_variable_id[f"{source_endpoint}.{variable}"],
                        port_slot_by_variable_id[f"{target_endpoint}.{variable}"],
                    )
                )
    volume_stage = add_stage("rhs.dynamic_volume", IRStageKind.DYNAMIC_VOLUME, volume_operations)

    property_operations: list[Any] = []
    for component in dynamic_components:
        writes = tuple(
            port_slot_by_variable_id[f"{component.name}.{definition.name}.{variable.name}"]
            for definition in component.active_port_definitions
            if definition.kind == "physical" and definition.domain == "pneumatic"
            for variable in definition.variables
            if variable.name in {"p", "h_outflow"}
        )
        property_operations.append(
            IRKernelCallOperation(
                component_kernel_index[component.name],
                component_index[component.name],
                IRKernelPhase.PROPERTY,
                (
                    *local_state_by_component[component.name],
                    *parameter_slots_by_component[component.name],
                    *port_slots_by_component[component.name],
                ),
                writes,
            )
        )
    property_stage = add_stage("rhs.property_bundle", IRStageKind.PROPERTY, property_operations)

    all_unknown_slots = tuple(item.slot for item in algebraic_unknowns)
    # Reference pressure-flow evaluation still reaches through the object model
    # to parameters, modes, medium constants and every physical port value.  C-01
    # records a conservative dependency superset so cache invalidation remains
    # safe until C-02 replaces this system call with pure component signatures.
    pressure_flow_model_read_slots = _unique_slots(
        IRSlotRef(kind, index)
        for kind in (
            IRBufferKind.TIME,
            IRBufferKind.LOCAL_STATE,
            IRBufferKind.ALGEBRAIC,
            IRBufferKind.SIGNAL,
            IRBufferKind.PARAMETER,
            IRBufferKind.CONSTANT,
            IRBufferKind.MODE,
            IRBufferKind.RUNTIME_INPUT,
        )
        for index in range(len(slots.initial[kind]))
    )
    global_pressure_stage = add_stage(
        "rhs.pressure_flow.global",
        IRStageKind.PRESSURE_FLOW,
        (
            IRKernelCallOperation(
                system_kernel_index,
                None,
                IRKernelPhase.RESIDUAL,
                pressure_flow_model_read_slots,
                (*all_unknown_slots, *residual_slots),
                tuple(range(len(pressure_equations))),
            ),
            IRCheckFiniteOperation(all_unknown_slots, "PRESSURE_FLOW_NONFINITE"),
        ),
    )

    # Causal v1 structural stages are intentionally not Python bindings.
    causal_effort_refs: list[IRCausalEffortStageRef] = []
    causal_flow_stage_indices: list[int] = []
    if causal_program is not None:
        for source_stage in causal_program.effort_stages:
            operations = tuple(
                IREffortBroadcastOperation(
                    variable=operation.variable,
                    anchor_slot=port_slot_by_variable_id[
                        causal_program.compatibility_slots[
                            operation.anchor_compatibility_slot
                        ].id
                    ],
                    residual_slot=residual_slots[equation_index_by_id[operation.equation_id]],
                    result_slot=causal_canonical_slots[operation.result_slot],
                    scatter_slots=tuple(
                        port_slot_by_variable_id[
                            causal_program.compatibility_slots[index].id
                        ]
                        for index in operation.scatter_compatibility_slots
                    ),
                    equation_id=operation.equation_id,
                    lower_bound=0.0 if operation.variable == "p" else None,
                )
                for operation in source_stage.operations
            )
            stage_index = add_stage(
                f"causal_v1.effort.{source_stage.variable}",
                IRStageKind.PRESSURE_FLOW,
                operations,
            )
            causal_effort_refs.append(
                IRCausalEffortStageRef(source_stage.variable, stage_index)
            )
        for flow_index, source_stage in enumerate(causal_program.flow_stages):
            operations = tuple(
                IRFlowAssignmentOperation(
                    value_slot=residual_slots[equation_index_by_id[equation_id]],
                    result_slot=causal_canonical_slots[target],
                    scatter_slots=(
                        port_slot_by_variable_id[
                            causal_program.compatibility_slots[compatibility].id
                        ],
                    ),
                    equation_id=equation_id,
                )
                for target, compatibility, equation_id in zip(
                    source_stage.target_slots,
                    source_stage.scatter_compatibility_slots,
                    source_stage.equation_ids,
                    strict=True,
                )
            )
            causal_flow_stage_indices.append(
                add_stage(
                    f"causal_v1.flow.{flow_index}",
                    IRStageKind.PRESSURE_FLOW,
                    operations,
                )
            )
        causal_plans.append(
            IRCausalPlan(
                plan_id="pressure_flow.global.causal_v1",
                scope_component_indices=tuple(range(len(components))),
                source_schema_version=int(causal_program.schema_version),
                source_structural_signature=causal_program.structural_signature,
                fallback_reason=None,
                canonical_slots=tuple(causal_canonical_slots),
                compatibility_slots=tuple(
                    port_slot_by_variable_id[item.id]
                    for item in causal_program.compatibility_slots
                ),
                reset_slots=tuple(
                    port_slot_by_variable_id[
                        causal_program.compatibility_slots[index].id
                    ]
                    for index in causal_program.reset_compatibility_slots
                ),
                external_effort_slots=tuple(
                    port_slot_by_variable_id[
                        causal_program.compatibility_slots[index].id
                    ]
                    for index in causal_program.external_effort_compatibility_slots
                ),
                effort_stages=tuple(causal_effort_refs),
                flow_stage_indices=tuple(causal_flow_stage_indices),
            )
        )
    else:
        causal_plans.append(
            IRCausalPlan(
                plan_id="pressure_flow.global.causal_v1",
                scope_component_indices=tuple(range(len(components))),
                source_schema_version=1,
                source_structural_signature=None,
                fallback_reason=causal_compilation.fallback_reason,
                canonical_slots=(),
                compatibility_slots=all_unknown_slots,
                reset_slots=(),
                external_effort_slots=(),
                effort_stages=(),
                flow_stage_indices=(),
            )
        )
        capability_issues.append(
            IRCapabilityIssue(
                code="IR_CAUSAL_PRESSURE_FLOW_FALLBACK",
                severity=IRDiagnosticSeverity.WARNING,
                scope_id="pressure_flow.global",
                message=(
                    "Causal IR v1 is unavailable; the Python reference path may "
                    "fall back to finite-difference least_squares: "
                    f"{causal_compilation.fallback_reason or 'unspecified reason'}."
                ),
            )
        )

    # Stream operations mirror the current resolver: non-dynamic components
    # iterate to an enthalpy fixed point; dynamic stream hooks run afterwards.
    connected_h_slots_by_component: dict[str, list[IRSlotRef]] = {
        component.name: [] for component in components
    }
    for connection in connections:
        if connection.kind != "physical":
            continue
        first, second = str(connection.endpoint_a), str(connection.endpoint_b)
        first_h = port_slot_by_variable_id.get(f"{first}.h_outflow")
        second_h = port_slot_by_variable_id.get(f"{second}.h_outflow")
        if first_h is not None and second_h is not None:
            connected_h_slots_by_component[connection.endpoint_a.component].append(second_h)
            connected_h_slots_by_component[connection.endpoint_b.component].append(first_h)
    stream_node_slots = tuple(
        variable.slot
        for port in port_specs
        if port.kind is IRPortKind.PHYSICAL
        for variable in port.variables
        if variable.name == "h_outflow"
    )

    def stream_operations(selected: Sequence[Component]) -> list[Any]:
        operations: list[Any] = []
        for component in selected:
            writes = tuple(
                variable.slot
                for port_index_value in port_indices_by_component[component.name]
                for variable in port_specs[port_index_value].variables
                if variable.name == "h_outflow"
            )
            operations.append(
                IRKernelCallOperation(
                    component_kernel_index[component.name],
                    component_index[component.name],
                    IRKernelPhase.PRIMAL,
                    tuple(connected_h_slots_by_component[component.name]),
                    writes,
                )
            )
        return operations

    stream_nondynamic_stage = add_stage(
        "rhs.stream.reference_fixed_point",
        IRStageKind.STREAM,
        stream_operations(system.stream_resolver._non_dynamic_components),
    )
    stream_dynamic_stage = add_stage(
        "rhs.stream.dynamic_outflows",
        IRStageKind.STREAM,
        stream_operations(dynamic_components),
    )
    temperature_reference_stage = add_stage(
        "rhs.stream.temperature_reference",
        IRStageKind.TEMPERATURE_REFERENCE,
        tuple(
            IRKernelCallOperation(
                component_kernel_index[component.name],
                component_index[component.name],
                IRKernelPhase.PRIMAL,
                tuple(connected_h_slots_by_component[component.name]),
                tuple(port_slots_by_component[component.name]),
            )
            for component in system.stream_resolver._flow_temperature_reference_components
        ),
    )

    # Pressure scopes: global plus exact selected blocks used after stream.
    scopes: list[IRPressureFlowScope] = []
    scopes.append(
        IRPressureFlowScope(
            scope_id="pressure_flow.global",
            kind=IRPressureFlowScopeKind.NETWORK,
            component_indices=tuple(range(len(components))),
            unknown_indices=tuple(range(len(algebraic_unknowns))),
            equation_indices=tuple(range(len(pressure_equations))),
            block_indices=tuple(range(len(pressure_blocks))),
            causal_plan_index=0,
            residual_tolerance=float(pressure_solver.residual_tolerance),
            max_evaluations=int(pressure_solver.max_evaluations),
            sparse_pattern_trusted=bool(pressure_solver.jacobian_sparsity_is_trusted),
            sparse_fallback_reason=pressure_solver.jacobian_sparsity_fallback_reason,
        )
    )
    if not pressure_solver.jacobian_sparsity_is_trusted:
        capability_issues.append(
            IRCapabilityIssue(
                code="IR_PRESSURE_FLOW_FINITE_DIFFERENCE_FALLBACK",
                severity=IRDiagnosticSeverity.WARNING,
                scope_id="pressure_flow.global",
                message=(
                    "The reference nonlinear solve cannot trust its structural "
                    "Jacobian and may use finite-difference least_squares: "
                    f"{pressure_solver.jacobian_sparsity_fallback_reason or 'unspecified reason'}."
                ),
            )
        )
    secondary_scope_indices: list[int] = []
    closure_plan = system._thermofluid_closure_plan
    if closure_plan.uses_conservative_global_solver:
        secondary_scope_indices.append(0)
        capability_issues.append(
            IRCapabilityIssue(
                code="IR_THERMOFLUID_CONSERVATIVE_GLOBAL_FALLBACK",
                severity=IRDiagnosticSeverity.WARNING,
                scope_id="thermofluid",
                message=(
                    "Stream-sensitive closure reuses the whole-network pressure-flow "
                    "solver: "
                    f"{closure_plan.conservative_fallback_reason or 'unspecified reason'}."
                ),
            )
        )
    else:
        for block_solver in closure_plan.secondary_block_solvers:
            if not block_solver.available:
                continue
            for block in block_solver.blocks:
                unknown_indices = tuple(
                    unknown_index_by_id[unknown.id] for unknown in block.unknowns
                )
                equation_indices = tuple(
                    equation_index_by_id[equation.id] for equation in block.equations
                )
                # block.jacobian_entries is already local row/column structure.
                rows: list[dict[int, float]] = [
                    {} for _ in range(len(equation_indices))
                ]
                for row, column in block.jacobian_entries:
                    rows[int(row)][int(column)] = 1.0
                local_pattern = _csr_matrix_from_rows(
                    len(equation_indices), len(unknown_indices), rows
                ).pattern
                block_index = len(pressure_blocks)
                key = (unknown_indices, equation_indices)
                existing = block_key_to_index.get(key)
                if existing is None:
                    block_key_to_index[key] = block_index
                    pressure_blocks.append(
                        IRPressureFlowBlock(
                            f"pressure_flow.secondary_block.{block_index}",
                            unknown_indices,
                            equation_indices,
                            local_pattern,
                        )
                    )
                else:
                    block_index = existing
                scope_index = len(scopes)
                secondary_scope_indices.append(scope_index)
                scopes.append(
                    IRPressureFlowScope(
                        scope_id=f"pressure_flow.secondary.{scope_index}",
                        kind=IRPressureFlowScopeKind.EQUATION_BLOCK,
                        component_indices=tuple(
                            component_index[name] for name in block.scope_components
                        ),
                        unknown_indices=unknown_indices,
                        equation_indices=equation_indices,
                        block_indices=(block_index,),
                        causal_plan_index=None,
                        residual_tolerance=float(
                            block_solver.pressure_flow_solver.residual_tolerance
                        ),
                        max_evaluations=int(
                            block_solver.pressure_flow_solver.max_evaluations
                        ),
                        sparse_pattern_trusted=True,
                        sparse_fallback_reason=None,
                    )
                )
        if not secondary_scope_indices:
            for pressure_scope_solver, names in zip(
                closure_plan.secondary_pressure_solvers,
                closure_plan.secondary_component_groups,
                strict=True,
            ):
                selected = frozenset(names)
                unknown_indices = tuple(
                    index
                    for index, unknown in enumerate(pressure_solver.unknowns)
                    if unknown.component in selected
                )
                equation_indices = tuple(
                    index
                    for index, equation in enumerate(pressure_solver.equation_templates)
                    if (
                        equation.owner == "component" and equation.owner_id in selected
                    )
                    or (
                        equation.owner == "connection"
                        and connection_specs[connection_index[equation.owner_id]].endpoint_a_port_index
                        in {
                            port_index
                            for name in selected
                            for port_index in port_indices_by_component[name]
                        }
                    )
                )
                scope_index = len(scopes)
                secondary_scope_indices.append(scope_index)
                scopes.append(
                    IRPressureFlowScope(
                        scope_id=f"pressure_flow.physical_island.{scope_index}",
                        kind=IRPressureFlowScopeKind.PHYSICAL_ISLAND,
                        component_indices=tuple(component_index[name] for name in names),
                        unknown_indices=unknown_indices,
                        equation_indices=equation_indices,
                        block_indices=(),
                        causal_plan_index=None,
                        residual_tolerance=float(pressure_scope_solver.residual_tolerance),
                        max_evaluations=int(pressure_scope_solver.max_evaluations),
                        sparse_pattern_trusted=bool(
                            pressure_scope_solver.jacobian_sparsity_is_trusted
                        ),
                        sparse_fallback_reason=(
                            pressure_scope_solver.jacobian_sparsity_fallback_reason
                        ),
                    )
                )
                if not pressure_scope_solver.jacobian_sparsity_is_trusted:
                    capability_issues.append(
                        IRCapabilityIssue(
                            code="IR_PRESSURE_FLOW_FINITE_DIFFERENCE_FALLBACK",
                            severity=IRDiagnosticSeverity.WARNING,
                            scope_id=f"pressure_flow.physical_island.{scope_index}",
                            message=(
                                "This secondary nonlinear scope cannot trust its "
                                "structural Jacobian and may use finite-difference "
                                "least_squares: "
                                f"{pressure_scope_solver.jacobian_sparsity_fallback_reason or 'unspecified reason'}."
                            ),
                        )
                    )

    secondary_pressure_stage = add_stage(
        "rhs.pressure_flow.stream_sensitive",
        IRStageKind.PRESSURE_FLOW,
        (
            IRKernelCallOperation(
                system_kernel_index,
                None,
                IRKernelPhase.RESIDUAL,
                pressure_flow_model_read_slots,
                (*all_unknown_slots, *residual_slots),
                tuple(range(len(pressure_equations))),
            ),
        )
        if secondary_scope_indices
        else (),
    )

    mechanical_acceleration_slots: dict[tuple[str, ...], IRSlotRef] = {}
    acceleration_operations: list[Any] = []
    for group in system.mechanical_state_reducer.groups:
        names = tuple(sorted(component.name for component in group.components))
        acceleration_slot = slots.add(
            IRBufferKind.WORK_FLOAT,
            "mechanical_group:" + ",".join(names) + ".acceleration",
            semantic="mechanical_acceleration",
            role="workspace",
            quantity="acceleration",
            unit="m/s2",
        )
        group_key = _mechanical_group_key(group)
        mechanical_acceleration_slots[group_key] = acceleration_slot
        reads = tuple(
            slot
            for name in names
            for slot in port_slots_by_component[name]
        )
        mode_slot = mode_slot_by_group.get(group_key)
        if mode_slot is not None:
            reads = (*reads, mode_slot)
        acceleration_operations.append(
            IRKernelCallOperation(
                system_kernel_index,
                None,
                IRKernelPhase.PRIMAL,
                reads,
                (acceleration_slot,),
            )
        )
    acceleration_stage = add_stage(
        "rhs.mechanical_acceleration",
        IRStageKind.MECHANICAL_ACCELERATION,
        acceleration_operations,
    )

    derivative_operations: list[Any] = []
    for entry in system.mechanical_state_reducer.state_entries:
        component = entry.representative if isinstance(entry, MechanicalConstraintGroup) else entry
        extra = (
            (mechanical_acceleration_slots[_mechanical_group_key(entry)],)
            if isinstance(entry, MechanicalConstraintGroup)
            else ()
        )
        derivative_operations.append(
            IRKernelCallOperation(
                component_kernel_index[component.name],
                component_index[component.name],
                IRKernelPhase.DERIVATIVE,
                (
                    *local_state_by_component[component.name],
                    *port_slots_by_component[component.name],
                    *connected_h_slots_by_component[component.name],
                    *parameter_slots_by_component[component.name],
                    *extra,
                ),
                local_derivative_by_component[component.name],
            )
        )
    derivative_operations.append(
        IRStateMapOperation(
            IRStateMapKind.DERIVATIVE_GATHER,
            tuple(local_derivative_slots),
            tuple(derivative_output_slots),
        )
    )
    derivative_stage = add_stage(
        "rhs.derivative_gather",
        IRStageKind.DERIVATIVE_REDUCE,
        derivative_operations,
    )

    # Blocks encode the nested current implementation: stream converges first,
    # then the outer stream/temperature/pressure loop monitors mass flows.
    blocks: list[IRExecutionBlock] = []
    stream_block_index = len(blocks)
    blocks.append(
        IRExecutionBlock(
            block_id="stream.reference_global_fixed_point",
            kind=IRBlockKind.STREAM_SCC,
            steps=(_step(stream_nondynamic_stage),),
            convergence=IRConvergenceSpec(
                monitor_slots=stream_node_slots,
                absolute_tolerance=0.0,
                relative_tolerance=float(system.stream_resolver.relative_tolerance),
                max_iterations=int(system.stream_resolver.max_iterations),
                relaxation=1.0,
                rollback_slots=stream_node_slots,
                failure_policy=IRFailurePolicy.RETRY_SMALLER_STEP,
            ),
        )
    )
    flow_slots = tuple(
        unknown.slot for unknown in algebraic_unknowns if unknown.variable == "m_flow"
    )
    transaction_plan_object = system._thermofluid_transaction_plan
    transaction_snapshot_slots = _unique_slots(
        port_slot_by_variable_id[
            f"{binding.component_name}.{binding.port_name}.{binding.variable}"
        ]
        for binding in transaction_plan_object.port_value_bindings
    )
    thermofluid_block_index = len(blocks)
    blocks.append(
        IRExecutionBlock(
            block_id="thermofluid.reference_fixed_point",
            kind=IRBlockKind.FIXED_POINT,
            steps=(
                _step(stream_block_index, IRStepKind.BLOCK),
                _step(stream_dynamic_stage),
                _step(temperature_reference_stage),
                _step(secondary_pressure_stage),
            ),
            convergence=IRConvergenceSpec(
                monitor_slots=flow_slots,
                absolute_tolerance=0.0,
                relative_tolerance=1.0e-12,
                max_iterations=25,
                relaxation=1.0,
                rollback_slots=transaction_snapshot_slots,
                failure_policy=IRFailurePolicy.RETRY_SMALLER_STEP,
            ),
        )
    )
    stream_plans = (
        IRStreamPlan(
            plan_id="stream.reference_global",
            node_slots=stream_node_slots,
            strongly_connected_components=(
                IRStreamSCC(
                    "stream.reference_global.scc",
                    stream_node_slots,
                    stream_block_index,
                ),
            )
            if stream_node_slots
            else (),
            condensed_edges=(),
            topological_scc_indices=(0,) if stream_node_slots else (),
        ),
    )

    sensitive_names = tuple(
        dict.fromkeys(
            name
            for block_solver in closure_plan.secondary_block_solvers
            for name in block_solver.sensitive_components
        )
    )
    thermofluid = IRThermofluidPlan(
        physical_port_indices=tuple(
            index
            for index, port in enumerate(port_specs)
            if port.kind is IRPortKind.PHYSICAL
        ),
        global_component_indices=tuple(
            component_index[name] for name in closure_plan.global_component_group
        ),
        stream_plan_index=0,
        secondary_pressure_scope_indices=tuple(secondary_scope_indices),
        sensitive_component_indices=tuple(
            component_index[name] for name in sensitive_names
        ),
        maximum_iterations=25,
        flow_relative_tolerance=1.0e-12,
        uses_conservative_global_solver=bool(
            closure_plan.uses_conservative_global_solver
        ),
        conservative_fallback_reason=closure_plan.conservative_fallback_reason,
    )

    pressure_flow = IRPressureFlowPlan(
        unknowns=tuple(algebraic_unknowns),
        equations=tuple(pressure_equations),
        blocks=tuple(pressure_blocks),
        scopes=tuple(scopes),
        global_scope_index=0,
        secondary_scope_indices=tuple(secondary_scope_indices),
        pressure_lower_bound=0.0,
    )

    rhs_steps = (
        _step(state_stage),
        _step(signal_stage),
        _step(mechanical_equivalence_stage),
        _step(volume_stage),
        _step(property_stage),
        _step(global_pressure_stage),
        _step(thermofluid_block_index, IRStepKind.BLOCK),
        _step(acceleration_stage),
        _step(derivative_stage),
    )

    # Event roots and reset plans. Signal event entries are event-source
    # families because concrete breakpoint times belong to the run plan.
    events: list[IREventSpec] = []
    event_operations: list[Any] = []
    pending_mechanical_events: list[
        tuple[str, MechanicalConstraintGroup, str, IRSlotRef]
    ] = []
    for group in event_groups:
        representative = group.representative
        position_slot = local_state_by_component[representative.name][1]
        owners = tuple(component_index[item.name] for item in group.components)
        for side, bound, direction in (
            ("lower", group.lower_bound, IREventDirection.DECREASING),
            ("upper", group.upper_bound, IREventDirection.INCREASING),
        ):
            if bound is None:
                continue
            event_id = (
                "mechanical_group:"
                + ",".join(sorted(item.name for item in group.components))
                + f".{side}_impact"
            )
            root_slot = slots.add(
                IRBufferKind.EVENT_OUTPUT,
                event_id,
                semantic="event_root",
                role="event",
                quantity="length",
                unit="m",
            )
            event_operations.append(
                IRLinearCombinationOperation((position_slot,), (1.0,), root_slot, -float(bound))
            )
            pending_mechanical_events.append((event_id, group, side, root_slot))
    for component_name, _event_times in system.signal_resolver._event_sources:
        event_id = f"{component_name}.signal_breakpoint_source"
        root_slot = slots.add(
            IRBufferKind.EVENT_OUTPUT,
            event_id,
            semantic="time_breakpoint_family",
            role="event",
            quantity="time",
            unit="s",
            owner_component_index=component_index[component_name],
        )
        event_operations.append(
            IRKernelCallOperation(
                component_kernel_index[component_name],
                component_index[component_name],
                IRKernelPhase.EVENT,
                (time_slot, *parameter_slots_by_component[component_name]),
                (root_slot,),
            )
        )
        events.append(
            IREventSpec(
                event_id=event_id,
                event_kind="known_time_breakpoint_source",
                owner_component_indices=(component_index[component_name],),
                root_slot=root_slot,
                direction=IREventDirection.ANY,
                terminal=False,
                priority=100,
                mode_guards=(),
                reset_steps=(),
                invalidated_caches=(
                    IRCacheKind.PRESSURE_FLOW,
                    IRCacheKind.STREAM,
                    IRCacheKind.JACOBIAN,
                    IRCacheKind.OUTPUT,
                ),
                restarts_integrator=True,
            )
        )
    event_stage = add_stage("events.evaluate", IRStageKind.EVENT, event_operations)
    for event_id, group, side, root_slot in pending_mechanical_events:
        group_key = _mechanical_group_key(group)
        state_offset = group_solver_offset[group_key]
        mode_slot = mode_slot_by_group[group_key]
        reset_stage = add_stage(
            f"reset.{event_id}",
            IRStageKind.RESET,
            (
                IRKernelCallOperation(
                    system_kernel_index,
                    None,
                    IRKernelPhase.RESET,
                    (
                        state_input_slots[state_offset],
                        state_input_slots[state_offset + 1],
                        mode_slot,
                    ),
                    (
                        state_input_slots[state_offset],
                        state_input_slots[state_offset + 1],
                        mode_slot,
                    ),
                ),
            ),
        )
        events.append(
            IREventSpec(
                event_id=event_id,
                event_kind=f"mechanical_{side}_impact",
                owner_component_indices=tuple(
                    component_index[item.name] for item in group.components
                ),
                root_slot=root_slot,
                direction=(
                    IREventDirection.DECREASING
                    if side == "lower"
                    else IREventDirection.INCREASING
                ),
                terminal=False,
                priority=10,
                mode_guards=(),
                reset_steps=(_step(reset_stage),),
                invalidated_caches=(
                    IRCacheKind.PRESSURE_FLOW,
                    IRCacheKind.STREAM,
                    IRCacheKind.JACOBIAN,
                    IRCacheKind.OUTPUT,
                ),
                restarts_integrator=True,
            )
        )

    # Fixed CSR state Jacobian and deterministic seed-0 coloring.
    jacobian_pattern = _csr_pattern(system.jacobian_sparsity())
    jacobian_value_slots = tuple(
        slots.add(
            IRBufferKind.JACOBIAN_VALUE,
            f"jacobian.value.{index}",
            semantic="jacobian_value",
            role="jacobian",
        )
        for index in range(jacobian_pattern.nonzero_count)
    )
    try:
        from scipy.optimize._numdiff import group_columns

        groups = group_columns(system.jacobian_sparsity(), order=0)
        color_count = int(groups.max(initial=-1)) + 1
        color_groups = tuple(
            tuple(int(index) for index, color in enumerate(groups) if int(color) == group)
            for group in range(color_count)
        )
    except (ImportError, AttributeError, TypeError):
        color_groups = tuple((index,) for index in range(len(initial_state)))
        capability_issues.append(
            IRCapabilityIssue(
                code="IR_JACOBIAN_COLORING_UNAVAILABLE",
                severity=IRDiagnosticSeverity.WARNING,
                scope_id="jacobian",
                message="Deterministic SciPy column coloring was unavailable.",
            )
        )
    positions_by_column: list[list[int]] = [
        [] for _ in range(jacobian_pattern.column_count)
    ]
    for row in range(jacobian_pattern.row_count):
        for position in range(
            jacobian_pattern.row_pointers[row],
            jacobian_pattern.row_pointers[row + 1],
        ):
            positions_by_column[jacobian_pattern.column_indices[position]].append(position)
    exact_rows = system._exact_ode_jacobian_rows()
    analytic_positions: list[int] = []
    for row, columns in exact_rows.items():
        for position in range(
            jacobian_pattern.row_pointers[row],
            jacobian_pattern.row_pointers[row + 1],
        ):
            if jacobian_pattern.column_indices[position] in columns:
                analytic_positions.append(position)
    analytic_position_set = set(analytic_positions)
    finite_difference_columns = tuple(
        IRFiniteDifferenceColumn(
            column_index=column,
            value_indices=tuple(
                position for position in positions if position not in analytic_position_set
            ),
            relative_step=1.4901161193847656e-08,
        )
        for column, positions in enumerate(positions_by_column)
        if any(position not in analytic_position_set for position in positions)
    )
    jacobian_stage = add_stage(
        "jacobian.fill",
        IRStageKind.JACOBIAN,
        (
            IRKernelCallOperation(
                system_kernel_index,
                None,
                IRKernelPhase.JACOBIAN,
                (time_slot, *state_input_slots),
                jacobian_value_slots,
            ),
        ),
    )
    jacobian = IRJacobianPlan(
        pattern=jacobian_pattern,
        value_slots=jacobian_value_slots,
        color_groups=color_groups,
        fill_steps=(*rhs_steps, _step(jacobian_stage)),
        analytic_value_indices=tuple(sorted(set(analytic_positions))),
        local_finite_difference_columns=finite_difference_columns,
    )

    output_operations: list[Any] = []
    for component in components:
        derived_sources = tuple(
            output_source_slots[index]
            for index in output_indices_by_component[component.name]
            if result_metadata[index].scope == "component"
        )
        if derived_sources:
            output_operations.append(
                IRKernelCallOperation(
                    component_kernel_index[component.name],
                    component_index[component.name],
                    IRKernelPhase.OUTPUT,
                    (
                        *local_state_by_component.get(component.name, ()),
                        *port_slots_by_component[component.name],
                        *parameter_slots_by_component[component.name],
                    ),
                    derived_sources,
                )
            )
    output_operations.extend(
        IRCopyOperation(source, target)
        for source, target in zip(output_source_slots, output_slots, strict=True)
    )
    output_stage = add_stage("outputs.project", IRStageKind.OUTPUT, output_operations)
    output_specs = tuple(
        IROutputSpec(
            output_id=metadata.key,
            component_index=component_index[metadata.component_id],
            scope=metadata.scope,
            port_name=metadata.port_name,
            name=metadata.name,
            label=metadata.label,
            category=metadata.category,
            quantity=metadata.quantity,
            unit=metadata.unit,
            order=int(metadata.order),
            source_slot=output_source_slots[index],
            output_slot=output_slots[index],
        )
        for index, metadata in enumerate(result_metadata)
    )

    # C-01 can only describe the current object kernels as reference calls.  Use
    # conservative, explicit input supersets so dependency slicing and cache
    # invalidation never omit parameters, modes, medium constants, or port
    # values that those Python methods may reach indirectly.  C-02 will replace
    # these supersets with independently declared, fixed kernel signatures.
    medium_slots_by_component: dict[int, list[IRSlotRef]] = {
        index: [] for index in range(len(components))
    }
    for medium in medium_specs:
        for owner_index in medium.component_indices:
            medium_slots_by_component[owner_index].extend(medium.parameter_slots)
    component_reference_reads = tuple(
        _unique_slots(
            (
                *component.parameter_slots,
                *component.state_slots,
                *component.mode_slots,
                *component.port_slots,
                *component.workspace_float_slots,
                *component.workspace_int_slots,
                *medium_slots_by_component[owner_index],
            )
        )
        for owner_index, component in enumerate(component_instances)
    )
    system_reference_reads = _unique_slots(
        IRSlotRef(kind, index)
        for kind in (
            IRBufferKind.TIME,
            IRBufferKind.STATE_INPUT,
            IRBufferKind.LOCAL_STATE,
            IRBufferKind.ALGEBRAIC,
            IRBufferKind.SIGNAL,
            IRBufferKind.PARAMETER,
            IRBufferKind.CONSTANT,
            IRBufferKind.MODE,
            IRBufferKind.RUNTIME_INPUT,
        )
        for index in range(len(slots.initial[kind]))
    )
    completed_stages: list[IRStage] = []
    for stage in stages:
        completed_operations: list[Any] = []
        for operation in stage.operations:
            if isinstance(operation, IRKernelCallOperation):
                extra_reads = (
                    system_reference_reads
                    if operation.component_index is None
                    else component_reference_reads[operation.component_index]
                )
                operation = replace(
                    operation,
                    read_slots=_unique_slots((*operation.read_slots, *extra_reads)),
                )
            completed_operations.append(operation)
        completed_stages.append(
            _stage(stage.stage_id, stage.kind, completed_operations)
        )
    stages = completed_stages

    # The Python PortState object has an ``m_flow`` attribute even for
    # mechanical ports, and the reference transaction helper consequently
    # carries inert bindings for those ports.  They are always hidden zeroes:
    # mechanical active metadata exposes ``f``, not ``m_flow``.  The numeric IR
    # records only the physical thermofluid closure quantities that can change
    # (active pneumatic ``m_flow`` slots); inventing hidden slots or treating
    # force as mass flow would give the native contract different semantics.
    transaction_flow_slots = _unique_slots(
        slot
        for item in transaction_plan_object.flow_bindings
        for slot in (
            port_slot_by_variable_id.get(
                f"{item.component_name}.{item.port_name}.m_flow"
            ),
        )
        if slot is not None
    )
    transaction = IRTransactionPlan(
        snapshot_slots=transaction_snapshot_slots,
        flow_slots=transaction_flow_slots,
        cache_component_indices=tuple(
            component_index[item.component.name]
            for item in transaction_plan_object.component_cache_bindings
        ),
        cache_attribute_ids=tuple(
            f"{item.component.name}.{name}"
            for item in transaction_plan_object.component_cache_bindings
            for name in item.attribute_names
        ),
        diagnostic_owner_ids=tuple(
            f"{type(owner).__module__}.{type(owner).__qualname__}:{index}"
            for index, owner in enumerate(transaction_plan_object.diagnostic_owners)
        ),
        restores_on_recoverable_failure=True,
        restores_on_fatal_failure=True,
    )

    component_capabilities = tuple(
        IRComponentCapability(
            component_index=index,
            level=IRCapabilityLevel.REFERENCE_ONLY,
            supported_phases=_component_phases(component),
            missing_features=("c02_pure_numeric_kernel", "native_implementation_version"),
        )
        for index, component in enumerate(components)
    )
    capability_issues.insert(
        0,
        IRCapabilityIssue(
            code="IR_NATIVE_KERNELS_NOT_DECLARED",
            severity=IRDiagnosticSeverity.WARNING,
            scope_id=network.name,
            message=(
                "System structure and conservative reference-kernel dependencies "
                "are described, but Python kernel effects remain opaque until C-02 "
                "declares pure numeric call signatures."
            ),
        ),
    )
    capabilities = IRCapabilityReport(
        system_level=(
            IRCapabilityLevel.UNSUPPORTED
            if any(issue.severity is IRDiagnosticSeverity.ERROR for issue in capability_issues)
            else IRCapabilityLevel.REFERENCE_ONLY
        ),
        components=component_capabilities,
        issues=tuple(capability_issues),
    )

    entry_points = (
        IREntryPoint(
            IREntryPointKind.RHS,
            rhs_steps,
            (time_slot, *state_input_slots),
            tuple(derivative_output_slots),
        ),
        IREntryPoint(
            IREntryPointKind.EVENTS,
            (_step(state_stage), _step(event_stage)),
            (time_slot, *state_input_slots),
            tuple(event.root_slot for event in events),
        ),
        IREntryPoint(
            IREntryPointKind.JACOBIAN,
            (*rhs_steps, _step(jacobian_stage)),
            (time_slot, *state_input_slots),
            jacobian_value_slots,
        ),
        IREntryPoint(
            IREntryPointKind.OUTPUTS,
            (
                _step(state_stage),
                _step(signal_stage),
                _step(mechanical_equivalence_stage),
                _step(volume_stage),
                _step(property_stage),
                _step(global_pressure_stage),
                _step(thermofluid_block_index, IRStepKind.BLOCK),
                _step(acceleration_stage),
                _step(output_stage),
            ),
            (time_slot, *state_input_slots),
            tuple(output_slots),
        ),
    )

    required_features = (
        "callback_free",
        "independent_entry_points",
        "transactional_closure",
        "fixed_csr_jacobian",
        "reference_kernel_dispatch",
    )

    program = SystemIR(
        version=CURRENT_SYSTEM_IR_VERSION,
        model_id=network.name,
        model_version=str(model_version),
        compiler_id=SYSTEM_NUMERIC_IR_COMPILER_ID,
        compiler_version=SYSTEM_NUMERIC_IR_COMPILER_VERSION,
        numeric_dtype=IRDType.FLOAT64,
        buffers=slots.buffers(),
        values=tuple(slots.values),
        kernels=tuple(kernels),
        components=component_instances,
        mediums=tuple(medium_specs),
        ports=tuple(port_specs),
        connections=tuple(connection_specs),
        state_reducer=state_reducer,
        causal_plans=tuple(causal_plans),
        pressure_flow=pressure_flow,
        stream_plans=stream_plans,
        thermofluid=thermofluid,
        stages=tuple(stages),
        blocks=tuple(blocks),
        entry_points=entry_points,
        transaction=transaction,
        modes=tuple(modes),
        jacobian=jacobian,
        events=tuple(events),
        outputs=output_specs,
        capabilities=capabilities,
        required_features=required_features,
    )

    # Validation is mandatory at the compiler boundary.  Validation owns the
    # exception type and complete invariant list; callers never receive an
    # unchecked program.
    from .validation import require_valid_system_ir

    return require_valid_system_ir(program)


__all__ = ["compile_system_ir"]