Merge branch 'develop' into verify-pr-2587-clean

Author: Shubham shukla

.gitmodules

@@ -30,3 +30,6 @@
 [submodule "externals/FADO"]
 	path = externals/FADO
 	url = https://github.com/pcarruscag/FADO.git
+[submodule "externals/eigen"]
+	path = externals/eigen
+	url = https://gitlab.com/libeigen/eigen.git

Common/include/CConfig.hpp

@@ -430,7 +430,7 @@ class CConfig {
   unsigned short nQuasiNewtonSamples;  /*!< \brief Number of samples used in quasi-Newton solution methods. */
   bool UseVectorization;       /*!< \brief Whether to use vectorized numerics schemes. */
   bool NewtonKrylov;           /*!< \brief Use a coupled Newton method to solve the flow equations. */
-  array<unsigned short,3> NK_IntParam{{20, 3, 2}}; /*!< \brief Integer parameters for NK method. */
+  array<unsigned short,4> NK_IntParam{{20, 3, 2, 0}}; /*!< \brief Integer parameters for NK method. */
   array<su2double,5> NK_DblParam{{-2.0, 0.1, -3.0, 1e-4, 1.0}}; /*!< \brief Floating-point parameters for NK method. */
   su2double NK_Relaxation = 1.0;
 
@@ -642,6 +642,7 @@ class CConfig {
   unsigned long Linear_Solver_Iter;              /*!< \brief Max iterations of the linear solver for the implicit formulation. */
   unsigned long Deform_Linear_Solver_Iter;       /*!< \brief Max iterations of the linear solver for the implicit formulation. */
   unsigned long Linear_Solver_Restart_Frequency; /*!< \brief Restart frequency of the linear solver for the implicit formulation. */
+  unsigned long Linear_Solver_Restart_Deflation; /*!< \brief Number of vectors used for deflated restarts. */
   unsigned long Linear_Solver_Prec_Threads;      /*!< \brief Number of threads per rank for ILU and LU_SGS preconditioners. */
   unsigned short Linear_Solver_ILU_n;            /*!< \brief ILU fill=in level. */
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
@@ -4342,6 +4343,11 @@ class CConfig {
    */
   unsigned long GetLinear_Solver_Restart_Frequency(void) const { return Linear_Solver_Restart_Frequency; }
 
+  /*!
+   * \brief Get the number of vectors used for deflated restarts.
+   */
+  unsigned long GetLinear_Solver_Restart_Deflation(void) const { return Linear_Solver_Restart_Deflation; }
+
   /*!
    * \brief Get the relaxation factor for iterative linear smoothers.
    * \return Relaxation factor.
@@ -4395,7 +4401,7 @@ class CConfig {
   /*!
    * \brief Get Newton-Krylov integer parameters.
    */
-  array<unsigned short,3> GetNewtonKrylovIntParam() const { return NK_IntParam; }
+  array<unsigned short,4> GetNewtonKrylovIntParam() const { return NK_IntParam; }
 
   /*!
    * \brief Get Newton-Krylov floating-point parameters.

Common/include/linear_algebra/CSysSolve.hpp

@@ -54,6 +54,15 @@ class CPreconditioner;
  *    Absolute tolerance, target residual is tol*||b||. ---*/
 enum class LinearToleranceType { RELATIVE, ABSOLUTE };
 
+/*!
+ * \brief Modes of using FGCRODR.
+ * \ingroup SpLinSys
+ */
+enum class FgcrodrMode {
+  NORMAL,   /*!< \brief Solve the linear system. */
+  SAME_MAT, /*!< \brief "NORMAL" but knowing the matrix did not change. */
+};
+
 /*!
  * \class CSysSolve
  * \ingroup SpLinSys
@@ -83,8 +92,8 @@ class CSysSolve {
   ScalarType Residual = 1e-20;  /*!< \brief Residual at the end of a call to Solve or Solve_b. */
   unsigned long Iterations = 0; /*!< \brief Iterations done in Solve or Solve_b. */
 
-  LINEAR_SOLVER_MODE
-  lin_sol_mode; /*!< \brief Type of operation for the linear system solver, changes the source of solver options. */
+  /*!< \brief Type of operation for the linear system solver, changes the source of solver options. */
+  LINEAR_SOLVER_MODE lin_sol_mode;
 
   mutable bool cg_ready;     /*!< \brief Indicate if memory used by CG is allocated. */
   mutable bool bcg_ready;    /*!< \brief Indicate if memory used by BCGSTAB is allocated. */
@@ -98,23 +107,24 @@ 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 std::vector<VectorType> W; /*!< \brief Large matrix used by FGMRES, w^i+1 = A * z^i. */
-  mutable std::vector<VectorType> Z; /*!< \brief Large matrix used by FGMRES, preconditioned W. */
-
-  VectorType
-      LinSysSol_tmp; /*!< \brief Temporary used when it is necessary to interface between active and passive types. */
-  VectorType
-      LinSysRes_tmp; /*!< \brief Temporary used when it is necessary to interface between active and passive types. */
-  VectorType*
-      LinSysSol_ptr; /*!< \brief Pointer to appropriate LinSysSol (set to original or temporary in call to Solve). */
-  const VectorType*
-      LinSysRes_ptr; /*!< \brief Pointer to appropriate LinSysRes (set to original or temporary in call to Solve). */
-
-  LinearToleranceType tol_type =
-      LinearToleranceType::ABSOLUTE; /*!< \brief How the linear solvers interpret the tolerance. */
-  bool xIsZero = false;              /*!< \brief If true assume the initial solution is always 0. */
-  bool recomputeRes = false;         /*!< \brief Recompute the residual after inner iterations, if monitoring. */
-  unsigned long monitorFreq = 10;    /*!< \brief Monitoring frequency. */
+  mutable unsigned long k = 0;
+  mutable std::vector<VectorType> Z, V; /*!< \brief Large matrices used by FGMRES, v^i+1 = A * z^i. */
+  mutable std::vector<VectorType> W, T; /*!< \brief Large matrices used by FGCRODR for deflation vectors. */
+
+  /*!< \brief Temporary used when it is necessary to interface between active and passive types. */
+  VectorType LinSysSol_tmp;
+  /*!< \brief Temporary used when it is necessary to interface between active and passive types. */
+  VectorType LinSysRes_tmp;
+  /*!< \brief Pointer to appropriate LinSysSol (set to original or temporary in call to Solve). */
+  VectorType* LinSysSol_ptr;
+  /*!< \brief Pointer to appropriate LinSysRes (set to original or temporary in call to Solve). */
+  const VectorType* LinSysRes_ptr;
+
+  /*!< \brief How the linear solvers interpret the tolerance. */
+  mutable LinearToleranceType tol_type = LinearToleranceType::ABSOLUTE;
+  mutable bool xIsZero = false;   /*!< \brief If true assume the initial solution is always 0. */
+  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;
@@ -225,72 +235,65 @@ class CSysSolve {
 
   /*!
    * \brief Used by Solve for compatibility between passive and active CSysVector.
-   * \note Same type specialization, temporary variables are not required.
-   * \param[in] LinSysRes - Linear system residual
-   * \param[in,out] LinSysSol - Linear system solution
-   */
-  template <class OtherType, su2enable_if<std::is_same<ScalarType, OtherType>::value> = 0>
-  void HandleTemporariesIn(const CSysVector<OtherType>& LinSysRes, CSysVector<OtherType>& LinSysSol) {
-    /*--- Set the pointers. ---*/
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      LinSysRes_ptr = &LinSysRes;
-      LinSysSol_ptr = &LinSysSol;
-    }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
-  }
-
-  /*!
-   * \brief Used by Solve for compatibility between passive and active CSysVector.
-   * \note Different type specialization, copy data into temporary solution and residual vectors.
    * \param[in] LinSysRes - Linear system residual
    * \param[in,out] LinSysSol - Linear system solution
    */
-  template <class OtherType, su2enable_if<!std::is_same<ScalarType, OtherType>::value> = 0>
+  template <class OtherType>
   void HandleTemporariesIn(const CSysVector<OtherType>& LinSysRes, CSysVector<OtherType>& LinSysSol) {
-    /*--- Copy data, the solution is also copied as it serves as initial condition. ---*/
-    LinSysRes_tmp.PassiveCopy(LinSysRes);
-    LinSysSol_tmp.PassiveCopy(LinSysSol);
-
-    /*--- Set the pointers. ---*/
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      LinSysRes_ptr = &LinSysRes_tmp;
-      LinSysSol_ptr = &LinSysSol_tmp;
+    if constexpr (std::is_same_v<ScalarType, OtherType>) {
+      /*--- Same type specialization, temporary variables are not required. ---*/
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+        LinSysRes_ptr = &LinSysRes;
+        LinSysSol_ptr = &LinSysSol;
+      }
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
+    } else {
+      /*--- Copy data, the solution is also copied as it serves as initial condition. ---*/
+      LinSysRes_tmp.PassiveCopy(LinSysRes);
+      LinSysSol_tmp.PassiveCopy(LinSysSol);
+
+      /*--- Set the pointers. ---*/
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+        LinSysRes_ptr = &LinSysRes_tmp;
+        LinSysSol_ptr = &LinSysSol_tmp;
+      }
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
     }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
   /*!
    * \brief Used by Solve for compatibility between passive and active CSysVector.
-   * \note Same type specialization, temporary variables are not required.
    * \param[out] LinSysSol - Linear system solution
    */
-  template <class OtherType, su2enable_if<std::is_same<ScalarType, OtherType>::value> = 0>
+  template <class OtherType>
   void HandleTemporariesOut(CSysVector<OtherType>& LinSysSol) {
-    /*--- Reset the pointers. ---*/
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      LinSysRes_ptr = nullptr;
-      LinSysSol_ptr = nullptr;
+    if constexpr (std::is_same_v<ScalarType, OtherType>) {
+      /*--- Same type specialization, temporary variables are not required. ---*/
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+        LinSysRes_ptr = nullptr;
+        LinSysSol_ptr = nullptr;
+      }
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
+    } else {
+      /*--- Copy data, only the temporary solution needs to be copied. ---*/
+      LinSysSol.PassiveCopy(LinSysSol_tmp);
+
+      /*--- Reset the pointers. ---*/
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+        LinSysRes_ptr = nullptr;
+        LinSysSol_ptr = nullptr;
+      }
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
     }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
-  /*!
-   * \brief Used by Solve for compatibility between passive and active CSysVector.
-   * \note Different type specialization, copy data from the temporary solution vector.
-   * \param[out] LinSysSol - Linear system solution
-   */
-  template <class OtherType, su2enable_if<!std::is_same<ScalarType, OtherType>::value> = 0>
-  void HandleTemporariesOut(CSysVector<OtherType>& LinSysSol) {
-    /*--- Copy data, only the temporary solution needs to be copied. ---*/
-    LinSysSol.PassiveCopy(LinSysSol_tmp);
-
-    /*--- Reset the pointers. ---*/
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      LinSysRes_ptr = nullptr;
-      LinSysSol_ptr = nullptr;
-    }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
-  }
+  /*--- TODO(pedro): The deflation part using Eigen does not compile in forward AD mode.
+   * So we need a dummy template to avoid instantiating this function for directdiff. ---*/
+  template <class Dummy = int>
+  unsigned long FGCRODR_LinSolverImpl(const VectorType& b, VectorType& x, const ProductType& mat_vec,
+                                      const PrecondType& precond, ScalarType tol, unsigned long max_iter,
+                                      ScalarType& residual, bool monitoring, const CConfig* config,
+                                      FgcrodrMode mode) const;
 
  public:
   /*!
@@ -335,7 +338,25 @@ class CSysSolve {
    */
   unsigned long RFGMRES_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,
                                   const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual,
