From 539b77acc0f3e1da7bbd8c86ecf8543d79ac153f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 5 Sep 2026 22:50:13 -0700 Subject: [PATCH 1/3] allow Q LU-SGS with CUDA, tweak the graph capture --- Common/include/linear_algebra/CSysMatrix.hpp | 26 ++-- Common/src/linear_algebra/CSysMatrix.cpp | 4 +- Common/src/linear_algebra/CSysMatrixGPU.cu | 141 ++++++++++++++----- 3 files changed, 119 insertions(+), 52 deletions(-) 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/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index e6292b4c10c..b49e1a4a4ab 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -224,9 +224,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..0c779e25334 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -27,6 +27,8 @@ #include #include +#include +#include #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" @@ -541,6 +543,54 @@ __global__ void QuantizedBlockLDU_SpMV_kernel( y[iRow * nVar + iVar] = sum; } +/*! + * \brief Report that a preconditioner CUDA graph had to be fully re-instantiated. + * \note This is the expensive path (cudaGraphInstantiate allocates and builds the executable + * graph) and it should only ever happen on the first call, since the topology is fixed by + * the level structure. If it shows up repeatedly the level structure is changing behind our + * back and the graphs are costing more than they save. Only a few occurrences are printed, + * from the master rank, to avoid flooding the output. + */ +inline void ReportGraphInstantiate(const char* what) { + constexpr int maxReports = 5; + static int nReports = 0; + if (nReports >= maxReports || SU2_MPI::GetRank() != MASTER_NODE) return; + ++nReports; + + std::cout << "Warning: re-instantiating the " << what << " CUDA graph, its topology changed." << std::endl; + if (nReports == maxReports) { + std::cout << "Further CUDA graph re-instantiation messages will not be printed." << std::endl; + } +} + +/*! + * \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(); + ReportGraphInstantiate(what); + gpuErrChk(cudaGraphExecDestroy(exec)); + exec = nullptr; + } + gpuErrChk(cudaGraphInstantiate(&exec, graph, nullptr, nullptr, 0)); +} + } // namespace template @@ -692,15 +742,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 +808,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 +829,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 +847,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 +912,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 +978,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; From e72a1d2904ae997a4fcfe18386a0b2a18dbeedc9 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 6 Sep 2026 20:07:20 -0700 Subject: [PATCH 2/3] Reject unsupported CUDA combinations, run the NK preconditioner on the GPU Three CUDA combinations that silently produced wrong results now fail at config time instead: - More than one MPI rank. The solver keeps its work vectors on the device for the whole Krylov iteration (CSysSolve::UploadSystem) while the halo exchange is host-side, so a partitioned run converges to the wrong answer. - The AD and direct differentiation solvers. nvcc cannot compile the CoDiPack types, so the kernels are only linked into the primal libraries (SU2_ENABLE_CUDA_KERNELS). This used to fail minutes into the run with a message claiming CUDA had not been compiled in, which is misleading for a build configured with -Denable-cuda=true; GPUNotAvailable now tells the two cases apart as well. Newton-Krylov was the third, and is fixed rather than rejected. Its outer matrix-vector product is matrix free and host resident, but the preconditioner is an ordinary CSysMatrix operation that the device can do. It went through neither CSysSolve::Solve (so the vectors were never uploaded) nor CSysMatrixVectorProduct (so the Jacobian was never uploaded), and every linear solve exited immediately having done nothing, with all preconditioners giving identical results. The Jacobian is now uploaded before Build(), and ApplyPreconditionerOnDevice (the mirror of the existing ApplyPreconditionerOnHost) transfers the vectors around the apply. Device expressions are enabled for the duration so a nested Krylov solve offloads too: verified with Tracy that all of its matrix-vector products and inner preconditioner applies run on the device, for the nested FGMRES, BCGSTAB and SMOOTHER inner solvers. On ONERA M6 RANS the GPU results track the CPU ones to 5-6 significant figures and are 1.3-1.9x faster per iteration. Also drops the SU2_OMP_MASTER around the force calculations, which #2870 parallelized with OpenMP. Co-Authored-By: Claude Opus 5 --- .../linear_algebra/CPreconditioner.hpp | 27 +++++++++++++++++++ Common/src/CConfig.cpp | 25 +++++++++++++++++ Common/src/linear_algebra/CSysMatrix.cpp | 5 +++- .../integration/CNewtonIntegration.hpp | 12 +++++++++ .../src/integration/CNewtonIntegration.cpp | 22 ++++++++++----- 5 files changed, 83 insertions(+), 8 deletions(-) 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/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 b49e1a4a4ab..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; " 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. ---*/ From 37c568947fc35d140df6180dde58aac486e6c3b6 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 6 Sep 2026 21:52:07 -0700 Subject: [PATCH 3/3] auto detect cuda [skip ci] --- Common/src/linear_algebra/CSysMatrixGPU.cu | 22 --------------- meson.build | 33 +++++++++++++++++++++- meson_options.txt | 1 + 3 files changed, 33 insertions(+), 23 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 0c779e25334..7e6cf83d1b8 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -27,7 +27,6 @@ #include #include -#include #include #include "../../include/linear_algebra/CMatrixInverse.hpp" @@ -543,26 +542,6 @@ __global__ void QuantizedBlockLDU_SpMV_kernel( y[iRow * nVar + iVar] = sum; } -/*! - * \brief Report that a preconditioner CUDA graph had to be fully re-instantiated. - * \note This is the expensive path (cudaGraphInstantiate allocates and builds the executable - * graph) and it should only ever happen on the first call, since the topology is fixed by - * the level structure. If it shows up repeatedly the level structure is changing behind our - * back and the graphs are costing more than they save. Only a few occurrences are printed, - * from the master rank, to avoid flooding the output. - */ -inline void ReportGraphInstantiate(const char* what) { - constexpr int maxReports = 5; - static int nReports = 0; - if (nReports >= maxReports || SU2_MPI::GetRank() != MASTER_NODE) return; - ++nReports; - - std::cout << "Warning: re-instantiating the " << what << " CUDA graph, its topology changed." << std::endl; - if (nReports == maxReports) { - std::cout << "Further CUDA graph re-instantiation messages will not be printed." << std::endl; - } -} - /*! * \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. @@ -584,7 +563,6 @@ inline void InstantiateOrUpdateGraph(cudaGraphExec_t& exec, cudaGraph_t graph, c /*--- 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(); - ReportGraphInstantiate(what); gpuErrChk(cudaGraphExecDestroy(exec)); exec = nullptr; } 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.')