Stop k-eigenvalue inner sweeps at round-off

CHANGELOG.md

@@ -6,6 +6,15 @@ All notable changes to ndiffusion are recorded here. The format follows
 
 ## [Unreleased]
 
+### Fixed
+
+- A k-eigenvalue solve with an `epsilon` tighter than round-off allows no longer
+  runs every inner iteration out to `max_inner`. The inner sweeps now stop once
+  their change stops shrinking. This mostly hit the CG inner solver on fine
+  meshes: a 247k-cell unstructured mesh at `epsilon=1e-12` took over 10 minutes
+  per outer iteration and now takes about 40 seconds. Converged results are
+  unchanged.
+
 ## [1.0.2] - 2026-10-12
 
 ### Fixed

cpp/include/ndiffusion/solver_detail.hpp

@@ -5,6 +5,7 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
+#include <limits>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -880,6 +881,42 @@ int cg_solve(int n,
     return it;
 }
 
+/// Ends an inner sweep loop once it has hit round-off.
+///
+/// The k-eigenvalue inner loops stop when the sweep-to-sweep change falls below
+/// `epsilon * 1e-3`, but that is unreachable when it sits below what round-off
+/// allows, which happens on fine meshes at a tight epsilon (the attainable
+/// residual grows with the condition number).  Each sweep then only shuffles
+/// the flux at round-off level and every outer iteration grinds through all
+/// max_inner sweeps.  A converging sweep loop keeps setting new lows in the
+/// change, so once `patience` sweeps have passed without beating the smallest
+/// change so far there is nothing left to converge.  Counting from the smallest
+/// change rather than the previous one matters: at round-off the change can
+/// alternate up and down indefinitely (it does on some platforms), and a run of
+/// consecutive increases may never come.  `patience` is 1 where each group solve
+/// is exact (Thomas, CG to tolerance) and more for the point and line
+/// Gauss-Seidel sweeps, whose change need not shrink monotonically.  Only sweeps
+/// under 1e-6 count, so an early transient cannot end the loop.
+class SweepFloor {
+public:
+    explicit SweepFloor(int patience = 1) : patience_(patience) {}
+
+    bool reached(double change) {
+        if (change < best_) {
+            best_ = change;
+            stalled_ = 0;
+        } else if (change < 1e-6) {
+            ++stalled_;
+        }
+        return stalled_ >= patience_;
+    }
+
+private:
+    int    patience_;
+    int    stalled_ = 0;
+    double best_    = std::numeric_limits<double>::infinity();
+};
+
 /// Read a boolean from the environment: true unless unset/empty or starting with
 /// one of `0 f F n N`.  Used to pick the inner linear solver at runtime
 /// (`NDIFFUSION_KEIG_CG=1`) without recompiling, for A/B benchmarking.

cpp/python/bindings.cpp

