完成求解器雅可比矩阵首轮优化，增加更新目录，整理了文档文件夹，增加了服务启动脚本

app/simulation/README.md

@@ -33,8 +33,8 @@ FastAPI 的 `GET /api/components/catalog` 会把注册表转换成前端组件
 公开临时库入口是 `components/amesim/library.py`。公开模型必须在
 模型类中声明 `MODEL_TYPE / MODEL_VERSION / PORTS / PARAMETERS /
 RESULT_VARIABLES / DISPLAY / create()`，再把类路径加入库清单。完整规范参见
-[`组件模型建模规范 v1`](../../docs/component-model-authoring-spec-v1.md)和
-[`组件库分类、发现与读取规范 v1`](../../docs/component-library-spec-v1.md)。
+[`组件模型建模规范 v1`](../../docs/standard/component-model-authoring-spec-v1.md)和
+[`组件库分类、发现与读取规范 v1`](../../docs/standard/component-library-spec-v1.md)。
 
 当前关键文件：
 
@@ -88,7 +88,7 @@ Jacobian 和 System XML XSD，完成后才开始接收请求。它不会运行
 `diagnostics.performance.propertyCache` 中返回。
 
 基准原始 JSON 默认放到已忽略的 `app/data/` 下。指标字段、实测结果和使用边界见
-[`仿真性能评估 2026-08-15`](../../docs/仿真性能评估-2026-08-15.md)。
+[`仿真性能评估 2026-08-15`](../../docs/other/仿真性能评估-2026-08-15.md)。
 
 ## 当前阶段进度
 
@@ -311,7 +311,7 @@ print(result.used_modelica_reference)
 ## 基线结果
 
 当前基线对比摘要来自：
-[`testmodel_modelica_comparison_summary.txt`](../../tests/baselines/simulation/testmodel/testmodel_modelica_comparison_summary.txt)
+[`testmodel_modelica_comparison_summary.txt`](../../tests/data/testmodel/testmodel_modelica_comparison_summary.txt)
 
 当前四个主变量的最大误差为：
 

app/simulation/components/amesim/flow/pipes.py

@@ -1,8 +1,9 @@
 from __future__ import annotations
 
-from collections.abc import Mapping
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
 from functools import lru_cache
-from math import isclose, log, log10, pi, sqrt, tanh
+from math import isclose, isfinite, log, log10, pi, sqrt, tanh
 
 from app.simulation.components.amesim.gases import (
     AMESIM_GAS_INDEX_PARAMETER,
@@ -22,11 +23,35 @@ from app.simulation.core.metadata import (
     ResultVariableDefinition,
     THERMODYNAMIC_VOLUME_RESULT_VARIABLES,
 )
-from app.simulation.core.medium import GasMedium, ThermodynamicProperties
+from app.simulation.core.medium import (
+    GasMedium,
+    ThermodynamicProperties,
+    ThermodynamicPropertiesLinearization,
+)
 from app.simulation.core.ports import PortDefinition
 from app.simulation.core.state import VolumeState
 
 
+@dataclass(frozen=True)
+class Pnl0001MassFlowLinearization:
+    value: float
+    partial_p_1: float
+    partial_p_2: float
+    partial_temperature: float
+    valid: bool = True
+    reason: str | None = None
+    direction: str = "forward"
+
+
+@dataclass(frozen=True)
+class Pnl0001DerivativeLinearization:
+    derivative: tuple[float, float]
+    tangents: tuple[tuple[float, ...], tuple[float, ...]]
+    properties: ThermodynamicPropertiesLinearization
+    valid: bool = True
+    reason: str | None = None
+
+
 _MAX_REPORTED_FRICTION_FACTOR = 64_000_000.0
 
 
@@ -829,6 +854,105 @@ class AmesimPnl0001(ThermodynamicVolumeComponent):
         )
         return magnitude if pressure_difference > 0.0 else -magnitude
 
+    def linearize_mass_flow(
+        self,
+        p_1: float,
+        p_2: float,
+        temperature: float,
+        *,
+        relative_step: float = 2.0 ** -26,
+        slope_relative_tolerance: float = 5.0e-3,
+    ) -> Pnl0001MassFlowLinearization:
+        """Audit local flow-law slopes without perturbing the full system RHS."""
+
+        p_1 = float(p_1)
+        p_2 = float(p_2)
+        temperature = float(temperature)
+        direction = "forward" if p_1 > p_2 else "reverse"
+        value = self.mass_flow(p_1, p_2, temperature)
+
+        def invalid(reason: str) -> Pnl0001MassFlowLinearization:
+            return Pnl0001MassFlowLinearization(
+                value=value,
+                partial_p_1=0.0,
+                partial_p_2=0.0,
+                partial_temperature=0.0,
+                valid=False,
+                reason=reason,
+                direction=direction,
+            )
+
+        if not all(isfinite(item) for item in (p_1, p_2, temperature, value)):
+            return invalid("non_finite_primal")
+        pressure_gap = abs(p_1 - p_2)
+        if pressure_gap <= 1.0e-8:
+            return invalid("flow_direction_boundary")
+        if temperature <= 1.0 * (1.0 + 1.0e-10):
+            return invalid("temperature_floor_boundary")
+        if relative_step <= 0.0 or slope_relative_tolerance <= 0.0:
+            raise ValueError("PNL0001 slope audit tolerances must be positive.")
+
+        pressure_step = min(
+            relative_step * max(abs(p_1), abs(p_2), 1.0),
+            0.25 * pressure_gap,
+        )
+        temperature_step = min(
+            relative_step * max(abs(temperature), 1.0),
+            0.25 * (temperature - 1.0),
+        )
+        if pressure_step <= 0.0 or temperature_step <= 0.0:
+            return invalid("unresolved_local_step")
+
+        arguments = (p_1, p_2, temperature)
+        argument_names = ("p_1", "p_2", "temperature")
+        steps = (pressure_step, pressure_step, temperature_step)
+        partials: list[float] = []
+        for argument_index, (argument, step) in enumerate(
+            zip(arguments, steps, strict=True)
+        ):
+            lower = list(arguments)
+            upper = list(arguments)
+            lower[argument_index] = argument - step
+            upper[argument_index] = argument + step
+            lower_value = self.mass_flow(*lower)
+            upper_value = self.mass_flow(*upper)
+            left_slope = (value - lower_value) / step
+            right_slope = (upper_value - value) / step
+            slope_scale = max(
+                abs(left_slope),
+                abs(right_slope),
+                abs(value) / max(abs(argument), 1.0),
+                1.0e-12,
+            )
+            if not all(
+                isfinite(item)
+                for item in (
+                    lower_value,
+                    upper_value,
+                    left_slope,
+                    right_slope,
+                )
+            ):
+                return invalid(
+                    f"non_finite_local_slope:{argument_names[argument_index]}"
+                )
+            if (
+                abs(left_slope - right_slope)
+                > slope_relative_tolerance * slope_scale
+            ):
+                return invalid(
+                    f"local_slope_disagreement:{argument_names[argument_index]}"
+                )
+            partials.append(0.5 * (left_slope + right_slope))
+
+        return Pnl0001MassFlowLinearization(
+            value=value,
+            partial_p_1=partials[0],
+            partial_p_2=partials[1],
+            partial_temperature=partials[2],
+            direction=direction,
+        )
+
     def component_result_values(self) -> Mapping[str, float]:
         props = self.properties()
         flow = self.mass_flow(self.port_1.p, props.p, props.T)
@@ -913,6 +1037,98 @@ class AmesimPnl0001(ThermodynamicVolumeComponent):
         )
         return derivative.as_vector()
 
+    def linearize_state_derivative(
+        self,
+        connected_h: Mapping[str, float],
+        *,
+        state_mass_tangent: Sequence[float],
+        state_energy_tangent: Sequence[float],
+        port_mass_flow_tangents: Mapping[str, Sequence[float]],
+        connected_h_tangents: Mapping[str, Sequence[float]],
+        property_linearization: ThermodynamicPropertiesLinearization | None = None,
+        flow_boundary_tolerance: float = 1.0e-12,
+    ) -> Pnl0001DerivativeLinearization:
+        """Linearize the pipe storage balance in a fixed stream mode."""
+
+        port_names = ("port_1", "port_2")
+        vectors = {
+            "state_mass": tuple(float(value) for value in state_mass_tangent),
+            "state_energy": tuple(float(value) for value in state_energy_tangent),
+        }
+        for port_name in port_names:
+            vectors[f"flow:{port_name}"] = tuple(
+                float(value) for value in port_mass_flow_tangents[port_name]
+            )
+            vectors[f"enthalpy:{port_name}"] = tuple(
+                float(value) for value in connected_h_tangents[port_name]
+            )
+        widths = {len(values) for values in vectors.values()}
+        if len(widths) != 1:
+            raise ValueError("PNL0001 tangent vectors must have equal lengths.")
+        width = len(vectors["state_mass"])
+        invalid_reason: str | None = None
+        if not all(isfinite(value) for values in vectors.values() for value in values):
+            invalid_reason = "non_finite_tangent_input"
+
+        properties = property_linearization or self.medium.linearize_properties_from_mU(
+            self.state.m,
+            self.state.U,
+            self.volume,
+            vectors["state_mass"],
+            vectors["state_energy"],
+            (0.0,) * width,
+        )
+        if properties.tangents.width != width:
+            raise ValueError(
+                "PNL0001 property tangent width must match balance tangents."
+            )
+        props = properties.properties
+        if not properties.valid:
+            invalid_reason = invalid_reason or properties.reason
+
+        mass_derivative = self.port_1.m_flow + self.port_2.m_flow
+        energy_derivative = self.thermal_energy_flow_w(props.T)
+        mass_tangent = [0.0] * width
+        thermal_coefficient = (
+            0.0 if self.mode == 1 else self.kth * self.exchange_area
+        )
+        energy_tangent = [
+            -thermal_coefficient * properties.tangents.T[index]
+            for index in range(width)
+        ]
+
+        for port_name in port_names:
+            port = self.get_port(port_name)
+            flow_tangent = vectors[f"flow:{port_name}"]
+            if (
+                abs(port.m_flow) <= flow_boundary_tolerance
+                and any(value != 0.0 for value in flow_tangent)
+            ):
+                invalid_reason = invalid_reason or (
+                    f"flow_direction_boundary:{port_name}"
+                )
+            if port.m_flow > 0.0:
+                inlet_h = connected_h[port_name]
+                inlet_h_tangent = vectors[f"enthalpy:{port_name}"]
+            else:
+                inlet_h = props.h
+                inlet_h_tangent = properties.tangents.h
+            energy_derivative += port.m_flow * inlet_h
+            for index in range(width):
+                mass_tangent[index] += flow_tangent[index]
+                energy_tangent[index] += (
+                    inlet_h * flow_tangent[index]
+                    + port.m_flow * inlet_h_tangent[index]
+                )
+
+        return Pnl0001DerivativeLinearization(
+            derivative=(mass_derivative, energy_derivative),
+            tangents=(tuple(mass_tangent), tuple(energy_tangent)),
+            properties=properties,
+            valid=invalid_reason is None,
+            reason=invalid_reason,
+        )
+
 
 class AmesimPnl0002(AmesimPnl0001):
     """AMESim PNL0002 R-C-R pneumatic pipe with one center compliance."""

app/simulation/components/amesim/mechanical/pistons.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
 
