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/include/ndiffusion/solver_detail.hpp b/cpp/include/ndiffusion/solver_detail.hpp index c9dba66..3768db6 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,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::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/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, 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 583b6d4..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; @@ -488,6 +490,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 +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; diff --git a/cpp/src/solver_2d_unstructured.cpp b/cpp/src/solver_2d_unstructured.cpp index 91d9dfa..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; @@ -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; diff --git a/tests/test_cg_ab.py b/tests/test_cg_ab.py index 61c1e00..34470d5 100644 --- a/tests/test_cg_ab.py +++ b/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) == []