Cleanup nested preconditioner implementation

AUTHORS.md

@@ -71,6 +71,7 @@ Eduardo Molina
 Edwin van der Weide
 Eitan Aberman
 Ethan Alan Hereth
+Ezgi Orbay Akcengiz
 Florian Dittmann
 Filip Hahs
 Francesco Poli

Common/include/CConfig.hpp

@@ -515,6 +515,7 @@ class CConfig {
   Kind_SlopeLimit_AdjFlow,      /*!< \brief Slope limiter for the adjoint equation.*/
   Kind_SlopeLimit_Heat,         /*!< \brief Slope limiter for the adjoint equation.*/
   Kind_SlopeLimit_Species;      /*!< \brief Slope limiter for the species equation.*/
+  LINEAR_SOLVER_INNER Kind_Linear_Solver_Inner; /*!< \brief Inner solver used in nested Krylov schemes. */
   unsigned short Kind_FluidModel,  /*!< \brief Kind of the Fluid Model: Ideal, van der Waals, etc. */
   Kind_InitOption,                 /*!< \brief Kind of Init option to choose if initializing with Reynolds number or with thermodynamic conditions   */
   Kind_GridMovement,               /*!< \brief Kind of the static mesh movement. */
@@ -4288,6 +4289,10 @@ class CConfig {
    */
   unsigned short GetKind_Linear_Solver(void) const { return Kind_Linear_Solver; }
 
+  /*!
+   * \brief Get the inner linear solver used in nested Krylov linear solvers.
+   */
+  LINEAR_SOLVER_INNER GetKind_Linear_Solver_Inner(void) const { return Kind_Linear_Solver_Inner; }
 
   /*!
    * \brief Get the kind of preconditioner for the implicit solver.

Common/include/linear_algebra/CPreconditioner.hpp

@@ -28,6 +28,7 @@
 
 #pragma once
 
+#include <functional>
 #include "../CConfig.hpp"
 #include "../geometry/CGeometry.hpp"
 #include "CSysVector.hpp"
@@ -305,6 +306,25 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
   inline void Build() override { sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact); }
 };
 
+/*!
+ * \class CAbstractPreconditioner
+ * \brief Applies a std::function as the preconditioning operation.
+ * \note This can be used to treat almost anything as a preconditioner.
+ */
+template <class ScalarType>
+class CAbstractPreconditioner final : public CPreconditioner<ScalarType> {
+ private:
+  std::function<void(const CSysVector<ScalarType>&, CSysVector<ScalarType>&)> impl;
+
+ public:
+  CAbstractPreconditioner() = delete;
+
+  template <class F>
+  explicit CAbstractPreconditioner(const F& function) : impl(function) {}
+
+  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override { impl(u, v); }
+};
+
 template <class ScalarType>
 CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOLVER_PREC kind,
                                                                  CSysMatrix<ScalarType>& jacobian, CGeometry* geometry,

Common/include/linear_algebra/CSysSolve.hpp

@@ -31,6 +31,7 @@
 #include "../containers/C2DContainer.hpp"
 
 #include <cmath>
+#include <memory>
 #include <vector>
 #include <iostream>
 #include <cstdlib>
@@ -115,6 +116,9 @@ class CSysSolve {
   bool recomputeRes = false;         /*!< \brief Recompute the residual after inner iterations, if monitoring. */
   unsigned long monitorFreq = 10;    /*!< \brief Monitoring frequency. */
 
+  /*!< \brief Inner solver for nested preconditioning. */
+  std::unique_ptr<CSysSolve<ScalarType>> inner_solver;
+
   /*!
    * \brief sign transfer function
    * \param[in] x - value having sign prescribed

Common/include/option_structure.hpp

@@ -2369,6 +2369,19 @@ static const MapType<std::string, ENUM_LINEAR_SOLVER> Linear_Solver_Map = {
   MakePair("PASTIX_LU", PASTIX_LU)
 };
 
+/*!
+ * \brief Inner solver for nested linear solver, only compatible with "flexible" linear solvers.
+ */
+enum class LINEAR_SOLVER_INNER {
+  NONE,     /*!< \brief Do not use a nested linear solver. */
+  BCGSTAB,  /*!< \brief Use BCGSTAB as the preconditioning linear solver. */
+};
+static const MapType<std::string, LINEAR_SOLVER_INNER> Inner_Linear_Solver_Map = {
+  MakePair("NONE", LINEAR_SOLVER_INNER::NONE)
+  MakePair("BCGSTAB", LINEAR_SOLVER_INNER::BCGSTAB)
+};
+
+
 /*!
  * \brief Types surface continuity at the intersection with the FFD
  */

Common/src/CConfig.cpp

@@ -1887,6 +1887,8 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("LINEAR_SOLVER_SMOOTHER_RELAXATION", Linear_Solver_Smoother_Relaxation, 1.0);
   /* DESCRIPTION: Custom number of threads used for additive domain decomposition for ILU and LU_SGS (0 is "auto"). */
   addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0);
+  /* DESCRIPTION: Use an inner linear solver. */
+  addEnumOption("LINEAR_SOLVER_INNER", Kind_Linear_Solver_Inner, Inner_Linear_Solver_Map, LINEAR_SOLVER_INNER::NONE);
   /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */
   addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0);
   /* DESCRIPTION: Relaxation of the CHT coupling */
@@ -1907,7 +1909,6 @@ void CConfig::SetConfig_Options() {
   addEnumOption("DISCADJ_LIN_SOLVER", Kind_DiscAdj_Linear_Solver, Linear_Solver_Map, FGMRES);
   /* DESCRIPTION: Preconditioner for the discrete adjoint Krylov linear solvers */
   addEnumOption("DISCADJ_LIN_PREC", Kind_DiscAdj_Linear_Prec, Linear_Solver_Prec_Map, ILU);
-  /* DESCRIPTION: Linear solver for the discete adjoint systems */
 
   /* DESCRIPTION: Maximum update ratio value for flow density and energy variables */
   addDoubleOption("MAX_UPDATE_FLOW", MaxUpdateFlow, 0.2);
@@ -7260,10 +7261,15 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
             case BCGSTAB:
             case FGMRES:
             case RESTARTED_FGMRES:
-              if (Kind_Linear_Solver == BCGSTAB)
+              if (Kind_Linear_Solver == BCGSTAB) {
                 cout << "BCGSTAB is used for solving the linear system." << endl;
-              else
-                cout << "FGMRES is used for solving the linear system." << endl;
+              } else {
+                if (Kind_Linear_Solver_Inner == LINEAR_SOLVER_INNER::BCGSTAB){
+                  cout << "Nested FGMRES (FGMRES with inner BiCGSTAB) is used for solving the linear system." << endl;
+                } else {
+                  cout << "FGMRES is used for solving the linear system." << endl;
+                }
+              }
               switch (Kind_Linear_Solver_Prec) {
                 case ILU: cout << "Using a ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break;
                 case LINELET: cout << "Using a linelet preconditioning."<< endl; break;
@@ -7310,6 +7316,8 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
               cout << "FGMRES is used for solving the linear system." << endl;
               cout << "Convergence criteria of the linear solver: "<< Linear_Solver_Error <<"."<< endl;
               cout << "Max number of iterations: "<< Linear_Solver_Iter <<"."<< endl;
+              if (Kind_Linear_Solver_Inner == LINEAR_SOLVER_INNER::BCGSTAB)
+                cout << "Nested BiCGSTAB is used as the inner solver." << endl;
               break;
             case CONJUGATE_GRADIENT:
               cout << "A Conjugate Gradient method is used for solving the linear system." << endl;

Common/src/linear_algebra/CSysSolve.cpp

@@ -35,6 +35,7 @@
 #include "../../include/linear_algebra/CPreconditioner.hpp"
 
 #include <limits>
+#include <memory>
 
 /*!
  * \brief Epsilon used in CSysSolve depending on datatype to
@@ -110,7 +111,7 @@ void CSysSolve<ScalarType>::SolveReduced(int n, const su2matrix<ScalarType>& Hsb
 
 template <class ScalarType>
 void CSysSolve<ScalarType>::ModGramSchmidt(bool shared_hsbg, int i, su2matrix<ScalarType>& Hsbg,
-                                           vector<CSysVector<ScalarType> >& w) const {
+                                           vector<CSysVector<ScalarType>>& w) const {
   const auto thread = omp_get_thread_num();
 
   /*--- If Hsbg is shared by multiple threads calling this function, only one
@@ -902,9 +903,8 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
       break;
     }
 
-    /*--- Normal mode
-     * assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', but does not enforce it to avoid compiler warning.
-     * ---*/
+    /*--- Normal mode assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD',
+     * but does not enforce it to avoid compiler warning.  ---*/
     default: {
       KindSolver = config->GetKind_Linear_Solver();
       KindPrecond = config->GetKind_Linear_Solver_Prec();
@@ -915,6 +915,17 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
     }
   }
 
+  const bool nested = (KindSolver == FGMRES || KindSolver == RESTARTED_FGMRES || KindSolver == SMOOTHER) &&
+                      config->GetKind_Linear_Solver_Inner() != LINEAR_SOLVER_INNER::NONE;
+
+  if (nested && !inner_solver) {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+      inner_solver = std::make_unique<CSysSolve<ScalarType>>(LINEAR_SOLVER_MODE::STANDARD);
+      inner_solver->SetxIsZero(true);
+    }
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+  }
+
   /*--- Stop the recording for the linear solver ---*/
   bool TapeActive = NO;
 
@@ -948,13 +959,25 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
 
     auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
 
-    const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
-
-    auto precond = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
-
     /*--- Build preconditioner. ---*/
 
-    precond->Build();
+    const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+    auto* normal_prec = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
+    normal_prec->Build();
+
+    CPreconditioner<ScalarType>* nested_prec = nullptr;
+    if (nested) {
+      nested_prec = new CAbstractPreconditioner<ScalarType>([&](const CSysVector<ScalarType>& u,
+                                                                CSysVector<ScalarType>& v) {
+        /*--- Initialize to 0 to be safe. ---*/
+        v = ScalarType{};
+        ScalarType unused{};
+        /*--- Handle other types here if desired but do not call Solve because
+         * that will create issues with the AD external function. ---*/
+        (void)inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, SolverTol, MaxIter, unused, false, config);
+      });
+    }
+    const auto* precond = nested ? nested_prec : normal_prec;
 
     /*--- Solve system. ---*/
 
