优化压力流量求解并达到四路性能门槛

app/simulation/solvers/algebraic.py

@@ -54,6 +54,23 @@ class ExplicitFlowAssignment:
 @dataclass(frozen=True)
 class ExplicitFlowStage:
     assignments: tuple[ExplicitFlowAssignment, ...]
+    direct_evaluations: tuple[tuple[int, Callable[[], float]], ...]
+    component_evaluations: tuple["ExplicitFlowComponentEvaluation", ...]
+
+
+@dataclass(frozen=True)
+class ExplicitFlowComponentEvaluation:
+    component: object
+    evaluate: Callable[[], tuple[float, ...]]
+    assignment_indices: tuple[int, ...]
+    equation_indices: tuple[int, ...]
+    equation_ids: tuple[str, ...]
+
+
+@dataclass(frozen=True)
+class EquationScalePlan:
+    variable_names: tuple[str, ...]
+    force_unknowns: tuple[AlgebraicUnknown, ...]
 
 
 @dataclass(frozen=True)
@@ -68,6 +85,13 @@ class ConnectionEquationEvaluation:
     evaluate: Callable[[], float]
 
 
+@dataclass(frozen=True)
+class ComponentEquationEvaluation:
+    component: object
+    evaluate: Callable[[], tuple[float, ...]]
+    templates: tuple[EquationResidual, ...]
+
+
 @dataclass(frozen=True)
 class PnorPnl0001SeriesBinding:
     orifice: AmesimPnor001
@@ -147,6 +171,18 @@ class PressureFlowSolver:
             for variable in ("p", "m_flow", "x", "v", "f")
         }
         self._component_equation_owners = tuple(network.components.values())
+        self._component_equation_plan = tuple(
+            ComponentEquationEvaluation(
+                component=component,
+                evaluate=component.pressure_flow_equation_values,
+                templates=component.pressure_flow_equation_residuals(),
+            )
+            for component in self._component_equation_owners
+        )
+        self._component_equation_plans_by_id = {
+            item.component.name: item
+            for item in self._component_equation_plan
+        }
         self._estimated_flow_components = tuple(
             component
             for component in self._component_equation_owners
@@ -164,6 +200,11 @@ class PressureFlowSolver:
             )
             for equation in network.connection_equation_residuals()
         )
+        self._equation_templates = tuple(
+            equation
+            for item in self._component_equation_plan
+            for equation in item.templates
+        ) + tuple(item.template for item in self._connection_equation_plan)
         self._effort_groups = {
             variable: self._build_effort_equality_groups(variable)
             for variable in ("p", "x", "v")
@@ -176,12 +217,19 @@ class PressureFlowSolver:
         )
         self._unilateral_contact_plan = self._build_unilateral_contact_plan()
         self._explicit_flow_plan = self._build_explicit_flow_plan()
+        self._explicit_flow_plans_by_variables = {
+            selected: self._filter_explicit_flow_plan(selected)
+            for selected in (frozenset(("f", "m_flow")), frozenset(("f",)))
+        }
+        self._equation_scale_plans = self._build_equation_scale_plans(
+            self._equation_templates
+        )
         self.last_diagnostics: AlgebraicSolveDiagnostics | None = None
 
     def _build_unknowns(self) -> tuple[AlgebraicUnknown, ...]:
         unknowns: list[AlgebraicUnknown] = []
         for component in self.network.components.values():
-            for definition in component.port_definitions:
+            for definition in component.active_port_definitions:
                 if definition.kind != "physical":
                     continue
                 state = component.get_port(definition.name)
@@ -277,8 +325,8 @@ class PressureFlowSolver:
                 union(first, second)
 
         component_equations = {
-            component.name: component.pressure_flow_equation_residuals()
-            for component in self.network.components.values()
+            item.component.name: item.templates
+            for item in self._component_equation_plan
         }
         for equations in component_equations.values():
             for equation in equations:
@@ -720,12 +768,10 @@ class PressureFlowSolver:
                 f"Unsupported connection equation relation: {equation.relation}."
             )
 
-        component = self.network.components[equation.owner_id]
+        evaluation = self._component_equation_plans_by_id[equation.owner_id]
+        component = evaluation.component
         equation_id = equation.id
