diff --git a/AUTHORS.md b/AUTHORS.md index eaea9531b71..fea157beba2 100644 --- a/AUTHORS.md +++ b/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 diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index ae30c779e7a..41d3c747cf8 100644 --- a/Common/include/code_config.hpp +++ b/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 diff --git a/Common/include/linear_algebra/CMatrixVectorProduct.hpp b/Common/include/linear_algebra/CMatrixVectorProduct.hpp index a0cecaa63d7..4069ff2fd00 100644 --- a/Common/include/linear_algebra/CMatrixVectorProduct.hpp +++ b/Common/include/linear_algebra/CMatrixVectorProduct.hpp @@ -72,7 +72,6 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct { const CSysMatrix& 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 { */ inline CSysMatrixVectorProduct(const CSysMatrix& 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) { + 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 { */ inline void operator()(const CSysVector& u, CSysVector& 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) { + 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 " diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index 532e0b53838..d28c729c4f0 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/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 +inline void ApplyPreconditionerOnHost(const CSysVector& u, CSysVector& v, Apply&& apply) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + 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 CPreconditioner::~CPreconditioner() {} +/*! + * \class CIdentityPreconditioner + * \brief No-op preconditioner used when Krylov solvers run without preconditioning. + */ +template +class CIdentityPreconditioner final : public CPreconditioner { + public: + inline void operator()(const CSysVector& u, CSysVector& 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 { * \param[out] v - CSysVector that is the result of the preconditioning. */ inline void operator()(const CSysVector& u, CSysVector& 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 { * \param[out] v - CSysVector that is the result of the preconditioning. */ inline void operator()(const CSysVector& u, CSysVector& 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 { CQuantizedLUSGSPreconditioner() = delete; inline void operator()(const CSysVector& u, CSysVector& 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 { * \param[out] v - CSysVector that is the result of the preconditioning. */ inline void operator()(const CSysVector& u, CSysVector& 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 { * \param[out] v - CSysVector that is the result of the preconditioning. */ inline void operator()(const CSysVector& u, CSysVector& 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* CPreconditioner::Create(ENUM_LINEAR_SOL CPreconditioner* prec = nullptr; switch (kind) { + case IDENTITY: + prec = new CIdentityPreconditioner(); + break; case JACOBI: prec = new CJacobiPreconditioner(jacobian, geometry, config); break; diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index d31c48fb8a0..3eddf892a7b 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/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& vec, CSysVector& 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& vec, CSysVector& prod, + CGeometry* geometry, const CConfig* config) const; + /*! * \brief Build the ILU preconditioner. */ diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index 2ea3cbf7df3..1f9bc851b92 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/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) { + 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) { + 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 - void HandleTemporariesIn(const CSysVector& LinSysRes, CSysVector& LinSysSol) { + void HandleTemporariesIn(const CSysVector& LinSysRes, CSysVector& LinSysSol, bool useCuda) { SU2_ZONE_SCOPED if constexpr (std::is_same_v) { /*--- 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 - void HandleTemporariesOut(CSysVector& LinSysSol) { + void HandleTemporariesOut(CSysVector& LinSysSol, bool useCuda) { SU2_ZONE_SCOPED + DownloadSolution(useCuda); if constexpr (std::is_same_v) { /*--- Same type specialization, temporary variables are not required. ---*/ BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { @@ -432,9 +473,9 @@ class CSysSolve { template > = 0> unsigned long Solve_b(MatrixType& Jacobian, const CSysVector& LinSysRes, CSysVector& 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; } diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 1498b549bbb..2ffbb9a7967 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -37,6 +37,35 @@ #include "vector_expressions.hpp" #include "../../include/CConfig.hpp" +#ifdef __CUDACC__ +#include "GPUComms.cuh" +#define SU2_CUDA_HOST_DEVICE __host__ __device__ +#else +#define SU2_CUDA_HOST_DEVICE +#endif + +template +class CSysVector; + +/*! + * \brief True for the plain floating-point scalar types the GPU vector kernels + * (CSysVectorGPU.cu) are instantiated for, in builds where those kernels + * exist at all. AD active types are never dispatched to the device: the + * tape/expression machinery those types pull in is not compatible with + * nvcc's device-code compilation, and device-resident autodiff is not + * supported by this GPU path. + * \note In an AD build this is false even for su2mixedfloat, because the .cu + * translation units are not linked into the AD libraries (see + * SU2_ENABLE_CUDA_KERNELS in code_config.hpp). + */ +#ifdef SU2_ENABLE_CUDA_KERNELS +template +inline constexpr bool su2_gpu_capable_v = std::is_floating_point_v; +#else +template +inline constexpr bool su2_gpu_capable_v = false; +#endif + /*! * \brief OpenMP worksharing construct used in CSysVector for loops. * \note The loop will only run in parallel if methods are called from a @@ -59,6 +88,106 @@ #define END_CSYSVEC_PARFOR #endif +/*! + * \brief Brackets device work so that it is issued by a single thread with the whole team + * synchronized before and after. + * \note The GPU is one shared resource and the device path does not use OpenMP worksharing. + * Issuing from one thread keeps kernel launches ordered on the default stream and, above + * all, stops part of a team from entering a worksharing construct that the rest skipped. + * Correctness relies on all threads reaching the same vector operations in the same order, + * which is the assumption the "nowait" clause on CSYSVEC_PARFOR already makes. These + * regions must not be nested; they are used by the operations of this class and by the + * matrix-vector product and preconditioner wrappers, and by nothing above those. + */ +#define SU2_DEVICE_REGION(...) SU2_OMP_SAFE_GLOBAL_ACCESS(__VA_ARGS__) +#define BEGIN_SU2_DEVICE_REGION BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS +#define END_SU2_DEVICE_REGION END_SU2_OMP_SAFE_GLOBAL_ACCESS + +namespace VecExpr { + +enum class DeviceAssignOp { Assign, Add, Subtract, Multiply, Divide }; + +/*! + * \brief Whether vector expressions are currently evaluated on the device. + * \note Defined in CSysVectorGPU.cu. This is a plain global, not per thread: every thread + * of a team has to agree on it or they would split over the worksharing constructs below. + * It is switched by CSysSolve at the same boundary that uploads and downloads the vectors + * (HandleTemporariesIn/Out), and nowhere else, so it does not change while a solve runs. + */ +#ifdef SU2_ENABLE_CUDA_KERNELS +bool UseDeviceExpressions(); +void SetUseDeviceExpressions(bool use); +#else +inline bool UseDeviceExpressions() { return false; } +inline void SetUseDeviceExpressions(bool) {} +#endif + +/*! + * \brief How a CSysVector is captured inside an expression: a bare pointer to whichever + * storage the expression is going to be evaluated from. + * \note Capturing by value (rather than a reference to the vector) is what makes an + * arbitrary expression tree trivially copyable, and therefore passable by value to the + * assignment kernel. The choice of storage is fixed when the expression is built, which is + * sound because there is no fallback: while UseDeviceExpressions() holds, every expression + * is evaluated by a kernel. + */ +template +class CVectorView : public CVecExpr, Scalar> { + private: + const Scalar* data = nullptr; + + public: + static constexpr bool StoreAsRef = false; + + CVectorView(const CSysVector& vector); + + SU2_CUDA_HOST_DEVICE FORCEINLINE const Scalar& operator[](size_t i) const { return data[i]; } +}; + +template +struct store_type> { + using type = CVectorView; +}; + +template +struct store_type> { + using type = CVectorView; +}; + +template +void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr& expr); + +#ifdef __CUDACC__ +template +__global__ void DeviceAssignKernel(Scalar* data, unsigned long size, T expr) { + const unsigned long i = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= size) return; + + if constexpr (Op == DeviceAssignOp::Assign) { + data[i] = expr[i]; + } else if constexpr (Op == DeviceAssignOp::Add) { + data[i] += expr[i]; + } else if constexpr (Op == DeviceAssignOp::Subtract) { + data[i] -= expr[i]; + } else if constexpr (Op == DeviceAssignOp::Multiply) { + data[i] *= expr[i]; + } else { + data[i] /= expr[i]; + } +} + +template +inline void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr& expr) { + if (size == 0) return; + constexpr unsigned block_size = 256; + const auto grid_size = static_cast((size + block_size - 1) / block_size); + DeviceAssignKernel<<>>(data, size, expr.derived()); + gpuErrChk(cudaPeekAtLastError()); +} +#endif + +} // namespace VecExpr + /*! * \class CSysVector * \ingroup SpLinSys @@ -110,6 +239,37 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } } + /*! + * \brief Evaluates an expression into the device storage of this vector. + * \note The kernel has to be instantiated for the expression type in CSysVectorGPU.cu, + * a shape that is not in that list is an undefined symbol at link time. + */ + template + CSysVector& AssignDevice(const VecExpr::CVecExpr& expr) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + BEGIN_SU2_DEVICE_REGION { + VecExpr::store_t stored_expr(expr.derived()); + VecExpr::AssignDeviceExpression(d_vec_val, nElm, stored_expr); + } + END_SU2_DEVICE_REGION + } +#endif + return *this; + } + + template + CSysVector& AssignDevice(ScalarType val) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + BEGIN_SU2_DEVICE_REGION + VecExpr::AssignDeviceExpression(d_vec_val, nElm, VecExpr::Bcast(val)); + END_SU2_DEVICE_REGION + } +#endif + return *this; + } + public: static constexpr bool StoreAsRef = true; /*! \brief Required by CVecExpr. */ @@ -235,16 +395,31 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> 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. + * \brief Dot product between this vector and another vector on the device. + * \note Explicit GPU helper for solver-side reductions. + * \param[in] other - Input vector. + * \return Dot product result. + */ + ScalarType GPUDot(const CSysVector& other) const; + + /*! + * \brief L2 norm of this vector on the device. + * \note Explicit GPU helper for solver-side reductions. + * \return L2 norm result. */ - void GPUSetVal(ScalarType val, bool trigger = true) const; + ScalarType GPUNorm() const; /*! * \brief return device pointer that points to the CSysVector values in GPU memory */ inline ScalarType* GetDevicePointer() const { return d_vec_val; } + /*! + * \brief return host pointer that points to the CSysVector values, counterpart of + * GetDevicePointer + */ + inline const ScalarType* GetHostPointer() const { return vec_val; } + /*! * \brief return the number of local elements in the CSysVector */ @@ -301,6 +476,11 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] other - Another vector. */ CSysVector& operator=(const CSysVector& other) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + if (VecExpr::UseDeviceExpressions()) return AssignDevice(other); + } +#endif CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR @@ -311,25 +491,31 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief Compound assignement operations with scalars and expressions. * \param[in] val/expr - Scalar value or expression. */ -#define MAKE_COMPOUND(OP) \ - CSysVector& operator OP(ScalarType val) { \ - CSYSVEC_PARFOR \ - for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ - END_CSYSVEC_PARFOR \ - return *this; \ - } \ - template \ - CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ - CSYSVEC_PARFOR \ - for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ - END_CSYSVEC_PARFOR \ - return *this; \ +#define MAKE_COMPOUND(OP, ASSIGN_OP) \ + CSysVector& operator OP(ScalarType val) { \ + if constexpr (su2_gpu_capable_v) { \ + if (VecExpr::UseDeviceExpressions()) return AssignDevice(val); \ + } \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ + END_CSYSVEC_PARFOR \ + return *this; \ + } \ + template \ + CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ + if constexpr (su2_gpu_capable_v) { \ + if (VecExpr::UseDeviceExpressions()) return AssignDevice(expr); \ + } \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ + END_CSYSVEC_PARFOR \ + return *this; \ } - MAKE_COMPOUND(=) - MAKE_COMPOUND(+=) - MAKE_COMPOUND(-=) - MAKE_COMPOUND(*=) - MAKE_COMPOUND(/=) + MAKE_COMPOUND(=, VecExpr::DeviceAssignOp::Assign) + MAKE_COMPOUND(+=, VecExpr::DeviceAssignOp::Add) + MAKE_COMPOUND(-=, VecExpr::DeviceAssignOp::Subtract) + MAKE_COMPOUND(*=, VecExpr::DeviceAssignOp::Multiply) + MAKE_COMPOUND(/=, VecExpr::DeviceAssignOp::Divide) #undef MAKE_COMPOUND /*! @@ -344,6 +530,21 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ template ScalarType dot(const VecExpr::CVecExpr& expr) const { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + using DeviceExpr = std::remove_cv_t>; + static_assert(std::is_same_v, + "On the device the dot product is a cuBLAS call, so it only takes vectors. " + "Assign the expression to a vector first."); + if (VecExpr::UseDeviceExpressions()) { + /*--- GPUDot reduces over MPI, which has to happen once for the team, so the result + * is published through the same scratch slot the host reduction below uses. ---*/ + SU2_DEVICE_REGION(dot_scratch[0] = GPUDot(expr.derived());) + return dot_scratch[0]; + } + } +#endif + /*--- All threads get the same "view" of the vectors. ---*/ SU2_OMP_BARRIER @@ -502,5 +703,14 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } }; +namespace VecExpr { + +template +CVectorView::CVectorView(const CSysVector& vector) + : data(UseDeviceExpressions() ? vector.GetDevicePointer() : vector.GetHostPointer()) {} + +} // namespace VecExpr + #undef CSYSVEC_PARFOR #undef END_CSYSVEC_PARFOR +#undef SU2_CUDA_HOST_DEVICE diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh index 13854381877..eb727477f21 100644 --- a/Common/include/linear_algebra/GPUComms.cuh +++ b/Common/include/linear_algebra/GPUComms.cuh @@ -25,6 +25,11 @@ * License along with SU2. If not, see . */ +#pragma once + +#ifndef SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH +#define SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH + #include #include @@ -51,3 +56,5 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t } #define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); } + +#endif // SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index a0d0ce28901..f9b941a0bdb 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -39,6 +39,12 @@ namespace VecExpr { /// \addtogroup VecExpr /// @{ +#ifdef __CUDACC__ +#define SU2_CUDA_HOST_DEVICE __host__ __device__ +#else +#define SU2_CUDA_HOST_DEVICE +#endif + /*! * \brief Base vector expression class. * \ingroup BLAS @@ -59,7 +65,7 @@ class CVecExpr { /*! * \brief Cast the expression to Derived, usually to allow evaluation via operator[]. */ - FORCEINLINE const Derived& derived() const { return static_cast(*this); } + SU2_CUDA_HOST_DEVICE FORCEINLINE const Derived& derived() const { return static_cast(*this); } // Allowed from C++14, allows nested expression propagation without // manually calling derived() on the expression being evaluated. @@ -76,8 +82,8 @@ class Bcast : public CVecExpr, Scalar> { public: static constexpr bool StoreAsRef = false; - FORCEINLINE Bcast(const Scalar& x_) : x(x_) {} - FORCEINLINE const Scalar& operator[](size_t) const { return x; } + SU2_CUDA_HOST_DEVICE FORCEINLINE Bcast(const Scalar& x_) : x(x_) {} + SU2_CUDA_HOST_DEVICE FORCEINLINE const Scalar& operator[](size_t) const { return x; } }; /*! @@ -101,7 +107,11 @@ struct add_lref_if { using type = remove_reference_t&; }; template -using store_t = typename add_lref_if::type; +struct store_type { + using type = typename add_lref_if::type; +}; +template +using store_t = typename store_type::type; /*--- Namespace from which the math function implementations come. ---*/ @@ -120,19 +130,19 @@ namespace math = ::std; /*--- Macro to create expression classes (EXPR) and overloads (FUN) for unary * functions, based on their coefficient-wise implementation (IMPL). ---*/ -#define MAKE_UNARY_FUN(FUN, EXPR, IMPL) \ - /*!--- Expression class. ---*/ \ - template \ - class EXPR : public CVecExpr, Scalar> { \ - store_t u; \ - \ - public: \ - static constexpr bool StoreAsRef = false; \ - FORCEINLINE EXPR(const U& u_) : u(u_) {} \ - FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \ - }; \ - /*!--- Function overload, returns an expression object. ---*/ \ - template \ +#define MAKE_UNARY_FUN(FUN, EXPR, IMPL) \ + /*!--- Expression class. ---*/ \ + template \ + class EXPR : public CVecExpr, Scalar> { \ + store_t u; \ + \ + public: \ + static constexpr bool StoreAsRef = false; \ + FORCEINLINE EXPR(const U& u_) : u(u_) {} \ + SU2_CUDA_HOST_DEVICE FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \ + }; \ + /*!--- Function overload, returns an expression object. ---*/ \ + template \ FORCEINLINE auto FUN(const CVecExpr& u) RETURNS(EXPR(u.derived())) #define sign_impl(x) Scalar(1 - 2 * (x < 0)) @@ -158,7 +168,7 @@ MAKE_UNARY_FUN(sign, sign_, sign_impl) public: \ static constexpr bool StoreAsRef = false; \ FORCEINLINE EXPR(const U& u_, const V& v_) : u(u_), v(v_) {} \ - FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i])) \ + SU2_CUDA_HOST_DEVICE FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i])) \ }; \ /*!--- Vector with vector function overload. ---*/ \ template \ @@ -241,4 +251,5 @@ MAKE_BINARY_FUN(operator>, gt_, gt_impl) #undef MAKE_BINARY_FUN /// @} +#undef SU2_CUDA_HOST_DEVICE } // namespace VecExpr diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 2bc09b47ada..2dcd241f9de 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2523,6 +2523,7 @@ static const MapType Sens_Smoothing_Map = { * \brief Types of preconditioners for the linear solver */ enum ENUM_LINEAR_SOLVER_PREC { + IDENTITY, /*!< \brief No preconditioner. */ JACOBI, /*!< \brief Jacobi preconditioner. */ LU_SGS, /*!< \brief LU SGS preconditioner. */ LINELET, /*!< \brief Line implicit preconditioner. */ @@ -2533,6 +2534,7 @@ enum ENUM_LINEAR_SOLVER_PREC { PASTIX_LDLT_P, /*!< \brief PaStiX LDLT as preconditioner. */ }; static const MapType Linear_Solver_Prec_Map = { + MakePair("NONE", IDENTITY) MakePair("JACOBI", JACOBI) MakePair("LU_SGS", LU_SGS) MakePair("LINELET", LINELET) diff --git a/Common/include/parallelization/omp_structure.hpp b/Common/include/parallelization/omp_structure.hpp index af581bd1d76..10bf9f7bca0 100644 --- a/Common/include/parallelization/omp_structure.hpp +++ b/Common/include/parallelization/omp_structure.hpp @@ -283,7 +283,10 @@ inline void atomicAdd(T rhs, T& lhs) { #define ATOMIC_COMPARE_FALLBACK /*--- Atomic max, shared = max(shared, local). ---*/ -#ifdef _OPENMP +/*--- nvcc's host pass drops the clause from "#pragma omp atomic compare", which the host + * compiler then rejects. The .cu sources do not use these functions, so they simply get + * the critical section fallback below. ---*/ +#if defined(_OPENMP) && !defined(__CUDACC__) #if _OPENMP >= ATOMIC_COMPARE_SINCE /*--- Atomic min/max are supported for arithmetic types. ---*/ template ::value> = 0> @@ -305,7 +308,7 @@ inline void atomicMax(const T& local, T& shared) { } /*--- Atomic min, shared = min(shared, local). ---*/ -#ifdef _OPENMP +#if defined(_OPENMP) && !defined(__CUDACC__) #if _OPENMP >= ATOMIC_COMPARE_SINCE template ::value> = 0> inline void atomicMin(const T& local, T& shared) { diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index fa47816f60a..a0013e24300 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -93,6 +93,7 @@ CSysMatrix::CSysMatrix() : rank(SU2_MPI::GetRank()), size(SU2_MPI::G q_blocks_d = nullptr; invM = nullptr; + d_invM = nullptr; #ifdef USE_MKL MatrixMatrixProductJitter = nullptr; @@ -129,6 +130,7 @@ CSysMatrix::~CSysMatrix() { GPUMemoryAllocation::gpu_free(gpu.col_ind_l); GPUMemoryAllocation::gpu_free(gpu.row_ptr_u); GPUMemoryAllocation::gpu_free(gpu.col_ind_u); + GPUMemoryAllocation::gpu_free(d_invM); } #ifdef USE_MKL @@ -282,6 +284,10 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn); + if (useCuda && diag_needed) { + d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + } + /*--- Thread parallel initialization. ---*/ int num_threads = omp_get_max_threads(); @@ -804,11 +810,29 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector& v template void CSysMatrix::BuildJacobiPreconditioner() { SU2_ZONE_SCOPED + /*--- Build Jacobi preconditioner (M = D), compute and store the inverses of the diagonal blocks. ---*/ SU2_OMP_FOR_DYN(omp_heavy_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) InverseDiagonalBlock(iPoint, &(invM[iPoint * nVar * nVar])); END_SU2_OMP_FOR + + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + BEGIN_SU2_DEVICE_REGION + gpuErrChk(cudaMemcpy(d_invM, invM, nPointDomain * nVar * nVar * sizeof(ScalarType), cudaMemcpyHostToDevice)); + END_SU2_DEVICE_REGION + } else { + SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION); + } +#else + SU2_MPI::Error( + "\nError in building Jacobi preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options " + "enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } } template @@ -816,6 +840,23 @@ void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector& prod, CGeometry* geometry, const CConfig* config) const { SU2_ZONE_SCOPED + + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + SU2_DEVICE_REGION(ComputeJacobiPreconditionerGPU(vec, prod, geometry, config);) + return; + } else { + SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION); + } +#else + SU2_MPI::Error( + "\nError in applying Jacobi preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options " + "enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + /*--- Apply Jacobi preconditioner, y = D^{-1} * x, the inverse of the diagonal is already known. ---*/ SU2_OMP_BARRIER SU2_OMP_FOR_DYN(omp_heavy_size) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 1ebba30097c..381379a0d81 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -34,12 +34,12 @@ */ template __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, - const unsigned long* __restrict__ row_ptr_l, - const unsigned long* __restrict__ col_ind_l, + const su2uint* __restrict__ row_ptr_l, + const su2uint* __restrict__ col_ind_l, const ScalarType* __restrict__ mat_l, const ScalarType* __restrict__ mat_d, - const unsigned long* __restrict__ row_ptr_u, - const unsigned long* __restrict__ col_ind_u, + const su2uint* __restrict__ row_ptr_u, + const su2uint* __restrict__ col_ind_u, const ScalarType* __restrict__ mat_u, const ScalarType* __restrict__ x, ScalarType* __restrict__ y) { const unsigned long iRow = blockIdx.x; @@ -84,7 +84,6 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector ScalarType* d_vec = vec.GetDevicePointer(); ScalarType* d_prod = prod.GetDevicePointer(); - vec.HtDTransfer(); dim3 blockDim(static_cast(nVar), 1, 1); dim3 gridDim(static_cast(nPointDomain), 1, 1); @@ -92,8 +91,15 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, gpu.l, gpu.d, gpu.row_ptr_u, gpu.col_ind_u, gpu.u, d_vec, d_prod); gpuErrChk(cudaGetLastError()); - - prod.DtHTransfer(); } +template void CSysMatrix::HtDTransfer(bool trigger) const; +template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const; -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. +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template void CSysMatrix::HtDTransfer(bool trigger) const; +template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const; +#endif diff --git a/Common/src/linear_algebra/CSysPreconditionerGPU.cu b/Common/src/linear_algebra/CSysPreconditionerGPU.cu new file mode 100644 index 00000000000..794726fbce3 --- /dev/null +++ b/Common/src/linear_algebra/CSysPreconditionerGPU.cu @@ -0,0 +1,84 @@ +/*! + * \file CSysPreconditionerGPU.cu + * \brief CUDA/GPU skeleton implementations for matrix-based preconditioners. + * \author Jesse Li + * \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 . + */ + +#include "../../include/linear_algebra/CSysMatrix.inl" +#include "../../include/linear_algebra/GPUComms.cuh" + +namespace { + +template +__global__ void ApplyJacobiPreconditionerKernel(const ScalarType* invM, const ScalarType* vec, ScalarType* prod, + unsigned long nPointDomain, unsigned long nVar) { + const auto iPoint = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (iPoint >= nPointDomain) return; + + const auto block = &invM[iPoint * nVar * nVar]; + const auto rhs = &vec[iPoint * nVar]; + auto out = &prod[iPoint * nVar]; + + for (auto iVar = 0ul; iVar < nVar; ++iVar) { + auto sum = ScalarType(0); + for (auto jVar = 0ul; jVar < nVar; ++jVar) { + sum += block[iVar * nVar + jVar] * rhs[jVar]; + } + out[iVar] = sum; + } +} + +} // namespace + +template +void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const { + (void)geometry; + (void)config; + + SU2_ZONE_SCOPED + + if (d_invM == nullptr) { + SU2_MPI::Error("CUDA Jacobi preconditioner used before BuildJacobiPreconditionerGPU.", CURRENT_FUNCTION); + } + + constexpr unsigned threadsPerBlock = 128; + const auto blocks = static_cast((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + ApplyJacobiPreconditionerKernel<<>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), + nPointDomain, nVar); + gpuErrChk(cudaPeekAtLastError()); +} + +template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, + CSysVector& prod, + CGeometry* geometry, + const CConfig* config) const; + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, + CSysVector& prod, + CGeometry* geometry, + const CConfig* config) const; +#endif diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index a50ed059625..dcd88c452c5 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -427,7 +427,7 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector 1; + const bool nestedParallel = !omp_in_parallel() && omp_get_max_threads() > 1 && !VecExpr::UseDeviceExpressions(); /*--- Check the subspace size ---*/ @@ -664,7 +664,7 @@ unsigned long CSysSolve::FGCRODR_LinSolverImpl(const CSysVector 1; + const bool nestedParallel = !omp_in_parallel() && omp_get_max_threads() > 1 && !VecExpr::UseDeviceExpressions(); /*--- Check the subspace size. ---*/ @@ -1464,7 +1464,7 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con auto externalFunction = [&]() { /*--- Create matrix-vector product, preconditioner, and solve the linear system ---*/ - HandleTemporariesIn(LinSysRes, LinSysSol); + HandleTemporariesIn(LinSysRes, LinSysSol, config->GetCUDA()); auto mat_vec = CSysMatrixVectorProduct(Jacobian, geometry, config); @@ -1539,7 +1539,7 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con } END_SU2_OMP_MASTER - HandleTemporariesOut(LinSysSol); + HandleTemporariesOut(LinSysSol, config->GetCUDA()); delete normal_prec; delete nested_prec; diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index f7df34633a0..54e47157d55 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -52,6 +52,9 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); + /*--- Device storage mirrors the host allocation; free first so that re-initializing a + * vector does not leak it. ---*/ + GPUMemoryAllocation::gpu_free(d_vec_val); d_vec_val = GPUMemoryAllocation::gpu_alloc(nElm * sizeof(ScalarType)); #ifdef HAVE_OMP @@ -78,6 +81,23 @@ const su2matrix& CSysVector::multiDot(const std::vector< if (n == 0 || m == 0) return shared; +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + if (VecExpr::UseDeviceExpressions()) { + BEGIN_SU2_DEVICE_REGION { + shared.resize(n, m); + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < m; ++j) { + shared(i, j) = V[i0 + i].GPUDot(W[j]); + } + } + } + END_SU2_DEVICE_REGION + return shared; + } + } +#endif + SU2_OMP_BARRIER const size_t size = V[0].nElmDomain; diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 94ec17bb88f..2be1215a7bd 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -27,23 +27,192 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" +#include +#include +#include +#include -template -void CSysVector::HtDTransfer(bool trigger) const -{ - if(trigger) gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType)*nElm), cudaMemcpyHostToDevice)); +namespace { + +/*--- cuBLAS handle for the reductions. Created on first use and kept for the lifetime of + * the program, matching the fact that CUDA is either on or off for the whole run. ---*/ +cublasHandle_t solver_blas_handle = nullptr; + +cublasHandle_t GetBlasHandle() { + if (solver_blas_handle == nullptr) { + if (cublasCreate(&solver_blas_handle) != CUBLAS_STATUS_SUCCESS) { + SU2_MPI::Error("cuBLAS handle creation failed for the GPU linear algebra.", CURRENT_FUNCTION); + } + } + return solver_blas_handle; } -template -void CSysVector::DtHTransfer(bool trigger) const -{ - if(trigger) gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType)*nElm), cudaMemcpyDeviceToHost)); +} // namespace + +namespace VecExpr { + +namespace { +/*--- Deliberately not thread local: every thread of an OpenMP team has to agree on this, + * otherwise the team splits over the worksharing constructs in CSysVector. It is only + * written by CSysSolve, outside any parallel region over the linear system. ---*/ +bool use_device_expressions = false; +} // namespace + +bool UseDeviceExpressions() { return use_device_expressions; } + +void SetUseDeviceExpressions(bool use) { use_device_expressions = use; } + +} // namespace VecExpr + +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 +ScalarType CSysVector::GPUDot(const CSysVector& other) const { + /*--- Both operands are already on the device, the caller owns the transfers. This + * reduces over MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). ---*/ + cublasHandle_t handle = GetBlasHandle(); + cublasStatus_t status = CUBLAS_STATUS_SUCCESS; + + ScalarType local_dot = ScalarType(0); + + if constexpr (std::is_same_v) { + status = cublasSdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, + &local_dot); + } else if constexpr (std::is_same_v) { + status = cublasDdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, + &local_dot); + } else { + SU2_MPI::Error("Unsupported ScalarType in CSysVector::GPUDot.", CURRENT_FUNCTION); + return ScalarType(0); + } + + if (status != CUBLAS_STATUS_SUCCESS) { + SU2_MPI::Error("cuBLAS dot failed in CSysVector::GPUDot.", CURRENT_FUNCTION); + return ScalarType(0); + } + + ScalarType global_dot = ScalarType(0); + const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE; + SelectMPIWrapper::W::Allreduce(&local_dot, &global_dot, 1, mpi_type, MPI_SUM, SU2_MPI::GetComm()); + + return global_dot; } -template -void CSysVector::GPUSetVal(ScalarType val, bool trigger) const -{ - if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm))); +template +ScalarType CSysVector::GPUNorm() const { + return sqrt(GPUDot(*this)); } -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. +/*--- Every expression the solvers assign to a CSysVector needs its assignment kernel + * instantiated here; the host compiler cannot emit one. A shape that is missing shows up + * as an undefined reference to VecExpr::AssignDeviceExpression at link time, and is fixed + * by adding a line to DEVICE_EXPRESSION_SHAPES below. The aliases use CSysVector (not + * CVectorView) because that is how the operator overloads name their operands; store_t + * turns it into a view when the node is built. ---*/ +namespace { + +template +using Vec = CSysVector; +template +using Sca = VecExpr::Bcast; + +/*--- Leaves and the shapes of the FGMRES/GMRES basis updates. ---*/ +template +using DeviceBcast = Sca; +template +using DeviceView = VecExpr::CVectorView; +template +using DeviceNeg = VecExpr::minus_, S>; + +/*--- vector * scalar and scalar * vector are distinct types, both are used. ---*/ +template +using DeviceScale = VecExpr::mul_, Sca, S>; +template +using DeviceLScale = VecExpr::mul_, Vec, S>; +template +using DeviceDivScale = VecExpr::div_, Sca, S>; + +/*--- Linear combinations, CSysSolve unrolls them up to four terms. ---*/ +template +using DeviceScale2 = VecExpr::add_, DeviceScale, S>; +template +using DeviceScale3 = VecExpr::add_, DeviceScale, S>; +template +using DeviceScale4 = VecExpr::add_, DeviceScale, S>; + +/*--- r = b - A_x, in CG, BCGSTAB, Smoother and FGCRODR. ---*/ +template +using DeviceSub = VecExpr::sub_, Vec, S>; + +/*--- p = beta * p + z, in CG. ---*/ +template +using DeviceLScalePlus = VecExpr::add_, Vec, S>; + +/*--- p = beta * (p - omega * v) + r, in BCGSTAB. ---*/ +template +using DeviceSubLScale = VecExpr::sub_, DeviceLScale, S>; +template +using DeviceLScaleSub = VecExpr::mul_, DeviceSubLScale, S>; +template +using DeviceBcgsDir = VecExpr::add_, Vec, S>; + +} // namespace + +#define DEVICE_EXPRESSION_SHAPES(SCALAR) \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceBcast); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceView); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceNeg); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceLScale); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceDivScale); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale2); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale3); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale4); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceSub); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceLScalePlus); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceSubLScale); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceLScaleSub); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceBcgsDir) + +#define INSTANTIATE_DEVICE_ASSIGN(SCALAR, OP, EXPR) \ + template void VecExpr::AssignDeviceExpression>( \ + SCALAR*, unsigned long, const VecExpr::CVecExpr, SCALAR>&) + +#define INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, EXPR) \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Assign, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Add, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Subtract, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Multiply, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Divide, EXPR) + +DEVICE_EXPRESSION_SHAPES(su2mixedfloat); + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +DEVICE_EXPRESSION_SHAPES(passivedouble); +#endif + +#undef DEVICE_EXPRESSION_SHAPES +#undef INSTANTIATE_DEVICE_ASSIGN_EXPR +#undef INSTANTIATE_DEVICE_ASSIGN + +template void CSysVector::HtDTransfer(bool trigger) const; +template void CSysVector::DtHTransfer(bool trigger) const; +template su2mixedfloat CSysVector::GPUDot(const CSysVector& other) const; +template su2mixedfloat CSysVector::GPUNorm() const; + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template void CSysVector::HtDTransfer(bool trigger) const; +template void CSysVector::DtHTransfer(bool trigger) const; +template passivedouble CSysVector::GPUDot(const CSysVector& other) const; +template passivedouble CSysVector::GPUNorm() const; +#endif diff --git a/Common/src/linear_algebra/meson.build b/Common/src/linear_algebra/meson.build index 7b880b29c1e..48ef65cb8db 100644 --- a/Common/src/linear_algebra/meson.build +++ b/Common/src/linear_algebra/meson.build @@ -5,6 +5,8 @@ common_src += files(['CSysSolve_b.cpp', 'CPastixWrapper.cpp', 'blas_structure.cpp']) - if get_option('enable-cuda') - common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu',]) +if get_option('enable-cuda') + # Kept apart from common_src: these are compiled without the CoDiPack defines and so + # must only go into the primal library, see common_cuda_src in Common/src/meson.build. + common_cuda_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu', 'CSysPreconditionerGPU.cu']) endif diff --git a/Common/src/meson.build b/Common/src/meson.build index f385c4a32ed..c34b014ab42 100644 --- a/Common/src/meson.build +++ b/Common/src/meson.build @@ -6,6 +6,25 @@ common_src =files(['graph_coloring_structure.cpp', '../include/parallelization/mpi_structure.cpp', '../include/parallelization/omp_structure.cpp']) +# nvcc cannot parse CoDiPack's tape machinery in the device pass, so the CUDA sources are +# compiled without the CODI_REVERSE_TYPE/CODI_FORWARD_TYPE defines. That gives them a +# different su2double (and hence a different layout for CConfig, CGeometry, ...) than the +# AD libraries, so they are collected separately and only linked into the primal library. +# The AD builds compile the device dispatch out entirely, see SU2_ENABLE_CUDA_KERNELS. +common_cuda_src = [] +common_cuda_cpp_args = [] +foreach arg : su2_cpp_args + if arg.startswith('-D') or arg.startswith('-U') + common_cuda_cpp_args += arg + endif +endforeach + +# Note that the OpenMP flags do reach nvcc (through omp_dep), so the .cu objects agree +# with the rest of the library about HAVE_OMP. They must still not instantiate anything +# containing OpenMP directives, because nvcc's host pass rewrites some of them; the +# device path is synchronized in the headers (see SU2_DEVICE_REGION in CSysVector.hpp), +# which only the .cpp sources instantiate. + subdir('linear_algebra') subdir('toolboxes') subdir('geometry') @@ -21,10 +40,11 @@ subdir('adt') if get_option('enable-normal') common = static_library('SU2Common', - common_src, + common_src, common_cuda_src, install : false, dependencies : su2_deps, - cpp_args: [default_warning_flags, su2_cpp_args]) + cpp_args: [default_warning_flags, su2_cpp_args], + cuda_args: common_cuda_cpp_args) common_dep = declare_dependency(link_with: common, include_directories : common_include) diff --git a/meson.build b/meson.build index e76354a6e86..a852cf7f113 100644 --- a/meson.build +++ b/meson.build @@ -19,11 +19,24 @@ python = pymod.find_installation() if get_option('enable-cuda') add_languages('cuda') - add_global_arguments('-arch=sm_86', language : 'cuda') + add_global_arguments('-arch=sm_89', language : 'cuda') + # nvcc's frontend does not recognize the AMX-tile builtins pulled in by + # newer glibc/gcc ; SU2 does not use AMX, so skip the header. + add_global_arguments('-D_AMXTILEINTRIN_H_INCLUDED', language : 'cuda') + # Any target linking CUDA object code is linked via nvcc instead of the C++ + # linker, which otherwise defaults to a plain host gcc and silently drops + # the MPI link flags that -Dcustom-mpi=true relies on mpicxx to provide. + add_global_arguments('-ccbin=' + meson.get_compiler('cpp').cmd_array()[0], language : 'cuda') + add_global_link_arguments('-ccbin=' + meson.get_compiler('cpp').cmd_array()[0], language : 'cuda') + cuda_deps = [ + meson.get_compiler('cuda').find_library('cublas', required : true), + ] +else + cuda_deps = [] endif su2_cpp_args = [] -su2_deps = [declare_dependency(include_directories: 'externals/CLI11')] +su2_deps = [declare_dependency(include_directories: 'externals/CLI11')] + cuda_deps default_warning_flags = [] if build_machine.system() != 'windows' @@ -231,8 +244,7 @@ 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 + su2_deps += meson.get_compiler('cuda').find_library('cudart', required : true) endif # blas-type dependencies diff --git a/su2omp.syntax.json b/su2omp.syntax.json index eb3f79653e2..382986c03d2 100644 --- a/su2omp.syntax.json +++ b/su2omp.syntax.json @@ -39,6 +39,7 @@ "CSYSVEC_PARFOR": "END_CSYSVEC_PARFOR", "CNEWTON_PARFOR": "END_CNEWTON_PARFOR", "BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS": "END_SU2_OMP_SAFE_GLOBAL_ACCESS", + "BEGIN_SU2_DEVICE_REGION": "END_SU2_DEVICE_REGION", "CPHYSGEO_PARFOR": "END_CPHYSGEO_PARFOR" } }