diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 1b2a45dc3ffe..446a158000d5 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -646,6 +646,9 @@ class CConfig { unsigned long Linear_Solver_Prec_Threads; /*!< \brief Number of threads per rank for ILU and LU_SGS preconditioners. */ unsigned short Linear_Solver_ILU_n; /*!< \brief ILU fill=in level. */ bool Linear_Solver_ILU_levels; /*!< \brief Use level scheduling for OMP parallelization of ILU. */ + /*!< \brief Colored-iterative sweep counts for the GPU ILU preconditioner: [0] builds the + * factorization (Gauss-Seidel), [1]/[2] apply it (Jacobi, forward/backward triangular solve). */ + array Linear_Solver_ILU_GPU_Sweeps{{1, 2, 2}}; su2double SemiSpan; /*!< \brief Wing Semi span. */ su2double MSW_Alpha; /*!< \brief Coefficient for blending states in the MSW scheme. */ su2double Roe_Kappa; /*!< \brief Relaxation of the Roe scheme. */ @@ -4400,6 +4403,12 @@ class CConfig { */ bool GetLinear_Solver_ILU_levels(void) const { return Linear_Solver_ILU_levels; } + /*! + * \brief Get the [build, forward, backward] colored-iterative sweep counts for the GPU ILU + * preconditioner, see Linear_Solver_ILU_GPU_Sweeps. + */ + array GetLinear_Solver_ILU_GPU_Sweeps(void) const { return Linear_Solver_ILU_GPU_Sweeps; } + /*! * \brief Get restart frequency of the linear solver for the implicit formulation. * \return Restart frequency of the linear solver for the implicit formulation. diff --git a/Common/include/geometry/meshreader/CSU2BinaryMeshReaderBase.hpp b/Common/include/geometry/meshreader/CSU2BinaryMeshReaderBase.hpp index f4b71d94e64e..ca66e6aa80c8 100644 --- a/Common/include/geometry/meshreader/CSU2BinaryMeshReaderBase.hpp +++ b/Common/include/geometry/meshreader/CSU2BinaryMeshReaderBase.hpp @@ -10,7 +10,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFEM.hpp b/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFEM.hpp index bec54f75fcb7..d4112d09e589 100644 --- a/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFEM.hpp +++ b/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFEM.hpp @@ -10,7 +10,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFVM.hpp b/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFVM.hpp index 03ba796289bb..894bd41486d8 100644 --- a/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFVM.hpp +++ b/Common/include/geometry/meshreader/CSU2BinaryMeshReaderFVM.hpp @@ -10,7 +10,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/geometry/meshreader/CSU2MeshReaderBase.hpp b/Common/include/geometry/meshreader/CSU2MeshReaderBase.hpp index 9e5217f5c2b3..0f80b73b23d9 100644 --- a/Common/include/geometry/meshreader/CSU2MeshReaderBase.hpp +++ b/Common/include/geometry/meshreader/CSU2MeshReaderBase.hpp @@ -10,7 +10,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/grid_movement/CVolumetricMovementFactory.hpp b/Common/include/grid_movement/CVolumetricMovementFactory.hpp index 702decb25e97..970d780676f8 100644 --- a/Common/include/grid_movement/CVolumetricMovementFactory.hpp +++ b/Common/include/grid_movement/CVolumetricMovementFactory.hpp @@ -8,7 +8,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/linear_algebra/CMatrixInverse.hpp b/Common/include/linear_algebra/CMatrixInverse.hpp new file mode 100644 index 000000000000..f3d64630349d --- /dev/null +++ b/Common/include/linear_algebra/CMatrixInverse.hpp @@ -0,0 +1,99 @@ +/*! + * \file CMatrixInverse.hpp + * \brief Dense small-matrix inversion via Gauss-Jordan elimination, shared between the host + * (CSysMatrix::MatrixInverse) and device (CSysPreconditionerGPU.cu) implementations. + * \author F. Palacios, A. Bueno, T. Economon, P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include + +#ifdef __CUDACC__ +#define SU2_CUDA_HOST_DEVICE __host__ __device__ +#else +#define SU2_CUDA_HOST_DEVICE +#endif + +namespace SU2_LinAlg { + +/*! + * \brief Regularize a pivot that is too small to prevent divide-by-zero, on host and device + * this needs to clamp to the same value so that the two produce the same factors. + */ +template +SU2_CUDA_HOST_DEVICE inline void RegularizePivot(ScalarType& pivot) { + const float eps = 1e-12; +#ifdef __CUDA_ARCH__ + if (fabs(pivot) < eps) pivot = copysign(ScalarType(eps), pivot); +#else + if (std::abs(pivot) < eps) pivot = std::copysign(ScalarType(eps), pivot); +#endif +} + +/*! + * \brief Invert the \p nVar by \p nVar dense matrix \p matrix into \p inverse via Gauss-Jordan + * elimination with partial pivoting on the diagonal. + * \note \p matrix is used as scratch space and destroyed, \p inverse must not alias it. + */ +template +SU2_CUDA_HOST_DEVICE inline void MatrixInverse(unsigned long nVar, ScalarType* matrix, ScalarType* inverse) { +#define A(I, J) matrix[(I)*nVar + (J)] +#define M(I, J) inverse[(I)*nVar + (J)] + + /*--- Initialize the inverse with the identity. ---*/ + for (auto iVar = 0ul; iVar < nVar; iVar++) + for (auto jVar = 0ul; jVar < nVar; jVar++) M(iVar, jVar) = ScalarType(iVar == jVar); + + /*--- Transform system in Upper Matrix. ---*/ + for (auto iVar = 1ul; iVar < nVar; iVar++) { + for (auto jVar = 0ul; jVar < iVar; jVar++) { + RegularizePivot(A(jVar, jVar)); + + const ScalarType weight = A(iVar, jVar) / A(jVar, jVar); + for (auto kVar = jVar; kVar < nVar; kVar++) A(iVar, kVar) -= weight * A(jVar, kVar); + + /*--- At this stage M is lower triangular so not all cols need updating. ---*/ + for (auto kVar = 0ul; kVar <= jVar; kVar++) M(iVar, kVar) -= weight * M(jVar, kVar); + } + } + + /*--- Backwards substitution. ---*/ + for (auto iVar = nVar; iVar > 0ul;) { + iVar--; // unsigned type + for (auto jVar = iVar + 1; jVar < nVar; jVar++) + for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) -= A(iVar, jVar) * M(jVar, kVar); + + RegularizePivot(A(iVar, iVar)); + + for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) /= A(iVar, iVar); + } + +#undef A +#undef M +} + +} // namespace SU2_LinAlg + +#undef SU2_CUDA_HOST_DEVICE diff --git a/Common/include/linear_algebra/CMatrixVectorProduct.hpp b/Common/include/linear_algebra/CMatrixVectorProduct.hpp index 4069ff2fd006..52614a45770b 100644 --- a/Common/include/linear_algebra/CMatrixVectorProduct.hpp +++ b/Common/include/linear_algebra/CMatrixVectorProduct.hpp @@ -105,28 +105,6 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct { * \param[out] v - CSysVector that is the result of the product */ inline void operator()(const CSysVector& u, CSysVector& v) const override { - if (config->GetCUDA()) { -#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); - } -#elif defined(HAVE_CUDA) - SU2_MPI::Error( - "\nError in launching Matrix-Vector Product Function\nENABLE_CUDA is set to YES\nThe GPU kernels are not " - "part of the AD libraries, use the primal build for GPU acceleration", - CURRENT_FUNCTION); -#else - SU2_MPI::Error( - "\nError in launching Matrix-Vector Product Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA " - "options enabled in Meson to access GPU Functions", - CURRENT_FUNCTION); -#endif - } else { - matrix.MatrixVectorProduct(u, v, geometry, config); - } + matrix.MatrixVectorProduct(u, v, geometry, config); } }; diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index d28c729c4f0f..e23f5de381f2 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -40,7 +40,7 @@ /*! * \brief Applies a preconditioner that only has a host implementation to vectors that live * on the device: bring the input down, apply, put the result back. - * \note This is what keeps ILU, LU-SGS, Linelet and PaStiX usable on the GPU path. The + * \note This is what keeps LU-SGS, Linelet and PaStiX usable on the GPU path. The * transfers are issued by one thread with the team synchronized around them, the apply * itself is the normal OpenMP parallel host code. */ @@ -199,7 +199,8 @@ class CILUPreconditioner final : public CPreconditioner { * \param[out] v - CSysVector that is the result of the preconditioning. */ inline void operator()(const CSysVector& u, CSysVector& v) const override { - ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config); }); + /*--- No host bracket, ILU has a device implementation and ComputeILUPreconditioner dispatches to it. ---*/ + sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config); } /*! diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 3eddf892a7bf..cf25c264e246 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -250,6 +250,7 @@ class CSysMatrix { 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 gpu_ilu; /*!< \brief ILU factorization, device (values and pattern in GPU memory). */ ScalarType* d_invM = nullptr; /*!< \brief Device inverse diagonal blocks for the Jacobi preconditioner. */ /*--- Quantized off-diagonal storage (used when quantized_mode == true). ---*/ @@ -271,7 +272,12 @@ class CSysMatrix { * Populated by QuantizeDiagonalBlocks(). */ QuantType* q_blocks_d; /*!< \brief Same as q_blocks_l for the diagonal entries. */ - bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ + bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ + + /*!< \brief Whether the inverse diagonal blocks are only needed on the device. False for the + * Linelet preconditioner, which builds the Jacobi one but reads invM on the host. */ + bool jacobi_on_device = false; + const su2uint* l_to_u_transp; /*!< \brief L-entry index -> U-entry index of its transpose. */ const su2uint* u_to_l_transp; /*!< \brief U-entry index -> L-entry index of its transpose. */ @@ -284,9 +290,61 @@ class CSysMatrix { unsigned short ilu_fill_in; /*!< \brief Fill level for the ILU preconditioner. */ - /*!< \brief Level structure for alternative shared memory parallelization of ILU. */ + /*!< \brief Level structure for alternative shared memory parallelization of ILU. + * Rows within a level are independent, rows in level k only depend on rows in levels < k. + * The same table drives the forward (increasing level) and backward (decreasing level) + * sweeps, because the U pattern is the transpose of the L pattern. */ CCompressedSparsePatternUL levels_ilu; + /*!< \brief Coloring of the (domain-only) ILU dependency graph, used only by the GPU iterative + * factorization and triangular solves (see ilu_color_ptr / d_ilu_color_idx below); the + * host/OMP path is unaffected and keeps using levels_ilu exactly as before. */ + CCompressedSparsePatternUL color_ilu; + + /*!< \brief Number of colored Gauss-Seidel sweeps used to build and apply the ILU factorization + * on the device, see IluFactorColorKernel. Fixed per solve (not adaptive) so the result is + * reproducible; set from config in Initialize(). */ + array ilu_gpu_sweeps{{1, 2, 2}}; + + /*--- Coloring of the ILU dependency graph: unlike levels_ilu, a color is a true independent + * set (no dependency between same-colored rows in either direction), so far fewer, wider + * colors are needed than levels, but a color launch is only exact as one step of an iterative + * refinement (see BuildILUPreconditionerGPU and ComputeILUPreconditionerGPU), not a single + * pass — this does not change the elimination order/pattern, so the factorization and both + * triangular solves converge to the exact same result levels_ilu would give, just reached by + * iterating instead of substituting. ---*/ + vector ilu_color_ptr; /*!< \brief Start of each color in d_ilu_color_idx, size nColors+1. */ + su2uint* d_ilu_color_idx = nullptr; /*!< \brief Row indices, grouped by color. */ + + /*!< \brief Fixed right-hand side of the backward triangular solve (the forward-solve result), + * kept separate from the evolving prod buffer. A color visits every row once per sweep, so + * after the first sweep prod[iRow] holds a solution *estimate*, not the right-hand side + * anymore; reading the right-hand side back out of prod past the first sweep would silently + * solve the wrong equation. Size nPointDomain*nVar, allocated once in Initialize(). */ + ScalarType* d_ilu_backward_rhs = nullptr; + + /*--- The per-color kernel launch sequence (ilu_gpu_sweeps[0] passes for the factorization, + * ilu_gpu_sweeps[1] / ilu_gpu_sweeps[2] passes for the forward/backward triangular solves) + * is identical on every call: same grid/block sizes, same device pointers (all fixed members, + * allocated once). It is captured once into a CUDA graph and replayed, which removes + * host-side launch overhead without changing the parallelization (unlike a persistent + * cooperative-groups kernel, this does not cap per-color parallelism to the occupancy-resident + * block count). ---*/ + /*--- Types are forward-declared as opaque structs (matching the real cudaGraphExec_t / + * cudaStream_t typedefs) so this header does not need to include the CUDA runtime. ---*/ + mutable struct CUgraphExec_st* ilu_build_graph_exec = nullptr; + mutable struct CUgraphExec_st* ilu_apply_graph_exec = nullptr; + mutable const ScalarType* ilu_apply_graph_vec = nullptr; /*!< \brief Pointers the apply graph + * was captured with, to detect when + * it must be recaptured. */ + mutable ScalarType* ilu_apply_graph_prod = nullptr; + /*--- The legacy default stream cannot be captured into a graph, so the ILU graphs are + * captured and replayed on this dedicated stream instead; every launch on it is followed by + * a sync back to the host before control returns to the rest of the (single-stream) solver, + * so this does not change execution order relative to everything else, which stays on the + * default stream. ---*/ + mutable struct CUstream_st* ilu_stream = nullptr; + ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/ @@ -506,6 +564,37 @@ class CSysMatrix { * ScalarType buffer and delegates to the scalar GaussElimination overload. */ inline void QuantizedGaussElimination(unsigned long block_i, ScalarType* rhs) const; + /*--- Hooks for GPU versions (implemented is in CSysMatrixGPU.cu). ---*/ + + /*! + * \brief Performs the product of a sparse matrix by a CSysVector on the device. + */ + void MatrixVectorProductGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, + const CConfig* config) const; + + /*! + * \brief Build the Jacobi preconditioner on the device, from the device copy of the matrix. + * \note Requires the device matrix to be up to date, see HtDTransfer. + */ + void BuildJacobiPreconditionerGPU(); + + /*! + * \brief Apply the Jacobi preconditioner on the GPU/device side. + */ + void ComputeJacobiPreconditionerGPU(const CSysVector& vec, CSysVector& prod, + CGeometry* geometry, const CConfig* config) const; + + /*! + * \brief Build the ILU preconditioner on the device, from the device copy of the matrix. + * \note Requires the device matrix to be up to date, see HtDTransfer. + */ + void BuildILUPreconditionerGPU(); + + /*! + * \brief Apply the ILU preconditioner on the device. + */ + void ComputeILUPreconditionerGPU(const CSysVector& vec, CSysVector& prod) const; + public: /*! * \brief Constructor of the class. @@ -1043,49 +1132,6 @@ class CSysMatrix { void MatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; - /*! - * \brief Performs the product of a sparse matrix by a CSysVector. - * \param[in] vec - CSysVector to be multiplied by the sparse matrix A. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] prod - Result of the product. - */ - void GPUMatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, - const CConfig* config) const; - - /*! - * \brief Performs first step of the LU_SGS Preconditioner building - * \param[in] vec - CSysVector to be multiplied by the sparse matrix A. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] prod - Result of the product. - */ - void GPUFirstSymmetricIteration(ScalarType& vec, ScalarType& prod, CGeometry* geometry, const CConfig* config) const; - - /*! - * \brief Performs second step of the LU_SGS Preconditioner building - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] prod - Result of the product. - */ - void GPUSecondSymmetricIteration(ScalarType& prod, CGeometry* geometry, const CConfig* config) const; - - /*! - * \brief Performs Gaussian Elimination between diagional blocks of the matrix and the prod vector - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] prod - Result of the product. - */ - void GPUGaussElimination(ScalarType& prod, CGeometry* geometry, const CConfig* config) const; - - /*! - * \brief Multiply CSysVector by the preconditioner all of which are stored on the device - * \param[in] vec - CSysVector to be multiplied by the preconditioner. - * \param[out] prod - Result of the product A*vec. - */ - void GPUComputeLU_SGSPreconditioner(ScalarType& vec, ScalarType& prod, CGeometry* geometry, - const CConfig* config) const; - /*! * \brief Build the Jacobi preconditioner. */ @@ -1101,14 +1147,6 @@ 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/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh index eb727477f214..39268c10e5a3 100644 --- a/Common/include/linear_algebra/GPUComms.cuh +++ b/Common/include/linear_algebra/GPUComms.cuh @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * -* Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) +* Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/toolboxes/SwapBytes.hpp b/Common/include/toolboxes/SwapBytes.hpp index 033dc49dd8dd..c09a91e8b12a 100644 --- a/Common/include/toolboxes/SwapBytes.hpp +++ b/Common/include/toolboxes/SwapBytes.hpp @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/include/toolboxes/random_toolbox.hpp b/Common/include/toolboxes/random_toolbox.hpp index 6ad052d9e52d..f2ad39a7aa2d 100644 --- a/Common/include/toolboxes/random_toolbox.hpp +++ b/Common/include/toolboxes/random_toolbox.hpp @@ -8,7 +8,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a9ca92ab1454..67850acae3da 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1972,6 +1972,9 @@ void CConfig::SetConfig_Options() { addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", Linear_Solver_ILU_n, 0); /* DESCRIPTION: Use level scheduling for OMP parallelization of the ILU preconditioner */ addBoolOption("LINEAR_SOLVER_ILU_LEVEL_SCHEDULING", Linear_Solver_ILU_levels, false); + /* DESCRIPTION: Colored-iterative sweep counts for the GPU ILU preconditioner: (build, forward, backward) */ + addUShortArrayOption("LINEAR_SOLVER_ILU_GPU_SWEEPS", Linear_Solver_ILU_GPU_Sweeps.size(), false, + Linear_Solver_ILU_GPU_Sweeps.data()); /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */ addUnsignedLongOption("LINEAR_SOLVER_RESTART_FREQUENCY", Linear_Solver_Restart_Frequency, 10); /* DESCRIPTION: Number of vectors used for deflated restarts */ @@ -4139,6 +4142,15 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } + for (auto n : Linear_Solver_ILU_GPU_Sweeps) { + if (n == 0) { + SU2_MPI::Error( + "LINEAR_SOLVER_ILU_GPU_SWEEPS entries must all be at least 1 (build, forward, backward); " + "0 sweeps never factorizes/solves.", + CURRENT_FUNCTION); + } + } + Radiation = (Kind_Radiation != RADIATION_MODEL::NONE); /*--- Check for unsupported features. ---*/ diff --git a/Common/src/geometry/meshreader/CSU2BinaryMeshReaderBase.cpp b/Common/src/geometry/meshreader/CSU2BinaryMeshReaderBase.cpp index ce7d03cb3e3c..9e0e9c02a65f 100644 --- a/Common/src/geometry/meshreader/CSU2BinaryMeshReaderBase.cpp +++ b/Common/src/geometry/meshreader/CSU2BinaryMeshReaderBase.cpp @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFEM.cpp b/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFEM.cpp index e3d4ceed2e6c..cbbc4826e9d1 100644 --- a/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFEM.cpp +++ b/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFEM.cpp @@ -10,7 +10,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFVM.cpp b/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFVM.cpp index 96905c0512a2..4c235416ebdc 100644 --- a/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFVM.cpp +++ b/Common/src/geometry/meshreader/CSU2BinaryMeshReaderFVM.cpp @@ -10,7 +10,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/src/geometry/meshreader/CSU2MeshReaderBase.cpp b/Common/src/geometry/meshreader/CSU2MeshReaderBase.cpp index ef919ef724b1..752d9cbdfd0d 100644 --- a/Common/src/geometry/meshreader/CSU2MeshReaderBase.cpp +++ b/Common/src/geometry/meshreader/CSU2MeshReaderBase.cpp @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index a0013e243002..35eb70ea6411 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -28,6 +28,7 @@ #include "../../include/linear_algebra/CSysMatrix.inl" #include "../../include/geometry/CGeometry.hpp" +#include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/toolboxes/allocation_toolbox.hpp" #include @@ -131,6 +132,20 @@ CSysMatrix::~CSysMatrix() { GPUMemoryAllocation::gpu_free(gpu.row_ptr_u); GPUMemoryAllocation::gpu_free(gpu.col_ind_u); GPUMemoryAllocation::gpu_free(d_invM); + GPUMemoryAllocation::gpu_free(gpu_ilu.d); + GPUMemoryAllocation::gpu_free(gpu_ilu.l); + GPUMemoryAllocation::gpu_free(gpu_ilu.u); + GPUMemoryAllocation::gpu_free(gpu_ilu.row_ptr_l); + GPUMemoryAllocation::gpu_free(gpu_ilu.col_ind_l); + GPUMemoryAllocation::gpu_free(gpu_ilu.row_ptr_u); + GPUMemoryAllocation::gpu_free(gpu_ilu.col_ind_u); + GPUMemoryAllocation::gpu_free(d_ilu_color_idx); + GPUMemoryAllocation::gpu_free(d_ilu_backward_rhs); +#ifdef SU2_ENABLE_CUDA_KERNELS + if (ilu_build_graph_exec != nullptr) cudaGraphExecDestroy(ilu_build_graph_exec); + if (ilu_apply_graph_exec != nullptr) cudaGraphExecDestroy(ilu_apply_graph_exec); + if (ilu_stream != nullptr) cudaStreamDestroy(ilu_stream); +#endif } #ifdef USE_MKL @@ -182,6 +197,10 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi const bool ilu_needed = (prec == ILU); const bool diag_needed = (prec == JACOBI) || (prec == LINELET); + + /*--- Linelet also builds the Jacobi preconditioner but reads the inverse diagonal blocks on + * the host, so only plain Jacobi can keep them exclusively on the device. ---*/ + jacobi_on_device = useCuda && (prec == JACOBI); #ifndef CODI_REVERSE_TYPE const bool q_lus_needed = allow_quant && !useCuda && (prec == Q_LU_SGS); #else @@ -260,6 +279,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (ilu_needed) { ilu_fill_in = config->GetLinear_Solver_ILU_n(); + ilu_gpu_sweeps = config->GetLinear_Solver_ILU_GPU_Sweeps(); const auto& pat_ilu = geometry->GetSparsePattern(type, ilu_fill_in); ilu.row_ptr_l = pat_ilu.l.outerPtr(); @@ -272,6 +292,27 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (omp_get_max_threads() > 1 && config->GetLinear_Solver_ILU_levels()) { levels_ilu = computeLevels(pat_ilu.l); } + + /*--- Coloring for the GPU iterative factorization, see IluFactorColorKernel. Colors are + * true independent sets of the (domain-only, symmetric) dependency graph, computed the same + * way SU2 already colors edges/elements for OMP loops, just applied to the ILU pattern + * instead. This does not change the elimination order/pattern (nothing here affects L/U + * membership), only how the build is scheduled on the device. ---*/ + if (useCuda) { + std::vector adjPtr(nPointDomain + 1, 0); + std::vector adjIdx; + adjIdx.reserve(ilu.nnz_l + ilu.nnz_u); + for (auto i = 0ul; i < nPointDomain; ++i) { + adjPtr[i] = static_cast(adjIdx.size()); + for (auto k = ilu.row_ptr_l[i]; k < ilu.row_ptr_l[i + 1]; ++k) adjIdx.push_back(ilu.col_ind_l[k]); + for (auto k = ilu.row_ptr_u[i]; k < ilu.row_ptr_u[i + 1]; ++k) { + const auto j = ilu.col_ind_u[k]; + if (j < nPointDomain) adjIdx.push_back(static_cast(j)); + } + } + adjPtr[nPointDomain] = static_cast(adjIdx.size()); + color_ilu = colorSparsePattern(CCompressedSparsePatternUL(adjPtr, adjIdx), 1, true, false); + } } /*--- Preconditioners. ---*/ @@ -284,10 +325,40 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn); - if (useCuda && diag_needed) { + if (jacobi_on_device) { d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); } + if (useCuda && ilu_needed) { + /*--- The factors are built and used on the device, only the pattern and the level table + * are uploaded (once, here) because they do not change. ---*/ + gpu_ilu.nnz_l = ilu.nnz_l; + gpu_ilu.nnz_u = ilu.nnz_u; + gpu_ilu.d = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + gpu_ilu.l = GPUMemoryAllocation::gpu_alloc(ilu.nnz_l * nVar * nEqn * sizeof(ScalarType)); + gpu_ilu.u = GPUMemoryAllocation::gpu_alloc(ilu.nnz_u * nVar * nEqn * sizeof(ScalarType)); + gpu_ilu.row_ptr_l = GPUMemoryAllocation::gpu_alloc_cpy(ilu.row_ptr_l, (nPointDomain + 1) * sizeof(su2uint)); + gpu_ilu.col_ind_l = GPUMemoryAllocation::gpu_alloc_cpy(ilu.col_ind_l, ilu.nnz_l * sizeof(su2uint)); + gpu_ilu.row_ptr_u = GPUMemoryAllocation::gpu_alloc_cpy(ilu.row_ptr_u, (nPointDomain + 1) * sizeof(su2uint)); + gpu_ilu.col_ind_u = GPUMemoryAllocation::gpu_alloc_cpy(ilu.col_ind_u, ilu.nnz_u * sizeof(su2uint)); + + /*--- Flatten the coloring, the index type differs from the one of the pattern. It drives + * the factorization and both triangular solves on the device. ---*/ + std::vector color_idx; + color_idx.reserve(nPointDomain); + ilu_color_ptr.clear(); + ilu_color_ptr.push_back(0); + for (auto color = 0ul; color < color_ilu.getOuterSize(); ++color) { + for (auto k = 0ul; k < color_ilu.getNumNonZeros(color); ++k) { + color_idx.push_back(static_cast(color_ilu.getInnerIdx(color, k))); + } + ilu_color_ptr.push_back(static_cast(color_idx.size())); + } + d_ilu_color_idx = GPUMemoryAllocation::gpu_alloc_cpy(color_idx.data(), color_idx.size() * sizeof(su2uint)); + + d_ilu_backward_rhs = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * sizeof(ScalarType)); + } + /*--- Thread parallel initialization. ---*/ int num_threads = omp_get_max_threads(); @@ -688,15 +759,16 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver assert((matrix != inverse) && "Output cannot be the same as the input."); + /*--- Inversion ---*/ +#ifdef USE_MKL_LAPACK + // With MKL_DIRECT_CALL enabled, this is significantly faster than native code on Intel Architectures. #define M(I, J) inverse[(I)*nVar + (J)] - /*--- Initialize the inverse with the identity. ---*/ + /*--- Initialize the inverse with the identity, LAPACKE_?getrs solves for it as the rhs. ---*/ for (auto iVar = 0ul; iVar < nVar; iVar++) for (auto jVar = 0ul; jVar < nVar; jVar++) M(iVar, jVar) = ScalarType(iVar == jVar); +#undef M - /*--- Inversion ---*/ -#ifdef USE_MKL_LAPACK - // With MKL_DIRECT_CALL enabled, this is significantly faster than native code on Intel Architectures. lapack_int ipiv[MAXNVAR]; if constexpr (std::is_same_v) { LAPACKE_dgetrf(LAPACK_ROW_MAJOR, nVar, nVar, matrix, nVar, ipiv); @@ -707,38 +779,9 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver LAPACKE_sgetrs(LAPACK_ROW_MAJOR, 'N', nVar, nVar, matrix, nVar, ipiv, inverse, nVar); } #else -#define A(I, J) matrix[(I)*nVar + (J)] - - /*--- Transform system in Upper Matrix ---*/ - for (auto iVar = 1ul; iVar < nVar; iVar++) { - for (auto jVar = 0ul; jVar < iVar; jVar++) { - /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ - RegularizePivot(A(jVar, jVar), jVar, jVar, "MatrixInverse"); - - ScalarType weight = A(iVar, jVar) / A(jVar, jVar); - for (auto kVar = jVar; kVar < nVar; kVar++) A(iVar, kVar) -= weight * A(jVar, kVar); - - /*--- at this stage M is lower triangular so not all cols need updating ---*/ - for (auto kVar = 0ul; kVar <= jVar; kVar++) M(iVar, kVar) -= weight * M(jVar, kVar); - } - } - - /*--- Backwards substitution ---*/ - for (auto iVar = nVar; iVar > 0ul;) { - iVar--; // unsigned type - for (auto jVar = iVar + 1; jVar < nVar; jVar++) - for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) -= A(iVar, jVar) * M(jVar, kVar); - - /*--- Regularize diagonal if too small ---*/ - RegularizePivot(A(iVar, iVar), iVar, iVar, "DEBUG MatrixInverse backsubst"); - - for (auto kVar = 0ul; kVar < nVar; kVar++) { - M(iVar, kVar) /= A(iVar, iVar); - } - } -#undef A + /*--- Shared with the device implementation, see CMatrixInverse.hpp. ---*/ + SU2_LinAlg::MatrixInverse(nVar, matrix, inverse); #endif -#undef M } template @@ -771,6 +814,25 @@ template void CSysMatrix::MatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const { SU2_ZONE_SCOPED + + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + BEGIN_SU2_DEVICE_REGION + MatrixVectorProductGPU(vec, prod, geometry, config); + END_SU2_DEVICE_REGION + return; + } else { + SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION); + } +#else + SU2_MPI::Error( + "\nError in launching Matrix-Vector Product Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA " + "options enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ #ifndef NDEBUG if ((nEqn != vec.GetNVar()) || (nVar != prod.GetNVar())) { @@ -811,18 +873,11 @@ 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) { + if (jacobi_on_device) { #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 + SU2_DEVICE_REGION(BuildJacobiPreconditionerGPU();) + return; } else { SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION); } @@ -833,6 +888,12 @@ void CSysMatrix::BuildJacobiPreconditioner() { CURRENT_FUNCTION); #endif } + + /*--- 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 } template @@ -872,6 +933,23 @@ void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector void CSysMatrix::BuildILUPreconditioner() { SU2_ZONE_SCOPED + + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + SU2_DEVICE_REGION(BuildILUPreconditionerGPU();) + return; + } else { + SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION); + } +#else + SU2_MPI::Error( + "\nError in building ILU preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options " + "enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + const auto blockSize = nVar * nVar; ScalarType Lij[MAXNVAR * MAXNVAR], Lij_Ujk[MAXNVAR * MAXNVAR]; @@ -1002,6 +1080,23 @@ template void CSysMatrix::ComputeILUPreconditioner(const CSysVector& vec, 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(ComputeILUPreconditionerGPU(vec, prod);) + return; + } else { + SU2_MPI::Error("GPU acceleration is not supported for AD scalar types.", CURRENT_FUNCTION); + } +#else + SU2_MPI::Error( + "\nError in applying ILU preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options " + "enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + /*--- Coherent view of vectors. ---*/ SU2_OMP_BARRIER diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 381379a0d812..601ae385cf4e 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -1,7 +1,7 @@ /*! * \file CSysMatrixGPU.cu * \brief Implementations of Kernels and Functions for Matrix Operations on the GPU - * \author A. Raj + * \author A. Raj, Jesse Li, P. Gomes * \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public @@ -25,9 +25,309 @@ * License along with SU2. If not, see . */ -#include "../../include/linear_algebra/CSysMatrix.hpp" +#include + +#include "../../include/linear_algebra/CMatrixInverse.hpp" +#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; + } +} + +/*--- ILU. The factorization and both triangular solves are all scheduled by coloring: colors + * are true independent sets of the ILU dependency graph (no two same-colored rows depend on each + * other in either direction), so a color's rows can be processed with zero races in one kernel + * launch, but since a color is wider and less ordered than a level, one pass over all colors is + * only an approximation, not an exact result — several sweeps (repeated passes) are needed to + * converge to the same fixed point an exact level-scheduled algorithm would give in one pass. The + * rows of a color are scattered through the matrix, hence the indirection through the color + * table. Throughout, one CUDA block works on one row. ---*/ + +/*! + * \brief The pointers of an LDU-partitioned matrix, all in device memory. This mirrors the + * private CSysMatrix::LDU, which the kernels cannot name. + */ +template +struct DeviceLDU { + ScalarType* d; + ScalarType* l; + ScalarType* u; + const su2uint* row_ptr_l; + const su2uint* col_ind_l; + const su2uint* row_ptr_u; + const su2uint* col_ind_u; +}; + +/*! + * \brief Start of block (i,j), or nullptr if it is not a nonzero of the pattern. + */ +template +__device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, unsigned long nVar, + unsigned long block_i, unsigned long block_j) { + const auto blockSize = nVar * nVar; + if (block_i == block_j) return M.d + block_i * blockSize; + + const bool lower = block_j < block_i; + const auto* row_ptr = lower ? M.row_ptr_l : M.row_ptr_u; + const auto* col_ind = lower ? M.col_ind_l : M.col_ind_u; + auto* vals = lower ? M.l : M.u; + + for (auto k = row_ptr[block_i]; k < row_ptr[block_i + 1]; ++k) { + if (col_ind[k] == block_j) return vals + k * blockSize; + } + return nullptr; +} + +/*! + * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock. + * \note Grid: one block per row, blockDim.x == nVar*nVar (one thread per block entry, they + * stage the input in shared memory). Dynamic shared memory: nVar*nVar scalars. + */ +template +__global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, + const ScalarType* __restrict__ mat_d, ScalarType* __restrict__ invM) { + const unsigned long iRow = blockIdx.x; + if (iRow >= nRows) return; + + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + + /*--- The inversion destroys its input, so it cannot work on the matrix itself. ---*/ + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* work = reinterpret_cast(smem); + + work[tid] = mat_d[iRow * blockSize + tid]; + __syncthreads(); + + if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); +} + +/*! + * \brief Factorize the rows of one color, one sweep of an iterative (colored Gauss-Seidel) + * ILU factorization: same order/pattern as the exact level-scheduled algorithm (this + * does not change L/U membership, so it converges to the exact same fixed point), but + * colors are true independent sets (zero dependency between same-colored rows in either + * direction), so a color can be processed with zero races in far fewer, wider launches + * than the number of levels — at the cost of needing several sweeps (repeated passes + * over all colors) instead of one exact pass, because for a fixed order the level count + * is already the minimum number of race-free single-pass groups (Mirsky's theorem). + * \note Every visit of a row (there is one per sweep) resets it from the original matrix first + * (folding in the device version of the InitIluRow helper of BuildILUPreconditioner), + * because the elimination below is a re-evaluation of the row's defining equation using + * the current (possibly stale) values of other rows, not an incremental accumulation. + * Grid: one block per row of the color, blockDim.x == nVar*nVar (one thread per block + * entry, so that the small matrix products are one dot product per thread). Dynamic + * shared memory: 2*nVar*nVar scalars. + */ +template +__global__ void IluFactorColorKernel(const su2uint* __restrict__ color_idx, unsigned long color_begin, + unsigned long color_size, unsigned long nRows, unsigned long nVar, + DeviceLDU A, DeviceLDU M) { + if (blockIdx.x >= color_size) return; + + const unsigned long iRow = color_idx[color_begin + blockIdx.x]; + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + const auto iVar = tid / nVar, jVar = tid % nVar; + + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* Lij = reinterpret_cast(smem); + auto* work = Lij + blockSize; + + /*--- Reset this row to the raw matrix entries (device version of InitIluRow, but for one + * row instead of the whole matrix, since here it runs once per row per sweep). ---*/ + M.d[iRow * blockSize + tid] = A.d[iRow * blockSize + tid]; + auto scatter = [&](const su2uint* a_row_ptr, const su2uint* a_col_ind, const ScalarType* a_vals, + const su2uint* m_row_ptr, const su2uint* m_col_ind, ScalarType* m_vals) { + auto ka = a_row_ptr[iRow]; + const auto ka_end = a_row_ptr[iRow + 1]; + for (auto k = m_row_ptr[iRow]; k < m_row_ptr[iRow + 1]; ++k) { + const auto jPoint = m_col_ind[k]; + while (ka < ka_end && a_col_ind[ka] < jPoint) ++ka; + if (ka < ka_end && a_col_ind[ka] == jPoint) { + m_vals[k * blockSize + tid] = a_vals[ka * blockSize + tid]; + } else { + m_vals[k * blockSize + tid] = ScalarType(0); + } + } + }; + scatter(A.row_ptr_l, A.col_ind_l, A.l, M.row_ptr_l, M.col_ind_l, M.l); + scatter(A.row_ptr_u, A.col_ind_u, A.u, M.row_ptr_u, M.col_ind_u, M.u); + __syncthreads(); + + /*--- For this row (unknown), loop over its lower diagonal entries. ---*/ + for (auto kl = M.row_ptr_l[iRow]; kl < M.row_ptr_l[iRow + 1]; ++kl) { + /*--- All threads must be done with the previous entry: Lij is about to be overwritten, + * and the blocks of this row updated below are read here across threads. ---*/ + __syncthreads(); + + /*--- jPoint is the column index (jPoint < iRow). ---*/ + const unsigned long jPoint = M.col_ind_l[kl]; + + /*--- Multiply the block by the inverse of the corresponding diagonal block. ---*/ + auto* Block_ij = M.l + kl * blockSize; + const auto* invUjj = M.d + jPoint * blockSize; + + ScalarType sum = 0; + for (auto k = 0ul; k < nVar; ++k) sum += Block_ij[iVar * nVar + k] * invUjj[k * nVar + jVar]; + Lij[tid] = sum; + __syncthreads(); + + /*--- Lij holds Aij*inv(Ujj). Jump to the upper part of the jPoint row. ---*/ + for (auto ku = M.row_ptr_u[jPoint]; ku < M.row_ptr_u[jPoint + 1]; ++ku) { + /*--- Get the column index (kPoint > jPoint), halo columns are not factorized. ---*/ + const unsigned long kPoint = M.col_ind_u[ku]; + if (kPoint >= nRows) break; + + /*--- If Aik exists, update it: Aik -= Lij * Ujk ---*/ + auto* Block_ik = GetBlockILU(M, nVar, iRow, kPoint); + if (Block_ik == nullptr) continue; + + /*--- Block_ik cannot alias Block_ij because kPoint > jPoint. ---*/ + const auto* Ujk = M.u + ku * blockSize; + ScalarType prod = 0; + for (auto k = 0ul; k < nVar; ++k) prod += Lij[iVar * nVar + k] * Ujk[k * nVar + jVar]; + Block_ik[tid] -= prod; + } + + /*--- Store Lij in the lower triangular part, each thread only writes its own entry. ---*/ + Block_ij[tid] = Lij[tid]; + } + + /*--- Invert the diagonal entry, Uii, for the rows that depend on it. The loop above may have + * updated it (when kPoint == iRow), so the whole block has to be done first. ---*/ + __syncthreads(); + work[tid] = M.d[iRow * blockSize + tid]; + __syncthreads(); + if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, M.d + iRow * blockSize); +} + +/*! + * \brief One colored-Jacobi sweep of forward substitution for the rows of one color, + * (L+I).prod = vec. + * \note Same idea as IluFactorColorKernel/IluBackwardColorKernel: colors are true independent + * sets, so a color can be processed with zero races, but an L-neighbor in a different, + * not-yet-processed color this sweep is read at whatever value it currently holds (a + * previous sweep's, or the zero-initialized buffer on the very first sweep). Repeated full + * passes over all colors (see ComputeILUPreconditionerGPU) converge this to the same + * result the exact level-scheduled solve would give, because the update is a genuine + * residual reduction, not an incremental accumulation. One thread per block entry (like + * the factorization kernel), so the inner dot product over a neighbor block is spread + * across nVar threads instead of done serially by one; each thread accumulates its own + * (iVar,jVar) partial product across every neighbor with no synchronization at all, and + * only the final nVar-way reduction (summing over jVar for each iVar) needs one + * __syncthreads() — not one per neighbor. Grid: one block per row of the color, + * blockDim.x == nVar*nVar. Dynamic shared memory: nVar*nVar scalars. + */ +template +__global__ void IluForwardColorKernel(const su2uint* __restrict__ color_idx, unsigned long color_begin, + unsigned long color_size, unsigned long nVar, DeviceLDU M, + const ScalarType* __restrict__ vec, ScalarType* __restrict__ prod) { + if (blockIdx.x >= color_size) return; + + const unsigned long iRow = color_idx[color_begin + blockIdx.x]; + const unsigned long tid = threadIdx.x; + const auto iVar = tid / nVar, jVar = tid % nVar; + + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* partial = reinterpret_cast(smem); + + ScalarType acc = 0; + /*--- The columns of L are rows of previous levels, so prod is final for all of them. ---*/ + for (auto kl = M.row_ptr_l[iRow]; kl < M.row_ptr_l[iRow + 1]; ++kl) { + const unsigned long jPoint = M.col_ind_l[kl]; + const auto* blk = M.l + kl * nVar * nVar; + acc += blk[iVar * nVar + jVar] * prod[jPoint * nVar + jVar]; + } + partial[tid] = acc; + __syncthreads(); + + if (jVar == 0) { + ScalarType sum = vec[iRow * nVar + iVar]; + for (auto j = 0ul; j < nVar; ++j) sum -= partial[iVar * nVar + j]; + prod[iRow * nVar + iVar] = sum; + } +} + +/*! + * \brief One colored-Jacobi sweep of backward substitution for the rows of one color, + * U.prod = rhs. + * \note The right-hand side is read from a separate fixed array (\p rhs), not from \p prod + * itself: a color visits every row once per sweep, so after the first sweep prod[iRow] + * holds a solution *estimate*, not the right-hand side. Reading the right-hand side back + * out of prod past the first sweep would silently solve a different (and divergent) + * recurrence — see CSysMatrix::d_ilu_backward_rhs. Colors are true independent sets, so a color is race-free + * within a sweep; a U-neighbor in a different, not-yet-processed color this sweep is read + * at whatever value it currently holds (a previous sweep's, or the forward-solve result on + * the very first sweep, since prod is not reset before backward starts). Repeated full + * passes over all colors converge to the same result the exact level-scheduled solve would + * give (validated on the host: same fixed point, ~4.4x error reduction per sweep). Each + * thread accumulates its own (iVar,jVar) partial product across every neighbor with no + * synchronization, then one __syncthreads() to reduce over jVar, then a second + * __syncthreads() before the diagonal multiply (which needs every iVar's result). Grid: + * one block per row of the color, blockDim.x == nVar*nVar. Dynamic shared memory: + * nVar*nVar + nVar scalars. + */ +template +__global__ void IluBackwardColorKernel(const su2uint* __restrict__ color_idx, unsigned long color_begin, + unsigned long color_size, unsigned long nRows, unsigned long nVar, + DeviceLDU M, const ScalarType* __restrict__ rhs, + ScalarType* __restrict__ prod) { + if (blockIdx.x >= color_size) return; + + const unsigned long iRow = color_idx[color_begin + blockIdx.x]; + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + const auto iVar = tid / nVar, jVar = tid % nVar; + + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* partial = reinterpret_cast(smem); + auto* aux = partial + blockSize; + + ScalarType acc = 0; + for (auto ku = M.row_ptr_u[iRow]; ku < M.row_ptr_u[iRow + 1]; ++ku) { + const unsigned long jPoint = M.col_ind_u[ku]; + if (jPoint >= nRows) break; + const auto* blk = M.u + ku * blockSize; + acc += blk[iVar * nVar + jVar] * prod[jPoint * nVar + jVar]; + } + partial[tid] = acc; + __syncthreads(); + + if (jVar == 0) { + ScalarType sum = rhs[iRow * nVar + iVar]; + for (auto j = 0ul; j < nVar; ++j) sum -= partial[iVar * nVar + j]; + /*--- The diagonal blocks are stored inverted by the factorization. ---*/ + aux[iVar] = sum; + } + __syncthreads(); + + if (jVar == 0) { + const auto* invUii = M.d + iRow * blockSize; + ScalarType out = 0; + for (auto k = 0ul; k < nVar; ++k) out += invUii[iVar * nVar + k] * aux[k]; + prod[iRow * nVar + iVar] = out; + } +} + /*! * \brief Block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row. * One CUDA block per block-row; threadIdx.x indexes output variable (0..nVar-1). @@ -67,8 +367,211 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, y[iRow * nVar + 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); + /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); +} + +template +void CSysMatrix::BuildJacobiPreconditionerGPU() { + SU2_ZONE_SCOPED + + if (d_invM == nullptr) { + SU2_MPI::Error("CUDA Jacobi preconditioner used without device storage.", CURRENT_FUNCTION); + } + if (nVar != nEqn) { + SU2_MPI::Error("CUDA Jacobi preconditioner requires square blocks.", CURRENT_FUNCTION); + } + if (nVar * nVar > 1024) { + SU2_MPI::Error("CUDA Jacobi preconditioner uses one thread per block entry, nVar is too large.", CURRENT_FUNCTION); + } + if (nPointDomain == 0) return; + + /*--- The matrix is expected to be on the device already, it is uploaded once per solve by + * CSysMatrixVectorProduct, which is created before the preconditioner is built. ---*/ + const auto blockSize = static_cast(nVar * nVar); + InvertDiagonalBlocksKernel + <<(nPointDomain), blockSize, blockSize * sizeof(ScalarType)>>>(nPointDomain, nVar, gpu.d, + d_invM); + /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); +} + +template +void CSysMatrix::BuildILUPreconditionerGPU() { + SU2_ZONE_SCOPED + + if (gpu_ilu.d == nullptr) { + SU2_MPI::Error("CUDA ILU preconditioner used without device storage.", CURRENT_FUNCTION); + } + if (nVar != nEqn) { + SU2_MPI::Error("CUDA ILU factorization requires square blocks.", CURRENT_FUNCTION); + } + if (nVar * nVar > 1024) { + SU2_MPI::Error("CUDA ILU factorization uses one thread per block entry, nVar is too large.", CURRENT_FUNCTION); + } + if (nPointDomain == 0) return; + + /*--- The matrix is expected to be on the device already, it is uploaded once per solve by + * CSysMatrixVectorProduct, which is created before the preconditioner is built. ---*/ + const DeviceLDU A{gpu.d, gpu.l, gpu.u, gpu.row_ptr_l, + gpu.col_ind_l, gpu.row_ptr_u, gpu.col_ind_u}; + const DeviceLDU M{gpu_ilu.d, gpu_ilu.l, gpu_ilu.u, gpu_ilu.row_ptr_l, + gpu_ilu.col_ind_l, gpu_ilu.row_ptr_u, gpu_ilu.col_ind_u}; + + const auto blockSize = static_cast(nVar * nVar); + const auto shared = 2 * blockSize * sizeof(ScalarType); + + /*--- The legacy default stream cannot be captured, so the graph lives on its own stream, + * created once. Every launch below is followed by a sync back to the host, so this does not + * change execution order relative to the rest of the (single-stream) solver. ---*/ + if (ilu_stream == nullptr) gpuErrChk(cudaStreamCreate(&ilu_stream)); + + /*--- The launch sequence (ilu_gpu_build_sweeps passes over all colors) is identical on every + * call: the grid and block sizes only depend on the (fixed) sparsity pattern/coloring and the + * device pointers are fixed members, allocated once. Capture it into a CUDA graph the first + * time and replay that from then on, which removes the per-launch host-side overhead without + * touching the parallelization of any individual kernel (unlike a persistent cooperative- + * groups kernel, this does not cap per-color parallelism to an occupancy-resident block + * count). See IluFactorColorKernel for why several sweeps over the colors are needed. + * Note that factors are not reset between calls to BuildILUPreconditionerGPU, so with + * LINEAR_SOLVER_ILU_GPU_SWEEPS set low (even 1), each call refines the previous one's + * result rather than reconverging from scratch, relying on the matrix changing little + * between outer/pseudo-time iterations. ---*/ + if (ilu_build_graph_exec == nullptr) { + cudaGraph_t graph; + gpuErrChk(cudaStreamBeginCapture(ilu_stream, cudaStreamCaptureModeThreadLocal)); + + for (unsigned short sweep = 0; sweep < ilu_gpu_sweeps[0]; ++sweep) { + for (auto color = 0ul; color + 1 < ilu_color_ptr.size(); ++color) { + const auto begin = ilu_color_ptr[color]; + const auto size = ilu_color_ptr[color + 1] - begin; + if (size == 0) continue; + IluFactorColorKernel + <<>>(d_ilu_color_idx, begin, size, nPointDomain, nVar, A, M); + } + } + + gpuErrChk(cudaStreamEndCapture(ilu_stream, &graph)); + gpuErrChk(cudaGraphInstantiate(&ilu_build_graph_exec, graph, nullptr, nullptr, 0)); + gpuErrChk(cudaGraphDestroy(graph)); + } + + gpuErrChk(cudaGraphLaunch(ilu_build_graph_exec, ilu_stream)); + gpuErrChk(cudaStreamSynchronize(ilu_stream)); + gpuErrChk(cudaGetLastError()); +} + +template +void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector& vec, + CSysVector& prod) const { + SU2_ZONE_SCOPED + + if (gpu_ilu.d == nullptr) { + SU2_MPI::Error("CUDA ILU preconditioner used before BuildILUPreconditionerGPU.", CURRENT_FUNCTION); + } + if (nPointDomain == 0) return; + + const DeviceLDU M{gpu_ilu.d, gpu_ilu.l, gpu_ilu.u, gpu_ilu.row_ptr_l, + gpu_ilu.col_ind_l, gpu_ilu.row_ptr_u, gpu_ilu.col_ind_u}; + + auto* d_vec = vec.GetDevicePointer(); + auto* d_prod = prod.GetDevicePointer(); + + const auto nColors = ilu_color_ptr.size() - 1; + + /*--- One thread per block entry, like the factorization kernel: spreads each row's neighbor + * dot products over nVar*nVar threads instead of doing them serially in nVar threads, without + * changing the number of blocks (still one per row), so this does not trade away SM coverage + * the way batching several rows into a block did. ---*/ + const auto threads = static_cast(nVar * nVar); + const auto sharedForward = threads * sizeof(ScalarType); + const auto sharedBackward = (threads + nVar) * sizeof(ScalarType); + + if (ilu_stream == nullptr) gpuErrChk(cudaStreamCreate(&ilu_stream)); + + /*--- Same idea as BuildILUPreconditionerGPU: the launch sequence only depends on the (fixed) + * coloring, plus the vec/prod device pointers. Those normally are the same temporary buffers + * on every call (owned by CSysSolve / CSysVector, allocated once), so the graph is captured + * once and replayed; if the pointers ever do change the graph is recaptured, which is no worse + * than the un-graphed loop, just not free. ---*/ + if (ilu_apply_graph_exec == nullptr || ilu_apply_graph_vec != d_vec || ilu_apply_graph_prod != d_prod) { + if (ilu_apply_graph_exec != nullptr) { + gpuErrChk(cudaGraphExecDestroy(ilu_apply_graph_exec)); + ilu_apply_graph_exec = nullptr; + } + + cudaGraph_t graph; + gpuErrChk(cudaStreamBeginCapture(ilu_stream, cudaStreamCaptureModeThreadLocal)); + + /*--- Forward substitution: colored-iterative, see IluForwardColorKernel. Initialize prod to + * vec, not zero: (L+I)x = vec, so x = vec is the zeroth-order approximation that ignores the + * off-diagonal L contributions, which is generally closer to the converged answer than an + * all-zero guess (and no worse: on the very first call this buffer may hold unrelated + * leftover data, so it has to be set to something well-defined regardless). ---*/ + gpuErrChk(cudaMemcpyAsync(d_prod, d_vec, nPointDomain * nVar * sizeof(ScalarType), cudaMemcpyDeviceToDevice, + ilu_stream)); + for (unsigned short sweep = 0; sweep < ilu_gpu_sweeps[1]; ++sweep) { + for (auto color = 0ul; color < nColors; ++color) { + const auto begin = ilu_color_ptr[color]; + const auto size = ilu_color_ptr[color + 1] - begin; + if (size == 0) continue; + IluForwardColorKernel + <<>>(d_ilu_color_idx, begin, size, nVar, M, d_vec, d_prod); + } + } + + /*--- Backward substitution: also colored-iterative, see IluBackwardColorKernel for why the + * right-hand side must be copied out to a fixed buffer first (unlike forward, which already + * had one in d_vec). Starting prod as its own initial guess (the forward-solve result) is as + * good a starting point as zero and saves a second memset. ---*/ + gpuErrChk(cudaMemcpyAsync(d_ilu_backward_rhs, d_prod, nPointDomain * nVar * sizeof(ScalarType), + cudaMemcpyDeviceToDevice, ilu_stream)); + for (unsigned short sweep = 0; sweep < ilu_gpu_sweeps[2]; ++sweep) { + for (auto color = 0ul; color < nColors; ++color) { + const auto begin = ilu_color_ptr[color]; + const auto size = ilu_color_ptr[color + 1] - begin; + if (size == 0) continue; + IluBackwardColorKernel<<>>( + d_ilu_color_idx, begin, size, nPointDomain, nVar, M, d_ilu_backward_rhs, d_prod); + } + } + + gpuErrChk(cudaStreamEndCapture(ilu_stream, &graph)); + gpuErrChk(cudaGraphInstantiate(&ilu_apply_graph_exec, graph, nullptr, nullptr, 0)); + gpuErrChk(cudaGraphDestroy(graph)); + ilu_apply_graph_vec = d_vec; + ilu_apply_graph_prod = d_prod; + } + + gpuErrChk(cudaGraphLaunch(ilu_apply_graph_exec, ilu_stream)); + gpuErrChk(cudaStreamSynchronize(ilu_stream)); + gpuErrChk(cudaGetLastError()); +} + template void CSysMatrix::HtDTransfer(bool trigger) const { + SU2_ZONE_SCOPED if (!trigger) return; gpuErrChk(cudaMemcpy(gpu.d, mat.d, sizeof(ScalarType) * nPoint * nVar * nEqn, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpy(gpu.l, mat.l, sizeof(ScalarType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); @@ -76,8 +579,9 @@ void CSysMatrix::HtDTransfer(bool trigger) const { } template -void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, CSysVector& prod, +void CSysMatrix::MatrixVectorProductGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const { + SU2_ZONE_SCOPED if (nVar != nEqn) { SU2_MPI::Error("CUDA CSysMatrix block-LDU SpMV requires square blocks.", CURRENT_FUNCTION); } @@ -90,16 +594,27 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector BlockLDU_SpMV_kernel<<>>( 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); + /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); gpuErrChk(cudaGetLastError()); } -template void CSysMatrix::HtDTransfer(bool trigger) const; -template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, - CSysVector& prod, CGeometry* geometry, - const CConfig* config) const; + +#define INSTANTIATE_MATRIX(TYPE) \ +template void CSysMatrix::HtDTransfer(bool trigger) const; \ +template void CSysMatrix::MatrixVectorProductGPU(const CSysVector& vec, \ + CSysVector& prod, \ + CGeometry* geometry, \ + const CConfig* config) const; \ +template void CSysMatrix::BuildJacobiPreconditionerGPU(); \ +template void CSysMatrix::BuildILUPreconditionerGPU(); \ +template void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector& vec, \ + CSysVector& prod) const; \ +template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, \ + CSysVector& prod, \ + CGeometry* geometry, \ + const CConfig* config) const; +INSTANTIATE_MATRIX(su2mixedfloat) #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; +INSTANTIATE_MATRIX(passivedouble) #endif diff --git a/Common/src/linear_algebra/CSysPreconditionerGPU.cu b/Common/src/linear_algebra/CSysPreconditionerGPU.cu deleted file mode 100644 index 794726fbce35..000000000000 --- a/Common/src/linear_algebra/CSysPreconditionerGPU.cu +++ /dev/null @@ -1,84 +0,0 @@ -/*! - * \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/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 2be1215a7bd0..01b146baf762 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public @@ -66,18 +66,21 @@ void SetUseDeviceExpressions(bool use) { use_device_expressions = use; } template void CSysVector::HtDTransfer(bool trigger) const { + SU2_ZONE_SCOPED if (trigger) gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType) * nElm), cudaMemcpyHostToDevice)); } template void CSysVector::DtHTransfer(bool trigger) const { + SU2_ZONE_SCOPED if (trigger) gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType) * nElm), cudaMemcpyDeviceToHost)); } template ScalarType CSysVector::GPUDot(const CSysVector& other) const { + SU2_ZONE_SCOPED /*--- 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(); @@ -110,6 +113,7 @@ ScalarType CSysVector::GPUDot(const CSysVector& other) const { template ScalarType CSysVector::GPUNorm() const { + SU2_ZONE_SCOPED return sqrt(GPUDot(*this)); } diff --git a/Common/src/linear_algebra/meson.build b/Common/src/linear_algebra/meson.build index 48ef65cb8db2..3b84b2373a55 100644 --- a/Common/src/linear_algebra/meson.build +++ b/Common/src/linear_algebra/meson.build @@ -8,5 +8,5 @@ common_src += files(['CSysSolve_b.cpp', 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']) + common_cuda_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu']) endif diff --git a/Common/src/toolboxes/SwapBytes.cpp b/Common/src/toolboxes/SwapBytes.cpp index 6f4d43504570..d60567e09ed0 100644 --- a/Common/src/toolboxes/SwapBytes.cpp +++ b/Common/src/toolboxes/SwapBytes.cpp @@ -9,7 +9,7 @@ * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) * * SU2 is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 443d367a10a8..e2c375f831a8 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -59,7 +59,7 @@ def main(): naca0012.cfg_dir = "euler/naca0012" naca0012.cfg_file = "inv_NACA0012_Roe.cfg" naca0012.test_iter = 20 - naca0012.test_vals = [-4.387423, -3.845336, 0.296291, 0.025144] + naca0012.test_vals = [-4.387446, -3.845358, 0.296291, 0.025143] test_list.append(naca0012) # Supersonic wedge @@ -157,7 +157,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.190528, -5.335496, 0.383385, 0.077606, 0.000000] + rae2822_sa.test_vals = [-2.190528, -5.335496, 0.383385, 0.077605, 0.000000] test_list.append(rae2822_sa) # RAE2822 SST @@ -197,7 +197,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.038042, -16.332088, 1.080346, 0.018385, 20.000000, -2.873258, 0.000000, -14.250270, 0.000000] + turb_naca0012_sa.test_vals = [-12.038067, -16.332088, 1.080346, 0.018385, 20.000000, -2.873679, 0.000000, -14.250269, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.038091, -16.332090, 1.080346, 0.018385, 20.000000, -2.873236, 0.000000, -14.250271, 0.000000] test_list.append(turb_naca0012_sa) @@ -206,7 +206,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.093984, -15.250705, -5.906323, 1.070413, 0.015775, -2.855331, 0.000000] + turb_naca0012_sst.test_vals = [-12.093958, -15.250719, -5.906323, 1.070413, 0.015775, -2.855776, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.075928, -15.246732, -5.861249, 1.070036, 0.015841, -2.835263, 0] test_list.append(turb_naca0012_sst) @@ -215,7 +215,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.080774, -14.837176, -5.732907, 1.000893, 0.019109, -2.120168] + turb_naca0012_sst_sust.test_vals = [-12.080818, -14.837175, -5.732906, 1.000893, 0.019109, -2.119717] turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073210, -14.836724, -5.732627, 1.000050, 0.019144, -2.629689] test_list.append(turb_naca0012_sst_sust) @@ -224,7 +224,7 @@ def main(): turb_naca0012_sst_fixedvalues.cfg_dir = "rans/naca0012" turb_naca0012_sst_fixedvalues.cfg_file = "turb_NACA0012_sst_fixedvalues.cfg" turb_naca0012_sst_fixedvalues.test_iter = 10 - turb_naca0012_sst_fixedvalues.test_vals = [-5.192390, -10.448219, 0.773965, 1.022535, 0.040529, -2.383282] + turb_naca0012_sst_fixedvalues.test_vals = [-5.192390, -10.448218, 0.773965, 1.022535, 0.040529, -2.383435] test_list.append(turb_naca0012_sst_fixedvalues) # NACA0012 (SST, explicit Euler for flow and turbulence equations) @@ -252,7 +252,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.083070, -5.374688, -8.880089, -4.073519, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-11.083070, -5.374685, -8.880088, -4.073510, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000] test_list.append(axi_rans_air_nozzle_restart) @@ -420,7 +420,7 @@ def main(): inc_weakly_coupled.cfg_dir = "disc_adj_heat" inc_weakly_coupled.cfg_file = "primal.cfg" inc_weakly_coupled.test_iter = 10 - inc_weakly_coupled.test_vals = [-18.106234, -16.302995, -16.484896, -15.006575, -17.858050, -14.024855, 5.609100] + inc_weakly_coupled.test_vals = [-18.106270, -16.302974, -16.484882, -15.006575, -17.858050, -14.024836, 5.609100] test_list.append(inc_weakly_coupled) ###################################### @@ -501,7 +501,7 @@ def main(): ddes_flatplate.cfg_dir = "ddes/flatplate" ddes_flatplate.cfg_file = "ddes_flatplate.cfg" ddes_flatplate.test_iter = 10 - ddes_flatplate.test_vals = [-2.714713, -5.763299, -0.214960, 0.023758, 0.000000] + ddes_flatplate.test_vals = [-2.714713, -5.763300, -0.214960, 0.023758, 0.000000] ddes_flatplate.unsteady = True test_list.append(ddes_flatplate) @@ -686,7 +686,7 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-2.830455, -1.746659, -2.366087, 110350] + statbeam3d.test_vals = [-2.862081, -1.795143, -2.380136, 110350.000000] statbeam3d.test_vals_aarch64 = [-2.777602, -1.710976, -2.445072, 110350] test_list.append(statbeam3d) @@ -715,7 +715,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.330727, -4.152808, 0.000000, 103.000000] + dyn_fsi.test_vals = [-4.330727, -4.152808, 0.000000, 102.000000] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) @@ -725,7 +725,7 @@ def main(): fsi_cht_restart.cfg_dir = "fea_fsi/stat_fsi" fsi_cht_restart.cfg_file = "config_restart.cfg" fsi_cht_restart.test_iter = 0 - fsi_cht_restart.test_vals = [5.000000, 0.006352, -1.960362, -9.327033, -9.627867, -9.318971, 608.380000, -0.012974, 0.000000, 20.000000] + fsi_cht_restart.test_vals = [5.000000, 0.006352, -1.960362, -9.327033, -9.643180, -9.319159, 608.380000, -0.012974, 0.000000, 20.000000] fsi_cht_restart.multizone = True test_list.append(fsi_cht_restart) @@ -755,7 +755,7 @@ def main(): mms_fvm_inc_ns.cfg_dir = "mms/fvm_incomp_navierstokes" mms_fvm_inc_ns.cfg_file = "lam_mms_fds.cfg" mms_fvm_inc_ns.test_iter = 20 - mms_fvm_inc_ns.test_vals = [-7.414945, -7.631547, 0.000000, 0.000000] + mms_fvm_inc_ns.test_vals = [-7.414944, -7.631546, 0.000000, 0.000000] test_list.append(mms_fvm_inc_ns) ########################## diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py index 4db360930312..9440e408ab15 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -78,7 +78,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.987158, 0.533077, 0.000004, -0.000000, 5.000000, -2.939652, 5.000000, -5.502411] + discadj_rans_naca0012_sa.test_vals = [-2.987158, 0.533077, 0.000004, -0.000000, 5.000000, -2.939652, 5.000000, -5.503137] test_list.append(discadj_rans_naca0012_sa) # Adjoint turbulent NACA0012 SST @@ -86,7 +86,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.201555, -0.175211, 3.045200, -0.041846] + discadj_rans_naca0012_sst.test_vals = [-2.201551, -0.175213, 3.045300, -0.041846] discadj_rans_naca0012_sst.test_vals_aarch64 = [-2.201855, -0.172443, 3.043400, -0.041820] test_list.append(discadj_rans_naca0012_sst) @@ -111,7 +111,7 @@ def main(): discadj_incomp_cylinder.cfg_dir = "disc_adj_incomp_navierstokes/cylinder" discadj_incomp_cylinder.cfg_file = "heated_cylinder.cfg" discadj_incomp_cylinder.test_iter = 20 - discadj_incomp_cylinder.test_vals = [20.000000, -1.665652, -6.239091, 0.000000] + discadj_incomp_cylinder.test_vals = [20.000000, -1.665652, -6.239103, 0.000000] discadj_incomp_cylinder.test_vals_aarch64 = [20.000000, -1.671920, -6.254841, 0.000000] discadj_incomp_cylinder.tol_aarch64 = 2e-1 test_list.append(discadj_incomp_cylinder) @@ -125,7 +125,7 @@ def main(): discadj_incomp_turb_NACA0012_sa.cfg_dir = "disc_adj_incomp_rans/naca0012" discadj_incomp_turb_NACA0012_sa.cfg_file = "turb_naca0012_sa.cfg" discadj_incomp_turb_NACA0012_sa.test_iter = 10 - discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.845989, -1.023525, 0.000000] + discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.845989, -1.023526, 0.000000] test_list.append(discadj_incomp_turb_NACA0012_sa) # Adjoint Incompressible Turbulent NACA 0012 SST @@ -133,7 +133,7 @@ def main(): discadj_incomp_turb_NACA0012_sst.cfg_dir = "disc_adj_incomp_rans/naca0012" discadj_incomp_turb_NACA0012_sst.cfg_file = "turb_naca0012_sst.cfg" discadj_incomp_turb_NACA0012_sst.test_iter = 10 - discadj_incomp_turb_NACA0012_sst.test_vals = [-3.775275, -3.089116, -7.143654, 0.000000, -0.896764] + discadj_incomp_turb_NACA0012_sst.test_vals = [-3.775320, -3.089107, -7.143663, 0.000000, -0.896754] test_list.append(discadj_incomp_turb_NACA0012_sst) ####################################################### @@ -159,7 +159,7 @@ def main(): discadj_cylinder.cfg_dir = "disc_adj_rans/cylinder" discadj_cylinder.cfg_file = "cylinder_Windowing_AD.cfg" discadj_cylinder.test_iter = 9 - discadj_cylinder.test_vals = [2.183381] + discadj_cylinder.test_vals = [2.183380] discadj_cylinder.unsteady = True discadj_cylinder.enabled_with_tsan = False test_list.append(discadj_cylinder) @@ -173,7 +173,7 @@ def main(): discadj_DT_1ST_cylinder.cfg_dir = "disc_adj_rans/cylinder_DT_1ST" discadj_DT_1ST_cylinder.cfg_file = "cylinder.cfg" discadj_DT_1ST_cylinder.test_iter = 9 - discadj_DT_1ST_cylinder.test_vals = [1.196350, -3.339010, -0.006213, 0.000020] + discadj_DT_1ST_cylinder.test_vals = [1.196351, -3.339009, -0.006213, 0.000020] discadj_DT_1ST_cylinder.unsteady = True discadj_DT_1ST_cylinder.enabled_with_tsan = False test_list.append(discadj_DT_1ST_cylinder) @@ -202,7 +202,7 @@ def main(): discadj_fea.cfg_dir = "disc_adj_fea" discadj_fea.cfg_file = "configAD_fem.cfg" discadj_fea.test_iter = 4 - discadj_fea.test_vals = [2.149620, 2.014985, -0.000364, -8.767900] + discadj_fea.test_vals = [2.180053, 2.075836, -0.000367, -8.730500] discadj_fea.test_vals_aarch64 = [1.794371, 2.005865, -0.000365, -8.718100] test_list.append(discadj_fea) @@ -228,7 +228,7 @@ def main(): pywrapper_FEA_AD_FlowLoad.cfg_dir = "py_wrapper/disc_adj_fea/flow_load_sens" pywrapper_FEA_AD_FlowLoad.cfg_file = "configAD_fem.cfg" pywrapper_FEA_AD_FlowLoad.test_iter = 100 - pywrapper_FEA_AD_FlowLoad.test_vals = [-0.132010, -0.554418, -0.000364, -0.003101] + pywrapper_FEA_AD_FlowLoad.test_vals = [-0.131741, -0.553305, -0.000364, -0.003101] pywrapper_FEA_AD_FlowLoad.test_vals_aarch64 = [-0.131745, -0.553214, -0.000364, -0.003101] pywrapper_FEA_AD_FlowLoad.command = TestCase.Command(exec = "python", param = "run_adjoint.py --parallel -f") pywrapper_FEA_AD_FlowLoad.timeout = 1600 @@ -243,7 +243,7 @@ def main(): pywrapper_CFD_AD_MeshDisp.cfg_dir = "py_wrapper/disc_adj_flow/mesh_disp_sens" pywrapper_CFD_AD_MeshDisp.cfg_file = "configAD_flow.cfg" pywrapper_CFD_AD_MeshDisp.test_iter = 1000 - pywrapper_CFD_AD_MeshDisp.test_vals = [30.000000, -2.496158, 1.441842, 0.000000] + pywrapper_CFD_AD_MeshDisp.test_vals = [30.000000, -2.496139, 1.441872, 0.000000] pywrapper_CFD_AD_MeshDisp.test_vals_aarch64 = [30.000000, -2.499079, 1.440068, 0.000000] pywrapper_CFD_AD_MeshDisp.command = TestCase.Command(exec = "python", param = "run_adjoint.py --parallel -f") pywrapper_CFD_AD_MeshDisp.timeout = 1600 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 9c08ad013e5e..e6242b093b1c 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -138,7 +138,7 @@ def main(): invwedge_msw.cfg_dir = "nonequilibrium/invwedge" invwedge_msw.cfg_file = "invwedge_msw.cfg" invwedge_msw.test_iter = 10 - invwedge_msw.test_vals = [-1.212335, -1.737098, -18.301768, -18.629149, -18.575169, 2.106171, 1.651949, 5.143958, 0.704444] + invwedge_msw.test_vals = [-1.212335, -1.737098, -18.301825, -18.629206, -18.575226, 2.106171, 1.651949, 5.143958, 0.704444] invwedge_msw.test_vals_aarch64 = [-1.212335, -1.737098, -18.299279, -18.626656, -18.572683, 2.106171, 1.651949, 5.143958, 0.704444] test_list.append(invwedge_msw) @@ -147,7 +147,7 @@ def main(): invwedge_roe.cfg_dir = "nonequilibrium/invwedge" invwedge_roe.cfg_file = "invwedge_roe.cfg" invwedge_roe.test_iter = 10 - invwedge_roe.test_vals = [-1.023216, -1.547979, -17.814656, -18.143616, -18.087775, 2.295078, 1.885054, 5.338487, 0.926120] + invwedge_roe.test_vals = [-1.023283, -1.548046, -17.814403, -18.143369, -18.087522, 2.295025, 1.884804, 5.338440, 0.926068] invwedge_roe.test_vals_aarch64 = [-1.052398, -1.577160, -17.794015, -18.122997, -18.067131, 2.266042, 1.849686, 5.304700, 0.899584] test_list.append(invwedge_roe) @@ -183,7 +183,7 @@ def main(): visc_cone.cfg_dir = "nonequilibrium/visc_wedge" visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 - visc_cone.test_vals = [-5.215235, -5.739371, -20.559852, -20.509281, -20.408911, 1.262701, -3.205457, -0.015696, 0.093205, 32637.000000] + visc_cone.test_vals = [-5.215234, -5.739371, -20.559852, -20.509281, -20.408911, 1.262701, -3.205457, -0.015696, 0.093205, 32637.000000] visc_cone.test_vals_aarch64 = [-5.222270, -5.746525, -20.560286, -20.510152, -20.409101, 1.255758, -3.208382, -0.016014, 0.093462, 32619.000000] test_list.append(visc_cone) @@ -326,7 +326,7 @@ def main(): ramp_msw.cfg_dir = "euler/ramp" ramp_msw.cfg_file = "inv_ramp_msw.cfg" ramp_msw.test_iter = 100 - ramp_msw.test_vals = [-7.059306, -1.300966, -0.077507, 0.054419] + ramp_msw.test_vals = [-7.219257, -1.444776, -0.077507, 0.054419] ramp_msw.tol = [0.2, 0.2, 0.00001, 0.00001] test_list.append(ramp_msw) @@ -820,7 +820,7 @@ def main(): turbmod_sa_neg_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_neg_rae2822.cfg_file = "turb_SA_NEG_RAE2822.cfg" turbmod_sa_neg_rae2822.test_iter = 10 - turbmod_sa_neg_rae2822.test_vals = [1.527546, 1.303378, -1.699437, 1.523548, 0.601919, 0.000000] + turbmod_sa_neg_rae2822.test_vals = [1.611007, 1.331579, -1.264006, 1.360053, 0.508605, 0.000000] turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.345593, 1.448310, 1.208721, -0.846597, 1.248410, 0.489117, 0.000000] test_list.append(turbmod_sa_neg_rae2822) @@ -1324,7 +1324,7 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-6.047457, -5.730708, -5.926093, 110190] + statbeam3d.test_vals = [-6.020105, -5.749794, -5.931330, 110190.000000] statbeam3d.test_vals_aarch64 = [-6.062693, -5.769132, -5.891190, 110190] test_list.append(statbeam3d) @@ -1333,7 +1333,7 @@ def main(): thermal_beam_3d.cfg_dir = "fea_fsi/ThermalBeam_3d" thermal_beam_3d.cfg_file = "configBeam_3d.cfg" thermal_beam_3d.test_iter = 4 - thermal_beam_3d.test_vals = [-8.070340, -7.802437, -7.856284, -13.978110, 217.000000, -4.047750, 39.000000, -4.072613, 136760.000000, 75.000000] + thermal_beam_3d.test_vals = [-8.141292, -7.848655, -7.909552, -13.978110, 217.000000, -4.099949, 39.000000, -4.072613, 136760.000000, 75.000000] test_list.append(thermal_beam_3d) # Static beam, 3d with coupled temperature, nonlinear elasticity @@ -1341,7 +1341,7 @@ def main(): thermal_beam_nl_3d.cfg_dir = "fea_fsi/ThermalBeam_3d" thermal_beam_nl_3d.cfg_file = "configBeamNonlinear_3d.cfg" thermal_beam_nl_3d.test_iter = 8 - thermal_beam_nl_3d.test_vals = [-7.564308, -2.992893, -12.242503, -14.068322, 57.000000, -4.017672, 24.000000, -4.204804, 138710.000000, 75.233000] + thermal_beam_nl_3d.test_vals = [-7.564308, -2.992893, -12.242503, -14.068322, 57.000000, -4.017665, 24.000000, -4.204804, 138710.000000, 75.233000] test_list.append(thermal_beam_nl_3d) # Rotating cylinder, 3d @@ -1352,7 +1352,7 @@ def main(): # For a thin disk with the inner and outer radius of this geometry, from # "Formulas for Stress, Strain, and Structural Matrices", 2nd Edition, figure 19-4, # the maximum stress is 165.6MPa, we get a von Mises stress very close to that. - rotating_cylinder_fea.test_vals = [-6.760497, -6.689264, -6.739355, 37.000000, -8.178510, 165020000.000000] + rotating_cylinder_fea.test_vals = [-6.760497, -6.689266, -6.739355, 37.000000, -8.178510, 165020000.000000] rotating_cylinder_fea.test_vals_aarch64 = [-6.861939, -6.835539, -6.895498, 22, -8.313847, 1.6502e+08] test_list.append(rotating_cylinder_fea) @@ -1361,7 +1361,7 @@ def main(): linear_plane_strain.cfg_dir = "fea_fsi/VonMissesVerif" linear_plane_strain.cfg_file = "linear_plane_strain_2d.cfg" linear_plane_strain.test_iter = 0 - linear_plane_strain.test_vals = [-6.406458, -5.995503, 0, 120140, 144, -8.122248] + linear_plane_strain.test_vals = [-6.246772, -6.060575, 0.000000, 120140.000000, 148.000000, -8.140997] test_list.append(linear_plane_strain) # 2D beam in plain stress with thermal expansion. This tests fixes to the 2D von Mises stress calculation, @@ -1370,7 +1370,7 @@ def main(): nonlinear_plane_stress.cfg_dir = "fea_fsi/VonMissesVerif" nonlinear_plane_stress.cfg_file = "nonlinear_plane_stress_2d.cfg" nonlinear_plane_stress.test_iter = 16 - nonlinear_plane_stress.test_vals = [-6.230217, -2.228380, -11.698370, 162480.000000, 30.000000, -4.169761] + nonlinear_plane_stress.test_vals = [-6.235980, -2.229690, -11.701484, 162480.000000, 30.000000, -4.152654] nonlinear_plane_stress.tol = [2e-4, 2e-4, 2e-4, 1e-5, 1e-5, 4e-4] test_list.append(nonlinear_plane_stress) @@ -1520,7 +1520,7 @@ def main(): pywrapper_custom_fea_load.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_custom_fea_load.cfg_file = "config.cfg" pywrapper_custom_fea_load.test_iter = 13 - pywrapper_custom_fea_load.test_vals = [-7.262040, -4.945686, -14.163208, 27.000000, -6.282188, 362.230000] + pywrapper_custom_fea_load.test_vals = [-7.262040, -4.945686, -14.163208, 27.000000, -6.285429, 362.230000] pywrapper_custom_fea_load.command = TestCase.Command("mpirun -np 2", "python", "run.py") test_list.append(pywrapper_custom_fea_load) @@ -1633,7 +1633,7 @@ def main(): mms_fvm_inc_euler.cfg_dir = "mms/fvm_incomp_euler" mms_fvm_inc_euler.cfg_file = "inv_mms_jst.cfg" mms_fvm_inc_euler.test_iter = 20 - mms_fvm_inc_euler.test_vals = [-9.128735, -9.441756, 0.000000, 0.000000] + mms_fvm_inc_euler.test_vals = [-9.128735, -9.441757, 0.000000, 0.000000] mms_fvm_inc_euler.tol = 0.0001 test_list.append(mms_fvm_inc_euler) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 8e66a9716b10..e68c3f27648e 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -231,7 +231,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, 2.549036, 2.313067, 2.139713, 0.736742] + discadj_trans_stator.test_vals = [79.000000, 2.549037, 2.313067, 2.139716, 0.736741] discadj_trans_stator.test_vals_aarch64 = [79.000000, 0.696755, 0.485950, 0.569475, -0.990065] test_list.append(discadj_trans_stator) @@ -244,7 +244,7 @@ def main(): discadj_fea.cfg_dir = "disc_adj_fea" discadj_fea.cfg_file = "configAD_fem.cfg" discadj_fea.test_iter = 4 - discadj_fea.test_vals = [-2.849947, -3.238801, -0.000364, -8.708700] + discadj_fea.test_vals = [-2.849664, -3.238571, -0.000364, -8.708700] discadj_fea.test_vals_aarch64 = [-2.849646, -3.238577, -0.000364, -8.708700] #last 4 columns test_list.append(discadj_fea) @@ -286,7 +286,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-3.824634, 1.979533, -3.863368, 0.295450, 3.839800] + discadj_fsi2.test_vals = [-3.824641, 1.979547, -3.863368, 0.295450, 3.839800] discadj_fsi2.test_vals_aarch64 = [-3.824870, 1.979160, -3.863368, 0.295450, 3.839800] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2) @@ -308,7 +308,7 @@ def main(): da_sp_pinArray_cht_2d_dp_hf.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" da_sp_pinArray_cht_2d_dp_hf.cfg_file = "DA_configMaster.cfg" da_sp_pinArray_cht_2d_dp_hf.test_iter = 100 - da_sp_pinArray_cht_2d_dp_hf.test_vals = [-11.620815, -6.488915, -13.316507] + da_sp_pinArray_cht_2d_dp_hf.test_vals = [-11.620815, -6.488915, -13.316675] da_sp_pinArray_cht_2d_dp_hf.multizone = True test_list.append(da_sp_pinArray_cht_2d_dp_hf) @@ -326,7 +326,7 @@ def main(): da_unsteadyCHT_cylinder.cfg_dir = "coupled_cht/disc_adj_unsteadyCHT_cylinder" da_unsteadyCHT_cylinder.cfg_file = "chtMaster.cfg" da_unsteadyCHT_cylinder.test_iter = 2 - da_unsteadyCHT_cylinder.test_vals = [-8.479629, -9.239920, -9.234868, -15.934511, -13.662007, 0.000000, 10.627000, 0.295190] + da_unsteadyCHT_cylinder.test_vals = [-8.479629, -9.239920, -9.234868, -15.934511, -13.662005, 0.000000, 10.627000, 0.295190] da_unsteadyCHT_cylinder.test_vals_aarch64 = [-8.479629, -9.239920, -9.234868, -15.934511, -13.662012, 0.000000, 89.932000, 0.295190] da_unsteadyCHT_cylinder.unsteady = True da_unsteadyCHT_cylinder.multizone = True @@ -521,7 +521,7 @@ def main(): pywrapper_wavy_wall_steady.cfg_dir = "py_wrapper/wavy_wall" pywrapper_wavy_wall_steady.cfg_file = "run_steady.py" pywrapper_wavy_wall_steady.test_iter = 100 - pywrapper_wavy_wall_steady.test_vals = [-1.353007, 2.581052, -2.900574] + pywrapper_wavy_wall_steady.test_vals = [-1.353007, 2.581051, -2.900574] pywrapper_wavy_wall_steady.command = TestCase.Command("mpirun -n 2", "python", "run_steady.py") pywrapper_wavy_wall_steady.timeout = 1600 pywrapper_wavy_wall_steady.tol = 0.00001 diff --git a/TestCases/py_wrapper/custom_source_buoyancy/lam_buoyancy_cavity.cfg b/TestCases/py_wrapper/custom_source_buoyancy/lam_buoyancy_cavity.cfg index 6329eaf739c5..6f33ec59f608 100644 --- a/TestCases/py_wrapper/custom_source_buoyancy/lam_buoyancy_cavity.cfg +++ b/TestCases/py_wrapper/custom_source_buoyancy/lam_buoyancy_cavity.cfg @@ -4,7 +4,7 @@ % Case description: Buoyancy-driven flow inside a cavity % % Author: Thomas D. Economon % % Date: 2018.06.10 % -% File Version 8.1.0 "Harrier" % +% File Version 8.5.0 "Harrier" % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/TestCases/py_wrapper/custom_source_buoyancy/run.py b/TestCases/py_wrapper/custom_source_buoyancy/run.py index f07dfdc3570c..190c2db49133 100644 --- a/TestCases/py_wrapper/custom_source_buoyancy/run.py +++ b/TestCases/py_wrapper/custom_source_buoyancy/run.py @@ -2,14 +2,14 @@ ## \file run.py # \brief Buoyancy force using user defines source term -# \version 8.1.0 "Harrier" +# \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-2024, SU2 Contributors (cf. AUTHORS.md) +# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) # # SU2 is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/run.py b/TestCases/py_wrapper/turbulent_premixed_psi/run.py index 8a157ef0ae97..27021528fdf6 100644 --- a/TestCases/py_wrapper/turbulent_premixed_psi/run.py +++ b/TestCases/py_wrapper/turbulent_premixed_psi/run.py @@ -3,14 +3,14 @@ ## \file run.py # \brief turbulent premixed dump combustor simulation (PSI flame) # phi=0.5, methane-air, U=40 m/s -# \version 8.1.0 "Harrier" +# \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-2024, SU2 Contributors (cf. AUTHORS.md) +# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) # # SU2 is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public @@ -95,7 +95,7 @@ def update_temperature(SU2Driver, iPoint): iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW'] iENTH = 3 - #h = + #h = SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iENTH, cp_u*(T-Tref)) diff --git a/TestCases/py_wrapper/updated_moving_frame_NACA12/forces_0.csv.ref b/TestCases/py_wrapper/updated_moving_frame_NACA12/forces_0.csv.ref index 0a0e8b9b8a41..0f39252ba9ae 100644 --- a/TestCases/py_wrapper/updated_moving_frame_NACA12/forces_0.csv.ref +++ b/TestCases/py_wrapper/updated_moving_frame_NACA12/forces_0.csv.ref @@ -18,7 +18,7 @@ 16, -7.83, 58.91, 0.00 17, -7.39, 56.13, 0.00 18, -6.87, 52.74, 0.00 -19, -6.30, 48.94, 0.00 +19, -6.30, 48.95, 0.00 20, -5.67, 44.59, 0.00 21, -5.01, 39.84, 0.00 22, -4.30, 34.65, 0.00 @@ -52,7 +52,7 @@ 50, 8.45, -142.76, 0.00 51, 6.40, -115.54, 0.00 52, 9.86, -191.93, 0.00 -53, 24.75, -524.72, 0.00 +53, 24.75, -524.71, 0.00 54, 26.02, -608.43, 0.00 55, 21.64, -567.32, 0.00 56, 18.59, -558.24, 0.00 @@ -83,35 +83,35 @@ 81, -34.19, -153.83, 0.00 82, -32.26, -131.97, 0.00 83, -30.00, -111.54, 0.00 -84, -26.56, -89.60, 0.00 +84, -26.57, -89.60, 0.00 85, -22.26, -67.94, 0.00 86, -16.65, -45.84, 0.00 -87, -10.42, -25.77, 0.00 +87, -10.42, -25.78, 0.00 88, -1.99, -4.39, 0.00 -89, 7.81, 15.27, 0.00 -90, 17.05, 29.26, 0.00 +89, 7.81, 15.28, 0.00 +90, 17.06, 29.26, 0.00 91, 27.28, 40.55, 0.00 -92, 36.43, 46.13, 0.00 +92, 36.43, 46.14, 0.00 93, 50.49, 53.24, 0.00 94, 63.17, 53.66, 0.00 -95, 73.70, 48.01, 0.00 -96, 73.03, 33.47, 0.00 -97, 67.54, 18.28, 0.00 +95, 73.70, 48.00, 0.00 +96, 73.02, 33.47, 0.00 +97, 67.53, 18.28, 0.00 98, 54.75, 7.35, 0.00 -99, 17.21, -0.00, 0.00 -100, 52.61, -7.06, 0.00 -101, 93.60, -25.33, 0.00 -102, 62.71, -28.74, 0.00 +99, 17.20, -0.00, 0.00 +100, 52.60, -7.06, 0.00 +101, 93.59, -25.33, 0.00 +102, 62.72, -28.75, 0.00 103, 27.14, -17.68, 0.00 -104, 23.12, -19.64, 0.00 +104, 23.11, -19.63, 0.00 105, 5.43, -5.72, 0.00 -106, -0.11, 0.15, 0.00 -107, -12.46, 18.52, 0.00 +106, -0.12, 0.15, 0.00 +107, -12.46, 18.53, 0.00 108, -19.25, 33.03, 0.00 109, -25.52, 49.90, 0.00 110, -31.57, 69.68, 0.00 111, -35.49, 87.76, 0.00 -112, -39.31, 108.24, 0.00 +112, -39.31, 108.23, 0.00 113, -41.08, 125.40, 0.00 114, -43.50, 146.70, 0.00 115, -44.37, 164.95, 0.00 @@ -153,7 +153,7 @@ 151, 50.35, 723.19, 0.00 152, 54.49, 748.60, 0.00 153, 50.92, 671.51, 0.00 -154, 6.40, 81.24, 0.00 +154, 6.40, 81.25, 0.00 155, -0.48, -5.87, 0.00 156, 2.99, 35.56, 0.00 157, 3.29, 37.89, 0.00 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 8adb14300f37..7a5a57026a24 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1069,7 +1069,7 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-6.175086, -5.939313, -6.084188, 110190] + statbeam3d.test_vals = [-6.192310, -5.950395, -6.079363, 110190.000000] statbeam3d.test_vals_aarch64 = [-6.166287, -5.938291, -6.069768, 110190] #last 4 columns test_list.append(statbeam3d) @@ -1107,7 +1107,7 @@ def main(): fsi_cht.cfg_dir = "fea_fsi/stat_fsi" fsi_cht.cfg_file = "config.cfg" fsi_cht.test_iter = 20 - fsi_cht.test_vals = [5.000000, -5.077003, -5.379449, -9.247804, -9.319626, -9.184904, 608.350000, -0.012973, 0.000000, 30.000000] + fsi_cht.test_vals = [5.000000, -5.076991, -5.379442, -9.247793, -9.319193, -9.184753, 608.350000, -0.012973, 0.000000, 30.000000] fsi_cht.multizone = True test_list.append(fsi_cht) diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 15598919e70e..2c51de9f9dd3 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -180,7 +180,7 @@ def main(): discadj_fea.cfg_dir = "disc_adj_fea" discadj_fea.cfg_file = "configAD_fem.cfg" discadj_fea.test_iter = 4 - discadj_fea.test_vals = [-2.849719, -3.238637, -0.000364, -8.708700] + discadj_fea.test_vals = [-2.849715, -3.238627, -0.000364, -8.708700] discadj_fea.test_vals_aarch64 = [-2.849588, -3.238523, -0.000364, -8.708700] discadj_fea.tol = 0.00007 test_list.append(discadj_fea) diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 4e039eea09b9..aefbc66f815d 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -302,7 +302,7 @@ def main(): tutorial_unst_naca0012.cfg_dir = "../Tutorials/compressible_flow/Unsteady_NACA0012" tutorial_unst_naca0012.cfg_file = "unsteady_naca0012.cfg" tutorial_unst_naca0012.test_iter = 520 - tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.301440, 0.000000, 0.313631, 0.803638, 0.002198, 0.014969] + tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.293170, 0.000000, 0.301553, 0.773822, 0.001267, 0.007555] tutorial_unst_naca0012.test_vals_aarch64 = [520, 0, -5.292359, 0, 0.284720, 0.766329, 0.000954, 0.007565] tutorial_unst_naca0012.unsteady = True test_list.append(tutorial_unst_naca0012) diff --git a/config_template.cfg b/config_template.cfg index bc6be98103ae..b12ed3ccd0f8 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1661,6 +1661,15 @@ DISCADJ_LIN_PREC= ILU % Linear solver ILU preconditioner fill-in level (0 by default) LINEAR_SOLVER_ILU_FILL_IN= 0 % +% Colored-iterative sweep counts for the GPU ILU preconditioner: (build, forward, backward), +% default (1, 2, 2). The build (factorization) is not reset between calls, so with the matrix +% changing little between outer/pseudo-time iterations, each call refines the previous one's +% result instead of reconverging from scratch, for cases that do few outer iterations (e.g. +% elasticity problems) it may be useful to increase it to 3-5. The forward/backward triangular +% solves have no such warm start (a new right-hand side every Krylov iteration), so more than +% one sweep is necessary, harder problems (e.g. higher CFL) may benefit from more than 3 sweeps. +LINEAR_SOLVER_ILU_GPU_SWEEPS= (1, 2, 2) +% % Minimum error of the linear solver for implicit formulations LINEAR_SOLVER_ERROR= 1E-6 %