From 3308a81727015a8b5aa333dcd370e8bb9b6797b6 Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Tue, 13 Oct 2026 20:57:20 -0600 Subject: [PATCH 1/4] Stop the CG group sweeps once their change bottoms out at round-off. --- cpp/include/ndiffusion/solver_detail.hpp | 24 ++++++++++++ cpp/src/solver_2d_structured.cpp | 4 +- cpp/src/solver_2d_unstructured.cpp | 5 ++- tests/test_cg_ab.py | 50 ++++++++++++++++++++++++ 4 files changed, 81 insertions(+), 2 deletions(-) diff --git a/cpp/include/ndiffusion/solver_detail.hpp b/cpp/include/ndiffusion/solver_detail.hpp index c9dba66..43fee29 100644 --- a/cpp/include/ndiffusion/solver_detail.hpp +++ b/cpp/include/ndiffusion/solver_detail.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -880,6 +881,29 @@ int cg_solve(int n, return it; } +/// Ends a Gauss-Seidel sweep over energy groups once it has hit round-off. +/// +/// The CG sweep 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 CG +/// residual grows with the condition number). Each sweep then only shuffles x +/// at round-off level and every outer iteration grinds through all max_inner +/// sweeps. The group iteration (and the deferred non-orthogonal correction, +/// whose CG tolerance tracks the change) shrinks the change geometrically, so a +/// change that fails to drop below the previous one means there is nothing left +/// to converge. The 1e-6 guard keeps an early transient from ending the loop. +class SweepFloor { +public: + bool reached(double change) { + const bool stalled = change >= last_ && change < 1e-6; + last_ = change; + return stalled; + } + +private: + double last_ = std::numeric_limits::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. diff --git a/cpp/src/solver_2d_structured.cpp b/cpp/src/solver_2d_structured.cpp index 583b6d4..7243541 100644 --- a/cpp/src/solver_2d_structured.cpp +++ b/cpp/src/solver_2d_structured.cpp @@ -488,6 +488,7 @@ bool KEigenSolver2D::solve_A_cg( std::vector 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 +525,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; diff --git a/cpp/src/solver_2d_unstructured.cpp b/cpp/src/solver_2d_unstructured.cpp index 91d9dfa..8a87f04 100644 --- a/cpp/src/solver_2d_unstructured.cpp +++ b/cpp/src/solver_2d_unstructured.cpp @@ -1054,6 +1054,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 +1100,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; diff --git a/tests/test_cg_ab.py b/tests/test_cg_ab.py index 61c1e00..81f10e4 100644 --- a/tests/test_cg_ab.py +++ b/tests/test_cg_ab.py @@ -15,6 +15,8 @@ with the NDIFFUSION_KEIG_CG environment variable. """ +import warnings + import numpy as np import ndiffusion as nd @@ -201,3 +203,51 @@ 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 group sweeps: the sweep stops once + its change stops shrinking. The outer cap still warns, as it should.""" + + def test_structured(self): + 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=50, + ) + s.set_use_cg(True) + assert _inner_warnings(s) == [] + + def test_unstructured(self): + m = two_group_mat() + mesh = make_quad_mesh(14, 14, 80.0, 80.0) + s = nd.KEigenSolverUnstructured2D( + mats=m, mesh=mesh, bc=[vacuum()] * 2, epsilon=1e-16, + max_outer=3, max_inner=50, use_cg=True, + ) + assert _inner_warnings(s) == [] + + def test_unstructured_skewed(self): + m = two_group_mat() + mesh = skew(make_quad_mesh(14, 14, 60.0, 60.0), 60.0, 0.3) + s = nd.KEigenSolverUnstructured2D( + mats=m, mesh=mesh, bc=[vacuum()] * 2, epsilon=1e-16, + max_outer=3, max_inner=50, use_cg=True, + ) + assert _inner_warnings(s) == [] From 8cec2f6ddc4445177f1432c3ccbbfcb37a40043f Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Tue, 13 Oct 2026 21:11:30 -0600 Subject: [PATCH 2/4] Apply the round-off sweep stop to the Gauss-Seidel and 1-D k-eigenvalue inner loops too. --- cpp/include/ndiffusion/solver_detail.hpp | 30 ++++++++++++-------- cpp/src/solver_1d.cpp | 4 ++- cpp/src/solver_2d_structured.cpp | 4 ++- cpp/src/solver_2d_unstructured.cpp | 4 ++- tests/test_cg_ab.py | 35 +++++++++++++++--------- 5 files changed, 49 insertions(+), 28 deletions(-) diff --git a/cpp/include/ndiffusion/solver_detail.hpp b/cpp/include/ndiffusion/solver_detail.hpp index 43fee29..4b74eae 100644 --- a/cpp/include/ndiffusion/solver_detail.hpp +++ b/cpp/include/ndiffusion/solver_detail.hpp @@ -881,27 +881,33 @@ int cg_solve(int n, return it; } -/// Ends a Gauss-Seidel sweep over energy groups once it has hit round-off. +/// Ends an inner sweep loop once it has hit round-off. /// -/// The CG sweep loops stop when the sweep-to-sweep change falls below +/// 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 CG -/// residual grows with the condition number). Each sweep then only shuffles x -/// at round-off level and every outer iteration grinds through all max_inner -/// sweeps. The group iteration (and the deferred non-orthogonal correction, -/// whose CG tolerance tracks the change) shrinks the change geometrically, so a -/// change that fails to drop below the previous one means there is nothing left -/// to converge. The 1e-6 guard keeps an early transient from ending the loop. +/// 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 shrinks the change geometrically, +/// so a change that keeps failing to drop below the previous one means there is +/// nothing left to converge. `patience` is how many such sweeps in a row it +/// takes: 1 where each group solve is exact (Thomas, CG to tolerance), more for +/// the point and line Gauss-Seidel sweeps, whose change need not shrink +/// monotonically. The 1e-6 guard keeps an early transient from ending the loop. class SweepFloor { public: + explicit SweepFloor(int patience = 1) : patience_(patience) {} + bool reached(double change) { - const bool stalled = change >= last_ && change < 1e-6; + stalled_ = (change >= last_ && change < 1e-6) ? stalled_ + 1 : 0; last_ = change; - return stalled; + return stalled_ >= patience_; } private: - double last_ = std::numeric_limits::infinity(); + int patience_; + int stalled_ = 0; + double last_ = std::numeric_limits::infinity(); }; /// Read a boolean from the environment: true unless unset/empty or starting with diff --git a/cpp/src/solver_1d.cpp b/cpp/src/solver_1d.cpp index f890199..e9ef55f 100644 --- a/cpp/src/solver_1d.cpp +++ b/cpp/src/solver_1d.cpp @@ -222,6 +222,7 @@ bool KEigenSolver::solve_A( std::vector lower_g(N_), diag_g(N_), upper_g(N_); std::vector 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; diff --git a/cpp/src/solver_2d_structured.cpp b/cpp/src/solver_2d_structured.cpp index 7243541..7caa98f 100644 --- a/cpp/src/solver_2d_structured.cpp +++ b/cpp/src/solver_2d_structured.cpp @@ -397,6 +397,7 @@ bool KEigenSolver2D::solve_A_gs( std::vector lower(N_x), diag_g(N_x), upper_g(N_x); std::vector 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; diff --git a/cpp/src/solver_2d_unstructured.cpp b/cpp/src/solver_2d_unstructured.cpp index 8a87f04..689daca 100644 --- a/cpp/src/solver_2d_unstructured.cpp +++ b/cpp/src/solver_2d_unstructured.cpp @@ -963,6 +963,7 @@ bool KEigenSolverUnstructured2D::solve_A_gs( std::vector& phi ) const { std::vector 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; diff --git a/tests/test_cg_ab.py b/tests/test_cg_ab.py index 81f10e4..34470d5 100644 --- a/tests/test_cg_ab.py +++ b/tests/test_cg_ab.py @@ -18,6 +18,7 @@ import warnings import numpy as np +import pytest import ndiffusion as nd @@ -220,34 +221,42 @@ def _inner_warnings(solver): class TestRoundOffFloor: """An epsilon below what round-off allows must not make every outer - iteration churn through all max_inner group sweeps: the sweep stops once - its change stops shrinking. The outer cap still warns, as it should.""" + 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.""" - def test_structured(self): + 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=50, + epsilon=1e-16, max_outer=3, max_inner=self.MAX_INNER, ) - s.set_use_cg(True) + s.set_use_cg(use_cg) assert _inner_warnings(s) == [] - def test_unstructured(self): + @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, 80.0, 80.0) + 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=50, use_cg=True, + max_outer=3, max_inner=self.MAX_INNER, use_cg=use_cg, ) assert _inner_warnings(s) == [] - def test_unstructured_skewed(self): + def test_one_dimensional(self): m = two_group_mat() - mesh = skew(make_quad_mesh(14, 14, 60.0, 60.0), 60.0, 0.3) - s = nd.KEigenSolverUnstructured2D( - mats=m, mesh=mesh, bc=[vacuum()] * 2, epsilon=1e-16, - max_outer=3, max_inner=50, use_cg=True, + 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) == [] From 9a9fc67926e4e77643cb55bd8020b44736e53afb Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Tue, 13 Oct 2026 21:11:30 -0600 Subject: [PATCH 3/4] Note the round-off sweep stop in the changelog and the unstructured k solver docstring. --- CHANGELOG.md | 9 +++++++++ cpp/python/bindings.cpp | 6 +++++- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b2f891b..eaef83c 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/cpp/python/bindings.cpp b/cpp/python/bindings.cpp index 91b7429..7b0da5c 100644 --- a/cpp/python/bindings.cpp +++ b/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, From b909371e08caf83b001430787c35ce56e9cfcf2c Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Wed, 14 Oct 2026 23:11:17 -0600 Subject: [PATCH 4/4] Count stalled sweeps from the smallest change so far, since at round-off the change can alternate instead of rising. --- cpp/include/ndiffusion/solver_detail.hpp | 25 +++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/cpp/include/ndiffusion/solver_detail.hpp b/cpp/include/ndiffusion/solver_detail.hpp index 4b74eae..3768db6 100644 --- a/cpp/include/ndiffusion/solver_detail.hpp +++ b/cpp/include/ndiffusion/solver_detail.hpp @@ -888,26 +888,33 @@ int cg_solve(int n, /// 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 shrinks the change geometrically, -/// so a change that keeps failing to drop below the previous one means there is -/// nothing left to converge. `patience` is how many such sweeps in a row it -/// takes: 1 where each group solve is exact (Thomas, CG to tolerance), more for -/// the point and line Gauss-Seidel sweeps, whose change need not shrink -/// monotonically. The 1e-6 guard keeps an early transient from ending the loop. +/// 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) { - stalled_ = (change >= last_ && change < 1e-6) ? stalled_ + 1 : 0; - last_ = change; + if (change < best_) { + best_ = change; + stalled_ = 0; + } else if (change < 1e-6) { + ++stalled_; + } return stalled_ >= patience_; } private: int patience_; int stalled_ = 0; - double last_ = std::numeric_limits::infinity(); + double best_ = std::numeric_limits::infinity(); }; /// Read a boolean from the environment: true unless unset/empty or starting with