-        equation_ids = tuple(
-            current.id
-            for current in component.pressure_flow_equation_residuals()
-        )
+        equation_ids = tuple(current.id for current in evaluation.templates)
         try:
             equation_index = equation_ids.index(equation_id)
         except ValueError as exc:
@@ -734,38 +780,49 @@ class PressureFlowSolver:
             ) from exc
 
         def read_component_equation() -> float:
-            current_equations = component.pressure_flow_equation_residuals()
-            if (
-                equation_index >= len(current_equations)
-                or current_equations[equation_index].id != equation_id
-            ):
+            current_values = evaluation.evaluate()
+            if equation_index >= len(current_values):
                 raise RuntimeError(
-                    f"Compiled algebraic equation disappeared at runtime: {equation_id}."
+                    "Compiled algebraic equation disappeared at runtime: "
+                    f"{equation_id}."
                 )
-            return float(current_equations[equation_index].value)
+            return float(current_values[equation_index])
 
         return read_component_equation
-    def _pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
-        """Evaluate live values through a precompiled connector topology."""
-        component_residuals = tuple(
-            residual
-            for component in self._component_equation_owners
-            for residual in component.pressure_flow_equation_residuals()
-        )
-        connection_residuals = tuple(
-            EquationResidual(
-                id=item.template.id,
-                owner=item.template.owner,
-                owner_id=item.template.owner_id,
-                relation=item.template.relation,
-                variables=item.template.variables,
-                value=item.evaluate(),
-                role=item.template.role,
-            )
-            for item in self._connection_equation_plan
-        )
-        return component_residuals + connection_residuals
 
+    def _pressure_flow_equation_values(self) -> tuple[float, ...]:
+        """Evaluate live equation values through the compiled topology."""
+
+        values: list[float] = []
+        for item in self._component_equation_plan:
+            current_values = item.evaluate()
+            if len(current_values) != len(item.templates):
+                raise RuntimeError(
+                    "Compiled algebraic equation count changed at runtime for "
+                    f"{item.component.name}."
+                )
+            values.extend(current_values)
+        values.extend(item.evaluate() for item in self._connection_equation_plan)
+        return tuple(values)
+
+    def _pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
+        """Materialize public residual objects only when explicitly requested."""
+
+        return tuple(
+            EquationResidual(
+                id=template.id,
+                owner=template.owner,
+                owner_id=template.owner_id,
+                relation=template.relation,
+                variables=template.variables,
+                value=value,
+                role=template.role,
+            )
+            for template, value in zip(
+                self._equation_templates,
+                self._pressure_flow_equation_values(),
+            )
+        )
 
     def _explicit_flow_assignment(
         self,
@@ -779,8 +836,9 @@ class PressureFlowSolver:
                 evaluate=self._equation_value_reader(equation),
             )
 
-        component = self.network.components[equation.owner_id]
-        equations = component.pressure_flow_equation_residuals()
+        evaluation = self._component_equation_plans_by_id[equation.owner_id]
+        component = evaluation.component
+        equations = evaluation.templates
         equation_ids = tuple(current.id for current in equations)
         try:
             equation_index = equation_ids.index(equation.id)
@@ -803,8 +861,8 @@ class PressureFlowSolver:
         seeded_ids: set[str] = set()
 
         initial_assignments: list[ExplicitFlowAssignment] = []
-        for component in self.network.components.values():
-            for equation in component.pressure_flow_equation_residuals():
+        for evaluation in self._component_equation_plan:
+            for equation in evaluation.templates:
                 if equation.relation != "constitutive" or equation.role != "flow":
                     continue
                 flow_unknowns = self._flow_unknowns_for_equation(equation)
@@ -818,9 +876,9 @@ class PressureFlowSolver:
                 )
                 seeded_ids.add(unknown.id)
         if initial_assignments:
-            stages.append(ExplicitFlowStage(tuple(initial_assignments)))
+            stages.append(self._compile_explicit_flow_stage(tuple(initial_assignments)))
 