@@ -708,7 +708,11 @@ PYBIND11_MODULE(_core, m) {
         "Uses cell-centered finite-volume method with point Gauss-Seidel.\n\n"
         "result.flux has shape (n_cells, n_groups).\n\n"
         "bc has size n_bc_types * n_groups; bc[tag*G+g] is the BC for\n"
-        "tag 'tag', group g.  Boundary faces with no matching bc_tag use tag 0.")
+        "tag 'tag', group g.  Boundary faces with no matching bc_tag use tag 0.\n\n"
+        "An epsilon tighter than round-off allows on the mesh cannot be met; the\n"
+        "inner sweeps stop once their change stops shrinking, so such a solve\n"
+        "ends at max_outer with a ConvergenceWarning instead of grinding through\n"
+        "max_inner sweeps every outer iteration.")
         .def(py::init<Materials,
                       UnstructuredMesh2D,
                       std::vector<BoundaryCondition>,

cpp/src/solver_1d.cpp

@@ -222,6 +222,7 @@ bool KEigenSolver::solve_A(
     std::vector<double> lower_g(N_), diag_g(N_), upper_g(N_);
     std::vector<double> rhs(N_), phi_g(N_), tw_c, tw_d, phi_prev;
 
+    SweepFloor floor;
     for (int inner = 0; inner < max_inner_; ++inner) {
         phi_prev = phi;
 
@@ -248,7 +249,8 @@ bool KEigenSolver::solve_A(
                 phi[g * N_ + i] = phi_g[i];
         }
 
-        if (rel_l2_diff(phi, phi_prev) < epsilon_ * 1e-3)
+        const double change = rel_l2_diff(phi, phi_prev);
+        if (change < epsilon_ * 1e-3 || floor.reached(change))
             return true;
     }
     return false;

cpp/src/solver_2d_structured.cpp

@@ -397,6 +397,7 @@ bool KEigenSolver2D::solve_A_gs(
 
     std::vector<double> lower(N_x), diag_g(N_x), upper_g(N_x);
     std::vector<double> rhs(N_x), phi_x(N_x), tw_c, tw_d, phi_prev;
+    SweepFloor floor(3);
 
     for (int inner = 0; inner < max_inner_; ++inner) {
         check_interrupt();
@@ -445,7 +446,8 @@ bool KEigenSolver2D::solve_A_gs(
             }
         }
 
-        if (rel_l2_diff(phi, phi_prev) < epsilon_ * 1e-3)
+        const double change = rel_l2_diff(phi, phi_prev);
+        if (change < epsilon_ * 1e-3 || floor.reached(change))
             return true;
     }
     return false;
@@ -488,6 +490,7 @@ bool KEigenSolver2D::solve_A_cg(
     std::vector<double> rhs_g(cells), x_g(cells), phi_prev;
     const int  max_cg   = 2 * cells + 50;
     const double cg_tol = std::min(epsilon_ * 1e-2, 1e-9);
+    SweepFloor floor;
 
     for (int sweep = 0; sweep < max_inner_; ++sweep) {
         check_interrupt();
@@ -524,7 +527,8 @@ bool KEigenSolver2D::solve_A_cg(
 
         // Single-group problems converge in one sweep (no scatter coupling).
         // Report failure if a within-group CG stalled so the caller can warn.
-        if (groups_ == 1 || rel_l2_diff(phi, phi_prev) < epsilon_ * 1e-3)
+        const double change = rel_l2_diff(phi, phi_prev);
+        if (groups_ == 1 || change < epsilon_ * 1e-3 || floor.reached(change))
             return cg_all_ok;
     }
     return false;

cpp/src/solver_2d_unstructured.cpp

@@ -963,6 +963,7 @@ bool KEigenSolverUnstructured2D::solve_A_gs(
           std::vector<double>& phi
 ) const {
     std::vector<double> phi_prev, gx, gy;
+    SweepFloor floor(3);
     for (int inner = 0; inner < max_inner_; ++inner) {
         check_interrupt();
         phi_prev = phi;
@@ -1002,7 +1003,8 @@ bool KEigenSolverUnstructured2D::solve_A_gs(
             }
         }
 
-        if (rel_l2_diff(phi, phi_prev) < epsilon_ * 1e-3)
+        const double change = rel_l2_diff(phi, phi_prev);
+        if (change < epsilon_ * 1e-3 || floor.reached(change))
             return true;
     }
     return false;
@@ -1054,6 +1056,7 @@ bool KEigenSolverUnstructured2D::solve_A_cg(
     // already satisfies takes no CG steps and reports zero change without having
     // converged anything.
     double change = 1.0;
+    SweepFloor floor;
     for (int sweep = 0; sweep < max_inner_; ++sweep) {
         const double tol = orthogonal_ ? cg_tol : std::max(cg_tol, 1e-2 * change);
         check_interrupt();
@@ -1099,7 +1102,9 @@ bool KEigenSolverUnstructured2D::solve_A_cg(
         // be iterated to consistency like any deferred correction.
         // Report failure if a within-group CG stalled so the caller can warn.
         change = rel_l2_diff(phi, phi_prev);
-        if (one_sweep || (tol == cg_tol && change < epsilon_ * 1e-3))
+        // floor sees every sweep: off an orthogonal mesh the tolerance tracks the
+        // change, so below round-off no sweep would ever count as tight.
+        if (one_sweep || (tol == cg_tol && change < epsilon_ * 1e-3) || floor.reached(change))
             return cg_all_ok;
     }
     return false;

tests/test_cg_ab.py

@@ -15,7 +15,10 @@
 with the NDIFFUSION_KEIG_CG environment variable.
 """
 
+import warnings
+
 import numpy as np
+import pytest
 
 import ndiffusion as nd
 
@@ -201,3 +204,59 @@ def test_skewed_mesh_agrees(self):
             cg = _unstructured(m, n, n, 60.0, True, skew_amount=0.3)
             assert abs(gs.keff - cg.keff) < 1e-7, (gs.keff, cg.keff)
             assert flux_cosine(gs.flux, cg.flux) > 1 - 1e-8
+
+
+# ---------------------------------------------------------------------------
+# Tolerance below round-off
+# ---------------------------------------------------------------------------
+
+
+def _inner_warnings(solver):
+    """Run a solve that is capped on purpose; return only the inner-cap warnings."""
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        solver.solve()
+    return [str(w.message) for w in caught if "inner solve" in str(w.message)]
+
+
+class TestRoundOffFloor:
+    """An epsilon below what round-off allows must not make every outer
+    iteration churn through all max_inner sweeps, with either inner solver:
+    the sweep loop stops once its change stops shrinking.  max_inner is set
+    well past what Gauss-Seidel needs to reach round-off, and the outer cap
+    still warns, as it should."""
+
+    MAX_INNER = 5000
+
+    @pytest.mark.parametrize("use_cg", [False, True])
+    def test_structured(self, use_cg):
+        m = two_group_mat()
+        s = nd.KEigenSolver2D(
+            mats=m, medium_map=[0] * 400,
+            edges_x=linspace(0.0, 80.0, 21), edges_y=linspace(0.0, 80.0, 21),
+            geom=nd.Geometry2D.XY, bc_x=[vacuum()] * 2, bc_y=[vacuum()] * 2,
+            epsilon=1e-16, max_outer=3, max_inner=self.MAX_INNER,
+        )
+        s.set_use_cg(use_cg)
+        assert _inner_warnings(s) == []
+
+    @pytest.mark.parametrize("use_cg", [False, True])
+    @pytest.mark.parametrize("skew_amount", [0.0, 0.3])
+    def test_unstructured(self, use_cg, skew_amount):
+        m = two_group_mat()
+        mesh = make_quad_mesh(14, 14, 60.0, 60.0)
+        if skew_amount:
+            mesh = skew(mesh, 60.0, skew_amount)
+        s = nd.KEigenSolverUnstructured2D(
+            mats=m, mesh=mesh, bc=[vacuum()] * 2, epsilon=1e-16,
+            max_outer=3, max_inner=self.MAX_INNER, use_cg=use_cg,
+        )
+        assert _inner_warnings(s) == []
+
+    def test_one_dimensional(self):
+        m = two_group_mat()
+        s = nd.KEigenSolver(
+            m, [0] * 50, linspace(0.0, 80.0, 51), nd.Geometry.Slab,
+            [vacuum()] * 2, epsilon=1e-16, max_outer=3, max_inner=self.MAX_INNER,
+        )
+        assert _inner_warnings(s) == []