Port ILU to GPU

Common/include/CConfig.hpp

@@ -646,6 +646,9 @@ class CConfig {
   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. */
   bool Linear_Solver_ILU_levels;                 /*!< \brief Use level scheduling for OMP parallelization of ILU. */
+  /*!< \brief Colored-iterative sweep counts for the GPU ILU preconditioner: [0] builds the
+   * factorization (Gauss-Seidel), [1]/[2] apply it (Jacobi, forward/backward triangular solve). */
+  array<unsigned short, 3> Linear_Solver_ILU_GPU_Sweeps{{1, 2, 2}};
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
   su2double MSW_Alpha;                  /*!< \brief Coefficient for blending states in the MSW scheme. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
@@ -4400,6 +4403,12 @@ class CConfig {
    */
   bool GetLinear_Solver_ILU_levels(void) const { return Linear_Solver_ILU_levels; }
 
+  /*!
+   * \brief Get the [build, forward, backward] colored-iterative sweep counts for the GPU ILU
+   *        preconditioner, see Linear_Solver_ILU_GPU_Sweeps.
+   */
+  array<unsigned short, 3> GetLinear_Solver_ILU_GPU_Sweeps(void) const { return Linear_Solver_ILU_GPU_Sweeps; }
+
   /*!
    * \brief Get restart frequency of the linear solver for the implicit formulation.
    * \return Restart frequency of the linear solver for the implicit formulation.

Common/include/geometry/meshreader/CSU2BinaryMeshReaderBase.hpp

@@ -10,7 +10,7 @@
  * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
- * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
  *
  * SU2 is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public

Common/include/geometry/meshreader/CSU2BinaryMeshReaderFEM.hpp

@@ -10,7 +10,7 @@
  * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
- * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
  *
  * SU2 is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public

Common/include/geometry/meshreader/CSU2BinaryMeshReaderFVM.hpp

@@ -10,7 +10,7 @@
  * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
- * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
  *
  * SU2 is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public

Common/include/geometry/meshreader/CSU2MeshReaderBase.hpp

@@ -10,7 +10,7 @@
  * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
- * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
  *
  * SU2 is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public

Common/include/grid_movement/CVolumetricMovementFactory.hpp

@@ -8,7 +8,7 @@
  * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
- * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md)
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
  *
  * SU2 is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public

Common/include/linear_algebra/CMatrixInverse.hpp

@@ -0,0 +1,99 @@
+/*!
+ * \file CMatrixInverse.hpp
+ * \brief Dense small-matrix inversion via Gauss-Jordan elimination, shared between the host
+ *        (CSysMatrix::MatrixInverse) and device (CSysPreconditionerGPU.cu) implementations.
+ * \author F. Palacios, A. Bueno, T. Economon, P. Gomes
+ * \version 8.5.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <cmath>
+
+#ifdef __CUDACC__
+#define SU2_CUDA_HOST_DEVICE __host__ __device__
+#else
+#define SU2_CUDA_HOST_DEVICE
+#endif
+
+namespace SU2_LinAlg {
+
+/*!
+ * \brief Regularize a pivot that is too small to prevent divide-by-zero, on host and device
+ *        this needs to clamp to the same value so that the two produce the same factors.
+ */
+template <class ScalarType>
+SU2_CUDA_HOST_DEVICE inline void RegularizePivot(ScalarType& pivot) {
+  const float eps = 1e-12;
+#ifdef __CUDA_ARCH__
+  if (fabs(pivot) < eps) pivot = copysign(ScalarType(eps), pivot);
+#else
+  if (std::abs(pivot) < eps) pivot = std::copysign(ScalarType(eps), pivot);
+#endif
+}
+
+/*!
+ * \brief Invert the \p nVar by \p nVar dense matrix \p matrix into \p inverse via Gauss-Jordan
+ *        elimination with partial pivoting on the diagonal.
+ * \note \p matrix is used as scratch space and destroyed, \p inverse must not alias it.
+ */
+template <class ScalarType>
+SU2_CUDA_HOST_DEVICE inline void MatrixInverse(unsigned long nVar, ScalarType* matrix, ScalarType* inverse) {
+#define A(I, J) matrix[(I)*nVar + (J)]
+#define M(I, J) inverse[(I)*nVar + (J)]
+
+  /*--- Initialize the inverse with the identity. ---*/
+  for (auto iVar = 0ul; iVar < nVar; iVar++)
+    for (auto jVar = 0ul; jVar < nVar; jVar++) M(iVar, jVar) = ScalarType(iVar == jVar);
+
+  /*--- Transform system in Upper Matrix. ---*/
+  for (auto iVar = 1ul; iVar < nVar; iVar++) {
+    for (auto jVar = 0ul; jVar < iVar; jVar++) {
+      RegularizePivot(A(jVar, jVar));
+
+      const ScalarType weight = A(iVar, jVar) / A(jVar, jVar);
+      for (auto kVar = jVar; kVar < nVar; kVar++) A(iVar, kVar) -= weight * A(jVar, kVar);
+
+      /*--- At this stage M is lower triangular so not all cols need updating. ---*/
+      for (auto kVar = 0ul; kVar <= jVar; kVar++) M(iVar, kVar) -= weight * M(jVar, kVar);
+    }
+  }
+
+  /*--- Backwards substitution. ---*/
+  for (auto iVar = nVar; iVar > 0ul;) {
+    iVar--;  // unsigned type
+    for (auto jVar = iVar + 1; jVar < nVar; jVar++)
+      for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) -= A(iVar, jVar) * M(jVar, kVar);
+
+    RegularizePivot(A(iVar, iVar));
+
+    for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) /= A(iVar, iVar);
+  }
+
+#undef A
+#undef M
+}
+
+}  // namespace SU2_LinAlg
+
+#undef SU2_CUDA_HOST_DEVICE

Common/include/linear_algebra/CMatrixVectorProduct.hpp

@@ -105,28 +105,6 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
    * \param[out] v - CSysVector that is the result of the product
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    if (config->GetCUDA()) {
-#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);
-      }
-#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 "
-          "options enabled in Meson to access GPU Functions",
-          CURRENT_FUNCTION);
-#endif
-    } else {
-      matrix.MatrixVectorProduct(u, v, geometry, config);
-    }
+    matrix.MatrixVectorProduct(u, v, geometry, config);
   }
 };

Common/include/linear_algebra/CPreconditioner.hpp

@@ -40,7 +40,7 @@
 /*!
  * \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
+ * \note This is what keeps 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.
  */
@@ -199,7 +199,8 @@ 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 {
-    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config); });
+    /*--- No host bracket, ILU has a device implementation and ComputeILUPreconditioner dispatches to it. ---*/
+    sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config);
   }
 
   /*!