-        equations = self._pressure_flow_equation_residuals()
+        equations = self._equation_templates
         while True:
             stage_assignments: list[ExplicitFlowAssignment] = []
             stage_unknown_ids: set[str] = set()
@@ -850,7 +908,9 @@ class PressureFlowSolver:
                 )
                 stage_unknown_ids.add(unknown.id)
             if stage_assignments:
-                stages.append(ExplicitFlowStage(tuple(stage_assignments)))
+                stages.append(
+                    self._compile_explicit_flow_stage(tuple(stage_assignments))
+                )
                 seeded_ids.update(stage_unknown_ids)
                 continue
 
@@ -880,43 +940,105 @@ class PressureFlowSolver:
                 break
             if fallback_assignment is None:
                 break
-            stages.append(ExplicitFlowStage((fallback_assignment,)))
+            stages.append(self._compile_explicit_flow_stage((fallback_assignment,)))
 
         return tuple(stages)
 
+    @staticmethod
+    def _compile_explicit_flow_stage(
+        assignments: tuple[ExplicitFlowAssignment, ...],
+    ) -> ExplicitFlowStage:
+        direct_evaluations: list[tuple[int, Callable[[], float]]] = []
+        assignments_by_component: dict[
+            object,
+            list[tuple[int, ExplicitFlowAssignment]],
+        ] = {}
+        for assignment_index, assignment in enumerate(assignments):
+            if assignment.component is None:
+                if assignment.evaluate is None:
+                    raise RuntimeError(
+                        "Explicit flow assignment has no compiled evaluator: "
+                        f"{assignment.equation_id}."
+                    )
+                direct_evaluations.append((assignment_index, assignment.evaluate))
+                continue
+            assignments_by_component.setdefault(assignment.component, []).append(
+                (assignment_index, assignment)
+            )
+
+        component_evaluations: list[ExplicitFlowComponentEvaluation] = []
+        for component, component_assignments in assignments_by_component.items():
+            if any(
+                assignment.equation_index is None
+                for _assignment_index, assignment in component_assignments
+            ):
+                raise RuntimeError(
+                    "Explicit component flow assignment has no equation index."
+                )
+            component_evaluations.append(
+                ExplicitFlowComponentEvaluation(
+                    component=component,
+                    evaluate=component.pressure_flow_equation_values,
+                    assignment_indices=tuple(
+                        assignment_index
+                        for assignment_index, _assignment in component_assignments
+                    ),
+                    equation_indices=tuple(
+                        int(assignment.equation_index)
+                        for _assignment_index, assignment in component_assignments
+                    ),
+                    equation_ids=tuple(
+                        assignment.equation_id
+                        for _assignment_index, assignment in component_assignments
+                    ),
+                )
+            )
+        return ExplicitFlowStage(
+            assignments=assignments,
+            direct_evaluations=tuple(direct_evaluations),
+            component_evaluations=tuple(component_evaluations),
+        )
+
+    def _filter_explicit_flow_plan(
+        self,
+        selected: frozenset[str],
+    ) -> tuple[ExplicitFlowStage, ...]:
+        return tuple(
+            self._compile_explicit_flow_stage(
+                tuple(
+                    assignment
+                    for assignment in stage.assignments
+                    if assignment.unknown.variable in selected
+                )
+            )
+            for stage in self._explicit_flow_plan
+        )
 
     @staticmethod
     def _evaluate_explicit_flow_stage(
-        assignments: tuple[ExplicitFlowAssignment, ...],
-    ) -> dict[str, float]:
-        values: dict[str, float] = {}
-        assignments_by_component: dict[object, list[ExplicitFlowAssignment]] = {}
-        for assignment in assignments:
-            if assignment.component is None:
-                assert assignment.evaluate is not None
-                values[assignment.equation_id] = assignment.evaluate()
-                continue
-            assignments_by_component.setdefault(assignment.component, []).append(
-                assignment
-            )
+        stage: ExplicitFlowStage,
+    ) -> tuple[float, ...]:
+        assignments = stage.assignments
+        values: list[float | None] = [None] * len(assignments)
+        for assignment_index, evaluate in stage.direct_evaluations:
+            values[assignment_index] = evaluate()
 
