Merge pull request #2668 from su2code/pedro/more_inner_solvers

Add smoother to the options for nested preconditioning

Author: pcarruscag

Common/include/linear_algebra/CSysSolve.hpp

@@ -295,6 +295,12 @@ class CSysSolve {
                                       ScalarType& residual, bool monitoring, const CConfig* config,
                                       FgcrodrMode mode) const;
 
+  /*!
+   * \brief Creates the inner solver for nested preconditioning if the settings allow it.
+   * \returns True if the inner solver can be used.
+   */
+  bool SetupInnerSolver(unsigned short kind_solver, const CConfig* config);
+
  public:
   /*!
    * \brief default constructor of the class.

Common/include/option_structure.hpp

@@ -2389,10 +2389,12 @@ static const MapType<std::string, ENUM_LINEAR_SOLVER> Linear_Solver_Map = {
 enum class LINEAR_SOLVER_INNER {
   NONE,     /*!< \brief Do not use a nested linear solver. */
   BCGSTAB,  /*!< \brief Use BCGSTAB as the preconditioning linear solver. */
+  SMOOTHER, /*!< \brief Iterative smoother. */
 };
 static const MapType<std::string, LINEAR_SOLVER_INNER> Inner_Linear_Solver_Map = {
   MakePair("NONE", LINEAR_SOLVER_INNER::NONE)
   MakePair("BCGSTAB", LINEAR_SOLVER_INNER::BCGSTAB)
+  MakePair("SMOOTHER", LINEAR_SOLVER_INNER::SMOOTHER)
 };
 
 

Common/src/linear_algebra/CSysSolve.cpp

@@ -312,8 +312,8 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
   /*--- Only compute the residuals in full communication mode. ---*/
 
   if (config->GetComm_Level() == COMM_FULL) {
-    norm_r = r.norm();
     norm0 = b.norm();
+    norm_r = xIsZero ? norm0 : r.norm();
 
     /*--- Set the norm to the initial initial residual value ---*/
 
@@ -1032,8 +1032,8 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
   /*--- Only compute the residuals in full communication mode. ---*/
 
   if (config->GetComm_Level() == COMM_FULL) {
-    norm_r = r.norm();
     norm0 = b.norm();
+    norm_r = xIsZero ? norm0 : r.norm();
 
     /*--- Set the norm to the initial initial residual value ---*/
 
@@ -1205,8 +1205,8 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
   /*--- Only compute the residuals in full communication mode. ---*/
 
   if (config->GetComm_Level() == COMM_FULL) {
-    norm_r = r.norm();
     norm0 = b.norm();
+    norm_r = xIsZero ? norm0 : r.norm();
 
     /*--- Set the norm to the initial initial residual value ---*/
 
@@ -1246,15 +1246,19 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
      current residual, the system is linear so this saves some computation
      compared to re-evaluating r = b-A*x. ---*/
 
-    mat_vec(z, A_x);
+    if (!fix_iter_mode || i != m - 1) {
+      mat_vec(z, A_x);
+    }
 
     /*--- Update solution and residual with relaxation omega. Mathematically this
      is a modified Richardson iteration for the left-preconditioned system
      M^{-1}(b-A*x) which converges if ||I-w*M^{-1}*A|| < 1. Combining this method
      with a Gauss-Seidel preconditioner and w>1 is NOT equivalent to SOR. ---*/
 
     x += omega * z;
-    r -= omega * A_x;
+    if (!fix_iter_mode || i != m - 1) {
+      r -= omega * A_x;
+    }
 
     /*--- Only compute the residuals in full communication mode. ---*/
     /*--- Check if solution has converged, else output the relative residual if necessary. ---*/
@@ -1282,6 +1286,33 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
   return i;
 }
 
+template <class ScalarType>
+bool CSysSolve<ScalarType>::SetupInnerSolver(unsigned short kind_solver, const CConfig* config) {
+  bool flexible = false;
+  switch (kind_solver) {
+    case FGMRES:
+    case FGCRODR:
+    case RESTARTED_FGMRES:
+    case SMOOTHER:
+      flexible = true;
+      break;
+    default:
+      flexible = false;
+  }
+  const bool is_linear = config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::SMOOTHER;
+
+  if (config->GetKind_Linear_Solver_Inner() != LINEAR_SOLVER_INNER::NONE && (flexible || is_linear)) {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+    if (!inner_solver) {
+      inner_solver = std::make_unique<CSysSolve<ScalarType>>(LINEAR_SOLVER_MODE::STANDARD);
+      inner_solver->SetxIsZero(true);
+    }
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+    return true;
+  }
+  return false;
+}
+
 template <class ScalarType>
 unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, const CSysVector<su2double>& LinSysRes,
                                            CSysVector<su2double>& LinSysSol, CGeometry* geometry,
@@ -1334,16 +1365,7 @@ 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
-  }
+  const bool nested = SetupInnerSolver(KindSolver, config);
 
   /*--- Stop the recording for the linear solver ---*/
   bool TapeActive = NO;
@@ -1389,11 +1411,15 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
       auto f = [&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
         /*--- Initialize to 0 to be safe. ---*/
         v = ScalarType{};
-        ScalarType unused{};
+        ScalarType res{};
         /*--- 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, sqrt(SolverTol), MaxIter, unused, false,
-                                              config);
+        if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) {
+          inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config);
+        } else {
+          const auto smooth_iter = static_cast<unsigned long>(std::round(fmax(2, sqrt(MaxIter))));
+          inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config);
+        }
       };
       nested_prec = new CAbstractPreconditioner<ScalarType>(f);
     }
@@ -1545,16 +1571,7 @@ 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
-  }
+  const bool nested = SetupInnerSolver(KindSolver, config);
 
   /*--- Set up preconditioner and matrix-vector product ---*/
 
@@ -1572,13 +1589,15 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
   CPreconditioner<ScalarType>* nested_prec = nullptr;
   if (nested) {
     auto f = [&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
-      /*--- Initialize to 0 to be safe. ---*/
+      /*--- See "Solve". ---*/
       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, sqrt(SolverTol), MaxIter, unused, false,
-                                            config);
+      ScalarType res{};
+      if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) {
+        inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config);
+      } else {
+        const auto smooth_iter = static_cast<unsigned long>(std::round(fmax(2, sqrt(MaxIter))));
+        inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config);
+      }
     };
     nested_prec = new CAbstractPreconditioner<ScalarType>(f);
   }

