Add fine-grained per-block diagnostics for convergence debugging

milanofthe · milanofthe · commit b547c2f49c5e · 2026-03-21T07:51:06.000+01:00
diff --git a/src/pathsim/simulation.py b/src/pathsim/simulation.py
@@ -37,6 +37,7 @@
 from .utils.deprecation import deprecated
 from .utils.portreference import PortReference
 from .utils.progresstracker import ProgressTracker
+from .utils.diagnostics import Diagnostics
 from .utils.logger import LoggerManager
 
 from .solvers import SSPRK22, SteadyState
@@ -165,17 +166,18 @@ class Simulation:
     """
 
     def __init__(
-        self, 
-        blocks=None, 
-        connections=None, 
+        self,
+        blocks=None,
+        connections=None,
         events=None,
-        dt=SIM_TIMESTEP, 
-        dt_min=SIM_TIMESTEP_MIN, 
-        dt_max=SIM_TIMESTEP_MAX, 
-        Solver=SSPRK22, 
-        tolerance_fpi=SIM_TOLERANCE_FPI, 
-        iterations_max=SIM_ITERATIONS_MAX, 
+        dt=SIM_TIMESTEP,
+        dt_min=SIM_TIMESTEP_MIN,
+        dt_max=SIM_TIMESTEP_MAX,
+        Solver=SSPRK22,
+        tolerance_fpi=SIM_TOLERANCE_FPI,
+        iterations_max=SIM_ITERATIONS_MAX,
         log=LOG_ENABLE,
+        diagnostics=False,
         **solver_kwargs
         ):
 
@@ -226,6 +228,12 @@ def __init__(
         #flag for setting the simulation active
         self._active = True
 
+        #diagnostics snapshot (None when disabled)
+        self.diagnostics = Diagnostics() if diagnostics else None
+
+        #diagnostics history (list of snapshots per timestep)
+        self._diagnostics_history = [] if diagnostics == "history" else None
+
         #initialize logging 
         logger_mgr = LoggerManager(
             enabled=bool(self.log),
@@ -815,6 +823,12 @@ def reset(self, time=0.0):
         for event in self.events:
             event.reset()
 
+        #reset diagnostics
+        if self.diagnostics is not None:
+            self.diagnostics = Diagnostics()
+        if self._diagnostics_history is not None:
+            self._diagnostics_history.clear()
+
         #evaluate system function
         self._update(self.time)
 
@@ -1027,18 +1041,24 @@ def _loops(self, t):
 
             #step boosters of loop closing connections
             max_err = 0.0
+            self._loop_errors = {}
             for con_booster in self.boosters:
                 err = con_booster.update()
+                self._loop_errors[con_booster] = err
                 if err > max_err:
                     max_err = err
-                       
+
             #check convergence
             if max_err <= self.tolerance_fpi:
+                self._loop_iterations = iteration
                 return
 
-        #not converged -> error
-        _msg = "algebraic loop not converged (iters: {}, err: {})".format(
-            self.iterations_max, max_err
+        self._loop_iterations = self.iterations_max
+
+        #not converged -> error with per-connection details
+        details = [f"  {b.connection}: {e:.2e}" for b, e in self._loop_errors.items()]
+        _msg = "algebraic loop not converged (iters: {}, err: {:.2e})\n{}".format(
+            self.iterations_max, max_err, "\n".join(details)
             )
         self.logger.error(_msg)
         raise RuntimeError(_msg)
@@ -1075,31 +1095,36 @@ def _solve(self, t, dt):
         #total evaluations of system equation
         total_evals = 0
 
+        #per-block residuals (overwritten each iteration, only final kept)
+        self._solve_errors = {}
+
         #perform fixed-point iterations to solve implicit update equation
         for it in range(self.iterations_max):
 
             #evaluate system equation (this is a fixed point loop)
             self._update(t)
-            total_evals += 1            
+            total_evals += 1
 
             #advance solution of implicit solver
             max_error = 0.0
             for block in self._blocks_dyn:
 
                 #skip inactive blocks
-                if not block: 
+                if not block:
                     continue
-                
+
                 #advance solution (internal optimizer)
                 error = block.solve(t, dt)
+                self._solve_errors[block] = error
                 if error > max_error:
                     max_error = error
 
             #check for convergence (only error)
             if max_error <= self.tolerance_fpi:
+                self._solve_iterations = it + 1
                 return True, total_evals, it+1
 
-        #not converged in 'self.iterations_max' steps
+        self._solve_iterations = self.iterations_max
         return False, total_evals, self.iterations_max
 
 
@@ -1156,8 +1181,10 @@ def steadystate(self, reset=False):
 
         #catch non convergence
         if not success:
-            _msg = "STEADYSTATE -> FINISHED (success: {}, evals: {}, iters: {}, runtime: {})".format(
-                success, evals, iters, T)
+            details = [f"  {b.__class__.__name__}: {e:.2e}"
+                       for b, e in self._solve_errors.items()]
+            _msg = "STEADYSTATE -> FAILED (evals: {}, iters: {}, runtime: {})\n{}".format(
+                evals, iters, T, "\n".join(details))
             self.logger.error(_msg)
             raise RuntimeError(_msg)
 
@@ -1278,6 +1305,7 @@ def _step(self, t, dt):
 
         #initial timestep rescale and error estimate
         success, max_error_norm, min_scale = True, 0.0, None
+        self._step_errors = {}
 
         #step blocks and get error estimates if available
         for block in self._blocks_dyn:
@@ -1287,6 +1315,7 @@ def _step(self, t, dt):
 
             #step the block
             suc, err_norm, scl = block.step(t, dt)
+            self._step_errors[block] = (suc, err_norm, scl)
 
             #check solver stepping success
             if not suc:
@@ -1471,8 +1500,23 @@ def timestep(self, dt=None, adaptive=True):
 
                     #adaptive implicit: revert if solver didn't converge
                     if not success and is_adaptive:
+                        details = [f"  {b.__class__.__name__}: {e:.2e}"
+                                   for b, e in self._solve_errors.items()]
+                        self.logger.warning(
+                            "implicit solver not converged, reverting step at t={:.6f}\n{}".format(
+                                time_stage, "\n".join(details))
+                            )
                         self._revert(self.time)
                         return False, 0.0, 0.5, total_evals + 1, total_solver_its
+
+                    #fixed implicit: warn if solver didn't converge
+                    if not success and not is_adaptive:
+                        details = [f"  {b.__class__.__name__}: {e:.2e}"
+                                   for b, e in self._solve_errors.items()]
+                        self.logger.warning(
+                            "implicit solver not converged at t={:.6f} (iters: {})\n{}".format(
+                                time_stage, solver_its, "\n".join(details))
+                            )
                 else:
                     #explicit: evaluate system equation
                     self._update(time_stage)
@@ -1511,6 +1555,19 @@ def timestep(self, dt=None, adaptive=True):
                 self._update(time_dt)
                 total_evals += 1
 
+        #update diagnostics snapshot for this timestep
+        if self.diagnostics is not None:
+            self.diagnostics = Diagnostics(
+                time=time_dt,
+                loop_residuals=dict(getattr(self, '_loop_errors', {})),
+                loop_iterations=getattr(self, '_loop_iterations', 0),
+                solve_residuals=dict(getattr(self, '_solve_errors', {})),
+                solve_iterations=getattr(self, '_solve_iterations', 0),
+                step_errors=dict(getattr(self, '_step_errors', {})),
+            )
+            if self._diagnostics_history is not None:
+                self._diagnostics_history.append(self.diagnostics)
+
         #sample data after successful timestep
         self._sample(time_dt, dt)
 
diff --git a/src/pathsim/utils/diagnostics.py b/src/pathsim/utils/diagnostics.py
@@ -0,0 +1,123 @@
+#########################################################################################
+##
+##                          DIAGNOSTICS FOR SIMULATION
+##                            (utils/diagnostics.py)
+##
+##         Fine-grained per-block and per-connection convergence metrics
+##         for debugging non-converging simulations.
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+from dataclasses import dataclass, field
+
+
+# DIAGNOSTICS ===========================================================================
+
+@dataclass
+class Diagnostics:
+    """Per-timestep convergence diagnostics snapshot.
+
+    Populated by the simulation engine after each successful timestep.
+    Stores per-block and per-connection residuals from the three
+    convergence loops: algebraic loop solver, implicit ODE solver,
+    and adaptive error control.
+
+    Attributes
+    ----------
+    time : float
+        simulation time
+    loop_residuals : dict
+        per-booster algebraic loop residuals (booster -> residual)
+    loop_iterations : int
+        number of algebraic loop iterations taken
+    solve_residuals : dict
+        per-block implicit solver residuals (block -> residual)
+    solve_iterations : int
+        number of implicit solver iterations taken
+    step_errors : dict
+        per-block adaptive step data (block -> (success, err_norm, scale))
+    """
+    time: float = 0.0
+    loop_residuals: dict = field(default_factory=dict)
+    loop_iterations: int = 0
+    solve_residuals: dict = field(default_factory=dict)
+    solve_iterations: int = 0
+    step_errors: dict = field(default_factory=dict)
+
+
+    @staticmethod
+    def _label(obj):
+        """Human-readable label for a block or booster."""
+        if hasattr(obj, 'connection'):
+            return str(obj.connection)
+        return obj.__class__.__name__
+
+
+    def worst_block(self):
+        """Block with the highest residual across solve and step errors.
+
+        Returns
+        -------
+        tuple[str, float] or None
+            (label, error) or None if no data
+        """
+        worst, worst_err = None, -1.0
+
+        for obj, err in self.solve_residuals.items():
+            if err > worst_err:
+                worst, worst_err = obj, err
+
+        for obj, (_, err_norm, _) in self.step_errors.items():
+            if err_norm > worst_err:
+                worst, worst_err = obj, err_norm
+
+        if worst is None:
+            return None
+        return self._label(worst), worst_err
+
+
+    def worst_booster(self):
+        """Connection booster with the highest algebraic loop residual.
+
+        Returns
+        -------
+        tuple[str, float] or None
+            (label, residual) or None if no data
+        """
+        if not self.loop_residuals:
+            return None
+
+        worst = max(self.loop_residuals, key=self.loop_residuals.get)
+        return self._label(worst), self.loop_residuals[worst]
+
+
+    def summary(self):
+        """Formatted summary of this diagnostics snapshot.
+
+        Returns
+        -------
+        str
+            human-readable diagnostics summary
+        """
+        lines = [f"Diagnostics at t = {self.time:.6f}"]
+
+        if self.step_errors:
+            lines.append(f"\n  Adaptive step errors:")
+            for obj, (suc, err, scl) in self.step_errors.items():
+                status = "OK" if suc else "FAIL"
+                scl_str = f"{scl:.3f}" if scl is not None else "N/A"
+                lines.append(f"    {status}  {self._label(obj)}: err={err:.2e}, scale={scl_str}")
+
+        if self.solve_residuals:
+            lines.append(f"\n  Implicit solver residuals ({self.solve_iterations} iterations):")
+            for obj, err in self.solve_residuals.items():
+                lines.append(f"    {self._label(obj)}: {err:.2e}")
+
+        if self.loop_residuals:
+            lines.append(f"\n  Algebraic loop residuals ({self.loop_iterations} iterations):")
+            for obj, err in self.loop_residuals.items():
+                lines.append(f"    {self._label(obj)}: {err:.2e}")
+
+        return "\n".join(lines)
diff --git a/tests/pathsim/test_diagnostics.py b/tests/pathsim/test_diagnostics.py