-        for component, component_assignments in assignments_by_component.items():
-            equations = component.pressure_flow_equation_residuals()
-            for assignment in component_assignments:
-                assert assignment.equation_index is not None
-                equation_index = assignment.equation_index
-                if (
-                    equation_index >= len(equations)
-                    or equations[equation_index].id != assignment.equation_id
-                ):
+        for evaluation in stage.component_evaluations:
+            equation_values = evaluation.evaluate()
+            for assignment_index, equation_index, equation_id in zip(
+                evaluation.assignment_indices,
+                evaluation.equation_indices,
+                evaluation.equation_ids,
+            ):
+                if equation_index >= len(equation_values):
                     raise RuntimeError(
                         "Compiled algebraic equation disappeared at runtime: "
-                        f"{assignment.equation_id}."
+                        f"{equation_id}."
                     )
-                values[assignment.equation_id] = float(
-                    equations[equation_index].value
-                )
-        return values
+                values[assignment_index] = float(equation_values[equation_index])
+        if any(value is None for value in values):
+            raise RuntimeError("Explicit flow evaluation plan returned no value.")
+        return tuple(float(value) for value in values)
 
     def _solve_explicit_flow_unknowns(
         self,
@@ -930,19 +1052,18 @@ class PressureFlowSolver:
                 unknown.write(0.0)
 
         seeded_ids: set[str] = set()
-        for stage in self._explicit_flow_plan:
-            assignments = tuple(
-                assignment
-                for assignment in stage.assignments
-                if assignment.unknown.variable in selected
-            )
-            values = self._evaluate_explicit_flow_stage(assignments)
+        plan = self._explicit_flow_plans_by_variables.get(selected)
+        if plan is None:
+            plan = self._filter_explicit_flow_plan(selected)
+        for stage in plan:
+            assignments = stage.assignments
+            values = self._evaluate_explicit_flow_stage(stage)
             targets = tuple(
                 (
                     assignment,
-                    assignment.unknown.read() - values[assignment.equation_id],
+                    assignment.unknown.read() - value,
                 )
-                for assignment in assignments
+                for assignment, value in zip(assignments, values)
             )
             for assignment, target_value in targets:
                 if not isfinite(target_value):
@@ -1101,12 +1222,38 @@ class PressureFlowSolver:
             orifice.get_port(orifice_port).p = pressure
             pipe.get_port(pipe_port).p = pressure
 
+    def _equation_scale_plan(self, equation) -> EquationScalePlan:
+        return EquationScalePlan(
+            variable_names=tuple(
+                variable.rsplit(".", 1)[-1]
+                for variable in equation.variables
+            ),
+            force_unknowns=tuple(
+                unknown
+                for variable in equation.variables
+                if (unknown := self._unknowns_by_id.get(variable)) is not None
+                and unknown.variable == "f"
+            ),
+        )
+
+    def _build_equation_scale_plans(
+        self,
+        equations: tuple[EquationResidual, ...],
+    ) -> tuple[EquationScalePlan, ...]:
+        return tuple(
+            self._equation_scale_plan(equation)
+            for equation in equations
+        )
+
     def _scales(self) -> dict[str, float]:
+        pressure_values = [
+            unknown.read() for unknown in self._unknowns_by_variable["p"]
+        ]
         pressure_scale = max(
             [
-                abs(unknown.read())
-                for unknown in self._unknowns_by_variable["p"]
-                if unknown.read() > 0.0
+                abs(value)
+                for value in pressure_values
+                if value > 0.0
             ]
             + [1e5]
         )
@@ -1166,36 +1313,24 @@ class PressureFlowSolver:
         scales = self._scales()
         pressure_scale = scales["p"]
         flow_scale = scales["m_flow"]
-        unknown_scales = {
-            unknown.id: (
-                max(abs(unknown.read()), 1.0)
-                if unknown.variable == "f"
-                else scales.get(unknown.variable, max(abs(unknown.read()), 1.0))
-            )
-            for unknown in self.unknowns
-        }
-        def variable_scale(unknown: AlgebraicUnknown) -> float:
-            return unknown_scales[unknown.id]
+        seeded_values = self._pressure_flow_equation_values()
 
-        seeded_equations = self._pressure_flow_equation_residuals()
-
-        def initial_equation_scale(equation) -> float:
-            variable_names = [
-                variable.rsplit(".", 1)[-1]
-                for variable in equation.variables
-            ]
+        def initial_equation_scale(
+            equation: EquationResidual,
+            value: float,
+            scale_plan: EquationScalePlan,
+        ) -> float:
+            variable_names = scale_plan.variable_names
             if equation.role == "flow":
                 force_scales = [
-                    unknown_scales[variable]
-                    for variable in equation.variables
-                    if variable in self._unknowns_by_id
-                    and self._unknowns_by_id[variable].variable == "f"
+                    max(abs(unknown.read()), 1.0)
+                    for unknown in scale_plan.force_unknowns
                 ]
                 if force_scales:
                     # Freeze force scaling per equation.  A 1e17 N source must
                     # not hide an unrelated 40 N piston/contact imbalance in a
                     # different mechanical branch.
-                    return max(force_scales + [abs(float(equation.value)), 1.0])
+                    return max(force_scales + [abs(value), 1.0])
                 return flow_scale
             if equation.role == "effort":
                 if "x" in variable_names:
@@ -1205,20 +1340,18 @@ class PressureFlowSolver:
                 return pressure_scale
             return max([scales.get(name, 1.0) for name in variable_names] + [1.0])
 
-        equation_scales = {
-            equation.id: initial_equation_scale(equation)
-            for equation in seeded_equations
-        }
-
-        def equation_scale(equation) -> float:
-            cached = equation_scales.get(equation.id)
-            if cached is not None:
-                return cached
-            return initial_equation_scale(equation)
+        equation_scales = tuple(
+            initial_equation_scale(equation, value, scale_plan)
+            for equation, value, scale_plan in zip(
+                self._equation_templates,
+                seeded_values,
+                self._equation_scale_plans,
+            )
+        )
 
         seeded_scaled = [
-            abs(equation.value / equation_scale(equation))
-            for equation in seeded_equations
+            abs(value / scale)
+            for value, scale in zip(seeded_values, equation_scales)
         ]
         seeded_max_scaled_residual = max(seeded_scaled, default=0.0)
         seeded_unknown_values = [
@@ -1242,13 +1375,25 @@ class PressureFlowSolver:
                 flow_scale=flow_scale,
                 max_scaled_residual=seeded_max_scaled_residual,
                 max_raw_residual=max(
-                    (abs(item.value) for item in seeded_equations),
+                    (abs(value) for value in seeded_values),
                     default=0.0,
                 ),
             )
             self.last_diagnostics = diagnostics
             return diagnostics
 
+        unknown_scales = {
+            unknown.id: (
+                max(abs(unknown.read()), 1.0)
+                if unknown.variable == "f"
+                else scales[unknown.variable]
+            )
+            for unknown in self.unknowns
+        }
+
+        def variable_scale(unknown: AlgebraicUnknown) -> float:
+            return unknown_scales[unknown.id]
+
         positive_pressures = [
             unknown.read()
             for unknown in self._unknowns_by_variable["p"]
@@ -1293,11 +1438,11 @@ class PressureFlowSolver:
         def scaled_residuals(values):
             assign(values)
             self._refresh_unilateral_contacts(contact_bindings)
-            equations = self._pressure_flow_equation_residuals()
+            equation_values = self._pressure_flow_equation_values()
             return np.asarray(
                 [
-                    equation.value / equation_scale(equation)
-                    for equation in equations
+                    value / scale
+                    for value, scale in zip(equation_values, equation_scales)
                 ],
                 dtype=float,
             )
@@ -1314,12 +1459,10 @@ class PressureFlowSolver:
         )
         assign(result.x)
         self._refresh_unilateral_contacts(contact_bindings)
-        equations = self._pressure_flow_equation_residuals()
+        equation_values = self._pressure_flow_equation_values()
         scaled = [
-            abs(
-                equation.value / equation_scale(equation)
-            )
-            for equation in equations
+            abs(value / scale)
+            for value, scale in zip(equation_values, equation_scales)
         ]
         max_scaled_residual = max(scaled, default=0.0)
         residuals_converged = (
@@ -1335,7 +1478,10 @@ class PressureFlowSolver:
             pressure_scale=pressure_scale,
             flow_scale=flow_scale,
             max_scaled_residual=max_scaled_residual,
-            max_raw_residual=max((abs(item.value) for item in equations), default=0.0),
+            max_raw_residual=max(
+                (abs(value) for value in equation_values),
+                default=0.0,
+            ),
         )
         self.last_diagnostics = diagnostics
         if not success:

app/simulation/solvers/mechanical.py

@@ -279,7 +279,7 @@ class MechanicalStateReducer:
         mechanical_ports = {
             (component.name, definition.name)
             for component in self.network.components.values()
-            for definition in component.port_definitions
+            for definition in component.active_port_definitions
             if definition.kind == "physical" and definition.domain == "mechanical"
         }
         parents = {
@@ -339,7 +339,7 @@ class MechanicalStateReducer:
         for component in masses:
             ports = [
                 (component.name, definition.name)
-                for definition in component.port_definitions
+                for definition in component.active_port_definitions
                 if definition.kind == "physical"
                 and definition.domain == "mechanical"
             ]
@@ -354,7 +354,7 @@ class MechanicalStateReducer:
         for component in masses:
             first_port = next(
                 (component.name, definition.name)
-                for definition in component.port_definitions
+                for definition in component.active_port_definitions
                 if definition.kind == "physical"
                 and definition.domain == "mechanical"
             )

app/simulation/solvers/pneumatic_storage.py

@@ -99,7 +99,7 @@ def _pressure_storage_endpoint_groups(
     pneumatic_endpoints = {
         Endpoint(component.name, definition.name)
         for component in network.components.values()
-        for definition in component.port_definitions
+        for definition in component.active_port_definitions
         if definition.kind == "physical" and definition.domain == "pneumatic"
     }
     parent = {endpoint: endpoint for endpoint in pneumatic_endpoints}

app/simulation/solvers/pneumatic_volume.py

@@ -38,7 +38,7 @@ class PneumaticVolumeResolver:
 
     def solve(self) -> PneumaticVolumeDiagnostics:
         for component in self.network.components.values():
-            for definition in component.port_definitions:
+            for definition in component.active_port_definitions:
                 if definition.kind == "physical" and definition.domain == "pneumatic":
                     port = component.get_port(definition.name)
                     port.volume = 0.0

app/simulation/solvers/solver.py

@@ -44,6 +44,36 @@ class SolveIVPConfig:
     first_step: float | None = None
 
 
+@dataclass(frozen=True)
+class SolverSegmentDiagnostics:
+    """Work performed by implicit solver instances inside one event segment."""
+
+    start_time: float
+    requested_stop_time: float
+    simulated_until: float
+    nfev: int = 0
+    njev: int = 0
+    nlu: int = 0
+    accepted_step_count: int = 0
+    solver_start_count: int = 0
+    state_transition_count: int = 0
+    recoverable_retry_count: int = 0
+
+    def as_dict(self) -> dict[str, float | int]:
+        return {
+            "startTime": self.start_time,
+            "requestedStopTime": self.requested_stop_time,
+            "simulatedUntil": self.simulated_until,
+            "nfev": self.nfev,
+            "njev": self.njev,
+            "nlu": self.nlu,
+            "acceptedStepCount": self.accepted_step_count,
+            "solverStartCount": self.solver_start_count,
+            "stateTransitionCount": self.state_transition_count,
+            "recoverableRetryCount": self.recoverable_retry_count,
+        }
+
+
 @dataclass(frozen=True)
 class ODESolution:
     t: list[float]
@@ -52,6 +82,7 @@ class ODESolution:
     message: str
     status: IntegrationStatus = "completed"
     error: Exception | None = None
+    solver_segments: tuple[SolverSegmentDiagnostics, ...] = ()
 
 
 def _vector_add(a: list[float], b: list[float], scale: float = 1.0) -> list[float]:
@@ -610,6 +641,7 @@ def _integrate_scipy_stepwise(
     same_time_transition_count = 0
     integration_progressed = False
     last_reported_step: float | None = None
+    solver_segments: list[SolverSegmentDiagnostics] = []
 
     def cancellation_message() -> str:
         return (
@@ -639,9 +671,17 @@ def _integrate_scipy_stepwise(
             math.nextafter(segment_end, -math.inf) if is_breakpoint else segment_end
         )
         has_integration_interval = integration_end > last_accepted_time
+        segment_start_time = last_accepted_time
         segment_max_step = float(config.max_step)
         recoverable_retry_count = 0
         last_recoverable_error: RecoverableTrialStateError | None = None
+        segment_nfev = 0
+        segment_njev = 0
+        segment_nlu = 0
+        segment_accepted_steps = 0
+        segment_solver_starts = 0
+        segment_state_transitions = 0
+        segment_recoverable_retries = 0
 
         while has_integration_interval and last_accepted_time < integration_end:
             if cancel_check():
@@ -681,6 +721,7 @@ def _integrate_scipy_stepwise(
                 break
             except RecoverableTrialStateError as exc:
                 recoverable_retry_count += 1
+                segment_recoverable_retries += 1
                 last_recoverable_error = exc
                 next_step = 0.5 * segment_max_step
                 minimum_step = 64.0 * math.ulp(max(abs(last_accepted_time), 1.0))
@@ -697,6 +738,7 @@ def _integrate_scipy_stepwise(
                 error = exc
                 break
 
+            segment_solver_starts += 1
 
             restart_at_transition = False
             restart_after_recoverable = False
@@ -720,6 +762,7 @@ def _integrate_scipy_stepwise(
                     break
                 except RecoverableTrialStateError as exc:
                     recoverable_retry_count += 1
+                    segment_recoverable_retries += 1
                     last_recoverable_error = exc
                     attempted_step = segment_max_step
                     next_step = 0.5 * attempted_step
@@ -744,6 +787,7 @@ def _integrate_scipy_stepwise(
                 if solver.status == "failed":
                     if last_recoverable_error is not None:
                         recoverable_retry_count += 1
+                        segment_recoverable_retries += 1
                         next_step = 0.5 * segment_max_step
                         minimum_step = 64.0 * math.ulp(
                             max(abs(last_accepted_time), 1.0)
@@ -759,6 +803,7 @@ def _integrate_scipy_stepwise(
                     message = str(step_message or "Integration step failed.")
                     break
 
+                segment_accepted_steps += 1
                 step_end_time = float(solver.t)
                 step_end_state = [float(value) for value in solver.y]
                 dense_output = (
@@ -807,6 +852,7 @@ def _integrate_scipy_stepwise(
                         break
 
                 if transition is not None:
+                    segment_state_transitions += 1
                     try:
                         same_time_transition_count = (
                             _next_same_time_transition_count(
@@ -896,6 +942,9 @@ def _integrate_scipy_stepwise(
                     )
                 report_step(reported_time)
 
+            segment_nfev += int(getattr(solver, "nfev", 0))
+            segment_njev += int(getattr(solver, "njev", 0))
+            segment_nlu += int(getattr(solver, "nlu", 0))
             if status != "completed":
                 break
             if restart_after_recoverable:
@@ -903,6 +952,22 @@ def _integrate_scipy_stepwise(
             if not restart_at_transition:
                 break
 
+        solver_segments.append(
+            SolverSegmentDiagnostics(
+                start_time=float(segment_start_time),
+                requested_stop_time=float(segment_end),
+                simulated_until=float(
+                    segment_end if status == "completed" else last_accepted_time
+                ),
+                nfev=segment_nfev,
+                njev=segment_njev,
+                nlu=segment_nlu,
+                accepted_step_count=segment_accepted_steps,
+                solver_start_count=segment_solver_starts,
+                state_transition_count=segment_state_transitions,
+                recoverable_retry_count=segment_recoverable_retries,
+            )
+        )
         if status != "completed":
             break
 
@@ -949,6 +1014,7 @@ def _integrate_scipy_stepwise(
         message=message,
         status=status,
         error=error,
+        solver_segments=tuple(solver_segments),
     )