整合求解器活动监控与步长回归证据

同步远端 PNL0003 诊断和大采样网格能力，语义合并活动感知的 60 秒真停滞判定与旧后端 15 分钟兼容兜底。

纳管热路径优化、15 单元运行证据、浏览器与 API 报告，并补充北京时间更新日志和遗留问题。

app/main.py

@@ -3,7 +3,7 @@ from __future__ import annotations
 from collections.abc import AsyncIterator, Callable, Iterator, Mapping
 from contextlib import asynccontextmanager
 import csv
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 from datetime import datetime, timezone
 import io
 import json
@@ -24,6 +24,7 @@ from pydantic import BaseModel, ConfigDict, Field, ValidationError
 
 from app.simulation.performance import performance_span, profile_phase, profile_run
 from app.simulation.property_cache import property_cache_run
+from app.simulation.solvers.solver import SolverActivityTracker
 from app.system_xml import (
     SystemXmlDocument,
     SystemXmlValidationReport,
@@ -86,6 +87,9 @@ SimulationTaskStatus = Literal[
 class SimulationTaskRecord:
     simulation_id: str
     cancel_event: threading.Event
+    activity_tracker: SolverActivityTracker = field(
+        default_factory=SolverActivityTracker
+    )
     status: SimulationTaskStatus = "queued"
     cancel_reason: SimulationCancelReason | None = None
     result: dict[str, object] | None = None
@@ -574,6 +578,7 @@ def _simulation_task_snapshot(task: SimulationTaskRecord) -> dict[str, object]:
             "cancelReason": task.cancel_reason,
             "result": task.result,
             "error": task.error,
+            **task.activity_tracker.snapshot().as_dict(),
         }
 
 
@@ -681,6 +686,7 @@ def run_system_xml_simulation(
     xml_bytes: bytes,
     progress_callback: SimulationProgressEmitter | None = None,
     cancel_check: Callable[[], bool] | None = None,
+    activity_tracker: SolverActivityTracker | None = None,
 ) -> dict[str, object]:
     with property_cache_run() as property_cache:
         with profile_run() as trace:
@@ -688,6 +694,7 @@ def run_system_xml_simulation(
                 xml_bytes,
                 progress_callback,
                 cancel_check,
+                activity_tracker,
             )
 
     performance = trace.snapshot()
@@ -712,6 +719,7 @@ def _run_system_xml_simulation_profiled(
     xml_bytes: bytes,
     progress_callback: SimulationProgressEmitter | None = None,
     cancel_check: Callable[[], bool] | None = None,
+    activity_tracker: SolverActivityTracker | None = None,
 ) -> dict[str, object]:
     from app.simulation.solvers.algebraic import AlgebraicSolveError
     from app.simulation.solvers.solver import SolveIVPConfig
@@ -774,6 +782,7 @@ def _run_system_xml_simulation_profiled(
             sample_step=document.simulation.sample_step,
             progress_callback=report_system_progress,
             cancel_check=cancel_check,
+            activity_tracker=activity_tracker,
         )
     except SimulationPreparationError as exc:
         raise HTTPException(
@@ -862,6 +871,9 @@ def simulation_event_stream(
     latest_message = "正在等待仿真任务启动"
     latest_simulated_time: float | None = None
     latest_total_time: float | None = None
+    activity_tracker = (
+        task.activity_tracker if task is not None else SolverActivityTracker()
+    )
 
     def emit_progress(
         progress: int,
@@ -889,6 +901,7 @@ def simulation_event_stream(
             event["simulatedTime"] = latest_simulated_time
         if latest_total_time is not None:
             event["totalTime"] = latest_total_time
+        event.update(activity_tracker.snapshot().as_dict())
         events.put(event)
 
     def worker() -> None:
@@ -899,28 +912,44 @@ def simulation_event_stream(
                 xml_bytes,
                 emit_progress,
                 task.cancel_event.is_set if task is not None else None,
+                activity_tracker,
             )
             if task is not None:
                 result = _mark_simulation_task_result(task, result)
             result_status = str(result.get("status", "completed"))
+            final_simulated_time = result.get("simulatedUntil")
+            final_activity_kind = (
+                "complete" if result_status == "completed" else result_status
+            )
+            if activity_tracker.snapshot().activity_kind != final_activity_kind:
+                activity_tracker.record_phase(
+                    final_activity_kind,
+                    (
+                        float(final_simulated_time)
+                        if isinstance(final_simulated_time, (int, float))
+                        and isfinite(final_simulated_time)
+                        else None
+                    ),
+                )
             result_messages = {
                 "completed": "仿真完成",
                 "stopped": "仿真已由用户终止，已保留部分结果",
                 "stalled": "仿真因进度连接异常而终止，已保留部分结果",
                 "failed": "仿真异常终止，已保留可用的部分结果",
             }
-            events.put(
-                {
-                    "event": "result",
-                    "progress": 100 if result_status == "completed" else latest_progress,
-                    "phase": result_status,
-                    "message": result_messages.get(result_status, "仿真任务结束"),
-                    "simulatedTime": result.get("simulatedUntil"),
-                    "totalTime": result.get("requestedStopTime"),
-                    "result": result,
-                }
-            )
+            result_event = {
+                "event": "result",
+                "progress": 100 if result_status == "completed" else latest_progress,
+                "phase": result_status,
+                "message": result_messages.get(result_status, "仿真任务结束"),
+                "simulatedTime": result.get("simulatedUntil"),
+                "totalTime": result.get("requestedStopTime"),
+                "result": result,
+            }
+            result_event.update(activity_tracker.snapshot().as_dict())
+            events.put(result_event)
         except HTTPException as exc:
+            activity_tracker.record_phase("failed")
             detail = exc.detail
             message = (
                 str(detail.get("message", "仿真失败"))
@@ -935,10 +964,12 @@ def simulation_event_stream(
                 "message": message,
                 "detail": detail,
             }
+            error_event.update(activity_tracker.snapshot().as_dict())
             if task is not None:
                 _mark_simulation_task_error(task, error_event)
             events.put(error_event)
         except Exception as exc:  # pragma: no cover - last-resort stream guard
+            activity_tracker.record_phase("failed")
             error_event = {
                 "event": "error",
                 "progress": latest_progress,
@@ -947,6 +978,7 @@ def simulation_event_stream(
                 "message": "仿真服务发生未预期错误。",
                 "detail": str(exc),
             }
+            error_event.update(activity_tracker.snapshot().as_dict())
             if task is not None:
                 _mark_simulation_task_error(task, error_event)
             events.put(error_event)
@@ -964,16 +996,18 @@ def simulation_event_stream(
             try:
                 event = events.get(timeout=SIMULATION_STREAM_HEARTBEAT_SECONDS)
             except queue.Empty:
+                heartbeat_event = {
+                    "event": "progress",
+                    "progress": latest_progress,
+                    "phase": latest_phase,
+                    "message": latest_message,
+                    "heartbeat": True,
+                    "simulatedTime": latest_simulated_time,
+                    "totalTime": latest_total_time,
+                }
+                heartbeat_event.update(activity_tracker.snapshot().as_dict())
                 yield json.dumps(
-                    {
-                        "event": "progress",
-                        "progress": latest_progress,
-                        "phase": latest_phase,
-                        "message": latest_message,
-                        "heartbeat": True,
-                        "simulatedTime": latest_simulated_time,
-                        "totalTime": latest_total_time,
-                    },
+                    heartbeat_event,
                     ensure_ascii=False,
                     separators=(",", ":"),
                 ) + "\n"

app/simulation/benchmark_regression.py

@@ -1038,9 +1038,18 @@ def runtime_snapshot() -> dict[str, object]:
             "SIMULATION_CAUSAL_COORDINATE_KERNEL": os.getenv(
                 "SIMULATION_CAUSAL_COORDINATE_KERNEL", "1"
             ),
+            "SIMULATION_CAUSAL_DIRECT_SUM_ASSIGNMENTS": os.getenv(
+                "SIMULATION_CAUSAL_DIRECT_SUM_ASSIGNMENTS", "1"
+            ),
+            "SIMULATION_CAUSAL_DIRECT_EQUATION_READERS": os.getenv(
+                "SIMULATION_CAUSAL_DIRECT_EQUATION_READERS", "1"
+            ),
             "SIMULATION_CAUSAL_FAST_PATH": os.getenv(
                 "SIMULATION_CAUSAL_FAST_PATH", "1"
             ),
+            "SIMULATION_MECHANICAL_ATOL_MODE": os.getenv(
+                "SIMULATION_MECHANICAL_ATOL_MODE", "legacy"
+            ),
         },
     }
 

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

@@ -60,6 +60,9 @@ class AmesimPnor001(AlgebraicComponent):
     MODEL_TYPE = "amesim_pnor001"
     MODEL_VERSION = "0.3.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = True
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_1", nominal_role="bidirectional"),
         PortDefinition.pneumatic("port_2", nominal_role="bidirectional"),
@@ -465,6 +468,9 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent):
     MODEL_TYPE = "amesim_pnvo001_fixed"
     MODEL_VERSION = "0.2.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = True
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_2", nominal_role="bidirectional"),
         PortDefinition.pneumatic("port_3", nominal_role="bidirectional"),
@@ -900,6 +906,10 @@ class AmesimPnvo001SignalOpening(AmesimPnvo001FixedOpening):
     MODEL_TYPE = "amesim_pnvo001"
     MODEL_VERSION = "0.2.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = True
+    # Repeat the exact-sum promise on this concrete subclass deliberately.
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.signal("res", nominal_role="input"),
         PortDefinition.pneumatic("port_2", nominal_role="bidirectional"),

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

@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from collections.abc import Mapping, Sequence
+from collections.abc import Callable, Mapping, Sequence
 from dataclasses import dataclass
 from functools import lru_cache
 from math import isclose, isfinite, log, log10, pi, sqrt, tanh
@@ -71,6 +71,47 @@ _PN2PIPEFR_TRANSITION_SHARPNESS = 8.37293695
 _PN2PIPEFR_ANALYTIC_LAMINAR_MAX_REYNOLDS = 1000.0
 
 
+def _pn2pipefr_friction_factor_with_precomputed_fully_rough(
+    reynolds_number: float,
+    *,
+    relative_roughness: float,
+    fully_rough: float | None,
+) -> float:
+    """Evaluate the built-in friction law with its rr-only term supplied."""
+
+    if reynolds_number <= 0.0:
+        return 64_000_000.0
+    laminar = 64.0 / reynolds_number
+    if reynolds_number <= _PN2PIPEFR_TRANSITION_START_REYNOLDS:
+        return laminar
+
+    # Keep this arithmetic in the same order as AmesimPnl00r.friction_factor.
+    # The specialized fixed-point path only moves the rr-only logarithm out of
+    # the iteration; subclasses continue to use the public virtual method.
+    smooth_turbulent = 1.0 / (
+        -1.8 * log10(6.9 / reynolds_number)
+    ) ** 2
+    if relative_roughness <= 0.0:
+        turbulent = smooth_turbulent
+    else:
+        assert fully_rough is not None
+        roughness_reynolds = reynolds_number * relative_roughness
+        roughness_reynolds_squared = roughness_reynolds * roughness_reynolds
+        roughness_weight = roughness_reynolds_squared / (
+            roughness_reynolds_squared + 180.0 * 180.0
+        )
+        turbulent = smooth_turbulent + roughness_weight * (
+            fully_rough - smooth_turbulent
+        )
+    transition_coordinate = (
+        (reynolds_number - _PN2PIPEFR_TRANSITION_START_REYNOLDS)
+        / _PN2PIPEFR_TRANSITION_SCALE_REYNOLDS
+    )
+    transition_power = transition_coordinate**_PN2PIPEFR_TRANSITION_SHARPNESS
+    transition_weight = transition_power / (1.0 + transition_power)
+    return laminar + transition_weight * (turbulent - laminar)
+
+
 def _reported_friction_factor(value: float) -> float:
     return min(float(value), _MAX_REPORTED_FRICTION_FACTOR)
 
@@ -123,6 +164,9 @@ class AmesimPnl00r(AlgebraicComponent):
     MODEL_TYPE = "amesim_pnl00r"
     MODEL_VERSION = "0.3.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = True
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_1", nominal_role="bidirectional"),
         PortDefinition.pneumatic("port_2", nominal_role="bidirectional"),
@@ -876,11 +920,11 @@ class AmesimPnl0001(ThermodynamicVolumeComponent):
                 value *= tanh(max(smoothing_argument, 0.0))
             return value
 
+        pressure_ratio_flow_parameter = mass_flow_parameter(pressure_ratio)
         if analytic_laminar:
-            flow_parameter = mass_flow_parameter(pressure_ratio)
             viscosity = self._dynamic_viscosity(T_up)
             laminar_mass_flow = (
-                self.area * p_up * flow_parameter
+                self.area * p_up * pressure_ratio_flow_parameter
             ) ** 2 / (16.0 * pi * viscosity * resistance_length * T_up)
             if (
                 self.reynolds_number(laminar_mass_flow, T_up)
@@ -888,27 +932,58 @@ class AmesimPnl0001(ThermodynamicVolumeComponent):
             ):
                 return laminar_mass_flow
 
+        sqrt_temperature = sqrt(T_up)
+        compiled_reynolds = (
+            type(self).reynolds_number is AmesimPnl0001.reynolds_number
+            and type(self)._dynamic_viscosity is AmesimPnl0001._dynamic_viscosity
+        )
+        reynolds_denominator = (
+            pi * self.diam * self._dynamic_viscosity(T_up)
+            if max_iterations > 0 and compiled_reynolds
+            else None
+        )
+        compiled_friction_factor = type(self) in (
+            AmesimPnl00r,
+            AmesimPnl0001,
+            AmesimPnl0002,
+        )
+        fully_rough = (
+            1.0 / (-2.0 * log10(self.rr / 3.7)) ** 2
+            if max_iterations > 0
+            and compiled_friction_factor
+            and self.rr > 0.0
+            else None
+        )
+
         def target_flow(mass_flow: float) -> float:
-            reynolds = self.reynolds_number(mass_flow, T_up)
-            friction = self.friction_factor(reynolds)
-            flow_coefficient = sqrt(
-                self.diam / (resistance_length * friction)
+            reynolds = (
+                4.0 * abs(mass_flow) / reynolds_denominator
+                if reynolds_denominator is not None
+                else self.reynolds_number(mass_flow, T_up)
+            )
+            friction = (
+                _pn2pipefr_friction_factor_with_precomputed_fully_rough(
+                    reynolds,
+                    relative_roughness=self.rr,
+                    fully_rough=fully_rough,
+                )
+                if compiled_friction_factor
+                else self.friction_factor(reynolds)
             )
             return (
-                flow_coefficient
+                sqrt(self.diam / (resistance_length * friction))
                 * self.area
                 * p_up
-                * mass_flow_parameter(pressure_ratio)
-                / sqrt(T_up)
+                * pressure_ratio_flow_parameter
+                / sqrt_temperature
             )
 
-        flow_coefficient = sqrt(self.diam / (resistance_length * 0.02))
         magnitude = (
-            flow_coefficient
+            sqrt(self.diam / (resistance_length * 0.02))
             * self.area
             * p_up
-            * mass_flow_parameter(pressure_ratio)
-            / sqrt(T_up)
+            * pressure_ratio_flow_parameter
+            / sqrt_temperature
         )
         for _iteration in range(max_iterations):
             next_magnitude = target_flow(magnitude)
@@ -1069,6 +1144,23 @@ class AmesimPnl0001(ThermodynamicVolumeComponent):
             self.port_1.m_flow - self.mass_flow(self.port_1.p, props.p, props.T),
         )
 
+    def pressure_flow_equation_value_readers(
+        self,
+    ) -> Mapping[str, Callable[[], float]]:
+        """Expose the independently evaluable state-pressure residual."""
+
+        def pressure_state_residual() -> float:
+            props = self.medium.properties_from_mU(
+                self.state.m,
+                self.state.U,
+                self.volume,
+            )
+            return self.port_2.p - props.p
+
+        return {
+            f"{self.name}:port_2_pressure_state": pressure_state_residual,
+        }
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume)
         return (

app/simulation/components/amesim/junctions/nodes.py

@@ -163,6 +163,9 @@ class AmesimPn3Node2(_AmesimPneumaticNode):
     MODEL_TYPE = "amesim_pn3node2"
     MODEL_VERSION = "0.3.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = False
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_1", nominal_role="bidirectional"),
         PortDefinition.pneumatic("port_2", nominal_role="bidirectional"),
@@ -200,6 +203,9 @@ class AmesimP4Node2(_AmesimPneumaticNode):
     MODEL_TYPE = "amesim_p4node2"
     MODEL_VERSION = "0.3.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = False
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_1", nominal_role="bidirectional"),
         PortDefinition.pneumatic("port_2", nominal_role="bidirectional"),

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

@@ -1386,6 +1386,9 @@ class AmesimLmechn1(AlgebraicComponent):
 
     MODEL_TYPE = "amesim_lmechn1"
     MODEL_VERSION = "0.2.0"
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("force_balance",)
+    )
     PORTS = tuple(
         PortDefinition.mechanical_translational(f"port_{index}")
         for index in range(1, 22)

app/simulation/components/experimental/flow/orifice.py

@@ -17,6 +17,9 @@ class Orifice(AlgebraicComponent):
     MODEL_TYPE = "orifice"
     MODEL_VERSION = "1.0.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = False
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_a", nominal_role="inlet"),
         PortDefinition.pneumatic("port_b", nominal_role="outlet"),
@@ -116,4 +119,3 @@ class Orifice(AlgebraicComponent):
     def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None:
         self.port_a.h_outflow = connected_h["port_b"]
         self.port_b.h_outflow = connected_h["port_a"]
-

app/simulation/components/experimental/flow/resistive_pipe.py

@@ -17,6 +17,9 @@ class ResistivePipe(AlgebraicComponent):
     MODEL_TYPE = "pipe"
     MODEL_VERSION = "1.0.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = False
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_a", nominal_role="inlet"),
         PortDefinition.pneumatic("port_b", nominal_role="outlet"),

app/simulation/components/experimental/junctions/tee.py

@@ -15,6 +15,9 @@ class Tee(AlgebraicComponent):
     MODEL_TYPE = "tee"
     MODEL_VERSION = "1.0.0"
     PRESSURE_FLOW_DEPENDS_ON_STREAM = False
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES = frozenset(
+        ("mass_flow_balance",)
+    )
     PORTS = (
         PortDefinition.pneumatic("port_in", nominal_role="bidirectional"),
         PortDefinition.pneumatic("port_out1", nominal_role="bidirectional"),

app/simulation/core/base.py

@@ -1,7 +1,7 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
-from collections.abc import Mapping
+from collections.abc import Callable, Mapping
 from typing import TYPE_CHECKING, Any, ClassVar
 
 from app.simulation.core.catalog import ComponentDisplaySpec
@@ -28,6 +28,13 @@ class Component(ABC):
     # they override either stream hook, the closure planner retains the legacy
     # full-network thermofluid fixed point.
     PRESSURE_FLOW_DEPENDS_ON_STREAM: ClassVar[bool | None] = None
+    # Exact residual suffixes whose declared variables are summed, in order,
+    # to form a ``sumToZero`` flow equation.  The causal solver deliberately
+    # reads this capability from the concrete class ``__dict__``: subclasses
+    # must repeat the promise after changing any equation semantics.
+    PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES: ClassVar[
+        frozenset[str]
+    ] = frozenset()
     PORTS: ClassVar[tuple[PortDefinition, ...]] = ()
     PARAMETERS: ClassVar[tuple[ParameterDefinition, ...]] = ()
     RESULT_VARIABLES: ClassVar[tuple[ResultVariableDefinition, ...]] = ()
@@ -245,6 +252,18 @@ class Component(ABC):
             for equation in self.pressure_flow_equation_residuals()
         )
 
+    def pressure_flow_equation_value_readers(
+        self,
+    ) -> Mapping[str, Callable[[], float]]:
+        """Return explicitly separable scalar residual readers.
+
+        The solver consumes this optional capability only when the concrete
+        component class declares the method itself.  Subclasses therefore
+        cannot accidentally inherit an equation-purity promise.
+        """
+
+        return {}
+
     def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None:
         """Update connector outflow properties from current flow directions."""
 

app/simulation/max_step_matrix.py

@@ -399,6 +399,19 @@ def _case_observation(
         checkpoint_details.append(detail)
         if not detail["available"]:
             issues.append("checkpointUnavailable")
+        requested_time = checkpoint.get("requestedTime")
+        actual_time = checkpoint.get("actualTime")
+        if (
+            not isinstance(requested_time, (int, float))
+            or not isinstance(actual_time, (int, float))
+            or not math.isclose(
+                float(requested_time),
+                float(actual_time),
+                rel_tol=0.0,
+                abs_tol=tolerance.checkpoint_time_absolute_seconds,
+            )
+        ):
+            issues.append("checkpointTimeMismatch")
         if expected_projection_count is not None and len(state_values) != int(
             expected_projection_count
         ):
@@ -455,6 +468,7 @@ def _pair_checkpoints(
     right: Mapping[str, object],
     *,
     tolerance: MatrixTolerance,
+    strict_prefix_stop: float | None = None,
 ) -> dict[str, object]:
     left_checkpoints = _checkpoints(_summary(left))
     right_checkpoints = _checkpoints(_summary(right))
@@ -463,6 +477,11 @@ def _pair_checkpoints(
         left_time = left_checkpoint.get("requestedTime")
         if not isinstance(left_time, (int, float)):
             continue
+        if strict_prefix_stop is not None and not (
+            float(left_time)
+            < strict_prefix_stop - tolerance.checkpoint_time_absolute_seconds
+        ):
+            continue
         match = next(
             (
                 right_checkpoint
@@ -554,6 +573,7 @@ def _pair_checkpoints(
         "commonCheckpointTimes": [
             pair[0].get("requestedTime") for pair in paired
         ],
+        "strictPrefixStopTime": strict_prefix_stop,
         "valueMismatchCount": mismatch_count,
         "keySetMismatchCount": key_set_mismatch_count,
         "nonnumericValueCount": nonnumeric_count,
@@ -569,6 +589,7 @@ def _compare_time_sequences(
     *,
     prefix_stop: float,
     absolute_tolerance: float,
+    strict_prefix: bool,
 ) -> dict[str, object]:
     left_times = _finite_event_times(left)
     right_times = _finite_event_times(right)
@@ -579,12 +600,12 @@ def _compare_time_sequences(
             "left": left,
             "right": right,
         }
-    left_prefix = [
-        value for value in left_times if value <= prefix_stop + absolute_tolerance
-    ]
-    right_prefix = [
-        value for value in right_times if value <= prefix_stop + absolute_tolerance
-    ]
+    if strict_prefix:
+        include = lambda value: value < prefix_stop - absolute_tolerance
+    else:
+        include = lambda value: value <= prefix_stop + absolute_tolerance
+    left_prefix = [value for value in left_times if include(value)]
+    right_prefix = [value for value in right_times if include(value)]
     passed = len(left_prefix) == len(right_prefix) and all(
         math.isclose(
             left_value,
@@ -607,6 +628,7 @@ def _pair_events(
     right: Mapping[str, object],
     *,
     prefix_stop: float,
+    strict_prefix: bool,
     tolerance: MatrixTolerance,
 ) -> dict[str, object]:
     left_trace = _event_trace(_summary(left))
@@ -616,12 +638,14 @@ def _pair_events(
         right_trace.get("signalEventTimes"),
         prefix_stop=prefix_stop,
         absolute_tolerance=tolerance.signal_event_time_absolute_seconds,
+        strict_prefix=strict_prefix,
     )
     mechanical = _compare_time_sequences(
         left_trace.get("mechanicalTransitionTimes"),
         right_trace.get("mechanicalTransitionTimes"),
         prefix_stop=prefix_stop,
         absolute_tolerance=tolerance.event_time_absolute_seconds,
+        strict_prefix=strict_prefix,
     )
     availability = (
         left_trace.get("mechanicalTransitionTimesAvailable") is not False
@@ -631,6 +655,7 @@ def _pair_events(
         "evaluated": bool(left_trace) and bool(right_trace),
         "passed": bool(signal["passed"] and mechanical["passed"] and availability),
         "prefixStopTime": prefix_stop,
+        "strictPrefix": strict_prefix,
         "signalEventTimes": signal,
         "mechanicalTransitionTimes": mechanical,
         "mechanicalTransitionTimesAvailable": availability,
@@ -701,17 +726,29 @@ def compare_matrix_cases(
     left_stop = float(left["stopTime"])
     right_stop = float(right["stopTime"])
     prefix_stop = min(left_stop, right_stop)
+    strict_prefix = not math.isclose(
+        left_stop,
+        right_stop,
+        rel_tol=0.0,
+        abs_tol=tolerance.checkpoint_time_absolute_seconds,
+    )
     if not _completed(left, left_stop) or not _completed(right, right_stop):
         return {
             "evaluated": False,
             "passed": False,
             "reason": "oneOrBothCasesDidNotComplete",
         }
-    state = _pair_checkpoints(left, right, tolerance=tolerance)
+    state = _pair_checkpoints(
+        left,
+        right,
+        tolerance=tolerance,
+        strict_prefix_stop=prefix_stop if strict_prefix else None,
+    )
     events = _pair_events(
         left,
         right,
         prefix_stop=prefix_stop,
+        strict_prefix=strict_prefix,
         tolerance=tolerance,
     )
     diagnostics = _pair_diagnostics(
@@ -728,6 +765,7 @@ def compare_matrix_cases(
         "evaluated": True,
         "passed": passed,
         "commonPrefixStopTime": prefix_stop,
+        "strictPrefix": strict_prefix,
         "stateProjection": state,
         "events": events,
         "diagnostics": diagnostics,
@@ -874,6 +912,36 @@ def run_max_step_matrix(
         if sampling_mode == "source"
         else _finite_positive(lane_config.get("sampleStep"), field="lane.sampleStep")
     )
+    source_start = float(source_config["tStart"])
+    checkpoint_grid_tolerance = max(
+        1.0e-12,
+        8.0 * math.ulp(max(1.0, abs(source_start))),
+        8.0 * math.ulp(max(1.0, abs(sample_step))),
+    )
+    for horizon in selected_horizons:
+        for checkpoint_time in horizon.checkpoint_times:
+            if math.isclose(
+                checkpoint_time,
+                horizon.stop_time,
+                rel_tol=0.0,
+                abs_tol=checkpoint_grid_tolerance,
+            ):
+                continue
+            grid_index = round((checkpoint_time - source_start) / sample_step)
+            grid_time = source_start + grid_index * sample_step
+            if not math.isclose(
+                checkpoint_time,
+                grid_time,
+                rel_tol=0.0,
+                abs_tol=checkpoint_grid_tolerance,
+            ):
+                raise RegressionManifestError(
+                    "Matrix checkpoint "
+                    f"{checkpoint_time:.17g} is not represented by the "
+                    f"{sample_step:.17g} s output grid for "
+                    f"{horizon.case_id!r}; use the event trace for off-grid "
+                    "transition times."
+                )
     instrumentation_mode = str(
         lane_config.get("instrumentationMode", "standard")
     )

app/simulation/solvers/algebraic.py

@@ -35,6 +35,12 @@ CAUSAL_EXECUTOR_V2_ENVIRONMENT_VARIABLE = "SIMULATION_CAUSAL_EXECUTOR_V2"
 CAUSAL_COORDINATE_KERNEL_ENVIRONMENT_VARIABLE = (
     "SIMULATION_CAUSAL_COORDINATE_KERNEL"
 )
+CAUSAL_DIRECT_SUM_ASSIGNMENTS_ENVIRONMENT_VARIABLE = (
+    "SIMULATION_CAUSAL_DIRECT_SUM_ASSIGNMENTS"
+)
+CAUSAL_DIRECT_EQUATION_READERS_ENVIRONMENT_VARIABLE = (
+    "SIMULATION_CAUSAL_DIRECT_EQUATION_READERS"
+)
 CAUSAL_FAST_PATH_AUDIT_INTERVAL = 64
 
 
@@ -57,6 +63,20 @@ def _causal_coordinate_kernel_environment_enabled() -> bool:
     return value.strip().lower() not in {"0", "false", "no", "off"}
 
 
+def _causal_direct_sum_assignments_environment_enabled() -> bool:
+    """Return whether exact sum-to-zero targets bypass component tuples."""
+
+    value = os.getenv(CAUSAL_DIRECT_SUM_ASSIGNMENTS_ENVIRONMENT_VARIABLE, "1")
+    return value.strip().lower() not in {"0", "false", "no", "off"}
+
+
+def _causal_direct_equation_readers_environment_enabled() -> bool:
+    """Return whether exact-class scalar residual readers are enabled."""
+
+    value = os.getenv(CAUSAL_DIRECT_EQUATION_READERS_ENVIRONMENT_VARIABLE, "1")
+    return value.strip().lower() not in {"0", "false", "no", "off"}
+
+
 class AlgebraicSolveError(RuntimeError):
     def __init__(
         self,
@@ -127,6 +147,7 @@ class EffortAnchor:
     unknown: AlgebraicUnknown
     evaluate: Callable[[], float]
     equation_id: str
+    causal_evaluate: Callable[[], float] | None = None
 
 
 @dataclass(frozen=True)
@@ -326,6 +347,14 @@ class PressureFlowSolver:
         self.residual_tolerance = residual_tolerance
         self.max_evaluations = max_evaluations
         self.scope_kind = scope_kind
+        self._causal_direct_sum_assignments_environment_enabled = (
+            _causal_direct_sum_assignments_environment_enabled()
+        )
+        self._causal_direct_sum_flow_assignment_count = 0
+        self._causal_direct_equation_readers_environment_enabled = (
+            _causal_direct_equation_readers_environment_enabled()
+        )
+        self._causal_direct_effort_anchor_count = 0
         self.unknowns = self._build_unknowns()
         self._unknowns_by_id = {unknown.id: unknown for unknown in self.unknowns}
         self._unknowns_by_variable = {
@@ -371,6 +400,9 @@ class PressureFlowSolver:
             for item in self._component_equation_plan
             for equation in item.templates
         ) + tuple(item.template for item in self._connection_equation_plan)
+        self._causal_direct_equation_readers = (
+            self._build_causal_direct_equation_readers()
+        )
         self._effort_groups = {
             variable: self._build_effort_equality_groups(variable)
             for variable in ("p", "x", "v")
@@ -418,6 +450,11 @@ class PressureFlowSolver:
             self._causal_fast_path_eligible,
             self._causal_fast_path_fallback_reason,
         ) = self._build_causal_execution_plan()
+        self._causal_direct_effort_anchor_count = sum(
+            assignment.anchor.causal_evaluate is not None
+            for assignments in self._causal_effort_plan_by_variable.values()
+            for assignment in assignments
+        )
         self._causal_fast_path_environment_enabled = (
             _causal_fast_path_environment_enabled()
         )
@@ -599,6 +636,18 @@ class PressureFlowSolver:
             "coordinateKernelFastSolveCount": (
                 self._causal_coordinate_fast_solve_count
             ),
+            "directSumAssignmentsConfigured": (
+                self._causal_direct_sum_assignments_environment_enabled
+            ),
+            "directSumFlowAssignmentCount": (
+                self._causal_direct_sum_flow_assignment_count
+            ),
+            "directEquationReadersConfigured": (
+                self._causal_direct_equation_readers_environment_enabled
+            ),
+            "directEffortAnchorCount": (
+                self._causal_direct_effort_anchor_count
+            ),
             "compiledEffortUnknownCount": (
                 self._causal_compiled_effort_unknown_count
             ),
@@ -850,6 +899,9 @@ class PressureFlowSolver:
             ] = {}
             component_evaluators: dict[int, Callable[[], tuple[float, ...]]] = {}
             for coordinate_index, assignment in indexed_assignments:
+                if assignment.anchor.causal_evaluate is not None:
+                    direct_targets.append((coordinate_index, assignment))
+                    continue
                 equation_index = equation_index_by_id[
                     assignment.anchor.equation_id
                 ]
@@ -934,10 +986,22 @@ class PressureFlowSolver:
             if stage is None:
                 return False
             for coordinate_index, assignment in stage.direct_targets:
-                workspace[coordinate_index] = (
-                    self._read_effort_anchor(assignment)
-                    - assignment.anchor.evaluate()
+                evaluate = (
+                    assignment.anchor.causal_evaluate
+                    if assignment.anchor.causal_evaluate is not None
+                    else assignment.anchor.evaluate
                 )
+                try:
+                    workspace[coordinate_index] = (
+                        self._read_effort_anchor(assignment)
+                        - evaluate()
+                    )
+                except MemoryError:
+                    raise
+                except Exception:
+                    if assignment.anchor.causal_evaluate is None:
+                        raise
+                    return False
             for evaluation in stage.component_evaluations:
                 equation_values = evaluation.evaluate()
                 for target in evaluation.targets:
@@ -975,7 +1039,19 @@ class PressureFlowSolver:
                 return False
             for assignment in assignments:
                 anchor = assignment.anchor
-                target = anchor.unknown.read() - anchor.evaluate()
+                evaluate = (
+                    anchor.causal_evaluate
+                    if anchor.causal_evaluate is not None
+                    else anchor.evaluate
+                )
+                try:
+                    target = anchor.unknown.read() - evaluate()
+                except MemoryError:
+                    raise
+                except Exception:
+                    if anchor.causal_evaluate is None:
+                        raise
+                    return False
                 if not isfinite(target) or (
                     variable == "p" and target <= PRESSURE_LOWER_BOUND_PA
                 ):
@@ -1431,6 +1507,9 @@ class PressureFlowSolver:
                         unknown=unknown,
                         evaluate=self._equation_value_reader(equation),
                         equation_id=equation.id,
+                        causal_evaluate=(
+                            self._causal_direct_equation_readers.get(equation.id)
+                        ),
                     )
                 )
 
@@ -1851,6 +1930,45 @@ class PressureFlowSolver:
 
         return read_component_equation
 
+    def _build_causal_direct_equation_readers(
+        self,
+    ) -> dict[str, Callable[[], float]]:
+        """Compile exact-class scalar residual capabilities, failing closed."""
+
+        if not self._causal_direct_equation_readers_environment_enabled:
+            return {}
+
+        readers: dict[str, Callable[[], float]] = {}
+        for evaluation in self._component_equation_plan:
+            component = evaluation.component
+            # A subclass must repeat the declaration after changing equation
+            # semantics; inherited purity promises are deliberately ignored.
+            declared = type(component).__dict__.get(
+                "pressure_flow_equation_value_readers"
+            )
+            if not callable(declared):
+                continue
+            try:
+                component_readers = declared(component)
+            except MemoryError:
+                raise
+            except Exception:
+                continue
+            if not isinstance(component_readers, Mapping):
+                continue
+            template_ids = frozenset(
+                template.id for template in evaluation.templates
+            )
+            if any(
+                not isinstance(equation_id, str)
+                or equation_id not in template_ids
+                or not callable(reader)
+                for equation_id, reader in component_readers.items()
+            ):
+                continue
+            readers.update(component_readers)
+        return readers
+
     def _connection_flow_target_reader(
         self,
         equation,
@@ -1870,6 +1988,69 @@ class PressureFlowSolver:
             f"Connection flow equation {equation.id} does not contain {unknown.id}."
         )
 
+    def _sum_to_zero_flow_target_reader(
+        self,
+        equation,
+        unknown: AlgebraicUnknown,
+    ) -> Callable[[], float]:
+        """Compile the exact target of a declared flow sum without callbacks."""
+
+        if equation.relation != "sumToZero":
+            raise ValueError(
+                f"Equation {equation.id} is not a sum-to-zero relation."
+            )
+        if len(set(equation.variables)) != len(equation.variables):
+            raise ValueError(
+                f"Equation {equation.id} repeats a sum-to-zero variable."
+            )
+        equation_unknowns: list[AlgebraicUnknown] = []
+        for variable in equation.variables:
+            current = self._unknowns_by_id.get(variable)
+            if current is None or current.role != "flow":
+                raise ValueError(
+                    f"Equation {equation.id} has a non-flow variable."
+                )
+            equation_unknowns.append(current)
+        if not any(current is unknown for current in equation_unknowns):
+            raise ValueError(
+                f"Equation {equation.id} does not contain {unknown.id}."
+            )
+        compiled_unknowns = tuple(equation_unknowns)
+        if len(compiled_unknowns) == 2:
+            first, second = compiled_unknowns
+            return lambda: 0.0 - (first.read() + second.read())
+
+        def target() -> float:
+            # The target unknown is already zeroed by the stage executor.  Keep
+            # its slot in the sum so the floating-point operation order matches
+            # the declared built-in residual, including near cancellation.
+            return 0.0 - sum(current.read() for current in compiled_unknowns)
+
+        return target
+
+    @staticmethod
+    def _component_declares_exact_sum_to_zero_equation(
+        component: object,
+        equation_id: str,
+    ) -> bool:
+        """Accept only an exact concrete-class promise for callback bypass."""
+
+        declared_suffixes = type(component).__dict__.get(
+            "PRESSURE_FLOW_EXACT_SUM_TO_ZERO_EQUATION_SUFFIXES"
+        )
+        if not isinstance(declared_suffixes, frozenset) or any(
+            not isinstance(suffix, str) or not suffix or ":" in suffix
+            for suffix in declared_suffixes
+        ):
+            return False
+        component_name = getattr(component, "name", None)
+        if not isinstance(component_name, str):
+            return False
+        return any(
+            equation_id == f"{component_name}:{suffix}"
+            for suffix in declared_suffixes
+        )
+
     def _pressure_flow_equation_values(self) -> tuple[float, ...]:
         """Evaluate live equation values through the compiled topology."""
 
@@ -1918,6 +2099,29 @@ class PressureFlowSolver:
 
         evaluation = self._component_equation_plans_by_id[equation.owner_id]
         component = evaluation.component
+        if (
+            self._causal_direct_sum_assignments_environment_enabled
+            and equation.relation == "sumToZero"
+            and self._component_declares_exact_sum_to_zero_equation(
+                component,
+                equation.id,
+            )
+        ):
+            try:
+                evaluate = self._sum_to_zero_flow_target_reader(
+                    equation,
+                    unknown,
+                )
+            except ValueError:
+                pass
+            else:
+                self._causal_direct_sum_flow_assignment_count += 1
+                return ExplicitFlowAssignment(
+                    equation_id=equation.id,
+                    unknown=unknown,
+                    evaluate=evaluate,
+                )
+
         equations = evaluation.templates
         equation_ids = tuple(current.id for current in equations)
         try:

app/simulation/solvers/causal_ir.py

@@ -378,13 +378,7 @@ class CausalNumericIR:
                         workspace.effort_written[output] = True
                 except MemoryError:
                     raise
-                except (
-                    ArithmeticError,
-                    IndexError,
-                    RuntimeError,
-                    TypeError,
-                    ValueError,
-                ) as exc:
+                except Exception as exc:
                     return failed(f"effortEvaluationFailed:{type(exc).__name__}")
             if any(not bool(workspace.effort_written[index]) for index in range(width)):
                 return failed("effortEvaluationCoverageMismatch")
@@ -530,6 +524,7 @@ def _unique_slots(items: Iterable[int]) -> tuple[int, ...]:
 def _compile_effort_evaluations(
     operations: tuple[CausalIREffortOperation, ...],
     anchor_evaluators: tuple[Callable[[], float], ...],
+    direct_residuals: tuple[bool, ...],
     component_locations: dict[
         str, tuple[object, Callable[[], tuple[float, ...]], int]
     ],
@@ -538,10 +533,14 @@ def _compile_effort_evaluations(
     grouped: dict[int, list[tuple[int, int, str]]] = {}
     component_callbacks: dict[int, Callable[[], tuple[float, ...]]] = {}
     direct: list[tuple[int, Callable[[], float], str]] = []
-    for output, (operation, anchor_evaluate) in enumerate(
-        zip(operations, anchor_evaluators)
+    for output, (operation, anchor_evaluate, direct_residual) in enumerate(
+        zip(operations, anchor_evaluators, direct_residuals)
     ):
-        location = component_locations.get(operation.equation_id)
+        location = (
+            None
+            if direct_residual
+            else component_locations.get(operation.equation_id)
+        )
         if location is None:
             direct.append((output, anchor_evaluate, operation.equation_id))
             continue
@@ -638,6 +637,7 @@ def compile_causal_numeric_ir(solver: object) -> CausalIRCompilation:
         for variable in ("p", "x", "v"):
             operations: list[CausalIREffortOperation] = []
             anchors: list[Callable[[], float]] = []
+            direct_residuals: list[bool] = []
             for assignment in effort_plan[variable]:
                 result = len(canonical_slots)
                 equation_id = str(assignment.anchor.equation_id)
@@ -665,7 +665,17 @@ def compile_causal_numeric_ir(solver: object) -> CausalIRCompilation:
                         equation_id,
                     )
                 )
-                anchors.append(assignment.anchor.evaluate)
+                causal_evaluate = getattr(
+                    assignment.anchor,
+                    "causal_evaluate",
+                    None,
+                )
+                anchors.append(
+                    causal_evaluate
+                    if causal_evaluate is not None
+                    else assignment.anchor.evaluate
+                )
+                direct_residuals.append(causal_evaluate is not None)
             operation_tuple = tuple(operations)
             effort_stages.append(
                 CausalIREffortStage(
@@ -674,6 +684,7 @@ def compile_causal_numeric_ir(solver: object) -> CausalIRCompilation:
                     _compile_effort_evaluations(
                         operation_tuple,
                         tuple(anchors),
+                        tuple(direct_residuals),
                         component_locations,
                         evaluators,
                     ),

app/simulation/solvers/mechanical.py

@@ -1,9 +1,12 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
+from math import isfinite
+import os
 from typing import Callable, Literal, Mapping, Sequence
 
 from app.simulation.components.amesim.mechanical.translational import (
+    AmesimLstp00a,
     AmesimMecmas21,
 )
 from app.simulation.core.base import DynamicComponent
@@ -12,15 +15,114 @@ from app.simulation.systems.network import SimulationNetwork
 
 
 ConstraintMode = Literal["uninitialized", "free", "lower", "upper"]
+MechanicalAbsoluteToleranceMode = Literal["legacy", "contact-aware-v1"]
 DenseState = Callable[[float], Sequence[float]]
 
 
+MECHANICAL_ATOL_MODE_ENVIRONMENT_VARIABLE = (
+    "SIMULATION_MECHANICAL_ATOL_MODE"
+)
+
+
+def _requested_mechanical_absolute_tolerance_mode(
+) -> MechanicalAbsoluteToleranceMode:
+    value = os.environ.get(
+        MECHANICAL_ATOL_MODE_ENVIRONMENT_VARIABLE,
+        "legacy",
+    ).strip().lower()
+    if value == "legacy":
+        return "legacy"
+    if value in {"contact-aware-v1", "contact_aware_v1", "contact-aware"}:
+        return "contact-aware-v1"
+    raise ValueError(
+        f"{MECHANICAL_ATOL_MODE_ENVIRONMENT_VARIABLE} must be "
+        "'legacy' or 'contact-aware-v1'."
+    )
+
+
+@dataclass(frozen=True)
+class MechanicalToleranceGroupPlan:
+    """One rigid-coordinate group's state tolerances and proof result."""
+
+    components: tuple[str, ...]
+    contacts: tuple[str, ...]
+    eligible: bool
+    reason: str
+    velocity_atol: float
+    position_atol: float
+    minimum_dvel: float | None
+    minimum_contact_damping_length: float | None
+    minimum_damping_strength_ratio: float | None
+    minimum_force_limited_velocity_atol: float | None
+
+    def as_dict(self) -> dict[str, object]:
+        return {
+            "components": list(self.components),
+            "contacts": list(self.contacts),
+            "eligible": self.eligible,
+            "reason": self.reason,
+            "velocityAtol": self.velocity_atol,
+            "positionAtol": self.position_atol,
+            "minimumDvel": self.minimum_dvel,
+            "minimumContactDampingLength": (
+                self.minimum_contact_damping_length
+            ),
+            "minimumDampingStrengthRatio": (
+                self.minimum_damping_strength_ratio
+            ),
+            "minimumForceLimitedVelocityAtol": (
+                self.minimum_force_limited_velocity_atol
+            ),
+        }
+
+
+@dataclass(frozen=True)
+class MechanicalAbsoluteTolerancePlan:
+    """State-aligned absolute tolerances with auditable group proofs."""
+
+    mode: MechanicalAbsoluteToleranceMode
+    default_atol: float
+    legacy_mechanical_atol: float
+    values: tuple[float, ...]
+    groups: tuple[MechanicalToleranceGroupPlan, ...]
+
+    def as_dict(self) -> dict[str, object]:
+        legacy_value = min(
+            self.default_atol,
+            self.legacy_mechanical_atol,
+        )
+        relaxed_groups = tuple(
+            group
+            for group in self.groups
+            if group.velocity_atol > legacy_value
+        )
+        return {
+            "mode": self.mode,
+            "defaultAtol": self.default_atol,
+            "legacyMechanicalAtol": self.legacy_mechanical_atol,
+            "stateCount": len(self.values),
+            "groupCount": len(self.groups),
+            "eligibleGroupCount": sum(group.eligible for group in self.groups),
+            "relaxedVelocityGroupCount": len(relaxed_groups),
+            "relaxedVelocityStateCount": len(relaxed_groups),
+            "relaxedPositionStateCount": 0,
+            "minimumEffectiveAtol": (
+                min(self.values) if self.values else None
+            ),
+            "maximumEffectiveAtol": (
+                max(self.values) if self.values else None
+            ),
+            "groups": [group.as_dict() for group in self.groups],
+        }
+
+
 @dataclass
 class MechanicalConstraintGroup:
     """MECMAS21 inertias that share one rigid translational coordinate."""
 
     components: tuple[AmesimMecmas21, ...]
     mode: ConstraintMode = "uninitialized"
+    contact_components: tuple[AmesimLstp00a, ...] = ()
 
     @property
     def representative(self) -> AmesimMecmas21:
@@ -361,9 +463,52 @@ class MechanicalStateReducer:
             roots = (find("x", first_port), find("v", first_port))
             masses_by_roots.setdefault(roots, []).append(component)
 
+        contacts_by_roots: dict[
+            tuple[tuple[str, str], tuple[str, str]],
+            dict[str, AmesimLstp00a],
+        ] = {}
+        for component in self.network.components.values():
+            if not isinstance(component, AmesimLstp00a):
+                continue
+            contact_roots = tuple(
+                (
+                    find("x", (component.name, definition.name)),
+                    find("v", (component.name, definition.name)),
+                )
+                for definition in component.active_port_definitions
+                if (
+                    definition.kind == "physical"
+                    and definition.domain == "mechanical"
+                )
+            )
+            if (
+                len({roots[0] for roots in contact_roots}) < 2
+                or len({roots[1] for roots in contact_roots}) < 2
+            ):
+                # A compliant contact whose two ports resolve to the same
+                # rigid coordinate cannot damp that coordinate.  Treating the
+                # self-loop as proof would relax an unrelated velocity state.
+                continue
+            for definition in component.active_port_definitions:
+                if (
+                    definition.kind != "physical"
+                    or definition.domain != "mechanical"
+                ):
+                    continue
+                endpoint = (component.name, definition.name)
+                roots = (find("x", endpoint), find("v", endpoint))
+                contacts_by_roots.setdefault(roots, {})[
+                    component.name
+                ] = component
+
         return tuple(
-            MechanicalConstraintGroup(tuple(components))
-            for components in masses_by_roots.values()
+            MechanicalConstraintGroup(
+                components=tuple(components),
+                contact_components=tuple(
+                    contacts_by_roots.get(roots, {}).values()
+                ),
+            )
+            for roots, components in masses_by_roots.items()
         )
 
     def _build_state_entries(self) -> tuple[StateEntry, ...]:
@@ -391,26 +536,253 @@ class MechanicalStateReducer:
     def has_state_events(self) -> bool:
         return any(group.discrete_endstop_components for group in self.groups)
 
+    def absolute_tolerance_plan(
+        self,
+        default: float,
+        *,
+        mechanical: float = 1.0e-12,
+        mode: MechanicalAbsoluteToleranceMode | None = None,
+    ) -> MechanicalAbsoluteTolerancePlan:
+        """Compile state tolerances without weakening non-smooth coordinates.
+
+        A scalar ``1e-8`` absolute tolerance makes SciPy perturb a zero-valued
+        endstop position across the much smaller unilateral boundary band while
+        constructing finite-difference Jacobians. Positions and ideal endstop
+        states therefore retain the legacy machine-scale tolerance.
+
+        Strongly damped, compliant LSTP contact can instead drive a *free*
+        velocity close to zero for hundreds of accepted steps. Only a
+        compile-proven smooth-contact group may use the bounded velocity floor;
+        the contact position coordinate remains unchanged.
+        """
+        default_atol = float(default)
+        mechanical_atol = float(mechanical)
+        if not isfinite(default_atol) or default_atol <= 0.0:
+            raise ValueError(
+                "default absolute tolerance must be finite and positive."
+            )
+        if not isfinite(mechanical_atol) or mechanical_atol <= 0.0:
+            raise ValueError(
+                "mechanical absolute tolerance must be finite and positive."
+            )
+        selected_mode = mode or _requested_mechanical_absolute_tolerance_mode()
+        if selected_mode not in {"legacy", "contact-aware-v1"}:
+            raise ValueError(
+                "mechanical absolute tolerance mode must be 'legacy' or "
+                "'contact-aware-v1'."
+            )
+
+        legacy_atol = min(default_atol, mechanical_atol)
+        values: list[float] = []
+        group_plans: list[MechanicalToleranceGroupPlan] = []
+        for entry in self.state_entries:
+            if not isinstance(entry, MechanicalConstraintGroup):
+                values.extend([default_atol] * entry.state_size)
+                continue
+
+            components = entry.components
+            contacts = entry.contact_components
+            positive_dvel = tuple(
+                float(component.dvel)
+                for component in components
+                if isfinite(float(component.dvel))
+                and float(component.dvel) > 0.0
+            )
+            positive_pdis = tuple(
+                float(contact.Pdis)
+                for contact in contacts
+                if isfinite(float(contact.Pdis))
+                and float(contact.Pdis) > 0.0
+            )
+            minimum_dvel = min(positive_dvel, default=None)
+            minimum_pdis = min(positive_pdis, default=None)
+            contact_scale_valid = True
+            contact_force_velocity_limits_list: list[float] = []
+            damping_strength_ratios_list: list[float] = []
+            if selected_mode == "contact-aware-v1" and minimum_dvel is not None:
+                for contact in contacts:
+                    stiffness = float(contact.kcont)
+                    damping_length = float(contact.Pdis)
+                    damping = float(contact.rcont)
+                    if not (
+                        isfinite(stiffness)
+                        and stiffness > 0.0
+                        and isfinite(damping_length)
+                        and damping_length > 0.0
+                        and isfinite(damping)
+                        and damping > 0.0
+                    ):
+                        contact_scale_valid = False
+                        continue
+                    elastic_force_scale = stiffness * damping_length
+                    damping_force_scale = damping * minimum_dvel
+                    if not (
+                        isfinite(elastic_force_scale)
+                        and elastic_force_scale > 0.0
+                        and isfinite(damping_force_scale)
+                        and damping_force_scale > 0.0
+                    ):
+                        contact_scale_valid = False
+                        continue
+                    force_velocity_limit = (
+                        1.0e-3 * elastic_force_scale / damping
+                    )
+                    damping_strength_ratio = (
+                        damping_force_scale / elastic_force_scale
+                    )
+                    if not (
+                        isfinite(force_velocity_limit)
+                        and force_velocity_limit > 0.0
+                        and isfinite(damping_strength_ratio)
+                        and damping_strength_ratio > 0.0
+                    ):
+                        contact_scale_valid = False
+                        continue
+                    contact_force_velocity_limits_list.append(
+                        force_velocity_limit
+                    )
+                    damping_strength_ratios_list.append(
+                        damping_strength_ratio
+                    )
+            contact_force_velocity_limits = tuple(
+                contact_force_velocity_limits_list
+            )
+            minimum_force_velocity_atol = min(
+                contact_force_velocity_limits,
+                default=None,
+            )
+            damping_strength_ratios = tuple(
+                damping_strength_ratios_list
+            )
+            minimum_damping_strength_ratio = min(
+                damping_strength_ratios,
+                default=None,
+            )
+
+            if selected_mode == "legacy":
+                eligible = False
+                reason = "legacyMode"
+            elif entry.discrete_endstop_components:
+                eligible = False
+                reason = "discreteEndstop"
+            elif any(int(component.stoptype) != 4 for component in components):
+                eligible = False
+                reason = "unsupportedStopType"
+            elif any(
+                component.use_friction and float(component.fcoul) != 0.0
+                for component in components
+            ):
+                eligible = False
+                reason = "dryFriction"
+            elif not contacts:
+                eligible = False
+                reason = "noFlexibleContact"
+            elif any(
+                not isfinite(float(contact.Pdis))
+                or float(contact.Pdis) <= 0.0
+                for contact in contacts
+            ):
+                eligible = False
+                reason = "nonSmoothContactDampingLength"
+            elif any(
+                not isfinite(float(contact.rcont))
+                or float(contact.rcont) <= 0.0
+                for contact in contacts
+            ):
+                eligible = False
+                reason = "undampedContact"
+            elif any(
+                not isfinite(float(contact.kcont))
+                or float(contact.kcont) <= 0.0
+                for contact in contacts
+            ):
+                eligible = False
+                reason = "invalidContactStiffness"
+            elif not contact_scale_valid:
+                eligible = False
+                reason = "invalidContactScale"
+            elif any(
+                int(contact.discContactOption) != 1
+                for contact in contacts
+            ):
+                eligible = False
+                reason = "clampedContactForce"
+            elif len(positive_dvel) != len(components):
+                eligible = False
+                reason = "invalidVelocityScale"
+            elif (
+                len(damping_strength_ratios) != len(contacts)
+                or minimum_damping_strength_ratio is None
+                or minimum_damping_strength_ratio < 1.0
+            ):
+                eligible = False
+                reason = "weakContactDamping"
+            else:
+                eligible = True
+                reason = "eligibleFlexibleContact"
+
+            velocity_atol = legacy_atol
+            if eligible:
+                assert minimum_dvel is not None
+                assert minimum_force_velocity_atol is not None
+                velocity_atol = min(
+                    default_atol,
+                    max(
+                        mechanical_atol,
+                        min(
+                            1.0e-9,
+                            1.0e-3 * minimum_dvel,
+                            minimum_force_velocity_atol,
+                        ),
+                    ),
+                )
+            position_atol = legacy_atol
+            values.extend((velocity_atol, position_atol))
+            group_plans.append(
+                MechanicalToleranceGroupPlan(
+                    components=tuple(
+                        component.name for component in components
+                    ),
+                    contacts=tuple(contact.name for contact in contacts),
+                    eligible=eligible,
+                    reason=reason,
+                    velocity_atol=velocity_atol,
+                    position_atol=position_atol,
+                    minimum_dvel=minimum_dvel,
+                    minimum_contact_damping_length=minimum_pdis,
+                    minimum_damping_strength_ratio=(
+                        minimum_damping_strength_ratio
+                    ),
+                    minimum_force_limited_velocity_atol=(
+                        minimum_force_velocity_atol
+                    ),
+                )
+            )
+
+        return MechanicalAbsoluteTolerancePlan(
+            mode=selected_mode,
+            default_atol=default_atol,
+            legacy_mechanical_atol=mechanical_atol,
+            values=tuple(values),
+            groups=tuple(group_plans),
+        )
+
     def absolute_tolerances(
         self,
         default: float,
         *,
         mechanical: float = 1.0e-12,
+        mode: MechanicalAbsoluteToleranceMode | None = None,
     ) -> list[float]:
-        """Return state-aligned tolerances with machine-scale mechanics.
+        """Return state-aligned values from the auditable tolerance plan."""
 
-        A scalar ``1e-8`` absolute tolerance makes SciPy perturb a zero-valued
-        endstop position across the much smaller unilateral boundary band while
-        constructing finite-difference Jacobians. Mechanical coordinates need
-        a tighter floor; thermodynamic states retain the caller's tolerance.
-        """
-        values: list[float] = []
-        for entry in self.state_entries:
-            if isinstance(entry, MechanicalConstraintGroup):
-                values.extend([min(default, mechanical)] * 2)
-            else:
-                values.extend([default] * entry.state_size)
-        return values
+        return list(
+            self.absolute_tolerance_plan(
+                default,
+                mechanical=mechanical,
+                mode=mode,
+            ).values
+        )
 
     def reset_constraint_modes(self) -> None:
         for group in self.groups:

app/simulation/solvers/solver.py

@@ -15,6 +15,135 @@ DenseState = Callable[[float], list[float]]
 JacobianCallable = Callable[[float, object], object]
 
 
+@dataclass(frozen=True)
+class SolverActivitySnapshot:
+    """Low-cost, additive view of work inside an integration task.
+
+    ``accepted_time`` deliberately changes only after an accepted solver step.
+    Trial evaluations may continue to advance ``activity_sequence`` and
+    ``current_trial_time`` while that public progress value stays fixed.
+    """
+
+    activity_sequence: int
+    activity_kind: str
+    current_trial_time: float | None
+    rhs_call_count: int
+    accepted_step_sequence: int
+    accepted_time: float | None
+    solver_step_sequence: int
+    jacobian_evaluation_count: int
+    thermofluid_closure_count: int
+
+    def as_dict(self) -> dict[str, object]:
+        return {
+            "activitySequence": self.activity_sequence,
+            "activityKind": self.activity_kind,
+            "currentTrialTime": self.current_trial_time,
+            "rhsCallCount": self.rhs_call_count,
+            "acceptedStepSequence": self.accepted_step_sequence,
+            "acceptedTime": self.accepted_time,
+            "solverStepSequence": self.solver_step_sequence,
+            "jacobianEvaluationCount": self.jacobian_evaluation_count,
+            "thermofluidClosureCount": self.thermofluid_closure_count,
+        }
+
+
+class SolverActivityTracker:
+    """Single-writer activity telemetry for a solver worker.
+
+    The solver thread is the only writer and the stream thread only snapshots
+    scalar attributes. The sequence is published last, so a reader never
+    treats partially published fields as a newer completed activity update.
+    Passing no tracker to :func:`integrate_ode` is the zero-cost opt-out path.
+    """
+
+    __slots__ = (
+        "_accepted_step_sequence",
+        "_accepted_time",
+        "_activity_kind",
+        "_activity_sequence",
+        "_current_trial_time",
+        "_jacobian_evaluation_count",
+        "_rhs_call_count",
+        "_solver_step_sequence",
+        "_thermofluid_closure_count",
+    )
+
+    def __init__(self) -> None:
+        self._activity_sequence = 0
+        self._activity_kind = "idle"
+        self._current_trial_time: float | None = None
+        self._rhs_call_count = 0
+        self._accepted_step_sequence = 0
+        self._accepted_time: float | None = None
+        self._solver_step_sequence = 0
+        self._jacobian_evaluation_count = 0
+        self._thermofluid_closure_count = 0
+
+    def _publish(self, kind: str, time: float | None = None) -> None:
+        self._activity_kind = kind
+        if time is not None:
+            self._current_trial_time = float(time)
+        self._activity_sequence += 1
+
+    def start_integration(self, time: float) -> None:
+        self._accepted_time = float(time)
+        self._publish("solver_initialization", time)
+
+    def record_phase(self, kind: str, time: float | None = None) -> None:
+        self._publish(kind, time)
+
+    def record_solver_step(self, time: float) -> None:
+        self._solver_step_sequence += 1
+        self._publish("solver_step", time)
+
+    def record_rhs(self, time: float) -> None:
+        self._rhs_call_count += 1
+        self._publish("rhs", time)
+
+    def record_jacobian(self, time: float) -> None:
+        self._jacobian_evaluation_count += 1
+        self._publish("jacobian", time)
+
+    def record_thermofluid_closure(self, time: float) -> None:
+        self._thermofluid_closure_count += 1
+        self._publish("thermofluid_closure", time)
+
+    def record_accepted_step(self, time: float) -> None:
+        accepted_time = float(time)
+        if (
+            self._accepted_time is not None
+            and accepted_time <= self._accepted_time
+        ):
+            return
+        self._accepted_step_sequence += 1
+        self._accepted_time = accepted_time
+        self._publish("accepted_step", accepted_time)
+
+    def snapshot(self) -> SolverActivitySnapshot:
+        # ``activity_sequence`` is read last because writers publish it last.
+        activity_kind = self._activity_kind
+        current_trial_time = self._current_trial_time
+        rhs_call_count = self._rhs_call_count
+        accepted_step_sequence = self._accepted_step_sequence
+        accepted_time = self._accepted_time
+        solver_step_sequence = self._solver_step_sequence
+        jacobian_evaluation_count = self._jacobian_evaluation_count
+        thermofluid_closure_count = self._thermofluid_closure_count
+        activity_sequence = self._activity_sequence
+        return SolverActivitySnapshot(
+            activity_sequence=activity_sequence,
+            activity_kind=activity_kind,
+            current_trial_time=current_trial_time,
+            rhs_call_count=rhs_call_count,
+            accepted_step_sequence=accepted_step_sequence,
+            accepted_time=accepted_time,
+            solver_step_sequence=solver_step_sequence,
+            jacobian_evaluation_count=jacobian_evaluation_count,
+            thermofluid_closure_count=thermofluid_closure_count,
+        )
+
+
 @dataclass(frozen=True)
 class StateTransition:
     """A state reset located inside an accepted integration step."""
@@ -773,6 +902,7 @@ def _integrate_scipy_stepwise(
     state_transition_handler: StateTransitionHandler | None = None,
     jac_sparsity=None,
     jac: JacobianCallable | None = None,
+    activity_tracker: SolverActivityTracker | None = None,
 ) -> ODESolution:
     """Initial stepwise integration path for breakpoints and state resets.
 
@@ -796,6 +926,18 @@ def _integrate_scipy_stepwise(
     if solver_type is None:
         raise ValueError(f"Unsupported integration method: {config.method}")
     implicit_jac = jac if config.method in {"BDF", "Radau"} else None
+    solver_jac = implicit_jac
+    if implicit_jac is not None and activity_tracker is not None:
+        original_jacobian = implicit_jac
+
+        def activity_jacobian(time, state):
+            activity_tracker.record_jacobian(float(time))
+            try:
+                return original_jacobian(time, state)
+            finally:
+                activity_tracker.record_phase("solver_step", float(time))
+
+        solver_jac = activity_jacobian
 
     times = [float(config.t_start)]
     states = [[float(value)] for value in initial_state]
@@ -886,8 +1028,8 @@ def _integrate_scipy_stepwise(
                 "max_step": segment_max_step,
             }
             if config.method in {"BDF", "Radau"}:
-                if implicit_jac is not None:
-                    solver_options["jac"] = implicit_jac
+                if solver_jac is not None:
+                    solver_options["jac"] = solver_jac
                 elif jac_sparsity is not None:
                     solver_options["jac_sparsity"] = jac_sparsity
             requested_first_step = (
@@ -909,6 +1051,11 @@ def _integrate_scipy_stepwise(
                 start_segment = getattr(implicit_jac, "start_segment", None)
                 if start_segment is not None:
                     start_segment()
+                if activity_tracker is not None:
+                    activity_tracker.record_phase(
+                        "solver_initialization",
+                        last_accepted_time,
+                    )
                 solver = solver_type(
                     cancellable_rhs,
                     last_accepted_time,
@@ -980,6 +1127,10 @@ def _integrate_scipy_stepwise(
                             integration_end - last_accepted_time
                         ),
                     )
+                    if activity_tracker is not None:
+                        activity_tracker.record_solver_step(
+                            last_accepted_time
+                        )
                     step_message = solver.step()
                 except IntegrationCancelled:
                     status = "cancelled"
@@ -1365,6 +1516,7 @@ def integrate_ode(
     jac_sparsity=None,
     jac: JacobianCallable | None = None,
     recoverable_trial_retries: bool = False,
+    activity_tracker: SolverActivityTracker | None = None,
 ):
     """Integrate an ODE, optionally restarting at equation discontinuities.
 
@@ -1381,8 +1533,34 @@ def integrate_ode(
     into the stepwise path so a ``RecoverableTrialStateError`` can rebuild the
     solver from its last accepted state.  It defaults to false to preserve the
     direct ``solve_ivp`` path for ordinary callers.
+
+    ``activity_tracker`` is optional and additive. When omitted, the numerical
+    call path and callback behavior are unchanged.
     """
 
+    integration_rhs = rhs
+    integration_accepted_step_callback = accepted_step_callback
+    if activity_tracker is not None:
+        activity_tracker.start_integration(config.t_start)
+        original_rhs = rhs
+
+        def activity_rhs(time, state):
+            numeric_time = float(time)
+            activity_tracker.record_rhs(numeric_time)
+            try:
+                return original_rhs(time, state)
+            finally:
+                activity_tracker.record_phase("solver_step", numeric_time)
+
+        integration_rhs = activity_rhs
+
+        def activity_accepted_step(time: float) -> None:
+            activity_tracker.record_accepted_step(float(time))
+            if accepted_step_callback is not None:
+                accepted_step_callback(float(time))
+
+        integration_accepted_step_callback = activity_accepted_step
+
     if (
         state_transition_handler is not None
         and config.t_stop < config.t_start
@@ -1406,22 +1584,22 @@ def integrate_ode(
     except ImportError:
         if normalized_breakpoints:
             return _runge_kutta_4_segmented(
-                rhs,
+                integration_rhs,
                 initial_state,
                 config,
                 t_eval,
                 normalized_breakpoints,
                 cancel_check,
-                accepted_step_callback,
+                integration_accepted_step_callback,
                 state_transition_handler,
             )
         return _runge_kutta_4(
-            rhs,
+            integration_rhs,
             initial_state,
             config,
             t_eval,
             cancel_check,
-            accepted_step_callback,
+            integration_accepted_step_callback,
             state_transition_handler,
         )
 
@@ -1432,25 +1610,26 @@ def integrate_ode(
         or recoverable_trial_retries
     ):
         return _integrate_scipy_stepwise(
-            rhs,
+            integration_rhs,
             initial_state,
             config,
             t_eval,
             cancel_check or (lambda: False),
-            accepted_step_callback,
+            integration_accepted_step_callback,
             normalized_breakpoints,
             state_transition_handler,
             jac_sparsity,
             jac,
+            activity_tracker,
         )
 
     implicit_jac = jac if config.method in {"BDF", "Radau"} else None
-    solve_rhs = rhs
+    solve_rhs = integration_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)
+                derivative = integration_rhs(time, state)
                 observer(float(time), state, derivative)
                 return derivative
 
@@ -1459,6 +1638,19 @@ def integrate_ode(
         if start_segment is not None:
             start_segment()
 
+    solve_jac = implicit_jac
+    if implicit_jac is not None and activity_tracker is not None:
+        original_jacobian = implicit_jac
+
+        def activity_jacobian(time, state):
+            activity_tracker.record_jacobian(float(time))
+            try:
+                return original_jacobian(time, state)
+            finally:
+                activity_tracker.record_phase("solver_step", float(time))
+
+        solve_jac = activity_jacobian
+
     solve_options = {
         "fun": solve_rhs,
         "t_span": (config.t_start, config.t_stop),
@@ -1471,8 +1663,11 @@ def integrate_ode(
     }
     if config.first_step is not None:
         solve_options["first_step"] = config.first_step
-    if implicit_jac is not None:
-        solve_options["jac"] = implicit_jac
+    if solve_jac is not None:
+        solve_options["jac"] = solve_jac
     elif jac_sparsity is not None and config.method in {"BDF", "Radau"}:
         solve_options["jac_sparsity"] = jac_sparsity
-    return solve_ivp(**solve_options)
+    direct_solution = solve_ivp(**solve_options)
+    if activity_tracker is not None and len(direct_solution.t):
+        activity_tracker.record_accepted_step(float(direct_solution.t[-1]))
+    return direct_solution

app/simulation/systems/generic.py

@@ -31,6 +31,7 @@ from app.simulation.solvers.solver import (
     IntegrationCancelled,
     ODESolution,
     SolveIVPConfig,
+    SolverActivityTracker,
     integrate_ode,
 )
 from app.simulation.solvers.signal import SignalResolver
@@ -453,6 +454,7 @@ class GenericFluidSystem:
         self.pneumatic_volume_propagation_count = 0
         self._jacobian_sparsity = None
         self._ode_tangent_provider: ThreePistonTangentProvider | None = None
+        self._activity_tracker: SolverActivityTracker | None = None
 
     def _request_causal_residual_audit(self) -> None:
         """Make topology or mode boundaries verify the next causal closure."""
@@ -957,6 +959,8 @@ class GenericFluidSystem:
         *,
         record_rhs_outcome: bool = False,
     ) -> dict[str, dict[str, float]]:
+        if self._activity_tracker is not None:
+            self._activity_tracker.record_thermofluid_closure(time)
         transaction = self._thermofluid_transaction_plan.capture()
         last_algebraic_diagnostics = self._last_algebraic_diagnostics
         last_algebraic_scope = self._last_algebraic_scope
@@ -1156,6 +1160,8 @@ class GenericFluidSystem:
         self,
         connected_h: dict[str, dict[str, float]],
     ) -> list[float]:
+        if self._activity_tracker is not None:
+            self._activity_tracker.record_phase("state_derivatives")
         return self.pneumatic_storage_reducer.coupled_derivatives(
             self.mechanical_state_reducer.state_derivatives(connected_h)
         )
@@ -1194,7 +1200,32 @@ class GenericFluidSystem:
         sample_step: float,
         progress_callback: SimulationProgressCallback | None = None,
         cancel_check: SimulationCancellationCheck | None = None,
+        activity_tracker: SolverActivityTracker | None = None,
     ) -> GenericSimulationResult:
+        previous_activity_tracker = self._activity_tracker
+        self._activity_tracker = activity_tracker
+        try:
+            return self._simulate(
+                config,
+                sample_step=sample_step,
+                progress_callback=progress_callback,
+                cancel_check=cancel_check,
+                activity_tracker=activity_tracker,
+            )
+        finally:
+            self._activity_tracker = previous_activity_tracker
+
+    def _simulate(
+        self,
+        config: SolveIVPConfig,
+        *,
+        sample_step: float,
+        progress_callback: SimulationProgressCallback | None = None,
+        cancel_check: SimulationCancellationCheck | None = None,
+        activity_tracker: SolverActivityTracker | None = None,
+    ) -> GenericSimulationResult:
+        if activity_tracker is not None:
+            activity_tracker.record_phase("initializing", config.t_start)
         last_reported_progress = -1.0
         last_reported_phase = ""
 
@@ -1223,12 +1254,21 @@ class GenericFluidSystem:
         report_progress(0.0, "initializing", force=True)
         with performance_span("simulation.sample_initialization"):
             integration_config = config
+            mechanical_tolerance_plan = None
             if isinstance(config.atol, (int, float)):
+                mechanical_tolerance_plan = (
+                    self.mechanical_state_reducer.absolute_tolerance_plan(
+                        float(config.atol),
+                        mode=(
+                            None
+                            if config.method == "BDF"
+                            else "legacy"
+                        ),
+                    )
+                )
                 integration_config = replace(
                     config,
-                    atol=self.mechanical_state_reducer.absolute_tolerances(
-                        float(config.atol)
-                    ),
+                    atol=list(mechanical_tolerance_plan.values),
                 )
             t_eval = simulation_sample_times(config, sample_step)
             signal_event_times = self.signal_resolver.event_times(
@@ -1321,6 +1361,8 @@ class GenericFluidSystem:
                     def evaluate_jacobian_rhs(time, state):
                         if cancel_check is not None and cancel_check():
                             raise IntegrationCancelled
+                        if activity_tracker is not None:
+                            activity_tracker.record_rhs(float(time))
                         return monitored_rhs(
                             time,
                             [float(value) for value in state],
@@ -1372,6 +1414,7 @@ class GenericFluidSystem:
                 jac_sparsity=jac_sparsity,
                 jac=jacobian,
                 recoverable_trial_retries=True,
+                activity_tracker=activity_tracker,
             )
         finally:
             self._ode_tangent_provider = None
@@ -1388,6 +1431,11 @@ class GenericFluidSystem:
             else max(0.0, last_reported_progress)
         )
         report_progress(postprocess_progress, "postprocessing", force=True)
+        if activity_tracker is not None:
+            activity_tracker.record_phase(
+                "postprocessing",
+                furthest_solver_time,
+            )
         times = [float(value) for value in solution.t]
         if isinstance(solution, ODESolution):
             solver_segment_diagnostics = [
@@ -1601,6 +1649,11 @@ class GenericFluidSystem:
                     self._close_current_state(times[time_index])
                     self._append_current_state(series)
                     series["time"].append(times[time_index])
+                    if activity_tracker is not None:
+                        activity_tracker.record_phase(
+                            "postprocessing",
+                            times[time_index],
+                        )
                 except Exception as exc:
                     run_status = "failed"
                     result_message = str(exc)
@@ -1621,6 +1674,14 @@ class GenericFluidSystem:
             diagnostics = {
                 "integration": {
                     "method": integration_config.method,
+                    "mechanicalAbsoluteTolerance": (
+                        mechanical_tolerance_plan.as_dict()
+                        if mechanical_tolerance_plan is not None
+                        else {
+                            "mode": "callerVector",
+                            "stateCount": len(initial_state),
+                        }
+                    ),
                     "jacobianSparsity": jacobian_diagnostics,
                     "jacobian": runtime_jacobian_diagnostics,
                     "segmentCount": len(solver_segment_diagnostics),
@@ -1730,6 +1791,19 @@ class GenericFluidSystem:
                 for variable in self.network.result_variable_metadata()
                 if variable.key in series
             )
+            final_simulated_time = (
+                float(series["time"][-1])
+                if series["time"]
+                else float(config.t_start)
+            )
+            if activity_tracker is not None:
+                activity_tracker.record_phase(
+                    "complete" if run_status == "completed" else run_status,
+                    final_simulated_time,
+                )
+                diagnostics["activity"] = (
+                    activity_tracker.snapshot().as_dict()
+                )
             report_progress(
                 1.0 if run_status == "completed" else max(0.0, last_reported_progress),
                 "complete" if run_status == "completed" else run_status,
@@ -1740,9 +1814,7 @@ class GenericFluidSystem:
                 status=run_status,
                 message=result_message,
                 simulated_until=(
-                    float(series["time"][-1])
-                    if series["time"]
-                    else float(config.t_start)
+                    final_simulated_time
                 ),
                 requested_stop_time=float(config.t_stop),
                 variables=variables,