SU2_CFD/include/integration/CNewtonIntegration.hpp

@@ -152,17 +152,25 @@ class CNewtonIntegration final : public CIntegration {
   template<class T, su2enable_if<std::is_same<T,MixedScalar>::value> = 0>
   inline unsigned long Preconditioner_impl(const CSysVector<T>& u, CSysVector<T>& v,
                                            unsigned long iters, Scalar& eps) const {
-    if (iters == 0) {
+    const auto inner_solver = config->GetKind_Linear_Solver_Inner();
+
+    if (iters == 0 || (iters == 1 && inner_solver == LINEAR_SOLVER_INNER::SMOOTHER)) {
       (*preconditioner)(u, v);
       return 0;
     }
     auto product = CSysMatrixVectorProduct<MixedScalar>(solvers[FLOW_SOL]->Jacobian, geometry, config);
     v = MixedScalar(0.0);
     MixedScalar eps_t = eps;
-    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);
+    switch (inner_solver) {
+      case LINEAR_SOLVER_INNER::NONE:
+        iters = solvers[FLOW_SOL]->System.FGMRES_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+        break;
+      case LINEAR_SOLVER_INNER::BCGSTAB:
+        iters = solvers[FLOW_SOL]->System.BCGSTAB_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+        break;
+      case LINEAR_SOLVER_INNER::SMOOTHER:
+        iters = solvers[FLOW_SOL]->System.Smoother_LinSolver(u, v, product, *preconditioner, 0, iters, eps_t, false, config);
+        break;
     }
     eps = eps_t;
     return iters;

config_template.cfg

@@ -1553,8 +1553,10 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 % BCGSTAB, FGMRES, FGCRODR, 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
+% Inner solver used for nested preconditioning (BCGSTAB is only possible with flexible linear solvers).
+% Options: NONE, BCGSTAB, SMOOTHER
+% Inner BCGSTAB does up to the number of iterations of the outer solver or sqrt(LINEAR_SOLVER_ERROR).
+% Inner SMOOTHER does max(2, sqrt(LINEAR_SOLVER_ITER))
 LINEAR_SOLVER_INNER= NONE
 
 % Same for discrete adjoint (smoothers not supported), replaces LINEAR_SOLVER in SU2_*_AD codes.