@@ -1000,7 +1023,8 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
 
     HandleTemporariesOut(LinSysSol);
 
-    delete precond;
+    delete normal_prec;
+    delete nested_prec;
 
     if (TapeActive) {
       /*--- To keep the behavior of SU2_DOT, but not strictly required since jacobian is symmetric(?). ---*/
@@ -1085,9 +1109,8 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
       break;
     }
 
-    /*--- Normal mode
-     * assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', but does not enforce it to avoid compiler warning.
-     * ---*/
+    /*--- Normal mode assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD',
+     * but does not enforce it to avoid compiler warning. ---*/
     default: {
       KindSolver = config->GetKind_Linear_Solver();
       KindPrecond = config->GetKind_Linear_Solver_Prec();
@@ -1098,19 +1121,43 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
     }
   }
 
+  const bool nested = (KindSolver == FGMRES || KindSolver == RESTARTED_FGMRES || KindSolver == SMOOTHER) &&
+                      config->GetKind_Linear_Solver_Inner() != LINEAR_SOLVER_INNER::NONE;
+
+  if (nested && !inner_solver) {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+      inner_solver = std::make_unique<CSysSolve<ScalarType>>(LINEAR_SOLVER_MODE::STANDARD);
+      inner_solver->SetxIsZero(true);
+    }
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+  }
+
   /*--- Set up preconditioner and matrix-vector product ---*/
 
