diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index af7aa99033b5..1d98b92534e7 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -132,6 +132,7 @@ class CConfig {
Sens_Remove_Sharp, /*!< \brief Flag for removing or not the sharp edges from the sensitivity computation. */
Hold_GridFixed, /*!< \brief Flag hold fixed some part of the mesh during the deformation. */
Axisymmetric, /*!< \brief Flag for axisymmetric calculations */
+ Enable_Cuda, /*!< \brief Flag for switching GPU computing*/
Integrated_HeatFlux; /*!< \brief Flag for heat flux BC whether it deals with integrated values.*/
su2double Buffet_k; /*!< \brief Sharpness coefficient for buffet sensor.*/
su2double Buffet_lambda; /*!< \brief Offset parameter for buffet sensor.*/
@@ -6192,6 +6193,12 @@ class CConfig {
*/
bool GetAxisymmetric(void) const { return Axisymmetric; }
+ /*!
+ * \brief Get information about GPU support.
+ * \return TRUE if cuda is enabled; otherwise FALSE.
+ */
+ bool GetCUDA(void) const { return Enable_Cuda; }
+
/*!
* \brief Subtract one to the index of the finest grid (full multigrid strategy).
* \return Change the index of the finest grid.
diff --git a/Common/include/linear_algebra/CMatrixVectorProduct.hpp b/Common/include/linear_algebra/CMatrixVectorProduct.hpp
index 586e8abe410b..83dc8e9c7025 100644
--- a/Common/include/linear_algebra/CMatrixVectorProduct.hpp
+++ b/Common/include/linear_algebra/CMatrixVectorProduct.hpp
@@ -51,6 +51,7 @@
* passed to a single implementation of the Krylov solvers.
* This abstraction may also be used to define matrix-free products.
*/
+
template
class CMatrixVectorProduct {
public:
@@ -94,6 +95,17 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct {
* \param[out] v - CSysVector that is the result of the product
*/
inline void operator()(const CSysVector& u, CSysVector& v) const override {
- matrix.MatrixVectorProduct(u, v, geometry, config);
+ if (config->GetCUDA()) {
+#ifdef HAVE_CUDA
+ matrix.GPUMatrixVectorProduct(u, v, geometry, config);
+#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);
+ }
}
};
diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp
index b026165e9e60..7d0e87b92e35 100644
--- a/Common/include/linear_algebra/CSysMatrix.hpp
+++ b/Common/include/linear_algebra/CSysMatrix.hpp
@@ -145,6 +145,12 @@ class CSysMatrix {
const unsigned long* col_ind; /*!< \brief Column index for each of the elements in val(). */
const unsigned long* col_ptr; /*!< \brief The transpose of col_ind, pointer to blocks with the same column index. */
+ ScalarType* d_matrix; /*!< \brief Device Pointer to store the matrix values on the GPU. */
+ const unsigned long* d_row_ptr; /*!< \brief Device Pointers to the first element in each row. */
+ const unsigned long* d_col_ind; /*!< \brief Device Column index for each of the elements in val(). */
+ bool useCuda; /*!< \brief Boolean that indicates whether user has enabled CUDA or not.
+ Mainly used to conditionally free GPU memory in the class destructor. */
+
ScalarType* ILU_matrix; /*!< \brief Entries of the ILU sparse matrix. */
unsigned long nnz_ilu; /*!< \brief Number of possible nonzero entries in the matrix (ILU). */
const unsigned long* row_ptr_ilu; /*!< \brief Pointers to the first element in each row (ILU). */
@@ -391,6 +397,12 @@ class CSysMatrix {
*/
void SetValDiagonalZero(void);
+ /*!
+ * \brief Performs the memory copy from host to device.
+ * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default.
+ */
+ void HtDTransfer(bool trigger = true) const;
+
/*!
* \brief Get a pointer to the start of block "ij"
* \param[in] block_i - Row index.
@@ -838,6 +850,49 @@ class CSysMatrix {
void MatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry,
const CConfig* config) const;
+ /*!
+ * \brief Performs the product of a sparse matrix by a CSysVector.
+ * \param[in] vec - CSysVector to be multiplied by the sparse matrix A.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[out] prod - Result of the product.
+ */
+ void GPUMatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry,
+ const CConfig* config) const;
+
+ /*!
+ * \brief Performs first step of the LU_SGS Preconditioner building
+ * \param[in] vec - CSysVector to be multiplied by the sparse matrix A.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[out] prod - Result of the product.
+ */
+ void GPUFirstSymmetricIteration(ScalarType& vec, ScalarType& prod, CGeometry* geometry, const CConfig* config) const;
+
+ /*!
+ * \brief Performs second step of the LU_SGS Preconditioner building
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[out] prod - Result of the product.
+ */
+ void GPUSecondSymmetricIteration(ScalarType& prod, CGeometry* geometry, const CConfig* config) const;
+
+ /*!
+ * \brief Performs Gaussian Elimination between diagional blocks of the matrix and the prod vector
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[out] prod - Result of the product.
+ */
+ void GPUGaussElimination(ScalarType& prod, CGeometry* geometry, const CConfig* config) const;
+
+ /*!
+ * \brief Multiply CSysVector by the preconditioner all of which are stored on the device
+ * \param[in] vec - CSysVector to be multiplied by the preconditioner.
+ * \param[out] prod - Result of the product A*vec.
+ */
+ void GPUComputeLU_SGSPreconditioner(ScalarType& vec, ScalarType& prod, CGeometry* geometry,
+ const CConfig* config) const;
+
/*!
* \brief Build the Jacobi preconditioner.
*/
diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp
index 7b0248b28b6b..7cec5b1b4e48 100644
--- a/Common/include/linear_algebra/CSysVector.hpp
+++ b/Common/include/linear_algebra/CSysVector.hpp
@@ -32,6 +32,7 @@
#include "../parallelization/omp_structure.hpp"
#include "../parallelization/vectorization.hpp"
#include "vector_expressions.hpp"
+#include "../../include/CConfig.hpp"
/*!
* \brief OpenMP worksharing construct used in CSysVector for loops.
@@ -71,6 +72,8 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType>
unsigned long nElmDomain = 0; /*!< \brief Total number of elements without Ghost cells. */
unsigned long nVar = 1; /*!< \brief Number of elements in a block. */
+ ScalarType* d_vec_val = nullptr; /*!< \brief Device Pointer to store the vector values on the GPU. */
+
/*!
* \brief Generic initialization from a scalar or array.
* \note If val==nullptr vec_val is not initialized, only allocated.
@@ -195,6 +198,29 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType>
END_CSYSVEC_PARFOR
}
+ /*!
+ * \brief Performs the memory copy from host to device.
+ * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default.
+ */
+ void HtDTransfer(bool trigger = true) const;
+
+ /*!
+ * \brief Performs the memory copy from device to host.
+ * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default.
+ */
+ void DtHTransfer(bool trigger = true) const;
+
+ /*!
+ * \brief Sets all the elements of the GPU vector to a certain value
+ * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default.
+ */
+ void GPUSetVal(ScalarType val, bool trigger = true) const;
+
+ /*!
+ * \brief return device pointer that points to the CSysVector values in GPU memory
+ */
+ inline ScalarType* GetDevicePointer() const { return d_vec_val; }
+
/*!
* \brief return the number of local elements in the CSysVector
*/
diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh
new file mode 100644
index 000000000000..05f40ee4d9b0
--- /dev/null
+++ b/Common/include/linear_algebra/GPUComms.cuh
@@ -0,0 +1,53 @@
+/*!
+\file GPUComms.cuh
+* \brief Header file containing universal functions that provide basic and essential utilities for other GPU processes
+* \author A. Raj
+* \version 8.2.0 "Harrier"
+*
+* SU2 Project Website: https://su2code.github.io
+*
+* The SU2 Project is maintained by the SU2 Foundation
+* (http://su2foundation.org)
+*
+* Copyright 2012-2024, 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 .
+*/
+
+#include
+#include
+
+namespace KernelParameters{
+
+inline constexpr int round_up_division(const int multiple, int x) { return ((x + multiple - 1) / multiple); }
+
+static constexpr int MVP_BLOCK_SIZE = 1024;
+static constexpr int MVP_WARP_SIZE = 32;
+}
+/*!
+* \brief assert style function that reads return codes after intercepting CUDA API calls.
+* It returns the result code and its location if the call is unsuccessful.
+* \param[in] code - result code of CUDA function
+* \param[in] file - name of file holding the function
+* \param[in] line - line containing the function
+*/
+
+inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true){
+ if (code != cudaSuccess){
+ fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
+ if (abort) exit(code);
+ }
+}
+
+#define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
diff --git a/Common/include/toolboxes/allocation_toolbox.hpp b/Common/include/toolboxes/allocation_toolbox.hpp
index 9caafbe151ec..74158670615d 100644
--- a/Common/include/toolboxes/allocation_toolbox.hpp
+++ b/Common/include/toolboxes/allocation_toolbox.hpp
@@ -36,6 +36,10 @@
#include
#endif
+#ifdef HAVE_CUDA
+#include "../linear_algebra/GPUComms.cuh"
+#endif
+
#include
#include
@@ -90,3 +94,52 @@ inline void aligned_free(T* ptr) noexcept {
}
} // namespace MemoryAllocation
+
+namespace GPUMemoryAllocation {
+/*!
+ * \brief Memory allocation for variables on the GPU.
+ * \param[in] size in bytes.
+ * \tparam ZeroInit, initialize memory to 0.
+ * \return Pointer to memory, always use gpu_free to deallocate.
+ */
+template
+inline T* gpu_alloc(size_t size) noexcept {
+ void* ptr = nullptr;
+
+#if defined(HAVE_CUDA)
+ gpuErrChk(cudaMalloc((void**)(&ptr), size));
+ if (ZeroInit) gpuErrChk(cudaMemset((void*)(ptr), 0.0, size));
+#else
+ return 0;
+#endif
+
+ return static_cast(ptr);
+}
+
+/*!
+ * \brief Free memory allocated on the GPU with gpu_alloc.
+ * \param[in] ptr, pointer to memory we want to release.
+ */
+template
+inline void gpu_free(T* ptr) noexcept {
+#ifdef HAVE_CUDA
+ gpuErrChk(cudaFree((void*)ptr));
+#endif
+}
+/*!
+ * \brief Memory allocation for variables on the GPU along with initialization from a source host array.
+ * \param[in] size in bytes.
+ * \return Pointer to memory, always use gpu_free to deallocate.
+ */
+template
+inline T* gpu_alloc_cpy(const T* src_ptr, size_t size) noexcept {
+ void* ptr = nullptr;
+
+#ifdef HAVE_CUDA
+ gpuErrChk(cudaMalloc((void**)(&ptr), size));
+ gpuErrChk(cudaMemcpy((void*)(ptr), (void*)src_ptr, size, cudaMemcpyHostToDevice));
+#endif
+
+ return static_cast(ptr);
+}
+} // namespace GPUMemoryAllocation
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 5fb7220e94e4..9222ed7f8928 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1149,6 +1149,8 @@ void CConfig::SetConfig_Options() {
/*\brief AXISYMMETRIC \n DESCRIPTION: Axisymmetric simulation \n DEFAULT: false \ingroup Config */
addBoolOption("AXISYMMETRIC", Axisymmetric, false);
+ /*\brief ENABLE_CUDA \n DESCRIPTION: GPU Acceleration \n DEFAULT: false \ingroup Config */
+ addBoolOption("ENABLE_CUDA", Enable_Cuda, false);
/* DESCRIPTION: Add the gravity force */
addBoolOption("GRAVITY_FORCE", GravityForce, false);
/* DESCRIPTION: Add the Vorticity Confinement term*/
diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp
index 92de9e3d704b..c0a6aa81ee04 100644
--- a/Common/src/linear_algebra/CSysMatrix.cpp
+++ b/Common/src/linear_algebra/CSysMatrix.cpp
@@ -68,6 +68,12 @@ CSysMatrix::~CSysMatrix() {
MemoryAllocation::aligned_free(matrix);
MemoryAllocation::aligned_free(invM);
+ if (useCuda) {
+ GPUMemoryAllocation::gpu_free(d_matrix);
+ GPUMemoryAllocation::gpu_free(d_row_ptr);
+ GPUMemoryAllocation::gpu_free(d_col_ind);
+ }
+
#ifdef USE_MKL
mkl_jit_destroy(MatrixMatrixProductJitter);
mkl_jit_destroy(MatrixVectorProductJitterBetaZero);
@@ -131,6 +137,30 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi
col_ind = csr.innerIdx();
dia_ptr = csr.diagPtr();
+ /*--- Allocate data. ---*/
+ auto allocAndInit = [](ScalarType*& ptr, unsigned long num) {
+ ptr = MemoryAllocation::aligned_alloc(64, num * sizeof(ScalarType));
+ };
+
+ allocAndInit(matrix, nnz * nVar * nEqn);
+
+ useCuda = config->GetCUDA();
+
+ if (useCuda) {
+ /*--- Allocate GPU data. ---*/
+ auto GPUAllocAndInit = [](ScalarType*& ptr, unsigned long num) {
+ ptr = GPUMemoryAllocation::gpu_alloc(num * sizeof(ScalarType));
+ };
+
+ auto GPUAllocAndCopy = [](const unsigned long*& ptr, const unsigned long*& src_ptr, unsigned long num) {
+ ptr = GPUMemoryAllocation::gpu_alloc_cpy(src_ptr, num * sizeof(const unsigned long));
+ };
+
+ GPUAllocAndInit(d_matrix, nnz * nVar * nEqn);
+ GPUAllocAndCopy(d_row_ptr, row_ptr, (nPointDomain + 1.0));
+ GPUAllocAndCopy(d_col_ind, col_ind, nnz);
+ }
+
if (needTranspPtr) col_ptr = geometry->GetTransposeSparsePatternMap(type).data();
if (type == ConnectivityType::FiniteVolume) {
@@ -151,13 +181,6 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi
nnz_ilu = csr_ilu.getNumNonZeros();
}
- /*--- Allocate data. ---*/
- auto allocAndInit = [](ScalarType*& ptr, unsigned long num) {
- ptr = MemoryAllocation::aligned_alloc(64, num * sizeof(ScalarType));
- };
-
- allocAndInit(matrix, nnz * nVar * nEqn);
-
/*--- Preconditioners. ---*/
if (ilu_needed) allocAndInit(ILU_matrix, nnz_ilu * nVar * nEqn);
diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu
new file mode 100644
index 000000000000..d63430c96f96
--- /dev/null
+++ b/Common/src/linear_algebra/CSysMatrixGPU.cu
@@ -0,0 +1,88 @@
+/*!
+ * \file CSysMatrixGPU.cu
+ * \brief Implementations of Kernels and Functions for Matrix Operations on the GPU
+ * \author A. Raj
+ * \version 8.2.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2024, 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 .
+ */
+
+#include "../../include/linear_algebra/CSysMatrix.hpp"
+#include "../../include/linear_algebra/GPUComms.cuh"
+
+template
+__global__ void GPUMatrixVectorProductAdd(matrixType* matrix, vectorType* vec, vectorType* prod, const unsigned long* d_row_ptr, const unsigned long* d_col_ind, unsigned long nPointDomain, unsigned long nVar, unsigned long nEqn)
+{
+ int row = (blockIdx.x * blockDim.x + threadIdx.x)/32;
+ int threadNo = threadIdx.x%32;
+ int activeThreads = nVar * (32/nVar);
+
+ int blockRow = (threadNo/nVar)%nVar;
+
+ if(row
+void CSysMatrix::HtDTransfer(bool trigger) const
+{
+ if(trigger) gpuErrChk(cudaMemcpy((void*)(d_matrix), (void*)&matrix[0], (sizeof(ScalarType)*nnz*nVar*nEqn), cudaMemcpyHostToDevice));
+}
+
+template
+void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, CSysVector& prod,
+ CGeometry* geometry, const CConfig* config) const
+ {
+
+ ScalarType* d_vec = vec.GetDevicePointer();
+ ScalarType* d_prod = prod.GetDevicePointer();
+
+ HtDTransfer();
+ vec.HtDTransfer();
+ prod.GPUSetVal(0.0);
+
+ dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1);
+ int gridx = KernelParameters::round_up_division(KernelParameters::MVP_WARP_SIZE, nPointDomain);
+ dim3 gridDim(gridx, 1, 1);
+
+ GPUMatrixVectorProductAdd<<>>(d_matrix, d_vec, d_prod, d_row_ptr, d_col_ind, nPointDomain, nVar, nEqn);
+ gpuErrChk( cudaPeekAtLastError() );
+
+ prod.DtHTransfer();
+
+}
+
+template class CSysMatrix; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.
diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp
index a2853a23cd6f..b99f03d4f499 100644
--- a/Common/src/linear_algebra/CSysVector.cpp
+++ b/Common/src/linear_algebra/CSysVector.cpp
@@ -52,6 +52,8 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB
if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType));
+ d_vec_val = GPUMemoryAllocation::gpu_alloc(nElm * sizeof(ScalarType));
+
if (val != nullptr) {
if (!valIsArray) {
for (auto i = 0ul; i < nElm; i++) vec_val[i] = *val;
@@ -66,6 +68,8 @@ CSysVector::~CSysVector() {
if (!std::is_trivial::value)
for (auto i = 0ul; i < nElm; i++) vec_val[i].~ScalarType();
MemoryAllocation::aligned_free(vec_val);
+
+ GPUMemoryAllocation::gpu_free(d_vec_val);
}
/*--- Explicit instantiations ---*/
diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu
new file mode 100644
index 000000000000..c07f84c3544f
--- /dev/null
+++ b/Common/src/linear_algebra/CSysVectorGPU.cu
@@ -0,0 +1,49 @@
+/*!
+ * \file CSysVectorGPU.cu
+ * \brief Implementations of Kernels and Functions for Vector Operations on the GPU
+ * \author A. Raj
+ * \version 8.2.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2024, 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 .
+ */
+
+#include "../../include/linear_algebra/CSysVector.hpp"
+#include "../../include/linear_algebra/GPUComms.cuh"
+
+template
+void CSysVector::HtDTransfer(bool trigger) const
+{
+ if(trigger) gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType)*nElm), cudaMemcpyHostToDevice));
+}
+
+template
+void CSysVector::DtHTransfer(bool trigger) const
+{
+ if(trigger) gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType)*nElm), cudaMemcpyDeviceToHost));
+}
+
+template
+void CSysVector::GPUSetVal(ScalarType val, bool trigger) const
+{
+ if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm)));
+}
+
+template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.
diff --git a/Common/src/linear_algebra/meson.build b/Common/src/linear_algebra/meson.build
index f8904b7f0f8a..7b880b29c1e3 100644
--- a/Common/src/linear_algebra/meson.build
+++ b/Common/src/linear_algebra/meson.build
@@ -4,3 +4,7 @@ common_src += files(['CSysSolve_b.cpp',
'CSysMatrix.cpp',
'CPastixWrapper.cpp',
'blas_structure.cpp'])
+
+ if get_option('enable-cuda')
+ common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu',])
+endif
diff --git a/config_template.cfg b/config_template.cfg
index 6b3d5ab396aa..32e6369e6e33 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1518,6 +1518,9 @@ DISCADJ_LIN_SOLVER= FGMRES
% Linear solver for the turbulent adjoint systems
ADJTURB_LIN_SOLVER= FGMRES
%
+% Use CUDA GPU Acceleration for FGMRES Linear Solver Only
+ENABLE_CUDA=NO
+%
% Preconditioner for the turbulent adjoint Krylov linear solvers
ADJTURB_LIN_PREC= ILU
%
diff --git a/doxyfile b/doxyfile
index b249d12d1b18..d7f30a3c83fb 100755
--- a/doxyfile
+++ b/doxyfile
@@ -291,7 +291,8 @@ OPTIMIZE_OUTPUT_VHDL = NO
# Note that for custom extensions you also need to set FILE_PATTERNS otherwise
# the files are not read by doxygen.
-EXTENSION_MAPPING =
+EXTENSION_MAPPING = cu=c++ \
+ cuh=c++
# If the MARKDOWN_SUPPORT tag is enabled then doxygen pre-processes all comments
# according to the Markdown format, which allows for more readable
@@ -818,7 +819,9 @@ INPUT_ENCODING = UTF-8
FILE_PATTERNS = *.cpp \
*.inl \
*.hpp \
- *.md
+ *.md \
+ *.cu \
+ *.cuh
# The RECURSIVE tag can be used to specify whether or not subdirectories should
# be searched for input files as well.
diff --git a/meson.build b/meson.build
index cdd905310798..443cb1f4e543 100644
--- a/meson.build
+++ b/meson.build
@@ -1,3 +1,4 @@
+
project('SU2', 'c', 'cpp',
version: '8.2.0 "Harrier"',
meson_version: '>=0.61.1',
@@ -15,6 +16,11 @@ endif
pymod = import('python')
python = pymod.find_installation()
+if get_option('enable-cuda')
+ add_languages('cuda')
+ add_global_arguments('-arch=sm_86', language : 'cuda')
+endif
+
su2_cpp_args = []
su2_deps = [declare_dependency(include_directories: 'externals/CLI11')]
@@ -196,6 +202,13 @@ if get_option('enable-pastix')
su2_deps += pastix_dep
endif
+# CUDA dependencies
+if get_option('enable-cuda')
+ su2_cpp_args += '-DHAVE_CUDA'
+ gpu_dep = dependency('cuda', version : '>=10', modules : ['cudart'])
+ su2_deps += gpu_dep
+endif
+
# blas-type dependencies
if get_option('enable-mkl')
diff --git a/meson_options.txt b/meson_options.txt
index 55782c62ff80..830b5f6574ee 100644
--- a/meson_options.txt
+++ b/meson_options.txt
@@ -9,6 +9,7 @@ option('enable-normal', type : 'boolean', value : true, description: 'enable no
option('enable-mkl', type : 'boolean', value : false, description: 'enable Intel-MKL support')
option('mkl_root', type : 'string', value : '/opt/intel/mkl', description: 'root of Intel-MKL installation (only for non-intel compilers)')
option('enable-openblas', type : 'boolean', value : false, description: 'enable BLAS and LAPACK support via OpenBLAS')
+option('enable-cuda', type : 'boolean', value : false, description: 'enable GPU acceleration using CUDA')
option('blas-name', type : 'string', value : 'openblas', description: 'name of the BLAS/LAPACK dependency')
option('enable-pastix', type : 'boolean', value : false, description: 'enable PaStiX support')
option('pastix_root', type : 'string', value : 'externals/pastix/', description: 'PaStiX base directory')