-                                  bool monitoring, const CConfig* config);
+                                  bool monitoring, const CConfig* config) const;
+
+  /*!
+   * \brief Flexible Generalized Conjugate Residual Method with Inner Orthogonalization and Deflated Restarting.
+   * \param[in] b - the right hand size vector
+   * \param[in,out] x - on entry the intial guess, on exit the solution
+   * \param[in] mat_vec - object that defines matrix-vector product
+   * \param[in] precond - object that defines preconditioner
+   * \param[in] tol - tolerance with which to solve the system
+   * \param[in] max_iter - maximum number of iterations
+   * \param[out] residual - final normalized residual
+   * \param[in] monitoring - turn on priting residuals from solver to screen.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] mode - See FgcrodrMode.
+   */
+  unsigned long FGCRODR_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,
+                                  const PrecondType& precond, ScalarType tol, unsigned long max_iter,
+                                  ScalarType& residual, bool monitoring, const CConfig* config,
+                                  FgcrodrMode mode = FgcrodrMode::NORMAL) const;
 
   /*!
    * \brief Biconjugate Gradient Stabilized Method (BCGSTAB)
@@ -390,7 +411,7 @@ class CSysSolve {
    * \param[in] directCall - If this method is called directly, or in AD context.
    */
   unsigned long Solve_b(MatrixType& Jacobian, const CSysVector<su2double>& LinSysRes, CSysVector<su2double>& LinSysSol,
-                        CGeometry* geometry, const CConfig* config, const bool directCall = true);
+                        CGeometry* geometry, const CConfig* config, bool directCall = true);
 
   /*!
    * \brief Get the number of iterations.
@@ -423,4 +444,9 @@ class CSysSolve {
    * \brief Set the screen output frequency during monitoring.
    */
   inline void SetMonitoringFrequency(bool frequency) { monitorFreq = frequency; }