-  const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+  auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
 
-  auto precond = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
+  const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+  auto* normal_prec = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
 
   /*--- If there was no call to solve first the preconditioner needs to be built here. ---*/
   if (directCall) {
     Jacobian.TransposeInPlace();
-    precond->Build();
+    normal_prec->Build();
   }
 
-  auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
+  CPreconditioner<ScalarType>* nested_prec = nullptr;
+  if (nested) {
+    nested_prec =
+        new CAbstractPreconditioner<ScalarType>([&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
+          /*--- Initialize to 0 to be safe. ---*/
+          v = ScalarType{};
+          ScalarType unused{};
+          /*--- Handle other types here if desired but do not call Solve because
+           * that will create issues with the AD external function. ---*/
+          (void)inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, SolverTol, MaxIter, unused, false, config);
+        });
+  }
+  const auto* precond = nested ? nested_prec : normal_prec;
 
   /*--- Solve the system ---*/
 
@@ -1154,7 +1201,8 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
 
   HandleTemporariesOut(LinSysSol);
 
-  delete precond;
+  delete normal_prec;
+  delete nested_prec;
 
   SU2_OMP_MASTER {
     Residual = residual;

SU2_CFD/include/integration/CNewtonIntegration.hpp

@@ -158,7 +158,11 @@ class CNewtonIntegration final : public CIntegration {
     auto product = CSysMatrixVectorProduct<MixedScalar>(solvers[FLOW_SOL]->Jacobian, geometry, config);
     v = MixedScalar(0.0);
     MixedScalar eps_t = eps;
-    iters = solvers[FLOW_SOL]->System.FGMRES_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+    if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::NONE) {
+      iters = solvers[FLOW_SOL]->System.FGMRES_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+    } else {
+      iters = solvers[FLOW_SOL]->System.BCGSTAB_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+    }
     eps = eps_t;
     return iters;
   }

