Safeguard reused FGMRES guesses against the current true residual

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<double>::quiet_NaN();
+                    if (!std::isfinite(guessResidualNorm) || guessResidualNorm > rightHandSideNorm) {
+                        correction = 0.0;
+                    }
+                }
                 m_solver.Mult(m_rightHandSide, correction);
                 const double solveSeconds =
                     MaximumRankValue(std::chrono::duration<double>(Clock::now() - start).count(), m_communicator);