Add initial end-to-end CUDA FGMRES solver path (#2825)

## Proposed Changes

This PR adds an end-to-end CUDA linear solve path: the Krylov solvers
keep their host control flow, and `CSysVector` operations, the SpMV and
the Jacobi preconditioner are dispatched to CUDA kernels when
`ENABLE_CUDA=YES`.

Transfers are explicit and owned by the object responsible for the data,
with no coherency or dirty-flag tracking in `CSysVector` / `CSysMatrix`:
- `CSysMatrixVectorProduct` uploads the matrix on construction
- the preconditioner uploads its data in `Build()`
- `CSysSolve` uploads `b` and `x` and downloads `x` in
`HandleTemporariesIn/Out`, which is also where device evaluation is
switched on and off
- preconditioners without a device implementation (ILU, LU-SGS, Linelet,
PaStiX) download the input, apply on the host, and upload the result

A solve is therefore fully device resident for Identity and Jacobi
preconditioners: **2 uploads and 1 download per linear system**,
independent of the Krylov subspace size.

Implementation notes:
- `cuBLAS` for `dot` / `norm`, custom kernels for the block-LDU SpMV,
the Jacobi apply, and `CSysVector` expression assignment
- `CSysVector` operands are captured in expressions by value, so an
arbitrary expression tree is trivially copyable into the assignment
kernel; the required expression shapes are explicitly instantiated in
`CSysVectorGPU.cu`, consistent with how `CSysMatrix` is instantiated
- device work is issued by a single thread with the OpenMP team
synchronized around it, so the GPU path is usable from inside the
existing parallel regions; this is internal to the linear algebra layer
- the CUDA translation units are only linked into the primal libraries,
since they cannot be compiled with the CoDiPack defines; device dispatch
is compiled out of the AD builds
- adds `LINEAR_SOLVER_PREC= NONE` (identity)

## Related Work

Follows the review direction in #2822 (show a working end-to-end GPU
linear solve before splitting out infrastructure) and the implementation
preferences in #2816.

## Validation

- CPU vs GPU on inviscid NACA0012, 25 iterations, RTX 4070 Ti SUPER
(sm_89), for FGMRES + `JACOBI`, FGMRES + `NONE`, FGMRES + `ILU` (host
preconditioner path) and BCGSTAB. Results agree to the printed precision
in double, and to ~6 significant figures in mixed precision. BCGSTAB
agrees exactly once the linear system is converged
(`LINEAR_SOLVER_ERROR=1e-10`); at loose tolerances the two paths diverge
through BCGSTAB's own sensitivity, not a difference in the algebra.
- Mixed, normal and single precision builds all compile and run; device
dispatch confirmed live in each (kernel launch counts track the subspace
size).
- `OMP_NUM_THREADS=1` and `4` give bit-identical results on the GPU
path.
- Builds verified: primal, primal + AD + directdiff, `enable-cuda` +
`with-omp`, mixed / normal / single precision.
- Transfer counts measured per solve: 2 H2D + 1 D2H + 1 matrix upload
for Jacobi and `NONE`, plus one download/upload pair per preconditioner
application for ILU.

Earlier validation of the original design (6 representative cases,
`nsys` / `ncu` profiling) predates the rework of the transfer and
dispatch model and should be repeated.

## PR Checklist