TestCases/hybrid_regression.py

@@ -404,7 +404,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.788405, -11.040877, 0.000008, 0.309505]
+    inc_turb_naca0012.test_vals = [-4.758062, -10.974496, -0.000005, -0.028654, 4, -5.397415, 2, -6.426845]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST

TestCases/incomp_rans/naca0012/naca0012.cfg

@@ -50,8 +50,8 @@ MARKER_MONITORING= ( airfoil )
 %
 NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
 CFL_NUMBER= 10.0
-CFL_ADAPT= NO
-CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+CFL_ADAPT= YES
+CFL_ADAPT_PARAM= ( 0.8, 1.1, 1.0, 100.0 )
 ITER= 2500
 
 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
@@ -64,6 +64,7 @@ SENS_REMOVE_SHARP= NO
 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
 %
 LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_INNER= BCGSTAB
 LINEAR_SOLVER_PREC= LU_SGS
 LINEAR_SOLVER_ERROR= 1E-4
 LINEAR_SOLVER_ITER= 5
@@ -107,4 +108,4 @@ GRAD_OBJFUNC_FILENAME= of_grad
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
 OUTPUT_WRT_FREQ= 100
-SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_NU_TILDE, LIFT, DRAG, TOTAL_HEATFLUX)
+SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_NU_TILDE, LIFT, DRAG, LINSOL)

TestCases/parallel_regression.py

@@ -647,7 +647,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.788595, -11.040942, -0.000002, 0.309519]
+    inc_turb_naca0012.test_vals = [-4.758062, -10.974496, -0.000006, -0.028654, 3, -5.740142, 2, -6.615162]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST

TestCases/serial_regression.py

@@ -442,7 +442,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.788495, -11.040895, 0.000023, 0.309502]
+    inc_turb_naca0012.test_vals = [-4.758063, -10.974497, -0.000005, -0.028655, 4, -6.041638, 2, -6.473951]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST

config_template.cfg

@@ -1543,7 +1543,11 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 % Linear solver or smoother for implicit formulations:
 % BCGSTAB, FGMRES, RESTARTED_FGMRES, CONJUGATE_GRADIENT (self-adjoint problems only), SMOOTHER.
 LINEAR_SOLVER= FGMRES
-%
+
+% Inner solver used for nested preconditioning (only possible with flexible linear solvers).
+% Options: NONE, BCGSTAB
+LINEAR_SOLVER_INNER= NONE
+
 % Same for discrete adjoint (smoothers not supported), replaces LINEAR_SOLVER in SU2_*_AD codes.
 DISCADJ_LIN_SOLVER= FGMRES
 %