From d89a5f86d078602c7dd70da65ced45b91ec5a9f5 Mon Sep 17 00:00:00 2001
From: bigfooted <bigfootedrockmidget@hotmail.com>
Date: Sun, 27 Sep 2026 14:22:01 +0200
Subject: [PATCH 1/2] BCGSTAB: return the iterate with the smallest residual

BCGSTAB does not reduce the residual monotonically. With a small iteration
limit (e.g. LINEAR_SOLVER_ITER= 3) the last iterate can have a larger
residual than the initial guess, and returning it gives a bad update. On
the RAE2822 SA case with multigrid at CFL 200-1500 this caused erratic
divergence; keeping the best iterate makes every CFL converge.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 Common/include/linear_algebra/CSysSolve.hpp |  1 +
 Common/src/linear_algebra/CSysSolve.cpp     | 14 ++++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp
index 1f9bc851b92e..4d5aea756dd5 100644
--- a/Common/include/linear_algebra/CSysSolve.hpp
+++ b/Common/include/linear_algebra/CSysSolve.hpp
@@ -111,6 +111,7 @@ class CSysSolve {
 
   mutable VectorType r_0; /*!< \brief The "arbitrary" vector in BCGSTAB. */
   mutable VectorType v;   /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */
+  mutable VectorType x_best; /*!< \brief BCGSTAB iterate with the smallest residual. */
 
   mutable bool ritz_failed = false;
   mutable unsigned long k = 0, k_new = 0;
diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp
index 8ec1b19b23b3..b480dcacdd95 100644
--- a/Common/src/linear_algebra/CSysSolve.cpp
+++ b/Common/src/linear_algebra/CSysSolve.cpp
@@ -1087,6 +1087,7 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
       p.Initialize(nBlk, nBlkDomain, nVar, nullptr);
       v.Initialize(nBlk, nBlkDomain, nVar, nullptr);
       z.Initialize(nBlk, nBlkDomain, nVar, nullptr);
+      x_best.Initialize(nBlk, nBlkDomain, nVar, nullptr);
 
       bcg_ready = true;
     }
@@ -1135,6 +1136,8 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
   /*--- Initialization ---*/
 
   ScalarType alpha = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0;
+  ScalarType norm_best = norm_r;
+  x_best = x;
   p = ScalarType(0.0);
   v = ScalarType(0.0);
   r_0 = r;
@@ -1195,6 +1198,10 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
       /*--- Check if solution has converged, else output the relative residual if necessary ---*/
 
       norm_r = r.norm();
+      if (norm_r < norm_best) {
+        norm_best = norm_r;
+        x_best = x;
+      }
       if (norm_r < tol * norm0) break;
       if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
         SU2_OMP_MASTER
@@ -1204,6 +1211,13 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
     }
   }
 
+  /*--- BCGSTAB does not reduce the residual monotonically, return the best iterate. ---*/
+
+  if ((config->GetComm_Level() == COMM_FULL) && (norm_r > norm_best)) {
+    x = x_best;
+    norm_r = norm_best;
+  }
+
   /*--- Recalculate final residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {

From e08ba03788791c07b9afea6f8b8aed4e6dff279c Mon Sep 17 00:00:00 2001
From: bigfooted <bigfootedrockmidget@hotmail.com>
Date: Sun, 27 Sep 2026 17:02:34 +0200
Subject: [PATCH 2/2] clang

---
 Common/include/linear_algebra/CSysSolve.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp
index 4d5aea756dd5..98ab059804ce 100644
--- a/Common/include/linear_algebra/CSysSolve.hpp
+++ b/Common/include/linear_algebra/CSysSolve.hpp
@@ -109,8 +109,8 @@ class CSysSolve {
   mutable VectorType p;   /*!< \brief Direction in CG and BCGSTAB. */
   mutable VectorType z;   /*!< \brief Preconditioned residual/direction in CG/BCGSTAB. */
 
-  mutable VectorType r_0; /*!< \brief The "arbitrary" vector in BCGSTAB. */
-  mutable VectorType v;   /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */
+  mutable VectorType r_0;    /*!< \brief The "arbitrary" vector in BCGSTAB. */
+  mutable VectorType v;      /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */
   mutable VectorType x_best; /*!< \brief BCGSTAB iterate with the smallest residual. */
 
   mutable bool ritz_failed = false;