- [x] I am submitting my contribution to the develop branch.
- [x] My contribution generates no new compiler warnings (try with
--warnlevel=3 when using meson).
- [x] My contribution is commented and consistent with SU2 style
(https://su2code.github.io/docs_v7/Style-Guide/).
- [x] I used the pre-commit hook to prevent dirty commits and used
`pre-commit run --all` to format old commits.
- [ ] I have added a test case that demonstrates my contribution, if
necessary.
- [ ] I have updated appropriate documentation (Tutorials, Docs Page,
config_template.cpp), if necessary.

---------

Co-authored-by: Pedro Gomes <pcarruscag@gmail.com>

Author: LwhJesse

AUTHORS.md

@@ -10,18 +10,18 @@ Thomas D. Economon    (formerly Stanford University)
 Juan J. Alonso                 (Stanford University)
 ```
 
-## Current Maintainer ##
+## Current Maintainers ##
 
 The SU2 project is maintained by members of the SU2 Foundation (https://su2foundation.org)
 
 ```
-Thomas D. Economon - Executive Director - tom@su2foundation.org
-Tim Albring        - Director           - tim@su2foundation.org
-Juan J. Alonso     - Director           - juan@su2foundation.org
-Eran Arad          - Director           - eran@su2foundation.org
-Piero Colonna      - Director           - piero@su2foundation.org
-Pedro Gomes        - Director           - pedro@su2foundation.org
-Daniel Mayer       - Director           - daniel@su2foundation.org
+Thomas D. Economon  - Chairperson
+Matteo Pini         - Vice Chairperson
+Nijso Beishuizen    - Treasurer
+Pedro Gomes         - Development Officer
+Giulio Gori         - Secretary
+Nitish Anand        - Editorial Officer
+Edwin van der Weide - Events Officer
 ```
 
 in collaboration with the following main contributors and research teams:
@@ -93,6 +93,7 @@ Jairo Paes Cavalcante Filho
 Jason Howison
 Jayant Mukhopadhaya
 Jeffrey van Oostrom
+Jesse Li
 Jessie Lauzon
 João Loureiro
 Johannes Blühdorn

Common/include/code_config.hpp

@@ -96,6 +96,16 @@ FORCEINLINE Out su2staticcast_p(In ptr) {
 #define HAVE_OMP
 #endif
 
+/*--- Detect whether the CUDA kernels are part of this build. The .cu translation units
+ * cannot be compiled with the CoDiPack defines (nvcc's device pass cannot parse the tape
+ * machinery), and an object compiled with a different definition of su2double must not be
+ * linked into an AD library. They are therefore only built into the primal libraries, and
+ * all device dispatch has to be compiled out of the AD builds, which HAVE_CUDA alone does
+ * not do because su2mixedfloat is a passive type there as well. ---*/
+#if defined(HAVE_CUDA) && !defined(CODI_REVERSE_TYPE) && !defined(CODI_FORWARD_TYPE)
+#define SU2_ENABLE_CUDA_KERNELS
+#endif
+
 /*--- No full single precision for AD builds. ---*/
 #if (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) && defined(USE_SINGLE_PRECISION)
 #undef USE_SINGLE_PRECISION

Common/include/linear_algebra/CMatrixVectorProduct.hpp

@@ -72,7 +72,6 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
   const CSysMatrix<ScalarType>& matrix; /*!< \brief pointer to matrix that defines the product. */
   CGeometry* geometry;                  /*!< \brief geometry associated with the matrix. */
   const CConfig* config;                /*!< \brief config of the problem. */
-  mutable bool matrix_uploaded = false; /*!< \brief Upload the matrix lazily on the first actual GPU matvec. */
 
  public:
   /*!
@@ -83,7 +82,17 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
    */
   inline CSysMatrixVectorProduct(const CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref,
                                  const CConfig* config_ref)
-      : matrix(matrix_ref), geometry(geometry_ref), config(config_ref) {}
+      : matrix(matrix_ref), geometry(geometry_ref), config(config_ref) {
+    /*--- The matrix does not change while this object lives, so it crosses the bus once,
+     * here. The vectors are uploaded by CSysSolve, see HandleTemporariesIn. ---*/
+#ifdef SU2_ENABLE_CUDA_KERNELS
+    if constexpr (su2_gpu_capable_v<ScalarType>) {
+      if (config->GetCUDA()) {
+        SU2_DEVICE_REGION(matrix.HtDTransfer();)
+      }
+    }
+#endif
+  }
 
   /*!
    * \note This class cannot be default constructed as that would leave us with invalid pointers.
@@ -97,12 +106,19 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
     if (config->GetCUDA()) {
-#ifdef HAVE_CUDA
-      if (!matrix_uploaded) {
-        matrix.HtDTransfer();
-        matrix_uploaded = true;
+#ifdef SU2_ENABLE_CUDA_KERNELS
+      if constexpr (su2_gpu_capable_v<ScalarType>) {
+        BEGIN_SU2_DEVICE_REGION
+        matrix.GPUMatrixVectorProduct(u, v, geometry, config);
+        END_SU2_DEVICE_REGION
+      } else {
+        SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION);
       }
-      matrix.GPUMatrixVectorProduct(u, v, geometry, config);
+#elif defined(HAVE_CUDA)
+      SU2_MPI::Error(
+          "\nError in launching Matrix-Vector Product Function\nENABLE_CUDA is set to YES\nThe GPU kernels are not "
+          "part of the AD libraries, use the primal build for GPU acceleration",
+          CURRENT_FUNCTION);
 #else
       SU2_MPI::Error(
           "\nError in launching Matrix-Vector Product Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA "

Common/include/linear_algebra/CPreconditioner.hpp

@@ -37,6 +37,33 @@
 /// \addtogroup SpLinSys
 /// @{
 
+/*!
+ * \brief Applies a preconditioner that only has a host implementation to vectors that live
+ * on the device: bring the input down, apply, put the result back.
+ * \note This is what keeps ILU, LU-SGS, Linelet and PaStiX usable on the GPU path. The
+ * transfers are issued by one thread with the team synchronized around them, the apply
+ * itself is the normal OpenMP parallel host code.
+ */
+template <class ScalarType, class Apply>
+inline void ApplyPreconditionerOnHost(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v, Apply&& apply) {
+#ifdef SU2_ENABLE_CUDA_KERNELS
+  if constexpr (su2_gpu_capable_v<ScalarType>) {
+    if (VecExpr::UseDeviceExpressions()) {
+      /*--- The host code must not see the device pointers of any expression it builds, so
+       * the switch is flipped for the duration of the apply. It is written inside the
+       * regions, by one thread, and published to the team by the trailing barrier. ---*/
+      SU2_DEVICE_REGION(u.DtHTransfer(); VecExpr::SetUseDeviceExpressions(false);)
+
+      apply();
+
+      SU2_DEVICE_REGION(VecExpr::SetUseDeviceExpressions(true); v.HtDTransfer();)
+      return;
+    }
+  }
+#endif
+  apply();
+}
+
 /*!
  * \class CPreconditioner
  * \brief Abstract base class for defining a preconditioning operation.
@@ -77,6 +104,18 @@ class CPreconditioner {
 template <class ScalarType>
 CPreconditioner<ScalarType>::~CPreconditioner() {}
 
+/*!
+ * \class CIdentityPreconditioner
+ * \brief No-op preconditioner used when Krylov solvers run without preconditioning.
+ */
+template <class ScalarType>
+class CIdentityPreconditioner final : public CPreconditioner<ScalarType> {
+ public:
+  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override { v = u; }
+
+  inline bool IsIdentity() const override { return true; }
+};
+
 /*!
  * \class CJacobiPreconditioner
  * \brief Specialization of preconditioner that uses CSysMatrix class.
@@ -160,7 +199,7 @@ class CILUPreconditioner final : public CPreconditioner<ScalarType> {
    * \param[out] v - CSysVector that is the result of the preconditioning.
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config);
+    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config); });
   }
 
   /*!
@@ -206,7 +245,7 @@ class CLU_SGSPreconditioner final : public CPreconditioner<ScalarType> {
    * \param[out] v - CSysVector that is the result of the preconditioning.
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config);
+    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); });
   }
 };
 
@@ -234,7 +273,7 @@ class CQuantizedLUSGSPreconditioner final : public CPreconditioner<ScalarType> {
   CQuantizedLUSGSPreconditioner() = delete;
 
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config);
+    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); });
   }
 
   /*! \brief Quantize the diagonal blocks (off diagonals are quantized on the fly). */
@@ -278,7 +317,7 @@ class CLineletPreconditioner final : public CPreconditioner<ScalarType> {
    * \param[out] v - CSysVector that is the result of the preconditioning.
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    sparse_matrix.ComputeLineletPreconditioner(u, v, geometry, config);
+    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLineletPreconditioner(u, v, geometry, config); });
   }
 
   /*!
@@ -328,7 +367,7 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
    * \param[out] v - CSysVector that is the result of the preconditioning.
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    sparse_matrix.ComputePastixPreconditioner(u, v, geometry, config);
+    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputePastixPreconditioner(u, v, geometry, config); });
   }
 
   /*!
@@ -363,6 +402,9 @@ CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOL
   CPreconditioner<ScalarType>* prec = nullptr;
 
   switch (kind) {
+    case IDENTITY:
+      prec = new CIdentityPreconditioner<ScalarType>();
+      break;
     case JACOBI:
       prec = new CJacobiPreconditioner<ScalarType>(jacobian, geometry, config);
       break;

Common/include/linear_algebra/CSysMatrix.hpp

@@ -247,9 +247,10 @@ class CSysMatrix {
     unsigned long nnz_u = 0;            /*!< \brief Number of U nonzeros. */
   };
 
-  LDU mat; /*!< \brief Host matrix (values owned via aligned_alloc; pattern from geometry). */
-  LDU gpu; /*!< \brief Device matrix (all pointers to GPU memory). */
-  LDU ilu; /*!< \brief ILU factorization, host (values owned; pattern from geometry). */
+  LDU mat;                      /*!< \brief Host matrix (values owned via aligned_alloc; pattern from geometry). */
+  LDU gpu;                      /*!< \brief Device matrix (all pointers to GPU memory). */
+  LDU ilu;                      /*!< \brief ILU factorization, host (values owned; pattern from geometry). */
+  ScalarType* d_invM = nullptr; /*!< \brief Device inverse diagonal blocks for the Jacobi preconditioner. */
 
   /*--- Quantized off-diagonal storage (used when quantized_mode == true). ---*/
   using QuantType = int8_t;
@@ -1100,6 +1101,14 @@ class CSysMatrix {
   void ComputeJacobiPreconditioner(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod, CGeometry* geometry,
                                    const CConfig* config) const;
 
+  /*!
+   * \brief Apply the Jacobi preconditioner on the GPU/device side.
+   * \note This helper is intended as the implementation hook for GPU-resident Krylov solvers.
+   *       The actual implementation belongs in CSysMatrixGPU.cu.
+   */
+  void ComputeJacobiPreconditionerGPU(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod,
+                                      CGeometry* geometry, const CConfig* config) const;
+
   /*!
    * \brief Build the ILU preconditioner.
    */

Common/include/linear_algebra/CSysSolve.hpp

@@ -240,13 +240,50 @@ class CSysSolve {
    */
   void WriteWarning(ScalarType res_calc, ScalarType res_true, ScalarType tol) const;
 
+  /*!
+   * \brief Moves the linear system to the device, if the GPU path is in use.
+   * \note This and DownloadSolution are the only places where b and x cross the bus. The
+   * work vectors of the solvers never do, they are allocated on the device and read back
+   * only through the reductions.
+   */
+  void UploadSystem(bool useCuda) const {
+#ifdef SU2_ENABLE_CUDA_KERNELS
+    if constexpr (su2_gpu_capable_v<ScalarType>) {
+      if (!useCuda) return;
+      BEGIN_SU2_DEVICE_REGION {
+        LinSysRes_ptr->HtDTransfer();
+        LinSysSol_ptr->HtDTransfer();
+        VecExpr::SetUseDeviceExpressions(true);
+      }
+      END_SU2_DEVICE_REGION
+    }
+#endif
+  }
+
+  /*!
+   * \brief Brings the solution back from the device and returns to host evaluation.
+   */
+  void DownloadSolution(bool useCuda) const {
+#ifdef SU2_ENABLE_CUDA_KERNELS
+    if constexpr (su2_gpu_capable_v<ScalarType>) {
+      if (!useCuda) return;
+      BEGIN_SU2_DEVICE_REGION {
+        LinSysSol_ptr->DtHTransfer();
+        VecExpr::SetUseDeviceExpressions(false);
+      }
+      END_SU2_DEVICE_REGION
+    }
+#endif
+  }
+
   /*!
    * \brief Used by Solve for compatibility between passive and active CSysVector.
    * \param[in] LinSysRes - Linear system residual
    * \param[in,out] LinSysSol - Linear system solution
+   * \param[in] useCuda - Whether to move the system to the device for the solve.
    */
   template <class OtherType>
-  void HandleTemporariesIn(const CSysVector<OtherType>& LinSysRes, CSysVector<OtherType>& LinSysSol) {
+  void HandleTemporariesIn(const CSysVector<OtherType>& LinSysRes, CSysVector<OtherType>& LinSysSol, bool useCuda) {
     SU2_ZONE_SCOPED
     if constexpr (std::is_same_v<ScalarType, OtherType>) {
       /*--- Same type specialization, temporary variables are not required. ---*/
@@ -255,6 +292,7 @@ class CSysSolve {
         LinSysSol_ptr = &LinSysSol;
       }
       END_SU2_OMP_SAFE_GLOBAL_ACCESS
+      UploadSystem(useCuda);
     } else {
       /*--- Copy data, the solution is also copied as it serves as initial condition. ---*/
       LinSysRes_tmp.PassiveCopy(LinSysRes);
@@ -266,16 +304,19 @@ class CSysSolve {
         LinSysSol_ptr = &LinSysSol_tmp;
       }
       END_SU2_OMP_SAFE_GLOBAL_ACCESS
+      UploadSystem(useCuda);
     }
   }
 
   /*!
    * \brief Used by Solve for compatibility between passive and active CSysVector.
    * \param[out] LinSysSol - Linear system solution
+   * \param[in] useCuda - Whether the system was solved on the device.
    */
   template <class OtherType>
-  void HandleTemporariesOut(CSysVector<OtherType>& LinSysSol) {
+  void HandleTemporariesOut(CSysVector<OtherType>& LinSysSol, bool useCuda) {
     SU2_ZONE_SCOPED
+    DownloadSolution(useCuda);
     if constexpr (std::is_same_v<ScalarType, OtherType>) {
       /*--- Same type specialization, temporary variables are not required. ---*/
       BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
@@ -432,9 +473,9 @@ class CSysSolve {
   template <class OtherType, su2enable_if<!std::is_same_v<ScalarType, OtherType>> = 0>
   unsigned long Solve_b(MatrixType& Jacobian, const CSysVector<OtherType>& LinSysRes, CSysVector<OtherType>& LinSysSol,
                         CGeometry* geometry, const CConfig* config, bool directCall = true) {
-    HandleTemporariesIn(LinSysRes, LinSysSol);
+    HandleTemporariesIn(LinSysRes, LinSysSol, false);
     auto iter = Solve_b(Jacobian, *LinSysRes_ptr, *LinSysSol_ptr, geometry, config, directCall);
-    HandleTemporariesOut(LinSysSol);
+    HandleTemporariesOut(LinSysSol, false);
     return iter;
   }