From 61174cf19225a6b81e29bd83b939392399288069 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 16 Sep 2026 14:49:47 -0400 Subject: [PATCH] Safeguard reused FGMRES guesses against the current true residual --- experiments/KRYLOV_HYPOTHESIS_2026-09-16.md | 37 +++++++++++++++ .../interface/solver/linear_backend.cppm | 31 ++++++++++-- tests/solver/stellar_equilibrium_runtime.cpp | 47 +++++++++++++++++++ 3 files changed, 111 insertions(+), 4 deletions(-) create mode 100644 experiments/KRYLOV_HYPOTHESIS_2026-09-16.md diff --git a/experiments/KRYLOV_HYPOTHESIS_2026-09-16.md b/experiments/KRYLOV_HYPOTHESIS_2026-09-16.md new file mode 100644 index 0000000..329aaf8 --- /dev/null +++ b/experiments/KRYLOV_HYPOTHESIS_2026-09-16.md @@ -0,0 +1,37 @@ +# Krylov information retention and correction reuse + +Baseline: `80cc3d8`. This branch evaluates two distinct effects without adding +public options or changing physical specifications, normalization, preconditioner, +linear tolerances, or nonlinear acceptance. + +The first control changes the existing FGMRES restart setting from 40 to 160. +Frozen solves in the audit established linear benefits; live Newton runs must +establish whether the resulting corrections remain useful and physically accurate. + +The production candidate safeguards a supplied linear initial guess using its +true residual under the current operator. It retains the guess when +`||b-Jp_guess|| <= ||b||`, otherwise starts at zero. A nonzero guess costs one +extra operator action, included in reported counts and solve time. Zero guesses +do not incur an extra action. This is neither Krylov basis recycling nor a +guarantee of fewer iterations: a smaller initial residual can still excite a +harder set of modes. The candidate must earn its cost in matched live runs. + +The retained `(1-alpha)p` heuristic from damped Newton remains intact, as does +the existing backend API. The safeguard only declines a guess that increases +the current true defect. It also handles a nonfinite operator action on an +otherwise finite guess by choosing the finite zero guess, collectively across +the communicator. Existing input checks still reject nonfinite supplied vectors. + +Required comparisons: + +- Same no-phase n=1.5 seed at restart 160 with baseline and candidate; record + each initial linear residual relative to its RHS, accepted step lengths, + true defects, total time and physical endpoint diagnostics. +- n=1 controls with full steps, where the Newton heuristic already produces + zero guesses; their physical trajectory should remain unchanged. +- Existing beneficial-warm-start test plus the new adverse-guess regression. + +Do not promote this change merely because the small regression passes. A live +improvement, unchanged validation, and sensible overhead are required. Virial +success means both measured definitions satisfy 1e-5 where pressure is defined; +interior EOS/Bernoulli consistency and algebraic convergence remain separate. diff --git a/libmeanfield/interface/solver/linear_backend.cppm b/libmeanfield/interface/solver/linear_backend.cppm index 47a99a0..88f2db2 100644 --- a/libmeanfield/interface/solver/linear_backend.cppm +++ b/libmeanfield/interface/solver/linear_backend.cppm @@ -62,8 +62,9 @@ export namespace mean_field::solver { * Numerical termination is data, not an exception. Implementations throw * for invalid controls, configuration, dimensions, or violated lifetime * contracts. All reported norms are communicator-global Euclidean norms. - * Solve uses the incoming correction as its initial guess and overwrites - * it with the final correction, so initialResidualNorm is ||b - A x_0||. + * Solve may safeguard the incoming correction before using it as an + * initial guess, and overwrites it with the final correction. + * initialResidualNorm is ||b - A x_0|| for the effective starting guess. * Convergence remains relative to the right-hand side rather than the * quality of a particular initial guess: * @@ -88,13 +89,13 @@ export namespace mean_field::solver { double reportedResidualNorm{0.0}; double trueResidualNorm{0.0}; double relativeTrueResidualNorm{0.0}; - // Includes MFEM's initial-guess residual application and the + // Includes any initial-guess safeguard, MFEM's initial residual, and the // post-solve application used to verify the true residual. std::uint64_t operatorApplications{0}; std::uint64_t inversePreconditionerApplications{0}; double solveSeconds{0.0}; // These totals include only completed applications. Operator time also - // includes the post-solve true-residual verification application. + // includes initial-guess safeguarding and post-solve true-residual verification. double operatorSeconds{0.0}; double inversePreconditionerSeconds{0.0}; @@ -769,6 +770,28 @@ export namespace mean_field::solver::linear { m_solver.SetMaxIter(control.maximumIterations); const Clock::time_point start = Clock::now(); + // A reused Newton correction belongs to an earlier linearization. + // Keep it only when its true defect is no worse than a cold start + // for this operator and right-hand side. This is a residual safeguard, + // not a claim that residual size alone predicts Krylov convergence. + bool locallyZeroGuess = true; + for (int index = 0; index < correction.Size(); ++index) { + if (correction(index) != 0.0) { + locallyZeroGuess = false; + break; + } + } + if (!AllRanksAgree(locallyZeroGuess, m_communicator)) { + m_countedOperation.Mult(correction, m_operationAction); + m_trueResidual = m_rightHandSide; + m_trueResidual -= m_operationAction; + const double guessResidualNorm = AllRanksAgree(LocallyFinite(m_trueResidual), m_communicator) + ? GlobalNorm(m_trueResidual, m_communicator) + : std::numeric_limits::quiet_NaN(); + if (!std::isfinite(guessResidualNorm) || guessResidualNorm > rightHandSideNorm) { + correction = 0.0; + } + } m_solver.Mult(m_rightHandSide, correction); const double solveSeconds = MaximumRankValue(std::chrono::duration(Clock::now() - start).count(), m_communicator); diff --git a/tests/solver/stellar_equilibrium_runtime.cpp b/tests/solver/stellar_equilibrium_runtime.cpp index 0d6162f..8e49e0b 100644 --- a/tests/solver/stellar_equilibrium_runtime.cpp +++ b/tests/solver/stellar_equilibrium_runtime.cpp @@ -743,6 +743,53 @@ TEST_CASE( CHECK(correction(1) == Approx(0.6).margin(1.0e-11)); } +TEST_CASE( + "MFEM FGMRES Discards A Supplied Guess That Worsens The Current True Residual", + "[solver][linear][fgmres][warm-start][safeguard]" +) { + using namespace mean_field; + using namespace stellar_solver_architecture_test; + using Catch::Approx; + + TwoByTwoOperator operation; + IdentityInverse inverse; + auto backend = solver::linear::prepareLinearBackend( + solver::linear::FGMRES({.restartLength = 1, .printLevel = -1}), operation, inverse, MPI_COMM_WORLD + ); + mfem::Vector rhs(2), correction(2); + rhs(0) = 1.0; + rhs(1) = 2.0; + const solver::LinearSolveControl control{ + .relativeTolerance = 1.0e-12, .absoluteTolerance = 1.0e-14, .maximumIterations = 200 + }; + + correction = 0.0; + const auto cold = backend.Solve(rhs, correction, control); + REQUIRE(cold.Converged()); + int rank = 0; + REQUIRE(MPI_Comm_rank(MPI_COMM_WORLD, &rank) == MPI_SUCCESS); + correction = 0.0; + SECTION("A finite adverse guess on one rank falls back collectively") { + if (rank == 0) { + correction(0) = 100.0; + correction(1) = -100.0; + } + } + SECTION("Finite coefficients with an overflowing operator action fall back collectively") { + if (rank == 0) { + correction = std::numeric_limits::max() / 2.0; + } + } + const auto guarded = backend.Solve(rhs, correction, control); + REQUIRE(guarded.Converged()); + CHECK(guarded.initialResidualNorm == Approx(cold.rightHandSideNorm)); + CHECK(guarded.iterations == cold.iterations); + CHECK(guarded.operatorApplications == cold.operatorApplications + 1); + CHECK(correction(0) == Approx(0.1).margin(1.0e-11)); + CHECK(correction(1) == Approx(0.6).margin(1.0e-11)); + CHECK(guarded.trueResidualNorm <= control.ConvergenceThreshold(guarded.rightHandSideNorm)); +} + TEST_CASE( "A User-Owned Context Certifies Views Only Through Evaluation Reports", "[solver][architecture][ownership][report][view]"