diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index 8bd7d515710..d8cbfd9a806 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -64,6 +64,33 @@ inline void ApplyPreconditionerOnHost(const CSysVector& u, CSysVecto apply(); } +/*! + * \brief Mirror of ApplyPreconditionerOnHost: applies a device preconditioner to host vectors. + * \note For callers that drive the Krylov solvers themselves and so never went through + * CSysSolve::Solve, which is what normally leaves the vectors on the device (Newton-Krylov). + * Device expressions are on for the duration so that a nested solve also uses the device copies. + * Only \p u is uploaded, \p v is always overwritten by the apply. + */ +template +inline void ApplyPreconditionerOnDevice(const CSysVector& u, CSysVector& v, bool useCuda, + Apply&& apply) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + if (useCuda && !VecExpr::UseDeviceExpressions()) { + SU2_DEVICE_REGION(u.HtDTransfer(); VecExpr::SetUseDeviceExpressions(true);) + + apply(); + + SU2_DEVICE_REGION(VecExpr::SetUseDeviceExpressions(false); v.DtHTransfer();) + return; + } + } +#else + (void)useCuda; +#endif + apply(); +} + /*! * \class CPreconditioner * \brief Abstract base class for defining a preconditioning operation. diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 384b0bf6212..f35efcef415 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -385,20 +385,24 @@ class CSysMatrix { mutable struct CUgraphExec_st* precond_bwd_graph_exec = nullptr; // LU-SGS backward only mutable const ScalarType* precond_fwd_graph_vec = nullptr; /*!< \brief Pointers the apply graph * was captured with, to detect when - * it must be recaptured. */ + * it must be recaptured (the + * executable graph itself is then + * updated in place, not rebuilt, + * see InstantiateOrUpdateGraph). */ mutable ScalarType* precond_fwd_graph_prod = nullptr; mutable ScalarType* precond_bwd_graph_prod = nullptr; - /*--- Non-default stream, needed for two mutually exclusive uses that never overlap on a given - * matrix (quantized_mode and ILU are alternative preconditioner choices, decided once in - * Initialize()): (1) the ILU build/apply CUDA graphs below, since the legacy default stream - * cannot be captured into a graph; (2) HtDTransfer's async H2D transfer of the quantized L/U - * blocks, so that transfer can run concurrently (copy engine) with kernels issued on the - * default stream (e.g. QuantizeDiagonalBlocksGPU, on the SM) instead of queueing behind them on - * the same stream. Because the two uses are mutually exclusive, sharing one stream (rather than - * a dedicated one per use) needs no extra synchronization between them. htd_event marks the end - * of the H2D transfer specifically, so the default-stream kernel that first reads the result - * (the quantized SpMV) can wait on it without a host-side block. ---*/ + /*--- Non-default stream, needed for two uses: (1) the preconditioner build/apply CUDA graphs + * below, since the legacy default stream cannot be captured into a graph; (2) HtDTransfer's + * async H2D transfer of the quantized L/U blocks, so that transfer can run concurrently (copy + * engine) with kernels issued on the default stream (e.g. QuantizeDiagonalBlocksGPU, on the SM) + * instead of queueing behind them on the same stream. The two are mutually exclusive for ILU + * (never quantized) but not for Q_LU_SGS, which uses both; sharing one stream still needs no + * extra synchronization, and in fact gives the right answer for free: the apply graph is + * launched into aux_stream, hence ordered after the transfer of the quantized blocks its + * kernels read. htd_event marks the end of the H2D transfer specifically, so a *default*-stream + * kernel that reads the result (the quantized SpMV) can wait on it without a host-side + * block. ---*/ mutable struct CUstream_st* aux_stream = nullptr; mutable struct CUevent_st* htd_event = nullptr; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c21116c6f3e..4a2f292bd61 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3645,6 +3645,31 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Multizone_Problem = YES; } + /*--- The solver vectors stay on the device but the halo exchange is host-side, so more than + * one rank would use stale halos. Use OpenMP for the host parts instead. ---*/ + if (Enable_Cuda && size > 1) { + SU2_MPI::Error("ENABLE_CUDA= YES is not supported with more than one MPI rank,\n" + " the halo exchange only happens on the host.\n" + " Use a single rank with OpenMP threads, e.g. 'SU2_CFD -t config.cfg'.", + CURRENT_FUNCTION); + } + + /*--- nvcc cannot compile the CoDiPack types, so the kernels are only built into the primal + * solver (see SU2_ENABLE_CUDA_KERNELS). Catch it here, not minutes into the run. ---*/ + if (Enable_Cuda) { +#ifndef SU2_ENABLE_CUDA_KERNELS +#ifdef HAVE_CUDA + SU2_MPI::Error("ENABLE_CUDA= YES is not available in the AD and direct differentiation solvers,\n" + " the CUDA kernels are only built into SU2_CFD.", + CURRENT_FUNCTION); +#else + SU2_MPI::Error("ENABLE_CUDA= YES but SU2 was not compiled with CUDA support,\n" + " reconfigure the build with -Denable-cuda=true.", + CURRENT_FUNCTION); +#endif +#endif + } + /*--- Set the default output files ---*/ if (!OptionIsSet("OUTPUT_FILES")){ nVolumeOutputFiles = 3; diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index e6292b4c10c..a4d9fbb1612 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -52,8 +52,11 @@ FORCEINLINE void RegularizePivot(ScalarType& pivot, unsigned long row, unsigned /*--- Common failure path for a device dispatch that is not available in this build/scalar type * combination, called with CURRENT_FUNCTION so the error names the right caller. ---*/ void GPUNotAvailable(const char* caller) { -#ifdef SU2_ENABLE_CUDA_KERNELS +#if defined(SU2_ENABLE_CUDA_KERNELS) SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", caller); +#elif defined(HAVE_CUDA) + /*--- AD build, the kernels are compiled out; normally rejected by CConfig::SetPostprocessing. ---*/ + SU2_MPI::Error("GPU acceleration is not available in the AD and direct differentiation solvers.", caller); #else SU2_MPI::Error( "ENABLE_CUDA is set to YES but SU2 was not compiled with CUDA support; " @@ -224,9 +227,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi * the host, so only plain (or quantized) Jacobi can keep them exclusively on the device. ---*/ jacobi_on_device = useCuda && (prec == JACOBI || prec == Q_JACOBI); #ifndef CODI_REVERSE_TYPE - /*--- Q_LU_SGS is still host-only. ---*/ - const bool quantized_offdiag_needed = - allow_quant && (prec == Q_JACOBI || prec == Q_IDENTITY || (prec == Q_LU_SGS && !useCuda)); + const bool quantized_offdiag_needed = allow_quant && (prec == Q_JACOBI || prec == Q_IDENTITY || prec == Q_LU_SGS); #else /*--- No quantization in adjoint mode for now because TransposeInPlace would get complicated. ---*/ const bool quantized_offdiag_needed = false; diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 723c8fb8268..7e6cf83d1b8 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -27,6 +27,7 @@ #include #include +#include #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" @@ -541,6 +542,33 @@ __global__ void QuantizedBlockLDU_SpMV_kernel( y[iRow * nVar + iVar] = sum; } +/*! + * \brief Instantiate the freshly captured \p graph into \p exec, or, when \p exec already holds a + * graph with the same topology, push the new node parameters into it in place. + * \note Re-capturing the topology is cheap, instantiating it is not: cudaGraphInstantiate + * allocates and builds the executable graph, at a cost that grows with the node count (one + * node per level here), so doing it on every call would cost more than simply launching the + * kernels and would defeat the purpose of using graphs at all. cudaGraphExecUpdate keeps the + * executable graph and only rewrites the kernel arguments that changed, which is what makes + * the graphs worth having on the flexible-FGMRES path where the vectors change every call. + * The full instantiation stays as the fallback for the first call and for the (unexpected) + * case of the topology actually changing. + */ +inline void InstantiateOrUpdateGraph(cudaGraphExec_t& exec, cudaGraph_t graph, const char* what) { + SU2_ZONE_SCOPED_N("Graph instantiate or update") + if (exec != nullptr) { + cudaGraphExecUpdateResultInfo info{}; + if (cudaGraphExecUpdate(exec, graph, &info) == cudaSuccess) return; + + /*--- A failed update is recoverable (we just instantiate again), but the runtime holds on to + * the error, so consume it before the next gpuErrChk mistakes it for a real failure. ---*/ + cudaGetLastError(); + gpuErrChk(cudaGraphExecDestroy(exec)); + exec = nullptr; + } + gpuErrChk(cudaGraphInstantiate(&exec, graph, nullptr, nullptr, 0)); +} + } // namespace template @@ -692,15 +720,13 @@ void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector::ComputeILUPreconditionerGPU(const CSysVector::ComputeILUPreconditionerGPU(const CSysVector +template __global__ void LU_SGS_ForwardKernel(const su2uint* __restrict__ level_idx, unsigned long level_begin, unsigned long level_size, unsigned long nVar, DeviceLDU M, const QuantType* __restrict__ q_l, const QuantScaleType* __restrict__ q_scale_l, const ScalarType* __restrict__ invD, const ScalarType* __restrict__ vec, - ScalarType* __restrict__ prod, bool quantized_mode) { + ScalarType* __restrict__ prod) { if (blockIdx.x >= level_size) return; const unsigned long iRow = level_idx[level_begin + blockIdx.x]; @@ -760,7 +786,7 @@ __global__ void LU_SGS_ForwardKernel(const su2uint* __restrict__ level_idx, unsi auto* aux = partial + blockSize; // skip nVar * nVar threads, serves nVar threads // Compute L.x* - if (quantized_mode) { + if constexpr (Quantized) { partial[tid] = QuantizedDeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_l, M.col_ind_l, q_l, q_scale_l, prod); } else { partial[tid] = DeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_l, M.col_ind_l, M.l, prod); @@ -781,12 +807,12 @@ __global__ void LU_SGS_ForwardKernel(const su2uint* __restrict__ level_idx, unsi * \brief Exact backward substitution for the rows of one level, x* = D^{-1}.(D.x* - U.x) = x* - D^{-1}.U.x * \note See notes in IluBackwardKernel for more details */ -template +template __global__ void LU_SGS_BackwardKernel(const su2uint* __restrict__ level_idx, unsigned long level_begin, unsigned long level_size, unsigned long nRows, unsigned long nVar, DeviceLDU M, const QuantType* __restrict__ q_u, const QuantScaleType* __restrict__ q_scale_u, const ScalarType* __restrict__ invD, - ScalarType* __restrict__ prod, bool quantized_mode) { + ScalarType* __restrict__ prod) { if (blockIdx.x >= level_size) return; const unsigned long iRow = level_idx[level_begin + blockIdx.x]; @@ -799,7 +825,7 @@ __global__ void LU_SGS_BackwardKernel(const su2uint* __restrict__ level_idx, uns auto* aux = partial + blockSize; // skip nVar * nVar threads, serves nVar threads // Compute U.x - if (quantized_mode) { + if constexpr (Quantized) { partial[tid] = QuantizedDeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_u, M.col_ind_u, q_u, q_scale_u, prod, nRows); } else { partial[tid] = DeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_u, M.col_ind_u, M.u, prod, nRows); @@ -864,25 +890,35 @@ void CSysMatrix::ComputeLU_SGSForwardGPU(const CSysVector<<>>(d_precond_level_idx, begin, size, nVar, M, d_q_blocks.l, d_q_scale.l, d_invM, d_vec, d_prod, quantized_mode); + /*--- Forward substitution: compute x* = D^{-1}.(vec - L.x*). Whether the off-diagonal blocks + * are quantized is fixed for the lifetime of the matrix (Initialize decides it from the + * preconditioner type), so it selects the kernel instantiation here rather than being tested + * by every thread: inside the kernel it is a compile-time constant and the unused branch is + * not compiled at all. ---*/ + auto RecordSweep = [&](auto quantized) { + for (auto level = 0ul; level < nLevels; ++level) { + const auto begin = precond_level_ptr[level]; + const auto size = precond_level_ptr[level + 1] - begin; + if (size == 0) continue; + LU_SGS_ForwardKernel + <<>>(d_precond_level_idx, begin, size, nVar, M, d_q_blocks.l, + d_q_scale.l, d_invM, d_vec, d_prod); + } + }; + if (quantized_mode) { + RecordSweep(std::true_type{}); + } else { + RecordSweep(std::false_type{}); } gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); - gpuErrChk(cudaGraphInstantiate(&precond_fwd_graph_exec, graph, nullptr, nullptr, 0)); + InstantiateOrUpdateGraph(precond_fwd_graph_exec, graph, "LU-SGS forward"); gpuErrChk(cudaGraphDestroy(graph)); precond_fwd_graph_vec = d_vec; precond_fwd_graph_prod = d_prod; @@ -920,26 +956,33 @@ void CSysMatrix::ComputeLU_SGSBackwardGPU(CSysVector& pr /*--- Second part of the symmetric iteration: (D+U).x_(1) = D.x* ---*/ if (precond_bwd_graph_exec == nullptr || precond_bwd_graph_prod != d_prod) { - if (precond_bwd_graph_exec != nullptr) { - gpuErrChk(cudaGraphExecDestroy(precond_bwd_graph_exec)); - precond_bwd_graph_exec = nullptr; - } + SU2_ZONE_SCOPED_N("LU-SGS bwd graph recapture") cudaGraph_t graph; gpuErrChk(cudaStreamBeginCapture(aux_stream, cudaStreamCaptureModeThreadLocal)); const auto nLevels = precond_level_ptr.size() - 1; - /*--- Backward substitution: compute x* = D^{-1}.(D.x* - U.x) = x* - D^{-1}.U.x ---*/ - for (auto level = nLevels; level > 0;) { - --level; - const auto begin = precond_level_ptr[level]; - const auto size = precond_level_ptr[level + 1] - begin; - if (size == 0) continue; - LU_SGS_BackwardKernel<<>>(d_precond_level_idx, begin, size, nPointDomain, nVar, M, d_q_blocks.u, d_q_scale.u, d_invM, d_prod, quantized_mode); + /*--- Backward substitution: compute x* = D^{-1}.(D.x* - U.x) = x* - D^{-1}.U.x. Quantization + * selects the kernel instantiation, see the forward sweep. ---*/ + auto RecordSweep = [&](auto quantized) { + for (auto level = nLevels; level > 0;) { + --level; + const auto begin = precond_level_ptr[level]; + const auto size = precond_level_ptr[level + 1] - begin; + if (size == 0) continue; + LU_SGS_BackwardKernel + <<>>(d_precond_level_idx, begin, size, nPointDomain, nVar, M, + d_q_blocks.u, d_q_scale.u, d_invM, d_prod); + } + }; + if (quantized_mode) { + RecordSweep(std::true_type{}); + } else { + RecordSweep(std::false_type{}); } gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); - gpuErrChk(cudaGraphInstantiate(&precond_bwd_graph_exec, graph, nullptr, nullptr, 0)); + InstantiateOrUpdateGraph(precond_bwd_graph_exec, graph, "LU-SGS backward"); gpuErrChk(cudaGraphDestroy(graph)); precond_bwd_graph_prod = d_prod; diff --git a/SU2_CFD/include/integration/CNewtonIntegration.hpp b/SU2_CFD/include/integration/CNewtonIntegration.hpp index da8ff67930a..25bb3139777 100644 --- a/SU2_CFD/include/integration/CNewtonIntegration.hpp +++ b/SU2_CFD/include/integration/CNewtonIntegration.hpp @@ -152,6 +152,18 @@ class CNewtonIntegration final : public CIntegration { template::value> = 0> inline unsigned long Preconditioner_impl(const CSysVector& u, CSysVector& v, unsigned long iters, Scalar& eps) const { + /*--- Unlike the matrix-free outer product this is a CSysMatrix operation, it can run on the + * device. The outer Krylov vectors are host resident, hence the transfers. ---*/ + unsigned long nIters = 0; + ApplyPreconditionerOnDevice(u, v, config->GetCUDA(), [&] { nIters = PreconditionerApply(u, v, iters, eps); }); + return nIters; + } + + /*! + * \brief The preconditioner on its own, or a nested solve with the approximate Jacobian. + */ + inline unsigned long PreconditionerApply(const CSysVector& u, CSysVector& v, + unsigned long iters, Scalar& eps) const { const auto inner_solver = config->GetKind_Linear_Solver_Inner(); if (iters == 0 || (iters == 1 && inner_solver == LINEAR_SOLVER_INNER::SMOOTHER)) { diff --git a/SU2_CFD/src/integration/CNewtonIntegration.cpp b/SU2_CFD/src/integration/CNewtonIntegration.cpp index feb07924caf..e745112fcb4 100644 --- a/SU2_CFD/src/integration/CNewtonIntegration.cpp +++ b/SU2_CFD/src/integration/CNewtonIntegration.cpp @@ -211,7 +211,18 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** solvers[FLOW_SOL]->PrepareImplicitIteration(geometry, solvers, config); - if (preconditioner) preconditioner->Build(); + if (preconditioner) { + /*--- The Jacobian is normally uploaded by CSysMatrixVectorProduct, but the outer product + * here is matrix free, so nothing else would upload it for Build(). ---*/ +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + if (config->GetCUDA()) { + SU2_DEVICE_REGION(solvers[FLOW_SOL]->Jacobian.HtDTransfer();) + } + } +#endif + preconditioner->Build(); + } auto CopyLinSysRes = [&](int sign, auto& dst) { SU2_OMP_FOR_STAT(omp_chunk_size) @@ -307,12 +318,9 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** solvers[FLOW_SOL]->Postprocessing(geometry, solvers, config, MESH_0); - SU2_OMP_MASTER { - solvers[FLOW_SOL]->Pressure_Forces(geometry, config); - solvers[FLOW_SOL]->Momentum_Forces(geometry, config); - solvers[FLOW_SOL]->Friction_Forces(geometry, config); - } - END_SU2_OMP_MASTER + solvers[FLOW_SOL]->Pressure_Forces(geometry, config); + solvers[FLOW_SOL]->Momentum_Forces(geometry, config); + solvers[FLOW_SOL]->Friction_Forces(geometry, config); /*--- At the end of the startup period the CFL is reset to the initial value. ---*/ diff --git a/meson.build b/meson.build index 78379a45507..e7efa33d1a6 100644 --- a/meson.build +++ b/meson.build @@ -19,7 +19,38 @@ python = pymod.find_installation() if get_option('enable-cuda') add_languages('cuda') - add_global_arguments('-arch=sm_89', language : 'cuda') + + # Compute capability to generate code for. Anything other than 'auto' is passed to nvcc as it + # is, so 'native', 'all-major', 'compute_XY', etc. all work. + cuda_arch = get_option('cuda-arch') + if cuda_arch == 'auto' + cuda_cc = '' + # Ships with the toolkit and prints e.g. "89"; it is what nvcc's own -arch=native uses. + cuda_query = find_program('__nvcc_device_query', required : false) + if cuda_query.found() + cuda_probe = run_command(cuda_query, check : false) + if cuda_probe.returncode() == 0 + cuda_cc = cuda_probe.stdout().strip().split('\n')[0].strip() + endif + endif + # Fall back to the driver, which prints e.g. "8.9". + if cuda_cc == '' + cuda_query = find_program('nvidia-smi', required : false) + if cuda_query.found() + cuda_probe = run_command(cuda_query, '--query-gpu=compute_cap', '--format=csv,noheader', check : false) + if cuda_probe.returncode() == 0 + cuda_cc = cuda_probe.stdout().strip().split('\n')[0].strip().replace('.', '') + endif + endif + endif + if cuda_cc == '' + error('Could not detect the CUDA compute capability of this machine (no GPU visible?), ' + + 'set it explicitly, for example -Dcuda-arch=sm_89.') + endif + cuda_arch = 'sm_' + cuda_cc + message('Detected CUDA compute capability, building for ' + cuda_arch) + endif + add_global_arguments('-arch=' + cuda_arch, 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') diff --git a/meson_options.txt b/meson_options.txt index 6adad19386b..db3bb47f552 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -11,6 +11,7 @@ option('enable-mkl', type : 'boolean', value : false, description: 'enable Intel 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('cuda-arch', type : 'string', value : 'auto', description: 'CUDA compute capability to generate code for, "auto" detects the GPU in this machine, otherwise passed to nvcc verbatim (e.g. sm_89, native, all-major)') 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('custom-mpi', type : 'boolean', value : false, description: 'enable MPI assuming the compiler and/or env vars give the correct include dirs and linker args.')