-from collections.abc import Mapping
-from math import pi
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+from math import isfinite, pi
 
 from app.simulation.components.amesim.gases import (
     AMESIM_GAS_INDEX_PARAMETER,
@@ -18,6 +19,18 @@ from app.simulation.core.ports import PortDefinition
 AMESIM_REFERENCE_PRESSURE_PA = 101300.0
 
 
+@dataclass(frozen=True)
+class Pnrp17Linearization:
+    volume: float
+    volume_flow: float
+    pressure_force: float
+    volume_tangent: tuple[float, ...]
+    volume_flow_tangent: tuple[float, ...]
+    pressure_force_tangent: tuple[float, ...]
+    valid: bool = True
+    reason: str | None = None
+
+
 class AmesimPnrp17(AlgebraicComponent):
     """AMESim PNRP17 pneumatic piston with two mechanical faces.
 
@@ -230,6 +243,53 @@ class AmesimPnrp17(AlgebraicComponent):
     def pneumatic_volume_outputs(self) -> Mapping[str, tuple[float, float]]:
         return {"port_1": (self.chamber_volume, self.chamber_volume_flow)}
 
+    def linearize_geometry_and_force(
+        self,
+        port_4_x_tangent: Sequence[float],
+        port_5_x_tangent: Sequence[float],
+        port_4_v_tangent: Sequence[float],
+        port_5_v_tangent: Sequence[float],
+        port_1_pressure_tangent: Sequence[float],
+    ) -> Pnrp17Linearization:
+        """Return exact piston geometry and pressure-force tangents."""
+
+        vectors = tuple(
+            tuple(float(value) for value in values)
+            for values in (
+                port_4_x_tangent,
+                port_5_x_tangent,
+                port_4_v_tangent,
+                port_5_v_tangent,
+                port_1_pressure_tangent,
+            )
+        )
+        widths = {len(values) for values in vectors}
+        if len(widths) != 1:
+            raise ValueError("PNRP17 tangent vectors must have equal lengths.")
+        valid = all(isfinite(value) for values in vectors for value in values)
+        area = self.effective_area
+        volume_tangent = tuple(
+            area * (right - left)
+            for left, right in zip(vectors[0], vectors[1], strict=True)
+        )
+        volume_flow_tangent = tuple(
+            area * (right - left)
+            for left, right in zip(vectors[2], vectors[3], strict=True)
+        )
+        pressure_force_tangent = tuple(
+            area * value for value in vectors[4]
+        )
+        return Pnrp17Linearization(
+            volume=self.chamber_volume,
+            volume_flow=self.chamber_volume_flow,
+            pressure_force=self.pressure_force,
+            volume_tangent=volume_tangent,
+            volume_flow_tangent=volume_flow_tangent,
+            pressure_force_tangent=pressure_force_tangent,
+            valid=valid,
+            reason=None if valid else "non_finite_tangent_input",
+        )
+
     def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None:
         self.port_1.h_outflow = connected_h.get(
             "port_1",

app/simulation/components/amesim/mechanical/translational.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
 
-from collections.abc import Mapping
-from math import expm1
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+from math import expm1, isfinite
 
 from app.simulation.core.base import AlgebraicComponent, DynamicComponent
 from app.simulation.core.catalog import (
@@ -20,6 +21,24 @@ from app.simulation.core.medium import IdealGasMedium
 from app.simulation.core.ports import PortDefinition
 
 
+@dataclass(frozen=True)
+class Mecmas21DerivativeLinearization:
+    derivative: tuple[float, float]
+    tangents: tuple[tuple[float, ...], tuple[float, ...]]
+    mode: str
+    valid: bool = True
+    reason: str | None = None
+
+
+@dataclass(frozen=True)
+class LstpContactForceLinearization:
+    force: float
+    force_tangent: tuple[float, ...]
+    mode: str
+    valid: bool = True
+    reason: str | None = None
+
+
 _MECMAS21_FRICTION_ENABLED = ParameterCondition("useFriction", (2.0,))
 _MECMAS21_NON_RESTITUTION = ParameterCondition("stoptype", (1.0, 2.0, 4.0))
 _MECMAS21_LIMITS_ENABLED = ParameterCondition("stoptype", (1.0, 2.0, 3.0))
@@ -723,6 +742,204 @@ class AmesimMecmas21(DynamicComponent):
         )
         return [self.acceleration(), velocity]
 
+    def linearize_state_derivative(
+        self,
+        port_1_force_tangent: Sequence[float],
+        port_2_force_tangent: Sequence[float],
+        velocity_tangent: Sequence[float],
+        position_tangent: Sequence[float],
+        *,
+        constraint_mode: str = "current",
+        boundary_tolerance: float = 1.0e-12,
+    ) -> Mecmas21DerivativeLinearization:
+        """Linearize one inertia in a declared fixed mechanical mode."""
+
+        vectors = tuple(
+            tuple(float(value) for value in values)
+            for values in (
+                port_1_force_tangent,
+                port_2_force_tangent,
+                velocity_tangent,
+                position_tangent,
+            )
+        )
+        widths = {len(values) for values in vectors}
+        if len(widths) != 1:
+            raise ValueError("MECMAS21 tangent vectors must have equal lengths.")
+        width = len(vectors[0])
+        invalid_reason: str | None = None
+        if not all(isfinite(value) for values in vectors for value in values):
+            invalid_reason = "non_finite_tangent_input"
+
+        requested_mode = constraint_mode
+        if requested_mode == "current":
+            fixed = (
+                self._constraint_acceleration == 0.0
+                and self._constraint_velocity == 0.0
+            )
+            mode = "fixed" if fixed else "free"
+            if self._constraint_acceleration is not None and not fixed:
+                invalid_reason = invalid_reason or (
+                    "group_acceleration_requires_aggregate"
+                )
+        elif requested_mode == "free":
+            mode = "free"
+        elif requested_mode in {"lower", "upper"}:
+            mode = requested_mode
+            fixed = (
+                self._constraint_acceleration == 0.0
+                and self._constraint_velocity == 0.0
+            )
+            if not fixed:
+                invalid_reason = invalid_reason or (
+                    "constraint_mode_not_statically_fixed"
+                )
+        elif requested_mode == "uninitialized":
+            mode = requested_mode
+            invalid_reason = invalid_reason or "constraint_mode_uninitialized"
+        else:
+            raise ValueError(
+                "MECMAS21 constraint_mode must be current, free, lower, upper, "
+                "or uninitialized."
+            )
+
+        if mode in {"fixed", "lower", "upper"}:
+            return Mecmas21DerivativeLinearization(
+                derivative=(self.acceleration(), 0.0),
+                tangents=((0.0,) * width, (0.0,) * width),
+                mode=mode,
+                valid=invalid_reason is None,
+                reason=invalid_reason,
+            )
+
+        force_1_tangent, force_2_tangent, dv, dx = vectors
+        acceleration_tangent = [
+            force_1_tangent[index] + force_2_tangent[index]
+            for index in range(width)
+        ]
+        if self.use_friction:
+            for index in range(width):
+                acceleration_tangent[index] += (
+                    -self.rvisc * dv[index]
+                    - 2.0 * self.wind * abs(self.v) * dv[index]
+                )
+            if (
+                self.fcoul != 0.0
+                and abs(self.v) <= boundary_tolerance
+                and any(value != 0.0 for value in dv)
+            ):
+                invalid_reason = invalid_reason or "dry_friction_direction_boundary"
+
+        def add_limit_tangent(
+            *,
+            side: str,
+            stiffness: float,
+            damping: float,
+            damping_penetration: float,
+            bound: float,
+            damping_sign: float,
+            force_sign: float,
+        ) -> None:
+            nonlocal invalid_reason
+            if int(self.stoptype) != 2:
+                return
+            penetration = (
+                bound - self.x if side == "lower" else self.x - bound
+            )
+            penetration_tangent = tuple(
+                (-value if side == "lower" else value) for value in dx
+            )
+            scale = max(abs(bound), abs(self.x), 1.0)
+            if penetration <= 0.0:
+                if (
+                    abs(penetration) <= boundary_tolerance * scale
+                    and any(value != 0.0 for value in penetration_tangent)
+                ):
+                    invalid_reason = invalid_reason or (
+                        f"soft_endstop_mode_boundary:{side}"
+                    )
+                return
+            if damping_penetration > 0.0:
+                fraction = min(penetration / damping_penetration, 1.0)
+                if penetration < damping_penetration:
+                    fraction_tangent = tuple(
+                        value / damping_penetration
+                        for value in penetration_tangent
+                    )
+                else:
+                    fraction_tangent = (0.0,) * width
+                    if (
+                        abs(penetration - damping_penetration)
+                        <= boundary_tolerance
+                        * max(abs(damping_penetration), 1.0)
+                        and any(value != 0.0 for value in penetration_tangent)
+                    ):
+                        invalid_reason = invalid_reason or (
+                            f"soft_endstop_damping_boundary:{side}"
+                        )
+            else:
+                fraction = 1.0
+                fraction_tangent = (0.0,) * width
+
+            raw_force = (
+                stiffness * penetration
+                + damping_sign * fraction * damping * self.v
+            )
+            raw_tangent = tuple(
+                stiffness * penetration_tangent[index]
+                + damping_sign
+                * damping
+                * (
+                    fraction * dv[index]
+                    + self.v * fraction_tangent[index]
+                )
+                for index in range(width)
+            )
+            if int(self.discContactOption) != 1 and raw_force <= 0.0:
+                if (
+                    abs(raw_force)
+                    <= boundary_tolerance
+                    * max(abs(stiffness * penetration), 1.0)
+                    and any(value != 0.0 for value in raw_tangent)
+                ):
+                    invalid_reason = invalid_reason or (
+                        f"soft_endstop_force_boundary:{side}"
+                    )
+                return
+            for index in range(width):
+                acceleration_tangent[index] += (
+                    force_sign * raw_tangent[index]
+                )
+
+        add_limit_tangent(
+            side="lower",
+            stiffness=self.Kbmin,
+            damping=self.Dbmin,
+            damping_penetration=self.Pdmin,
+            bound=self.xmin,
+            damping_sign=-1.0,
+            force_sign=1.0,
+        )
+        add_limit_tangent(
+            side="upper",
+            stiffness=self.Kbmax,
+            damping=self.Dbmax,
+            damping_penetration=self.Pdmax,
+            bound=self.xmax,
+            damping_sign=1.0,
+            force_sign=-1.0,
+        )
+        acceleration_tangent = tuple(
+            value / self.mass for value in acceleration_tangent
+        )
+        return Mecmas21DerivativeLinearization(
+            derivative=(self.unconstrained_acceleration(), self.v),
+            tangents=(acceleration_tangent, tuple(dv)),
+            mode=mode,
+            valid=invalid_reason is None,
+            reason=invalid_reason,
+        )
+
     def component_result_values(self) -> Mapping[str, float]:
         return {
             "a": self.acceleration(),
@@ -979,6 +1196,112 @@ class AmesimLstp00a(AlgebraicComponent):
         )
         return force if int(self.discContactOption) == 1 else max(force, 0.0)
 
+    def linearize_contact_force(
+        self,
+        port_1_x_tangent: Sequence[float],
+        port_2_x_tangent: Sequence[float],
+        port_1_velocity_tangent: Sequence[float],
+        port_2_velocity_tangent: Sequence[float],
+        *,
+        boundary_tolerance: float = 1.0e-12,
+    ) -> LstpContactForceLinearization:
+        """Linearize the elastic contact in its current unilateral mode."""
+
+        vectors = tuple(
+            tuple(float(value) for value in values)
+            for values in (
+                port_1_x_tangent,
+                port_2_x_tangent,
+                port_1_velocity_tangent,
+                port_2_velocity_tangent,
+            )
+        )
+        widths = {len(values) for values in vectors}
+        if len(widths) != 1:
+            raise ValueError("LSTP00A tangent vectors must have equal lengths.")
+        width = len(vectors[0])
+        if not all(isfinite(value) for values in vectors for value in values):
+            return LstpContactForceLinearization(
+                force=self.contact_force,
+                force_tangent=(0.0,) * width,
+                mode="invalid",
+                valid=False,
+                reason="non_finite_tangent_input",
+            )
+
+        dx_1, dx_2, dv_1, dv_2 = vectors
+        penetration_tangent = tuple(
+            left - right for left, right in zip(dx_1, dx_2, strict=True)
+        )
+        velocity_tangent = tuple(
+            left - right for left, right in zip(dv_1, dv_2, strict=True)
+        )
+        overlap = -self.gap
+        force = self.contact_force
+        scale = max(abs(self.gap0), abs(self.port_1.x), abs(self.port_2.x), 1.0)
+        if overlap <= 0.0:
+            on_boundary = abs(overlap) <= boundary_tolerance * scale
+            crossing = any(value != 0.0 for value in penetration_tangent)
+            return LstpContactForceLinearization(
+                force=force,
+                force_tangent=(0.0,) * width,
+                mode="boundary" if on_boundary else "inactive",
+                valid=not (on_boundary and crossing),
+                reason=(
+                    "contact_mode_boundary"
+                    if on_boundary and crossing
+                    else None
+                ),
+            )
+
+        penetration = overlap
+        if self.Pdis > 0.0:
+            damping_fraction = -expm1(-penetration / self.Pdis)
+            damping_fraction_tangent = tuple(
+                (1.0 - damping_fraction) * value / self.Pdis
+                for value in penetration_tangent
+            )
+        else:
+            damping_fraction = 1.0
+            damping_fraction_tangent = (0.0,) * width
+        relative_velocity = self.penetration_velocity
+        raw_force = (
+            self.kcont * penetration
+            + damping_fraction * self.rcont * relative_velocity
+        )
+        raw_tangent = tuple(
+            self.kcont * penetration_tangent[index]
+            + self.rcont
+            * (
+                damping_fraction * velocity_tangent[index]
+                + relative_velocity * damping_fraction_tangent[index]
+            )
+            for index in range(width)
+        )
+        if int(self.discContactOption) != 1 and raw_force <= 0.0:
+            on_boundary = (
+                abs(raw_force)
+                <= boundary_tolerance
+                * max(abs(self.kcont * penetration), 1.0)
+            )
+            crossing = any(value != 0.0 for value in raw_tangent)
+            return LstpContactForceLinearization(
+                force=force,
+                force_tangent=(0.0,) * width,
+                mode="force_boundary" if on_boundary else "clamped",
+                valid=not (on_boundary and crossing),
+                reason=(
+                    "contact_force_boundary"
+                    if on_boundary and crossing
+                    else None
+                ),
+            )
+        return LstpContactForceLinearization(
+            force=force,
+            force_tangent=raw_tangent,
+            mode="active",
+        )
+
     def clear_causal_contact(self) -> None:
         self._causal_penetration = None
         self._causal_contact_force = None

app/simulation/components/amesim/media/mediums.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
 
-from collections.abc import Callable
+from collections.abc import Callable, Sequence
 from dataclasses import dataclass
+from math import isfinite
 from typing import ClassVar
 
 from app.simulation.core.errors import RecoverableTrialStateError
@@ -9,6 +10,8 @@ from app.simulation.core.medium import (
     GasMedium,
     IdealGasMedium,
     ThermodynamicProperties,
+    ThermodynamicPropertiesLinearization,
+    ThermodynamicPropertyTangents,
 )
 from app.simulation.core.peng_robinson import HELIUM_PR, PengRobinsonFluid
 from app.simulation.performance import profile_property, record_property_iterations
@@ -309,6 +312,153 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium):
             ),
         )
 
+    def linearize_properties_from_mU(
+        self,
+        m: float,
+        U: float,
+        V: float,
+        dm: Sequence[float],
+        dU: Sequence[float],
+        dV: Sequence[float],
+        *,
+        properties: ThermodynamicProperties | None = None,
+    ) -> ThermodynamicPropertiesLinearization:
+        """Implicitly differentiate the Peng-Robinson m/U/V recovery."""
+
+        dm_values = tuple(float(value) for value in dm)
+        dU_values = tuple(float(value) for value in dU)
+        dV_values = tuple(float(value) for value in dV)
+        if not (len(dm_values) == len(dU_values) == len(dV_values)):
+            raise ValueError("Thermodynamic tangent vectors must have equal lengths.")
+        props = properties or self.properties_from_mU(m, U, V)
+        width = len(dm_values)
+
+        def invalid(reason: str) -> ThermodynamicPropertiesLinearization:
+            return ThermodynamicPropertiesLinearization(
+                properties=props,
+                tangents=ThermodynamicPropertyTangents.zeros(width),
+                valid=False,
+                reason=reason,
+            )
+
+        expected_density = m / V
+        expected_internal_energy = U / m
+        if (
+            abs(props.rho - expected_density)
+            > 1.0e-12 * max(abs(expected_density), 1.0)
+            or abs(props.u - expected_internal_energy)
+            > 1.0e-12 * max(abs(expected_internal_energy), 1.0)
+        ):
+            return invalid("properties_primal_mismatch")
+        if not all(
+            isfinite(value)
+            for values in (dm_values, dU_values, dV_values)
+            for value in values
+        ):
+            return invalid("non_finite_tangent_input")
+        if props.T <= 2.2 * (1.0 + 1.0e-10):
+            return invalid("temperature_floor_boundary")
+
+        pressure_temperature_derivative = (
+            self.fluid.pressure_temperature_derivative_at_density(
+                props.T,
+                props.rho,
+            )
+        )
+        pressure_density_derivative = (
+            self.fluid.pressure_density_derivative_at_temperature(
+                props.T,
+                props.rho,
+            )
+        )
+        cv = (
+            self.cv_at_temperature(props.T)
+            + self.fluid.residual_isochoric_heat_capacity_at_density(
+                props.T,
+                props.rho,
+            )
+        )
+        recovered_internal_energy = (
+            self.specific_internal_energy(props.T)
+            + self.fluid.residual_specific_internal_energy_at_density(
+                props.T,
+                props.rho,
+            )
+        )
+        recovery_scale = max(
+            abs(props.u),
+            abs(cv * props.T) if isfinite(cv) else 0.0,
+            1.0,
+        )
+        if (
+            not all(
+                isfinite(value)
+                for value in (
+                    pressure_temperature_derivative,
+                    pressure_density_derivative,
+                    cv,
+                    recovered_internal_energy,
+                )
+            )
+            or cv <= 0.0
+        ):
+            return invalid("invalid_peng_robinson_derivative")
+        if abs(recovered_internal_energy - props.u) > 1.0e-8 * recovery_scale:
+            return invalid("properties_recovery_not_converged")
+
+        internal_energy_density_derivative = (
+            props.p - props.T * pressure_temperature_derivative
+        ) / (props.rho * props.rho)
+        drho: list[float] = []
+        du: list[float] = []
+        dT: list[float] = []
+        dp: list[float] = []
+        dh: list[float] = []
+        for mass_tangent, energy_tangent, volume_tangent in zip(
+            dm_values,
+            dU_values,
+            dV_values,
+            strict=True,
+        ):
+            density_tangent = (
+                mass_tangent / V - m * volume_tangent / (V * V)
+            )
+            internal_energy_tangent = (
+                energy_tangent / m - U * mass_tangent / (m * m)
+            )
+            temperature_tangent = (
+                internal_energy_tangent
+                - internal_energy_density_derivative * density_tangent
+            ) / cv
+            pressure_tangent = (
+                pressure_temperature_derivative * temperature_tangent
+                + pressure_density_derivative * density_tangent
+            )
+            enthalpy_tangent = (
+                internal_energy_tangent
+                + pressure_tangent / props.rho
+                - props.p * density_tangent / (props.rho * props.rho)
+            )
+            drho.append(density_tangent)
+            du.append(internal_energy_tangent)
+            dT.append(temperature_tangent)
+            dp.append(pressure_tangent)
+            dh.append(enthalpy_tangent)
+
+        tangent_values = (*drho, *du, *dT, *dp, *dh)
+        if not all(isfinite(value) for value in tangent_values):
+            return invalid("non_finite_property_tangent")
+        return ThermodynamicPropertiesLinearization(
+            properties=props,
+            tangents=ThermodynamicPropertyTangents(
+                p=tuple(dp),
+                T=tuple(dT),
+                rho=tuple(drho),
+                u=tuple(du),
+                h=tuple(dh),
+            ),
+        )
+
 
 @dataclass(frozen=True)
 class AmesimGasPropertyModelSpec:

app/simulation/components/amesim/storage/chambers.py

@@ -1,6 +1,8 @@
 from __future__ import annotations
 
-from collections.abc import Mapping
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+from math import isfinite
 
 from app.simulation.components.amesim.gases import (
     AMESIM_GAS_INDEX_PARAMETER,
@@ -14,11 +16,24 @@ from app.simulation.core.metadata import (
     ResultVariableDefinition,
     THERMODYNAMIC_VOLUME_RESULT_VARIABLES,
 )
-from app.simulation.core.medium import GasMedium, ThermodynamicProperties
+from app.simulation.core.medium import (
+    GasMedium,
+    ThermodynamicProperties,
+    ThermodynamicPropertiesLinearization,
+)
 from app.simulation.core.ports import PortDefinition
 from app.simulation.core.state import VolumeState
 
 
+@dataclass(frozen=True)
+class Pnch012DerivativeLinearization:
+    derivative: tuple[float, float]
+    tangents: tuple[tuple[float, ...], tuple[float, ...]]
+    properties: ThermodynamicPropertiesLinearization
+    valid: bool = True
+    reason: str | None = None
+
+
 class AmesimPnch023(ThermodynamicVolumeComponent):
     """AMESim PNCH023 simple pneumatic chamber with heat exchange.
 
@@ -518,6 +533,132 @@ class AmesimPnch012(ThermodynamicVolumeComponent):
         energy_derivative -= props.p * self.total_volume_rate()
         return VolumeState(m=mass_derivative, U=energy_derivative).as_vector()
 
+    def linearize_state_derivative(
+        self,
+        connected_h: Mapping[str, float],
+        *,
+        state_mass_tangent: Sequence[float],
+        state_energy_tangent: Sequence[float],
+        external_volume_tangent: Sequence[float],
+        external_volume_rate_tangent: Sequence[float],
+        port_mass_flow_tangents: Mapping[str, Sequence[float]],
+        connected_h_tangents: Mapping[str, Sequence[float]],
+        property_linearization: ThermodynamicPropertiesLinearization | None = None,
+        flow_boundary_tolerance: float = 1.0e-12,
+    ) -> Pnch012DerivativeLinearization:
+        """Linearize the chamber balance while keeping stream modes fixed."""
+
+        port_names = ("port_1", "port_2", "port_3", "port_4")
+        vectors = {
+            "state_mass": tuple(float(value) for value in state_mass_tangent),
+            "state_energy": tuple(float(value) for value in state_energy_tangent),
+            "volume": tuple(float(value) for value in external_volume_tangent),
+            "volume_rate": tuple(
+                float(value) for value in external_volume_rate_tangent
+            ),
+        }
+        for port_name in port_names:
+            vectors[f"flow:{port_name}"] = tuple(
+                float(value) for value in port_mass_flow_tangents[port_name]
+            )
+            vectors[f"enthalpy:{port_name}"] = tuple(
+                float(value) for value in connected_h_tangents[port_name]
+            )
+        widths = {len(values) for values in vectors.values()}
+        if len(widths) != 1:
+            raise ValueError("PNCH012 tangent vectors must have equal lengths.")
+        width = len(vectors["state_mass"])
+        invalid_reason: str | None = None
+        if not all(isfinite(value) for values in vectors.values() for value in values):
+            invalid_reason = "non_finite_tangent_input"
+
+        raw_volume = (
+            self.cvol0
+            + sum(self.external_volumes.values())
+            + self.connected_external_volume()
+        )
+        minimum_volume = self.cvol0 / 100.0
+        volume_scale = max(abs(raw_volume), abs(minimum_volume), 1.0e-18)
+        on_volume_boundary = (
+            abs(raw_volume - minimum_volume) <= 1.0e-12 * volume_scale
+        )
+        supplied_volume_tangent = vectors["volume"]
+        if raw_volume < minimum_volume or on_volume_boundary:
+            used_volume_tangent = (0.0,) * width
+            used_volume_rate_tangent = (0.0,) * width
+            if on_volume_boundary and any(
+                value != 0.0
+                for value in (
+                    *supplied_volume_tangent,
+                    *vectors["volume_rate"],
+                )
+            ):
+                invalid_reason = invalid_reason or "volume_floor_boundary"
+        else:
+            used_volume_tangent = supplied_volume_tangent
+            used_volume_rate_tangent = vectors["volume_rate"]
+
+        properties = property_linearization or self.medium.linearize_properties_from_mU(
+            self.state.m,
+            self.state.U,
+            self.total_volume(),
+            vectors["state_mass"],
+            vectors["state_energy"],
+            used_volume_tangent,
+        )
+        if properties.tangents.width != width:
+            raise ValueError(
+                "PNCH012 property tangent width must match balance tangents."
+            )
+        props = properties.properties
+        if not properties.valid:
+            invalid_reason = invalid_reason or properties.reason
+
+        mass_derivative = sum(
+            self.get_port(port_name).m_flow for port_name in port_names
+        )
+        volume_rate = self.total_volume_rate()
+        energy_derivative = self.thermal_energy_flow_w(props.T) - props.p * volume_rate
+        mass_tangent = [0.0] * width
+        energy_tangent = [
+            -self.kth * self.sth * properties.tangents.T[index]
+            - volume_rate * properties.tangents.p[index]
+            - props.p * used_volume_rate_tangent[index]
+            for index in range(width)
+        ]
+
+        for port_name in port_names:
+            port = self.get_port(port_name)
+            flow_tangent = vectors[f"flow:{port_name}"]
+            if (
+                abs(port.m_flow) <= flow_boundary_tolerance
+                and any(value != 0.0 for value in flow_tangent)
+            ):
+                invalid_reason = invalid_reason or (
+                    f"flow_direction_boundary:{port_name}"
+                )
+            if port.m_flow > 0.0:
+                inlet_h = connected_h[port_name]
+                inlet_h_tangent = vectors[f"enthalpy:{port_name}"]
+            else:
+                inlet_h = props.h
+                inlet_h_tangent = properties.tangents.h
+            energy_derivative += port.m_flow * inlet_h
+            for index in range(width):
+                mass_tangent[index] += flow_tangent[index]
+                energy_tangent[index] += (
+                    inlet_h * flow_tangent[index]
+                    + port.m_flow * inlet_h_tangent[index]
+                )
+
+        return Pnch012DerivativeLinearization(
+            derivative=(mass_derivative, energy_derivative),
+            tangents=(tuple(mass_tangent), tuple(energy_tangent)),
+            properties=properties,
+            valid=invalid_reason is None,
+            reason=invalid_reason,
+        )
+
     def pressure_flow_equation_values(self) -> tuple[float, ...]:
         pressure = self.medium.properties_from_mU(
             self.state.m,

app/simulation/components/example.md

@@ -1,7 +1,7 @@
 # 元件建模规范与示例
 
 规范的权威版本位于
-[`docs/component-model-authoring-spec-v1.md`](../../../docs/component-model-authoring-spec-v1.md)。
+[`docs/standard/component-model-authoring-spec-v1.md`](../../../docs/standard/component-model-authoring-spec-v1.md)。
 本文档保留在组件目录中，作为离模型源码最近的完整示例；若两者不一致，应在同一次
 修改中同步，不能让示例形成另一套规则。
 
@@ -279,4 +279,4 @@ models=(
 10. 是否补充参数边界、端口契约、目录输出、结果元数据和最小仿真的自动测试。
 
 组件库、分类和自动发现的完整规则参见
-[`组件库分类、发现与读取规范 v1`](../../../docs/component-library-spec-v1.md)。
+[`组件库分类、发现与读取规范 v1`](../../../docs/standard/component-library-spec-v1.md)。

app/simulation/core/medium.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
-from typing import Protocol
+from math import isfinite
+from typing import Protocol, Sequence
 
 from app.simulation.core.errors import RecoverableTrialStateError
 from app.simulation.performance import profile_property
@@ -16,6 +17,36 @@ class ThermodynamicProperties:
     h: float
 
 
+@dataclass(frozen=True)
+class ThermodynamicPropertyTangents:
+    """Directional derivatives of a recovered thermodynamic state."""
+
+    p: tuple[float, ...]
+    T: tuple[float, ...]
+    rho: tuple[float, ...]
+    u: tuple[float, ...]
+    h: tuple[float, ...]
+
+    @property
+    def width(self) -> int:
+        return len(self.p)
+
+    @classmethod
+    def zeros(cls, width: int) -> "ThermodynamicPropertyTangents":
+        values = (0.0,) * width
+        return cls(p=values, T=values, rho=values, u=values, h=values)
+
+
+@dataclass(frozen=True)
+class ThermodynamicPropertiesLinearization:
+    """Primal properties and a validity-checked directional linearization."""
+
+    properties: ThermodynamicProperties
+    tangents: ThermodynamicPropertyTangents
+    valid: bool = True
+    reason: str | None = None
+
+
 class GasMedium(Protocol):
     """Thermodynamic contract required by pneumatic components.
 
@@ -75,6 +106,18 @@ class GasMedium(Protocol):
         V: float,
     ) -> ThermodynamicProperties: ...
 
+    def linearize_properties_from_mU(
+        self,
+        m: float,
+        U: float,
+        V: float,
+        dm: Sequence[float],
+        dU: Sequence[float],
+        dV: Sequence[float],
+        *,
+        properties: ThermodynamicProperties | None = None,
+    ) -> ThermodynamicPropertiesLinearization: ...
+
 
 @dataclass(frozen=True)
 class IdealGasMedium:
@@ -224,3 +267,100 @@ class IdealGasMedium:
         u = U / m
         h = self.specific_enthalpy(T)
         return ThermodynamicProperties(p=p, T=T, rho=rho, u=u, h=h)
+
+    def linearize_properties_from_mU(
+        self,
+        m: float,
+        U: float,
+        V: float,
+        dm: Sequence[float],
+        dU: Sequence[float],
+        dV: Sequence[float],
+        *,
+        properties: ThermodynamicProperties | None = None,
+    ) -> ThermodynamicPropertiesLinearization:
+        """Linearize properties_from_mU for several seed directions."""
+
+        dm_values = tuple(float(value) for value in dm)
+        dU_values = tuple(float(value) for value in dU)
+        dV_values = tuple(float(value) for value in dV)
+        if not (len(dm_values) == len(dU_values) == len(dV_values)):
+            raise ValueError("Thermodynamic tangent vectors must have equal lengths.")
+        props = properties or self.properties_from_mU(m, U, V)
+        width = len(dm_values)
+        expected_density = m / V
+        expected_internal_energy = U / m
+        if (
+            abs(props.rho - expected_density)
+            > 1.0e-12 * max(abs(expected_density), 1.0)
+            or abs(props.u - expected_internal_energy)
+            > 1.0e-12 * max(abs(expected_internal_energy), 1.0)
+        ):
+            return ThermodynamicPropertiesLinearization(
+                properties=props,
+                tangents=ThermodynamicPropertyTangents.zeros(width),
+                valid=False,
+                reason="properties_primal_mismatch",
+            )
+        if not all(
+            isfinite(value)
+            for values in (dm_values, dU_values, dV_values)
+            for value in values
+        ):
+            return ThermodynamicPropertiesLinearization(
+                properties=props,
+                tangents=ThermodynamicPropertyTangents.zeros(width),
+                valid=False,
+                reason="non_finite_tangent_input",
+            )
+
+        cv = self.cv_at_temperature(props.T)
+        cp = self.cp_at_temperature(props.T)
+        if not isfinite(cv) or not isfinite(cp) or cv <= 0.0 or cp <= 0.0:
+            return ThermodynamicPropertiesLinearization(
+                properties=props,
+                tangents=ThermodynamicPropertyTangents.zeros(width),
+                valid=False,
+                reason="non_positive_heat_capacity",
+            )
+
+        drho: list[float] = []
+        du: list[float] = []
+        dT: list[float] = []
+        dp: list[float] = []
+        dh: list[float] = []
+        for mass_tangent, energy_tangent, volume_tangent in zip(
+            dm_values,
+            dU_values,
+            dV_values,
+            strict=True,
+        ):
+            density_tangent = mass_tangent / V - m * volume_tangent / (V * V)
+            internal_energy_tangent = (
+                energy_tangent / m - U * mass_tangent / (m * m)
+            )
+            temperature_tangent = internal_energy_tangent / cv
+            pressure_tangent = self.R_gas * (
+                props.T * density_tangent + props.rho * temperature_tangent
+            )
+            enthalpy_tangent = cp * temperature_tangent
+            drho.append(density_tangent)
+            du.append(internal_energy_tangent)
+            dT.append(temperature_tangent)
+            dp.append(pressure_tangent)
+            dh.append(enthalpy_tangent)
+
+        tangent_values = (*drho, *du, *dT, *dp, *dh)
+        valid = all(isfinite(value) for value in tangent_values)
+        return ThermodynamicPropertiesLinearization(
+            properties=props,
+            tangents=ThermodynamicPropertyTangents(
+                p=tuple(dp),
+                T=tuple(dT),
+                rho=tuple(drho),
+                u=tuple(du),
+                h=tuple(dh),
+            ),
+            valid=valid,
+            reason=None if valid else "non_finite_property_tangent",
+        )

app/simulation/solvers/solver.py

@@ -12,6 +12,7 @@ CancellationCheck = Callable[[], bool]
 AcceptedStepCallback = Callable[[float], None]
 IntegrationStatus = Literal["completed", "cancelled", "failed"]
 DenseState = Callable[[float], list[float]]
+JacobianCallable = Callable[[float, object], object]
 
 
 @dataclass(frozen=True)
@@ -30,7 +31,9 @@ StateTransitionHandler = Callable[
 _MAX_STATE_TRANSITIONS_AT_SAME_TIME = 64
 
 
-class _IntegrationCancelled(Exception):
+class IntegrationCancelled(Exception):
+    """Internal control-flow signal shared by RHS and Jacobian evaluation."""
+
     pass
 
 
@@ -59,9 +62,19 @@ class SolverSegmentDiagnostics:
     solver_start_count: int = 0
     state_transition_count: int = 0
     recoverable_retry_count: int = 0
+    jacobian_evaluation_count: int = 0
+    jacobian_full_build_count: int = 0
+    jacobian_secant_reuse_count: int = 0
+    jacobian_audit_failure_count: int = 0
+    finite_difference_rhs_evaluation_count: int = 0
+    jacobian_base_rhs_evaluation_count: int = 0
+    jacobian_jv_audit_rhs_evaluation_count: int = 0
+    exact_column_build_count: int = 0
+    exact_column_fallback_count: int = 0
+    jacobian_assembly_seconds: float = 0.0
 
     def as_dict(self) -> dict[str, float | int]:
-        return {
+        result: dict[str, float | int] = {
             "startTime": self.start_time,
             "requestedStopTime": self.requested_stop_time,
             "simulatedUntil": self.simulated_until,
@@ -73,6 +86,61 @@ class SolverSegmentDiagnostics:
             "stateTransitionCount": self.state_transition_count,
             "recoverableRetryCount": self.recoverable_retry_count,
         }
+        if (
+            self.jacobian_evaluation_count
+            or self.finite_difference_rhs_evaluation_count
+            or self.jacobian_assembly_seconds
+        ):
+            result.update(
+                {
+                    "jacobianEvaluationCount": self.jacobian_evaluation_count,
+                    "jacobianFullBuildCount": self.jacobian_full_build_count,
+                    "jacobianSecantReuseCount": self.jacobian_secant_reuse_count,
+                    "jacobianAuditFailureCount": self.jacobian_audit_failure_count,
+                    "finiteDifferenceRhsEvaluationCount": (
+                        self.finite_difference_rhs_evaluation_count
+                    ),
+                    "jacobianBaseRhsEvaluationCount": (
+                        self.jacobian_base_rhs_evaluation_count
+                    ),
+                    "jacobianJvAuditRhsEvaluationCount": (
+                        self.jacobian_jv_audit_rhs_evaluation_count
+                    ),
+                    "exactColumnBuildCount": self.exact_column_build_count,
+                    "exactColumnFallbackCount": (
+                        self.exact_column_fallback_count
+                    ),
+                    "jacobianAssemblySeconds": self.jacobian_assembly_seconds,
+                }
+            )
+        return result
+
+
+_JACOBIAN_DIAGNOSTIC_KEYS = (
+    "jacobianEvaluationCount",
+    "fullBuildCount",
+    "secantReuseCount",
+    "auditFailureCount",
+    "finiteDifferenceRhsEvaluationCount",
+    "baseRhsEvaluationCount",
+    "jvAuditEvaluationCount",
+    "exactColumnBuildCount",
+    "exactColumnFallbackCount",
+    "assemblySeconds",
+)
+
+
+def _jacobian_diagnostic_snapshot(
+    jac: JacobianCallable | None,
+) -> dict[str, float]:
+    diagnostics = getattr(jac, "diagnostics", None)
+    if diagnostics is None:
+        return {key: 0.0 for key in _JACOBIAN_DIAGNOSTIC_KEYS}
+    values = diagnostics()
+    return {
+        key: float(values.get(key, 0.0))
+        for key in _JACOBIAN_DIAGNOSTIC_KEYS
+    }
 
 
 @dataclass(frozen=True)
@@ -309,7 +377,7 @@ def _runge_kutta_4(
         for target_time in t_eval[1:]:
             while current_time < target_time:
                 if cancel_check is not None and cancel_check():
-                    raise _IntegrationCancelled
+                    raise IntegrationCancelled
                 dt = min(config.max_step, target_time - current_time)
                 k1 = rhs(current_time, state)
                 k2 = rhs(current_time + 0.5 * dt, _vector_add(state, k1, 0.5 * dt))
@@ -374,7 +442,7 @@ def _runge_kutta_4(
                 report_step(current_time)
 
             _append_solution_sample(times, states, target_time, state)
-    except _IntegrationCancelled:
+    except IntegrationCancelled:
         status = "cancelled"
         message = "Simulation was stopped before reaching the requested end time."
         _append_solution_sample(times, states, current_time, state)
@@ -451,7 +519,7 @@ def _runge_kutta_4_segmented(
         nonlocal current_time, last_transition, same_time_transition_count, state
         while current_time < target_time:
             if cancel_check is not None and cancel_check():
-                raise _IntegrationCancelled
+                raise IntegrationCancelled
             dt = min(config.max_step, target_time - current_time)
             k1 = rhs(current_time, state)
             k2 = rhs(
@@ -565,7 +633,7 @@ def _runge_kutta_4_segmented(
                     sample_time = float(sample_times[sample_index])
                     _append_solution_sample(times, states, sample_time, state)
                     sample_index += 1
-    except _IntegrationCancelled:
+    except IntegrationCancelled:
         status = "cancelled"
         message = "Simulation was stopped before reaching the requested end time."
         _append_solution_sample(times, states, current_time, state)
@@ -595,6 +663,7 @@ def _integrate_scipy_stepwise(
     breakpoints: Sequence[float] = (),
     state_transition_handler: StateTransitionHandler | None = None,
     jac_sparsity=None,
+    jac: JacobianCallable | None = None,
 ) -> ODESolution:
     """Initial stepwise integration path for breakpoints and state resets.
 
@@ -617,6 +686,7 @@ def _integrate_scipy_stepwise(
     solver_type = solver_types.get(config.method)
     if solver_type is None:
         raise ValueError(f"Unsupported integration method: {config.method}")
+    implicit_jac = jac if config.method in {"BDF", "Radau"} else None
 
     times = [float(config.t_start)]
     states = [[float(value)] for value in initial_state]
@@ -632,8 +702,13 @@ def _integrate_scipy_stepwise(
 
     def cancellable_rhs(time, state):
         if cancel_check():
-            raise _IntegrationCancelled
-        return rhs(float(time), [float(value) for value in state])
+            raise IntegrationCancelled
+        normalized_state = [float(value) for value in state]
+        derivative = rhs(float(time), normalized_state)
+        observer = getattr(implicit_jac, "observe", None)
+        if observer is not None:
+            observer(float(time), normalized_state, derivative)
+        return derivative
 
     status: IntegrationStatus = "completed"
     message = "The solver successfully reached the end of the integration interval."
@@ -683,6 +758,7 @@ def _integrate_scipy_stepwise(
         segment_solver_starts = 0
         segment_state_transitions = 0
         segment_recoverable_retries = 0
+        jacobian_work_start = _jacobian_diagnostic_snapshot(implicit_jac)
 
         while has_integration_interval and last_accepted_time < integration_end:
             if cancel_check():
@@ -695,8 +771,11 @@ def _integrate_scipy_stepwise(
                 "atol": config.atol,
                 "max_step": segment_max_step,
             }
-            if jac_sparsity is not None and config.method in {"BDF", "Radau"}:
-                solver_options["jac_sparsity"] = jac_sparsity
+            if config.method in {"BDF", "Radau"}:
+                if implicit_jac is not None:
+                    solver_options["jac"] = implicit_jac
+                elif jac_sparsity is not None:
+                    solver_options["jac_sparsity"] = jac_sparsity
             requested_first_step = (
                 0.1 * segment_max_step
                 if last_recoverable_error is not None
@@ -709,6 +788,9 @@ def _integrate_scipy_stepwise(
                 )
 
             try:
+                start_segment = getattr(implicit_jac, "start_segment", None)
+                if start_segment is not None:
+                    start_segment()
                 solver = solver_type(
                     cancellable_rhs,
                     last_accepted_time,
@@ -716,7 +798,7 @@ def _integrate_scipy_stepwise(
                     integration_end,
                     **solver_options,
                 )
-            except _IntegrationCancelled:
+            except IntegrationCancelled:
                 status = "cancelled"
                 message = cancellation_message()
                 break
@@ -755,7 +837,7 @@ def _integrate_scipy_stepwise(
                 step_start_state = list(last_accepted_state)
                 try:
                     step_message = solver.step()
-                except _IntegrationCancelled:
+                except IntegrationCancelled:
                     status = "cancelled"
                     message = (
                         "Simulation was stopped before reaching the requested end time."
@@ -957,6 +1039,11 @@ def _integrate_scipy_stepwise(
             if not restart_at_transition:
                 break
 
+        jacobian_work_end = _jacobian_diagnostic_snapshot(implicit_jac)
+        jacobian_work = {
+            key: jacobian_work_end[key] - jacobian_work_start[key]
+            for key in _JACOBIAN_DIAGNOSTIC_KEYS
+        }
         solver_segments.append(
             SolverSegmentDiagnostics(
                 start_time=float(segment_start_time),
@@ -971,6 +1058,34 @@ def _integrate_scipy_stepwise(
                 solver_start_count=segment_solver_starts,
                 state_transition_count=segment_state_transitions,
                 recoverable_retry_count=segment_recoverable_retries,
+                jacobian_evaluation_count=int(
+                    jacobian_work["jacobianEvaluationCount"]
+                ),
+                jacobian_full_build_count=int(
+                    jacobian_work["fullBuildCount"]
+                ),
+                jacobian_secant_reuse_count=int(
+                    jacobian_work["secantReuseCount"]
+                ),
+                jacobian_audit_failure_count=int(
+                    jacobian_work["auditFailureCount"]
+                ),
+                finite_difference_rhs_evaluation_count=int(
+                    jacobian_work["finiteDifferenceRhsEvaluationCount"]
+                ),
+                jacobian_base_rhs_evaluation_count=int(
+                    jacobian_work["baseRhsEvaluationCount"]
+                ),
+                jacobian_jv_audit_rhs_evaluation_count=int(
+                    jacobian_work["jvAuditEvaluationCount"]
+                ),
+                exact_column_build_count=int(
+                    jacobian_work["exactColumnBuildCount"]
+                ),
+                exact_column_fallback_count=int(
+                    jacobian_work["exactColumnFallbackCount"]
+                ),
+                jacobian_assembly_seconds=jacobian_work["assemblySeconds"],
             )
         )
         if status != "completed":
@@ -1034,6 +1149,7 @@ def integrate_ode(
     breakpoints: Sequence[float] | None = None,
     state_transition_handler: StateTransitionHandler | None = None,
     jac_sparsity=None,
+    jac: JacobianCallable | None = None,
 ):
     """Integrate an ODE, optionally restarting at equation discontinuities.
 
@@ -1104,10 +1220,26 @@ def integrate_ode(
             normalized_breakpoints,
             state_transition_handler,
             jac_sparsity,
+            jac,
         )
 
+    implicit_jac = jac if config.method in {"BDF", "Radau"} else None
+    solve_rhs = rhs
+    if implicit_jac is not None:
+        observer = getattr(implicit_jac, "observe", None)
+        if observer is not None:
+            def observed_rhs(time, state):
+                derivative = rhs(time, state)
+                observer(float(time), state, derivative)
+                return derivative
+
+            solve_rhs = observed_rhs
+        start_segment = getattr(implicit_jac, "start_segment", None)
+        if start_segment is not None:
+            start_segment()
+
     solve_options = {
-        "fun": rhs,
+        "fun": solve_rhs,
         "t_span": (config.t_start, config.t_stop),
         "y0": initial_state,
         "method": config.method,
@@ -1118,6 +1250,8 @@ def integrate_ode(
     }
     if config.first_step is not None:
         solve_options["first_step"] = config.first_step
-    if jac_sparsity is not None and config.method in {"BDF", "Radau"}:
+    if implicit_jac is not None:
+        solve_options["jac"] = implicit_jac
+    elif jac_sparsity is not None and config.method in {"BDF", "Radau"}:
         solve_options["jac_sparsity"] = jac_sparsity
     return solve_ivp(**solve_options)

app/simulation/solvers/tangent.py

@@ -0,0 +1,850 @@
+"""Proof-gated tangent columns for the three-piston reference network.
+
+This module is deliberately narrower than the generic algebraic solver.  It
+only compiles a tangent provider after proving the state layout, component
+types, physical connections, and causal execution plan used by the committed
+three-piston XML.  A failed proof leaves the ordinary seed-0 numerical
+Jacobian in control; a runtime mode boundary requests the same one-build
+fallback through :class:`ExactColumnsUnavailable`.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+from math import isfinite
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+from app.simulation.solvers.jacobian import ExactColumnsUnavailable
+from app.simulation.solvers.mechanical import MechanicalConstraintGroup
+from app.simulation.systems.network import Endpoint
+
+if TYPE_CHECKING:
+    from app.simulation.systems.generic import GenericFluidSystem
+
+
+_TARGET_BRANCH_NAMES = (
+    (
+        "mass_friction_endstops_10",
+        "pn_brp2_8",
+        "pn_c1_8",
+        "pneumatic_69",
+        "elasticendstop_8",
+    ),
+    (
+        "mass_friction_endstops_11",
+        "pn_brp2_9",
+        "pn_c1_9",
+        "pneumatic_68",
+        "elasticendstop_9",
+    ),
+    (
+        "mass_friction_endstops_12",
+        "pn_brp2_10",
+        "pn_c1_10",
+        "pneumatic_66",
+        "elasticendstop_10",
+    ),
+)
+
+
+@dataclass(frozen=True)
+class ThreePistonBranch:
+    mass: object
+    piston: object
+    chamber: object
+    pipe: object
+    contact: object
+    velocity_index: int
+    position_index: int
+
+
+@dataclass(frozen=True)
+class ThreePistonTangentCompilation:
+    eligible: bool
+    reason: str | None
+    columns: tuple[int, ...] = ()
+    provider: "ThreePistonTangentProvider | None" = None
+    reached_assignment_count: int = 0
+
+    def diagnostics(self) -> dict[str, object]:
+        return {
+            "eligible": self.eligible,
+            "fallbackReason": self.reason,
+            "columns": list(self.columns),
+            "columnCount": len(self.columns),
+            "reachedAssignmentCount": self.reached_assignment_count,
+        }
+
+
+@dataclass(frozen=True)
+class _PrimalContext:
+    time: float
+    state: np.ndarray
+    connected_h: dict[str, dict[str, float]]
+
+
+def _failed(reason: str) -> ThreePistonTangentCompilation:
+    return ThreePistonTangentCompilation(False, reason)
+
+
+class ThreePistonTangentProvider:
+    """Batched six-direction provider compiled for one system instance."""
+
+    def __init__(
+        self,
+        system: "GenericFluidSystem",
+        branches: tuple[ThreePistonBranch, ...],
+        columns: tuple[int, ...],
+        state_offsets: Mapping[str, tuple[int, int]],
+        reached_assignment_count: int,
+    ) -> None:
+        self.system = system
+        self.branches = branches
+        self.columns = columns
+        self.state_offsets = dict(state_offsets)
+        self.reached_assignment_count = reached_assignment_count
+        self._context: _PrimalContext | None = None
+        self._capture_requested = False
+
+    def request_primal_capture(self) -> None:
+        """Capture exactly the next completed RHS primal closure."""
+
+        self._capture_requested = True
+
+    def cancel_primal_capture(self) -> None:
+        """Discard a pending capture after an interrupted base RHS."""
+
+        self._capture_requested = False
+
+    def record_primal(
+        self,
+        time: float,
+        state: Sequence[float],
+        connected_h: Mapping[str, Mapping[str, float]],
+    ) -> None:
+        if not self._capture_requested:
+            return
+        self._capture_requested = False
+        required_components = {
+            item.name
+            for branch in self.branches
+            for item in (branch.chamber, branch.pipe)
+        }
+        self._context = _PrimalContext(
+            time=float(time),
+            state=np.asarray(state, dtype=float).copy(),
+            connected_h={
+                name: {port: float(value) for port, value in values.items()}
+                for name, values in connected_h.items()
+                if name in required_components
+            },
+        )
+
+    def __call__(
+        self,
+        time: float,
+        state: np.ndarray,
+        normalized_indexes: tuple[int, ...],
+    ) -> np.ndarray:
+        if normalized_indexes != self.columns:
+            raise ValueError("Compiled tangent columns were requested out of order.")
+        context = self._context
+        values = np.asarray(state, dtype=float)
+        if (
+            context is None
+            or float(time) != context.time
+            or values.shape != context.state.shape
+            or not np.array_equal(values, context.state)
+        ):
+            raise ExactColumnsUnavailable("stalePrimalContext")
+        return self._evaluate(context)
+
+    @staticmethod
+    def _zeros(width: int) -> tuple[float, ...]:
+        return (0.0,) * width
+
+    @staticmethod
+    def _has_signal(vector: Sequence[float]) -> bool:
+        return any(float(value) != 0.0 for value in vector)
+
+    @staticmethod
+    def _require_valid(value: object, prefix: str) -> None:
+        if not bool(getattr(value, "valid", False)):
+            reason = getattr(value, "reason", None) or "invalid"
+            raise ExactColumnsUnavailable(f"{prefix}:{reason}")
+
+    def _evaluate(self, context: _PrimalContext) -> np.ndarray:
+        # The primal RHS may disable the causal fast path after a residual
+        # audit.  Recheck after that closure and before replaying its compiled
+        # assignments so a dynamic downgrade uses the ordinary full FD build.
+        if not self.system.pressure_flow_solver.causal_fast_path_enabled:
+            raise ExactColumnsUnavailable("causalFastPathDisabled")
+        width = len(self.columns)
+        zero = self._zeros(width)
+        out = np.zeros((len(context.state), width), dtype=float)
+        seeds = {
+            column: tuple(float(index == seed_index) for index in range(width))
+            for seed_index, column in enumerate(self.columns)
+        }
+        solver = self.system.pressure_flow_solver
+        tangents: dict[str, tuple[float, ...]] = {}
+
+        def tangent(key: str) -> tuple[float, ...]:
+            return tangents.get(key, zero)
+
+        def negative_sum(
+            vectors: Sequence[Sequence[float]],
+        ) -> tuple[float, ...]:
+            return tuple(
+                -sum(float(vector[index]) for vector in vectors)
+                for index in range(width)
+            )
+
+        for variable in ("x", "v"):
+            plan = solver._causal_effort_plan_by_variable.get(variable)
+            if plan is None:
+                raise ExactColumnsUnavailable("causalEffortPlanUnavailable")
+            for assignment in plan:
+                matched = [
+                    branch
+                    for branch in self.branches
+                    if any(
+                        unknown.component == branch.mass.name
+                        for unknown in assignment.members
+                    )
+                ]
+                if len(matched) > 1:
+                    raise ExactColumnsUnavailable("coupledTargetMechanicalSeeds")
+                vector = zero
+                if matched:
+                    branch = matched[0]
+                    column = (
+                        branch.position_index
+                        if variable == "x"
+                        else branch.velocity_index
+                    )
+                    vector = seeds[column]
+                for member in assignment.members:
+                    tangents[member.id] = vector
+
+        geometry: dict[str, object] = {}
+        chamber_properties: dict[str, object] = {}
+        for branch in self.branches:
+            piston = branch.piston
+            linearization = piston.linearize_geometry_and_force(
+                tangent(f"{piston.name}.port_4.x"),
+                tangent(f"{piston.name}.port_5.x"),
+                tangent(f"{piston.name}.port_4.v"),
+                tangent(f"{piston.name}.port_5.v"),
+                zero,
+            )
+            self._require_valid(linearization, "pistonGeometry")
+            chamber = branch.chamber
+            volume = float(chamber.total_volume())
+            if volume <= float(chamber.cvol0) / 100.0:
+                raise ExactColumnsUnavailable("volumeFloor")
+            primal = chamber.medium.properties_from_mU(
+                chamber.state.m,
+                chamber.state.U,
+                volume,
+            )
+            properties = chamber.medium.linearize_properties_from_mU(
+                chamber.state.m,
+                chamber.state.U,
+                volume,
+                zero,
+                zero,
+                linearization.volume_tangent,
+                properties=primal,
+            )
+            self._require_valid(properties, "chamberProperties")
+            geometry[piston.name] = linearization
+            chamber_properties[chamber.name] = properties
+
+        pressure_plan = solver._causal_effort_plan_by_variable.get("p")
+        if pressure_plan is None:
+            raise ExactColumnsUnavailable("causalPressurePlanUnavailable")
+        for assignment in pressure_plan:
+            matched = [
+                branch
+                for branch in self.branches
+                if any(
+                    unknown.component == branch.chamber.name
+                    for unknown in assignment.members
+                )
+            ]
+            if len(matched) > 1:
+                raise ExactColumnsUnavailable("coupledTargetPressureSeeds")
+            vector = (
+                chamber_properties[matched[0].chamber.name].tangents.p
+                if matched
+                else zero
+            )
+            for member in assignment.members:
+                tangents[member.id] = tuple(vector)
+
+        contacts: dict[str, object] = {}
+        for branch in self.branches:
+            piston = branch.piston
+            linearization = piston.linearize_geometry_and_force(
+                tangent(f"{piston.name}.port_4.x"),
+                tangent(f"{piston.name}.port_5.x"),
+                tangent(f"{piston.name}.port_4.v"),
+                tangent(f"{piston.name}.port_5.v"),
+                tangent(f"{piston.name}.port_1.p"),
+            )
+            self._require_valid(linearization, "pistonPressureForce")
+            geometry[piston.name] = linearization
+            contact = branch.contact
+            contact_linearization = contact.linearize_contact_force(
+                tangent(f"{contact.name}.port_1.x"),
+                tangent(f"{contact.name}.port_2.x"),
+                tangent(f"{contact.name}.port_1.v"),
+                tangent(f"{contact.name}.port_2.v"),
+            )
+            self._require_valid(contact_linearization, "contactMode")
+            contacts[contact.name] = contact_linearization
+
+        equations = {
+            equation.id: equation for equation in solver.equation_templates
+        }
+        branch_component = {
+            item.name: branch
+            for branch in self.branches
+            for item in (
+                branch.mass,
+                branch.piston,
+                branch.chamber,
+                branch.pipe,
+                branch.contact,
+            )
+        }
+        pipe_flows: dict[str, object] = {}
+        for stage in solver._explicit_flow_plan:
+            pending: list[tuple[str, tuple[float, ...]]] = []
+            for assignment in stage.assignments:
+                equation = equations.get(assignment.equation_id)
+                if equation is None:
+                    raise ExactColumnsUnavailable("causalEquationMissing")
+                dependencies = [
+                    tangent(variable)
+                    for variable in equation.variables
+                    if variable != assignment.unknown.id
+                ]
+                if not any(
+                    self._has_signal(vector) for vector in dependencies
+                ):
+                    pending.append((assignment.unknown.id, zero))
+                    continue
+                if equation.relation == "sumToZero":
+                    vectors = [
+                        tangent(variable)
+                        for variable in equation.variables
+                        if variable != assignment.unknown.id
+                        and variable in solver._unknowns_by_id
+                        and solver._unknowns_by_id[variable].role == "flow"
+                    ]
+                    pending.append(
+                        (assignment.unknown.id, negative_sum(vectors))
+                    )
+                    continue
+
+                component = assignment.component
+                model = getattr(component, "MODEL_TYPE", None)
+                if model == "amesim_pnl0001":
+                    branch = branch_component.get(component.name)
+                    if branch is None or component is not branch.pipe:
+                        raise ExactColumnsUnavailable("unexpectedPipeReach")
+                    flow_linearization = pipe_flows.get(component.name)
+                    if flow_linearization is None:
+                        properties = component.medium.properties_from_mU(
+                            component.state.m,
+                            component.state.U,
+                            component.volume,
+                        )
+                        flow_linearization = component.linearize_mass_flow(
+                            component.port_1.p,
+                            component.port_2.p,
+                            properties.T,
+                        )
+                        self._require_valid(
+                            flow_linearization,
+                            "pipeMassFlow",
+                        )
+                        pipe_flows[component.name] = flow_linearization
+                    vector = tuple(
+                        flow_linearization.partial_p_1 * first
+                        + flow_linearization.partial_p_2 * second
+                        for first, second in zip(
+                            tangent(f"{component.name}.port_1.p"),
+                            tangent(f"{component.name}.port_2.p"),
+                            strict=True,
+                        )
+                    )
+                elif model == "amesim_lstp00a":
+                    contact_linearization = contacts.get(component.name)
+                    if contact_linearization is None:
+                        raise ExactColumnsUnavailable(
+                            f"unexpectedContactReach:{component.name}"
+                        )
+                    sign = (
+                        1.0
+                        if assignment.unknown.port == "port_1"
+                        else -1.0
+                    )
+                    vector = tuple(
+                        sign * value
+                        for value in contact_linearization.force_tangent
+                    )
+                elif model == "amesim_pnrp17":
+                    piston_geometry = geometry.get(component.name)
+                    if piston_geometry is None:
+                        raise ExactColumnsUnavailable(
+                            f"unexpectedPistonReach:{component.name}"
+                        )
+                    other = next(
+                        (
+                            tangent(variable)
+                            for variable in equation.variables
+                            if variable != assignment.unknown.id
+                            and variable.endswith(".f")
+                        ),
+                        zero,
+                    )
+                    sign = (
+                        -1.0
+                        if equation.id.endswith(
+                            "piston_side_force_balance"
+                        )
+                        else 1.0
+                    )
+                    vector = tuple(
+                        -force + sign * pressure
+                        for force, pressure in zip(
+                            other,
+                            piston_geometry.pressure_force_tangent,
+                            strict=True,
+                        )
+                    )
+                else:
+                    raise ExactColumnsUnavailable(
+                        f"unsupportedReach:{model or 'connection'}"
+                    )
+                pending.append((assignment.unknown.id, tuple(vector)))
+            for key, vector in pending:
+                tangents[key] = vector
+
+        outflow: dict[Endpoint, tuple[float, ...]] = {}
+        for branch in self.branches:
+            enthalpy_tangent = tuple(
+                chamber_properties[branch.chamber.name].tangents.h
+            )
+            for port_name in branch.chamber.ports:
+                outflow[Endpoint(branch.chamber.name, port_name)] = (
+                    enthalpy_tangent
+                )
+            # PNRP17 mirrors the connected chamber enthalpy on its sole
+            # pneumatic port at the already closed primal point.
+            outflow[Endpoint(branch.piston.name, "port_1")] = enthalpy_tangent
+
+        connected_tangent: dict[
+            str,
+            dict[str, tuple[float, ...]],
+        ] = {name: {} for name in self.system.network.components}
+        for connection in self.system.network.connections:
+            if connection.kind != "physical":
+                continue
+            first, second = connection.endpoints
+            connected_tangent[first.component][first.port] = outflow.get(
+                second,
+                zero,
+            )
+            connected_tangent[second.component][second.port] = outflow.get(
+                first,
+                zero,
+            )
+
+        for component_name, port_values in connected_tangent.items():
+            component = self.system.network.components[component_name]
+            if not getattr(component, "PRESSURE_FLOW_DEPENDS_ON_STREAM", False):
+                continue
+            if any(
+                self._has_signal(vector) for vector in port_values.values()
+            ):
+                raise ExactColumnsUnavailable(
+                    f"streamSensitiveReach:{component_name}"
+                )
+
+        for branch in self.branches:
+            chamber = branch.chamber
+            chamber_linearization = chamber.linearize_state_derivative(
+                context.connected_h[chamber.name],
+                state_mass_tangent=zero,
+                state_energy_tangent=zero,
+                external_volume_tangent=(
+                    geometry[branch.piston.name].volume_tangent
+                ),
+                external_volume_rate_tangent=(
+                    geometry[branch.piston.name].volume_flow_tangent
+                ),
+                port_mass_flow_tangents={
+                    name: tangent(f"{chamber.name}.{name}.m_flow")
+                    for name in chamber.ports
+                },
+                connected_h_tangents=connected_tangent[chamber.name],
+                property_linearization=chamber_properties[chamber.name],
+            )
+            self._require_valid(chamber_linearization, "chamberDerivative")
+            offset, size = self.state_offsets[chamber.name]
+            if size != 2:
+                raise ValueError("Target chamber state layout changed.")
+            out[offset, :], out[offset + 1, :] = (
+                chamber_linearization.tangents
+            )
+
+            pipe = branch.pipe
+            pipe_properties = pipe.medium.linearize_properties_from_mU(
+                pipe.state.m,
+                pipe.state.U,
+                pipe.volume,
+                zero,
+                zero,
+                zero,
+            )
+            self._require_valid(pipe_properties, "pipeProperties")
+            pipe_linearization = pipe.linearize_state_derivative(
+                context.connected_h[pipe.name],
+                state_mass_tangent=zero,
+                state_energy_tangent=zero,
+                port_mass_flow_tangents={
+                    name: tangent(f"{pipe.name}.{name}.m_flow")
+                    for name in pipe.ports
+                },
+                connected_h_tangents=connected_tangent[pipe.name],
+                property_linearization=pipe_properties,
+            )
+            self._require_valid(pipe_linearization, "pipeDerivative")
+            offset, size = self.state_offsets[pipe.name]
+            if size != 2:
+                raise ValueError("Target pipe state layout changed.")
+            out[offset, :], out[offset + 1, :] = pipe_linearization.tangents
+
+        target_pneumatic = {
+            item.name
+            for branch in self.branches
+            for item in (branch.chamber, branch.pipe)
+        }
+        for entry in self.system.mechanical_state_reducer.state_entries:
+            if (
+                isinstance(entry, MechanicalConstraintGroup)
+                or entry.name in target_pneumatic
+            ):
+                continue
+            for port_name, port in entry.ports.items():
+                definition = port.definition
+                if (
+                    definition is not None
+                    and definition.kind == "physical"
+                    and definition.domain == "pneumatic"
+                    and self._has_signal(
+                        tangent(f"{entry.name}.{port_name}.m_flow")
+                    )
+                ):
+                    raise ExactColumnsUnavailable(
+                        f"unsupportedDynamicReach:{entry.name}"
+                    )
+
+        mass_seeds = {
+            branch.mass.name: (
+                seeds[branch.velocity_index],
+                seeds[branch.position_index],
+            )
+            for branch in self.branches
+        }
+        for group in self.system.mechanical_state_reducer.groups:
+            if len(group.components) != 1:
+                raise ExactColumnsUnavailable("reachableRigidMassGroup")
+            mass = group.representative
+            force_1 = tangent(f"{mass.name}.port_1.f")
+            force_2 = tangent(f"{mass.name}.port_2.f")
+            velocity, position = mass_seeds.get(mass.name, (zero, zero))
+            if not any(
+                self._has_signal(vector)
+                for vector in (force_1, force_2, velocity, position)
+            ):
+                continue
+            fixed = (
+                mass._constraint_acceleration == 0.0
+                and mass._constraint_velocity == 0.0
+            )
+            if fixed and (
+                abs(group.total_unconstrained_force())
+                <= 1.0e-12
+                * max(
+                    abs(mass.port_1.f),
+                    abs(mass.port_2.f),
+                    1.0,
+                )
+            ):
+                raise ExactColumnsUnavailable("mechanicalReleaseBoundary")
+            mass_linearization = mass.linearize_state_derivative(
+                force_1,
+                force_2,
+                velocity,
+                position,
+                constraint_mode="current" if fixed else "free",
+            )
+            self._require_valid(
+                mass_linearization,
+                "mechanicalDerivative",
+            )
+            offset, size = self.state_offsets[mass.name]
+            if size != 2:
+                raise ValueError("Mechanical state layout changed.")
+            out[offset, :], out[offset + 1, :] = mass_linearization.tangents
+
+        if not np.all(np.isfinite(out)):
+            raise ExactColumnsUnavailable("nonFiniteTangentColumns")
+        return out
+
+
+def compile_three_piston_tangent_provider(
+    system: "GenericFluidSystem",
+) -> ThreePistonTangentCompilation:
+    """Compile the proof-gated target provider, or return a stable reason."""
+
+    solver = system.pressure_flow_solver
+    if not solver.causal_fast_path_eligible:
+        return _failed("causalFastPathIneligible")
+    if not solver.causal_fast_path_enabled:
+        return _failed("causalFastPathDisabled")
+    closure = system._thermofluid_closure_plan
+    if closure.uses_conservative_global_solver:
+        return _failed("conservativeThermofluidClosure")
+    if system.pneumatic_storage_reducer.groups:
+        return _failed("coupledPneumaticStorage")
+
+    state_offsets: dict[str, tuple[int, int]] = {}
+    cursor = 0
+    group_by_name: dict[str, MechanicalConstraintGroup] = {}
+    for entry in system.mechanical_state_reducer.state_entries:
+        if isinstance(entry, MechanicalConstraintGroup):
+            if len(entry.components) != 1:
+                cursor += 2
+                continue
+            component = entry.representative
+            state_offsets[component.name] = (cursor, 2)
+            group_by_name[component.name] = entry
+            cursor += 2
+        else:
+            state_offsets[entry.name] = (cursor, int(entry.state_size))
+            cursor += int(entry.state_size)
+
+    expected_types = (
+        "amesim_mecmas21",
+        "amesim_pnrp17",
+        "amesim_pnch012",
+        "amesim_pnl0001",
+        "amesim_lstp00a",
+    )
+    branches: list[ThreePistonBranch] = []
+    for names in _TARGET_BRANCH_NAMES:
+        try:
+            components = tuple(system.network.components[name] for name in names)
+        except KeyError:
+            return _failed("targetComponentMissing")
+        if tuple(getattr(item, "MODEL_TYPE", None) for item in components) != expected_types:
+            return _failed("targetComponentTypeMismatch")
+        mass, piston, chamber, pipe, contact = components
+        if mass.name not in group_by_name or mass.name not in state_offsets:
+            return _failed("targetMechanicalStateLayout")
+        offset, size = state_offsets[mass.name]
+        if size != 2:
+            return _failed("targetMechanicalStateLayout")
+        branches.append(
+            ThreePistonBranch(
+                mass=mass,
+                piston=piston,
+                chamber=chamber,
+                pipe=pipe,
+                contact=contact,
+                velocity_index=offset,
+                position_index=offset + 1,
+            )
+        )
+
+    required_pairs: set[frozenset[Endpoint]] = set()
+    for branch in branches:
+        required_pairs.update(
+            {
+                frozenset((Endpoint(branch.mass.name, "port_1"), Endpoint(branch.piston.name, "port_2"))),
+                frozenset((Endpoint(branch.piston.name, "port_1"), Endpoint(branch.chamber.name, "port_3"))),
+                frozenset((Endpoint(branch.chamber.name, "port_1"), Endpoint(branch.pipe.name, "port_1"))),
+                frozenset((Endpoint(branch.piston.name, "port_5"), Endpoint(branch.contact.name, "port_1"))),
+            }
+        )
+    actual_pairs = {
+        frozenset(connection.endpoints)
+        for connection in system.network.connections
+        if connection.kind == "physical"
+    }
+    if not required_pairs <= actual_pairs:
+        return _failed("targetTopologyMismatch")
+
+    required_methods = (
+        ("piston", "linearize_geometry_and_force"),
+        ("chamber", "linearize_state_derivative"),
+        ("pipe", "linearize_mass_flow"),
+        ("pipe", "linearize_state_derivative"),
+        ("contact", "linearize_contact_force"),
+        ("mass", "linearize_state_derivative"),
+    )
+    for branch in branches:
+        for owner, method in required_methods:
+            if not callable(getattr(getattr(branch, owner), method, None)):
+                return _failed(f"missingTangentPrimitive:{owner}.{method}")
+        if not callable(getattr(branch.chamber.medium, "linearize_properties_from_mU", None)):
+            return _failed("missingTangentPrimitive:medium.linearize_properties_from_mU")
+
+    if any(
+        len(group.components) != 1
+        for group in system.mechanical_state_reducer.groups
+    ):
+        return _failed("rigidMassAggregation")
+
+    target_mass_names = {branch.mass.name for branch in branches}
+    target_chamber_names = {branch.chamber.name for branch in branches}
+    target_pipe_names = {branch.pipe.name for branch in branches}
+    target_piston_names = {branch.piston.name for branch in branches}
+    target_contact_names = {branch.contact.name for branch in branches}
+    reached_ids: set[str] = set()
+    for variable in ("x", "v"):
+        for assignment in solver._causal_effort_plan_by_variable.get(
+            variable,
+            (),
+        ):
+            if any(
+                member.component in target_mass_names
+                for member in assignment.members
+            ):
+                reached_ids.update(member.id for member in assignment.members)
+    for assignment in solver._causal_effort_plan_by_variable.get("p", ()):
+        if any(
+            member.component in target_chamber_names
+            for member in assignment.members
+        ):
+            reached_ids.update(member.id for member in assignment.members)
+
+    equations = {
+        item.id: item for item in solver.equation_templates
+    }
+    reached_assignments = []
+    allowed_reached_models = {
+        "amesim_pnl0001",
+        "amesim_pnrp17",
+        "amesim_lstp00a",
+    }
+    for stage in solver._explicit_flow_plan:
+        stage_reached = []
+        for assignment in stage.assignments:
+            equation = equations.get(assignment.equation_id)
+            if equation is None:
+                return _failed("causalEquationMissing")
+            dependencies = set(equation.variables) - {
+                assignment.unknown.id
+            }
+            if not dependencies.intersection(reached_ids):
+                continue
+            component = assignment.component
+            if equation.relation != "sumToZero":
+                model_type = getattr(component, "MODEL_TYPE", None)
+                if model_type not in allowed_reached_models:
+                    return _failed(f"unsupportedReach:{model_type}")
+                expected_names = {
+                    "amesim_pnl0001": target_pipe_names,
+                    "amesim_pnrp17": target_piston_names,
+                    "amesim_lstp00a": target_contact_names,
+                }[model_type]
+                if component.name not in expected_names:
+                    return _failed(
+                        f"unexpectedReach:{model_type}:{component.name}"
+                    )
+            if bool(
+                getattr(
+                    component,
+                    "PRESSURE_FLOW_DEPENDS_ON_STREAM",
+                    False,
+                )
+            ):
+                return _failed(f"streamSensitiveReach:{component.name}")
+            stage_reached.append(assignment)
+        reached_assignments.extend(stage_reached)
+        reached_ids.update(item.unknown.id for item in stage_reached)
+
+    # The only non-zero h_outflow seeds are the target PNCH ports.  Each
+    # direct neighbour must consume it in a supported target balance, mirror
+    # it through the one-port piston, or terminate at a fixed PNPL01 cap.
+    # Secondary pressure blocks may contain the same causal flow coordinates,
+    # so membership alone is not evidence of a stream derivative.  The direct
+    # enthalpy reach proof below, plus the runtime dynamic-owner gate, is the
+    # relevant condition for this target-specific program.
+    neighbor_by_endpoint: dict[Endpoint, Endpoint] = {}
+    for connection in system.network.connections:
+        if connection.kind != "physical":
+            continue
+        first, second = connection.endpoints
+        neighbor_by_endpoint[first] = second
+        neighbor_by_endpoint[second] = first
+    permitted_h_neighbors = {
+        "amesim_pnl0001",
+        "amesim_pnrp17",
+        "amesim_pnpl01",
+    }
+    for branch in branches:
+        for port_name in branch.chamber.ports:
+            neighbor = neighbor_by_endpoint.get(
+                Endpoint(branch.chamber.name, port_name)
+            )
+            if neighbor is None:
+                return _failed("targetStreamBindingMissing")
+            component = system.network.components[neighbor.component]
+            if (
+                getattr(component, "MODEL_TYPE", None)
+                not in permitted_h_neighbors
+            ):
+                return _failed(
+                    f"unsupportedStreamReach:{component.name}"
+                )
+            if bool(
+                getattr(
+                    component,
+                    "PRESSURE_FLOW_DEPENDS_ON_STREAM",
+                    False,
+                )
+            ):
+                return _failed(f"streamSensitiveReach:{component.name}")
+
+    columns = tuple(
+        sorted(
+            index
+            for branch in branches
+            for index in (branch.velocity_index, branch.position_index)
+        )
+    )
+    provider = ThreePistonTangentProvider(
+        system,
+        tuple(branches),
+        columns,
+        state_offsets,
+        reached_assignment_count=len(reached_assignments),
+    )
+    return ThreePistonTangentCompilation(
+        True,
+        None,
+        columns,
+        provider,
+        reached_assignment_count=len(reached_assignments),
+    )

app/simulation/systems/generic.py

@@ -3,6 +3,7 @@ from __future__ import annotations
 from collections.abc import Callable
 from dataclasses import dataclass, replace
 from math import floor, isfinite
+import os
 from typing import Literal
 
 from app.simulation.core.base import Component, DynamicComponent
@@ -12,6 +13,10 @@ from app.simulation.performance import performance_span, profile_phase
 from app.simulation.property_cache import with_property_cache
 from app.simulation.solvers.algebraic import PressureFlowSolver
 from app.simulation.solvers.algebraic_blocks import StreamPressureBlockSolver
+from app.simulation.solvers.jacobian import (
+    SparseJacobianCompatibilityError,
+    SparseSecantJacobian,
+)
 from app.simulation.solvers.mechanical import (
     MechanicalConstraintGroup,
     MechanicalStateReducer,
@@ -21,15 +26,58 @@ from app.simulation.solvers.pneumatic_storage import (
     ideal_storage_group_is_reducible,
 )
 from app.simulation.solvers.pneumatic_volume import PneumaticVolumeResolver
-from app.simulation.solvers.solver import ODESolution, SolveIVPConfig, integrate_ode
+from app.simulation.solvers.solver import (
+    IntegrationCancelled,
+    ODESolution,
+    SolveIVPConfig,
+    integrate_ode,
+)
 from app.simulation.solvers.signal import SignalResolver
 from app.simulation.solvers.stream import StreamResolver
+from app.simulation.solvers.tangent import (
+    ThreePistonTangentCompilation,
+    ThreePistonTangentProvider,
+    compile_three_piston_tangent_provider,
+)
 from app.simulation.systems.network import Endpoint, SimulationNetwork
 
 
 SimulationProgressCallback = Callable[[float, str], None]
 SimulationCancellationCheck = Callable[[], bool]
 SimulationRunStatus = Literal["completed", "cancelled", "failed"]
+ODE_JACOBIAN_MODE_ENVIRONMENT_VARIABLE = "SIMULATION_ODE_JACOBIAN_MODE"
+
+
+def _requested_ode_jacobian_mode() -> Literal[
+    "optimized",
+    "hybrid",
+    "semi-analytic",
+    "scipy",
+]:
+    value = os.getenv(
+        ODE_JACOBIAN_MODE_ENVIRONMENT_VARIABLE,
+        "scipy",
+    ).strip().lower()
+    if value in {"optimized", "colored"}:
+        return "optimized"
+    if value in {"hybrid", "secant"}:
+        return "hybrid"
+    if value in {"semi-analytic", "semi_analytic", "analytic"}:
+        return "semi-analytic"
+    if value in {
+        "scipy",
+        "native",
+        "finite-difference",
+        "0",
+        "false",
+        "no",
+        "off",
+    }:
+        return "scipy"
+    raise ValueError(
+        f"{ODE_JACOBIAN_MODE_ENVIRONMENT_VARIABLE} must be "
+        "'optimized', 'hybrid', 'semi-analytic', or 'scipy'."
+    )
 
 
 @dataclass(frozen=True)
@@ -398,6 +446,7 @@ class GenericFluidSystem:
         self.signal_propagation_count = 0
         self.pneumatic_volume_propagation_count = 0
         self._jacobian_sparsity = None
+        self._ode_tangent_provider: ThreePistonTangentProvider | None = None
 
     def _request_causal_residual_audit(self) -> None:
         """Make topology or mode boundaries verify the next causal closure."""
@@ -874,6 +923,23 @@ class GenericFluidSystem:
             "colorGroupCount": group_count,
         }
 
+    def _exact_ode_jacobian_rows(self) -> dict[int, dict[int, float]]:
+        """Return mode-independent kinematic rows safe to evaluate exactly."""
+
+        rows: dict[int, dict[int, float]] = {}
+        cursor = 0
+        for entry in self.mechanical_state_reducer.state_entries:
+            if isinstance(entry, MechanicalConstraintGroup):
+                # A discrete endstop can replace x' = v with x' = 0 for the
+                # active constrained mode.  Keep those rows numerical; free
+                # mechanical groups always have d(x')/d(v) = 1.
+                if not entry.discrete_endstop_components:
+                    rows[cursor + 1] = {cursor: 1.0}
+                cursor += 2
+            else:
+                cursor += entry.state_size
+        return rows
+
     @profile_phase(
         "simulation.closure",
         minimum_mode="audit",
@@ -1065,7 +1131,11 @@ class GenericFluidSystem:
     def rhs(self, _time: float, state_vector: list[float]) -> list[float]:
         self.apply_state_vector(state_vector)
         connected_h = self._close_current_state(_time)
-        return self._state_derivatives(connected_h)
+        derivatives = self._state_derivatives(connected_h)
+        provider = self._ode_tangent_provider
+        if provider is not None:
+            provider.record_primal(_time, state_vector, connected_h)
+        return derivatives
 
     def _append_current_state(self, series: dict[str, list[float]]) -> None:
         for component in self.network.components.values():
@@ -1156,29 +1226,112 @@ class GenericFluidSystem:
                 report_solver_time(time)
             return self.rhs(time, state_vector)
 
+        jacobian = None
+        jacobian_fallback_reason: str | None = None
+        tangent_compilation: ThreePistonTangentCompilation | None = None
+        selected_tangent_provider: ThreePistonTangentProvider | None = None
+        self._ode_tangent_provider = None
+        requested_jacobian_mode = (
+            _requested_ode_jacobian_mode()
+            if jac_sparsity is not None
+            else "scipy"
+        )
+        if (
+            jac_sparsity is not None
+            and requested_jacobian_mode
+            in {"optimized", "hybrid", "semi-analytic"}
+        ):
+            state_count = int(jac_sparsity.shape[0])
+            if (
+                requested_jacobian_mode == "hybrid"
+                and not self.pressure_flow_solver.causal_fast_path_enabled
+            ):
+                jacobian_fallback_reason = "causalAlgebraicExecutionUnavailable"
+            elif int(jac_sparsity.nnz) >= state_count * state_count:
+                jacobian_fallback_reason = "denseStateDependencyPattern"
+            else:
+                exact_columns = None
+                if requested_jacobian_mode == "semi-analytic":
+                    tangent_compilation = (
+                        compile_three_piston_tangent_provider(self)
+                    )
+                    if tangent_compilation.eligible:
+                        provider = tangent_compilation.provider
+                        if provider is None:
+                            raise RuntimeError(
+                                "An eligible tangent compilation has no provider."
+                            )
+                        selected_tangent_provider = provider
+                        exact_columns = (
+                            tangent_compilation.columns,
+                            provider,
+                        )
+                    else:
+                        jacobian_fallback_reason = (
+                            f"semiAnalytic:{tangent_compilation.reason}"
+                        )
+
+                if (
+                    requested_jacobian_mode != "semi-analytic"
+                    or tangent_compilation is not None
+                    and tangent_compilation.eligible
+                ):
+                    def evaluate_jacobian_rhs(time, state):
+                        if cancel_check is not None and cancel_check():
+                            raise IntegrationCancelled
+                        return monitored_rhs(
+                            time,
+                            [float(value) for value in state],
+                        )
+
+                    try:
+                        jacobian = SparseSecantJacobian(
+                            evaluate_jacobian_rhs,
+                            jac_sparsity,
+                            integration_config.atol,
+                            exact_rows=self._exact_ode_jacobian_rows(),
+                            exact_columns=exact_columns,
+                            max_consecutive_reuses=(
+                                1
+                                if requested_jacobian_mode == "hybrid"
+                                else 0
+                            ),
+                        )
+                        self._ode_tangent_provider = (
+                            selected_tangent_provider
+                        )
+                    except SparseJacobianCompatibilityError as exc:
+                        jacobian_fallback_reason = (
+                            f"scipyCompatibility:{type(exc).__name__}"
+                        )
+
         def handle_state_transition(*args):
             transition = self.mechanical_state_reducer.state_transition(*args)
             if transition is not None:
                 self._request_causal_residual_audit()
             return transition
 
-        solution = integrate_ode(
-            rhs=monitored_rhs,
-            initial_state=initial_state,
-            config=integration_config,
-            t_eval=t_eval,
-            cancel_check=cancel_check,
-            accepted_step_callback=(
-                report_solver_time if cancel_check is not None else None
-            ),
-            breakpoints=signal_event_times,
-            state_transition_handler=(
-                handle_state_transition
-                if self.mechanical_state_reducer.has_state_events
-                else None
-            ),
-            jac_sparsity=jac_sparsity,
-        )
+        try:
+            solution = integrate_ode(
+                rhs=monitored_rhs,
+                initial_state=initial_state,
+                config=integration_config,
+                t_eval=t_eval,
+                cancel_check=cancel_check,
+                accepted_step_callback=(
+                    report_solver_time if cancel_check is not None else None
+                ),
+                breakpoints=signal_event_times,
+                state_transition_handler=(
+                    handle_state_transition
+                    if self.mechanical_state_reducer.has_state_events
+                    else None
+                ),
+                jac_sparsity=jac_sparsity,
+                jac=jacobian,
+            )
+        finally:
+            self._ode_tangent_provider = None
         if isinstance(solution, ODESolution):
             run_status: SimulationRunStatus = solution.status
             integration_error = solution.error
@@ -1212,6 +1365,44 @@ class GenericFluidSystem:
                     "recoverableRetryCount": 0,
                 }
             ]
+            if jacobian is not None:
+                direct_jacobian = jacobian.diagnostics()
+                solver_segment_diagnostics[0].update(
+                    {
+                        "jacobianEvaluationCount": int(
+                            direct_jacobian["jacobianEvaluationCount"]
+                        ),
+                        "jacobianFullBuildCount": int(
+                            direct_jacobian["fullBuildCount"]
+                        ),
+                        "jacobianSecantReuseCount": int(
+                            direct_jacobian["secantReuseCount"]
+                        ),
+                        "jacobianAuditFailureCount": int(
+                            direct_jacobian["auditFailureCount"]
+                        ),
+                        "finiteDifferenceRhsEvaluationCount": int(
+                            direct_jacobian[
+                                "finiteDifferenceRhsEvaluationCount"
+                            ]
+                        ),
+                        "jacobianBaseRhsEvaluationCount": int(
+                            direct_jacobian["baseRhsEvaluationCount"]
+                        ),
+                        "jacobianJvAuditRhsEvaluationCount": int(
+                            direct_jacobian["jvAuditEvaluationCount"]
+                        ),
+                        "exactColumnBuildCount": int(
+                            direct_jacobian["exactColumnBuildCount"]
+                        ),
+                        "exactColumnFallbackCount": int(
+                            direct_jacobian["exactColumnFallbackCount"]
+                        ),
+                        "jacobianAssemblySeconds": float(
+                            direct_jacobian["assemblySeconds"]
+                        ),
+                    }
+                )
         solver_total_keys = (
             "nfev",
             "njev",
@@ -1225,21 +1416,123 @@ class GenericFluidSystem:
             key: sum(int(segment[key]) for segment in solver_segment_diagnostics)
             for key in solver_total_keys
         }
+        jacobian_work_keys = (
+            "jacobianEvaluationCount",
+            "jacobianFullBuildCount",
+            "jacobianSecantReuseCount",
+            "jacobianAuditFailureCount",
+            "finiteDifferenceRhsEvaluationCount",
+            "jacobianBaseRhsEvaluationCount",
+            "jacobianJvAuditRhsEvaluationCount",
+            "exactColumnBuildCount",
+            "exactColumnFallbackCount",
+            "jacobianAssemblySeconds",
+        )
+        for key in jacobian_work_keys:
+            if any(key in segment for segment in solver_segment_diagnostics):
+                solver_totals[key] = sum(
+                    segment.get(key, 0)
+                    for segment in solver_segment_diagnostics
+                )
         jacobian_diagnostics = (
             self.jacobian_sparsity_diagnostics()
             if integration_config.method in {"BDF", "Radau"}
             else None
         )
+        runtime_jacobian_diagnostics: dict[str, object] | None = None
         if jacobian_diagnostics is not None:
             color_group_count = int(jacobian_diagnostics["colorGroupCount"])
-            for segment in solver_segment_diagnostics:
-                segment["finiteDifferenceRhsEstimate"] = (
-                    int(segment["njev"]) * color_group_count
+            if jacobian is None:
+                for segment in solver_segment_diagnostics:
+                    segment["finiteDifferenceRhsEstimate"] = (
+                        int(segment["njev"]) * color_group_count
+                    )
+                runtime_jacobian_diagnostics = {
+                    "mode": "scipySparseFiniteDifference",
+                    "fallbackReason": jacobian_fallback_reason,
+                    "jacobianEvaluationCount": int(solver_totals["njev"]),
+                    "fullBuildCount": int(solver_totals["njev"]),
+                    "finiteDifferenceRhsEstimateIsExact": False,
+                }
+            else:
+                for segment in solver_segment_diagnostics:
+                    segment["finiteDifferenceRhsEstimate"] = int(
+                        segment.get("finiteDifferenceRhsEvaluationCount", 0)
+                    ) + int(
+                        segment.get("jacobianJvAuditRhsEvaluationCount", 0)
+                    )
+                runtime_jacobian_diagnostics = dict(jacobian.diagnostics())
+                runtime_jacobian_diagnostics.update(
+                    {
+                        "fallbackReason": jacobian_fallback_reason,
+                        "finiteDifferenceRhsEstimateIsExact": True,
+                    }
                 )
+            if tangent_compilation is not None:
+                runtime_jacobian_diagnostics["tangentCompilation"] = (
+                    tangent_compilation.diagnostics()
+                )
+                if tangent_compilation.eligible and jacobian is not None:
+                    runtime_jacobian_diagnostics["mode"] = (
+                        "semiAnalyticExactColumns"
+                    )
+                    exact_builds = int(
+                        runtime_jacobian_diagnostics[
+                            "exactColumnBuildCount"
+                        ]
+                    )
+                    exact_fallbacks = int(
+                        runtime_jacobian_diagnostics[
+                            "exactColumnFallbackCount"
+                        ]
+                    )
+                    if exact_builds == 0 and exact_fallbacks == 0:
+                        effective_mode = "notEvaluated"
+                    elif exact_builds == 0:
+                        effective_mode = "numericalFallbackOnly"
+                    elif exact_fallbacks:
+                        effective_mode = "mixedExactAndNumericalFallback"
+                    else:
+                        effective_mode = "exactColumns"
+                    runtime_jacobian_diagnostics["effectiveMode"] = (
+                        effective_mode
+                    )
+                    if exact_fallbacks:
+                        runtime_jacobian_diagnostics[
+                            "runtimeFallbackReason"
+                        ] = runtime_jacobian_diagnostics[
+                            "lastExactColumnFallbackReason"
+                        ]
             solver_totals["finiteDifferenceRhsEstimate"] = sum(
                 int(segment["finiteDifferenceRhsEstimate"])
                 for segment in solver_segment_diagnostics
             )
+            if jacobian is None:
+                solver_totals["jacobianRhsEvaluationCountEstimate"] = (
+                    int(solver_totals["finiteDifferenceRhsEstimate"])
+                    + int(solver_totals["njev"])
+                )
+            else:
+                solver_totals["jacobianRhsEvaluationCount"] = (
+                    int(
+                        solver_totals.get(
+                            "finiteDifferenceRhsEvaluationCount",
+                            0,
+                        )
+                    )
+                    + int(
+                        solver_totals.get(
+                            "jacobianBaseRhsEvaluationCount",
+                            0,
+                        )
+                    )
+                    + int(
+                        solver_totals.get(
+                            "jacobianJvAuditRhsEvaluationCount",
+                            0,
+                        )
+                    )
+                )
 
         with performance_span("simulation.postprocessing"):
             series: dict[str, list[float]] = {"time": []}
@@ -1286,6 +1579,7 @@ class GenericFluidSystem:
                 "integration": {
                     "method": integration_config.method,
                     "jacobianSparsity": jacobian_diagnostics,
+                    "jacobian": runtime_jacobian_diagnostics,
                     "segmentCount": len(solver_segment_diagnostics),
                     "segments": solver_segment_diagnostics,
                     "totals": solver_totals,