+
+  /*!
+   * \brief Discard FGCRODR's deflation vectors for the next solve.
+   */
+  inline void ResetDeflation() const { k = 0; }
 };

Common/include/linear_algebra/CSysVector.hpp

@@ -66,6 +66,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
  private:
   enum { OMP_MAX_SIZE = 4096 }; /*!< \brief Maximum chunk size used in parallel for loops. */
 
+  /// NOTE: Update swap() if you add member variables.
   unsigned long omp_chunk_size = OMP_MAX_SIZE; /*!< \brief Static chunk size used in loops. */
   ScalarType* vec_val = nullptr;               /*!< \brief Storage, 64 byte aligned (do not use normal new/delete). */
   unsigned long nElm = 0;       /*!< \brief Total number of elements (or number elements on this processor). */
@@ -155,6 +156,18 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    */
   CSysVector(const CSysVector& u) { Initialize(u.GetNBlk(), u.GetNBlkDomain(), u.nVar, u.vec_val, true); }
 
+  /*!
+   * \brief Swap contents with another vector.
+   */
+  void swap(CSysVector& other) {
+    std::swap(omp_chunk_size, other.omp_chunk_size);
+    std::swap(vec_val, other.vec_val);
+    std::swap(d_vec_val, other.d_vec_val);
+    std::swap(nElm, other.nElm);
+    std::swap(nElmDomain, other.nElmDomain);
+    std::swap(nVar, other.nVar);
+  }
+
   /*!
    * \brief Initialize the class with a scalar.
    * \param[in] numBlk - number of blocks locally

Common/include/option_structure.hpp

@@ -2352,9 +2352,10 @@ static const MapType<std::string, ENUM_FFD_BLENDING> Blending_Map = {
  */
 enum ENUM_LINEAR_SOLVER {
   CONJUGATE_GRADIENT,   /*!< \brief Preconditionated conjugate gradient method for grid deformation. */
-  FGMRES,               /*!< \brief Flexible Generalized Minimal Residual method. */
   BCGSTAB,              /*!< \brief BCGSTAB - Biconjugate Gradient Stabilized Method (main solver). */
+  FGMRES,               /*!< \brief Flexible Generalized Minimal Residual method. */
   RESTARTED_FGMRES,     /*!< \brief Flexible Generalized Minimal Residual method with restart. */
+  FGCRODR,              /*!< \brief Flexible Generalized Conjugate Residual Method with Inner Orthogonalization and Deflated Restarting. */
   SMOOTHER,             /*!< \brief Iterative smoother. */
   PASTIX_LDLT,          /*!< \brief PaStiX LDLT (complete) factorization. */
   PASTIX_LU,            /*!< \brief PaStiX LU (complete) factorization. */
@@ -2364,6 +2365,7 @@ static const MapType<std::string, ENUM_LINEAR_SOLVER> Linear_Solver_Map = {
   MakePair("BCGSTAB", BCGSTAB)
   MakePair("FGMRES", FGMRES)
   MakePair("RESTARTED_FGMRES", RESTARTED_FGMRES)
+  MakePair("FGCRODR", FGCRODR)
   MakePair("SMOOTHER", SMOOTHER)
   MakePair("PASTIX_LDLT", PASTIX_LDLT)
   MakePair("PASTIX_LU", PASTIX_LU)

Common/src/CConfig.cpp

@@ -1883,6 +1883,8 @@ void CConfig::SetConfig_Options() {
   addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", Linear_Solver_ILU_n, 0);
   /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */
   addUnsignedLongOption("LINEAR_SOLVER_RESTART_FREQUENCY", Linear_Solver_Restart_Frequency, 10);
+  /* DESCRIPTION: Number of vectors used for deflated restarts */
+  addUnsignedLongOption("LINEAR_SOLVER_RESTART_DEFLATION", Linear_Solver_Restart_Deflation, 4);
   /* DESCRIPTION: Relaxation factor for iterative linear smoothers (SMOOTHER_ILU/JACOBI/LU-SGS/LINELET) */
   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"). */
@@ -7264,10 +7266,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
               if (Kind_Linear_Solver == BCGSTAB) {
                 cout << "BCGSTAB is used for solving the linear system." << endl;
               } else {
+                const std::string name = Kind_Linear_Solver == FGCRODR ? "FGCRODR" : "FGMRES";
                 if (Kind_Linear_Solver_Inner == LINEAR_SOLVER_INNER::BCGSTAB){
-                  cout << "Nested FGMRES (FGMRES with inner BiCGSTAB) is used for solving the linear system." << endl;
+                  cout << "Nested " << name << " (with inner BiCGSTAB) is used for solving the linear system." << endl;
                 } else {
-                  cout << "FGMRES is used for solving the linear system." << endl;
+                  cout << name << " is used for solving the linear system." << endl;
                 }
               }
               switch (Kind_Linear_Solver_Prec) {
@@ -7319,6 +7322,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
               if (Kind_Linear_Solver_Inner == LINEAR_SOLVER_INNER::BCGSTAB)
                 cout << "Nested BiCGSTAB is used as the inner solver." << endl;
               break;
+            case FGCRODR:
+              cout << "FGCRODR 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;
+              break;
             case CONJUGATE_GRADIENT:
               cout << "A Conjugate Gradient method is used for solving the linear system." << endl;
               cout << "Convergence criteria of the linear solver: "<< Linear_Solver_Error <<"."<< endl;