diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index c4e3b652b06..753a774d2b3 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -233,7 +233,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -282,7 +282,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} -a "--tapetests" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -330,7 +330,7 @@ jobs: PMIX_MCA_gds: hash with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054 with: @@ -375,7 +375,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: diff --git a/Common/include/basic_types/datatype_structure.hpp b/Common/include/basic_types/datatype_structure.hpp index 623c410e551..d9b6bb0b58f 100644 --- a/Common/include/basic_types/datatype_structure.hpp +++ b/Common/include/basic_types/datatype_structure.hpp @@ -135,6 +135,10 @@ template <> struct Passive { FORCEINLINE static passivedouble Value(const su2double& val) { return GetValue(val); } }; +template +FORCEINLINE auto PassiveValue(const T& val) { + return Passive::Value(val); +} /*! * \brief Casts the primitive value to int (uses GetValue, already implemented for each type). diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index bc02c2830aa..ae30c779e7a 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -26,6 +26,7 @@ */ #pragma once +#include #include #include @@ -100,6 +101,10 @@ FORCEINLINE Out su2staticcast_p(In ptr) { #undef USE_SINGLE_PRECISION #endif +/*--- Default integer types. Currently used for rank-local sparse patterns. ---*/ +using su2uint = uint32_t; +using su2int = int32_t; + /*--- This type can be used for (rare) compatibility cases or for * computations that are intended to be (always) passive. ---*/ #ifdef USE_SINGLE_PRECISION diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 5d1e7e5a83a..97d3f845537 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -205,10 +205,10 @@ class CGeometry { LDUSparsePattern finiteElementPatternFill0; /*!< \brief FEM sparsity with 0-fill (structural pattern). */ LDUSparsePattern finiteElementPatternFillN; /*!< \brief FEM sparsity with N-fill (e.g. for ILU-N). */ - su2vector finiteVolumeLToUTranspMap; /*!< \brief FVM L-entry -> U-entry of its transpose. */ - su2vector finiteVolumeUToLTranspMap; /*!< \brief FVM U-entry -> L-entry of its transpose. */ - su2vector finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */ - su2vector finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */ + su2vector finiteVolumeLToUTranspMap; /*!< \brief FVM L-entry -> U-entry of its transpose. */ + su2vector finiteVolumeUToLTranspMap; /*!< \brief FVM U-entry -> L-entry of its transpose. */ + su2vector finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */ + su2vector finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */ /*--- Edge and element colorings. ---*/ @@ -1892,7 +1892,7 @@ class CGeometry { * \param[in] type - Finite volume or finite element. * \return Reference to the l_to_u map. */ - const su2vector& GetLToUTransposeSparsePatternMap(ConnectivityType type); + const su2vector& GetLToUTransposeSparsePatternMap(ConnectivityType type); /*! * \brief Get the bijective map from U-entry indices to L-entry indices of their transposes. @@ -1900,7 +1900,7 @@ class CGeometry { * \param[in] type - Finite volume or finite element. * \return Reference to the u_to_l map. */ - const su2vector& GetUToLTransposeSparsePatternMap(ConnectivityType type); + const su2vector& GetUToLTransposeSparsePatternMap(ConnectivityType type); /*! * \brief Get the edge coloring. diff --git a/Common/include/linear_algebra/CPastixWrapper.hpp b/Common/include/linear_algebra/CPastixWrapper.hpp index a471a60eaff..44fb9d503db 100644 --- a/Common/include/linear_algebra/CPastixWrapper.hpp +++ b/Common/include/linear_algebra/CPastixWrapper.hpp @@ -28,6 +28,8 @@ #pragma once +#include "../code_config.hpp" + #ifdef HAVE_PASTIX #ifdef CODI_FORWARD_TYPE @@ -74,11 +76,11 @@ class CPastixWrapper { unsigned long nPointDomain = 0; unsigned long blkSz = 0; /*!< \brief Block size (nVar * nVar) for value assembly. */ - const unsigned long* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */ - const unsigned long* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */ - const ScalarType* d = nullptr; /*!< \brief Diagonal blocks (matrix-owned). */ - const ScalarType* l = nullptr; /*!< \brief Lower blocks (matrix-owned). */ - const ScalarType* u = nullptr; /*!< \brief Upper blocks (matrix-owned). */ + const su2uint* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */ + const su2uint* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */ + const ScalarType* d = nullptr; /*!< \brief Diagonal blocks (matrix-owned). */ + const ScalarType* l = nullptr; /*!< \brief Lower blocks (matrix-owned). */ + const ScalarType* u = nullptr; /*!< \brief Upper blocks (matrix-owned). */ unsigned long size_rhs() const { return nPointDomain * nVar; } } matrix; /*!< \brief Dimensions and LDU pointers captured from the owning CSysMatrix. */ @@ -150,9 +152,9 @@ class CPastixWrapper { * \param[in] col_ind_l/u - LDU lower/upper column indices (geometry-owned). * \param[in] d/l/u - LDU value blocks (matrix-owned, must outlive wrapper). */ - void SetLDU(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const unsigned long* row_ptr_l, - const unsigned long* col_ind_l, const unsigned long* row_ptr_u, const unsigned long* col_ind_u, - const ScalarType* d, const ScalarType* l, const ScalarType* u) { + void SetLDU(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const su2uint* row_ptr_l, + const su2uint* col_ind_l, const su2uint* row_ptr_u, const su2uint* col_ind_u, const ScalarType* d, + const ScalarType* l, const ScalarType* u) { if (issetup) return; matrix.nVar = nVar; matrix.nPoint = nPoint; diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index e4fc7cf159f..532e0b53838 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -210,6 +210,37 @@ class CLU_SGSPreconditioner final : public CPreconditioner { } }; +/*! + * \class CQuantizedLUSGSPreconditioner + * \brief Specialization of preconditioner that uses CSysMatrix class. + */ +template +class CQuantizedLUSGSPreconditioner final : public CPreconditioner { + private: + CSysMatrix& sparse_matrix; + CGeometry* geometry; + const CConfig* config; + + public: + inline CQuantizedLUSGSPreconditioner(CSysMatrix& matrix_ref, CGeometry* geometry_ref, + const CConfig* config_ref) + : sparse_matrix(matrix_ref) { + if ((geometry_ref == nullptr) || (config_ref == nullptr)) + SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION); + geometry = geometry_ref; + config = config_ref; + } + + CQuantizedLUSGSPreconditioner() = delete; + + inline void operator()(const CSysVector& u, CSysVector& v) const override { + sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); + } + + /*! \brief Quantize the diagonal blocks (off diagonals are quantized on the fly). */ + inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } +}; + /*! * \class CLineletPreconditioner * \brief Specialization of preconditioner that uses CSysMatrix class. @@ -341,6 +372,9 @@ CPreconditioner* CPreconditioner::Create(ENUM_LINEAR_SOL case LU_SGS: prec = new CLU_SGSPreconditioner(jacobian, geometry, config); break; + case Q_LU_SGS: + prec = new CQuantizedLUSGSPreconditioner(jacobian, geometry, config); + break; case ILU: prec = new CILUPreconditioner(jacobian, geometry, config); break; diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 11baf6edf8e..d31c48fb8a0 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -33,6 +33,7 @@ #include "CPastixWrapper.hpp" #include "../toolboxes/graph_toolbox.hpp" +#include #include #include #include @@ -108,6 +109,96 @@ struct CSysMatrixComms { MPI_QUANTITIES commType = MPI_QUANTITIES::SOLUTION_MATRIX); }; +/*! + * \brief Reconstruct the float row-scale from a stored int8 binary exponent. + * The exponent \p e was packed as (e + 127) into the IEEE 754 biased-exponent field + * with a zero mantissa, giving an exact power of two: 2^e. + * This is the inverse of the encoding in EncodeQuantBlock. + */ +FORCEINLINE float DecodeQuantScale(int8_t e) noexcept { + const uint32_t bits = static_cast(std::max(0, static_cast(e) + 127)) << 23; + float scale; + memcpy(&scale, &bits, sizeof(bits)); + return scale; +} + +/*! + * \brief Encode one nVar×nVar block into per-row int8 quantized storage. + * \p f(r,c) is called twice per entry (max-abs scan then encoding); it should be cheap. + * Stores a per-row scale exponent in \p qs and clamped int8 values in \p qv. + */ +template +FORCEINLINE void EncodeQuantBlock(const F& f, int8_t* __restrict qs, int8_t* __restrict qv, + unsigned long nVar) noexcept { + for (auto r = 0ul; r < nVar; ++r) { + constexpr uint32_t eps_bits = 0x34000000u; + uint32_t max_abs_bits = eps_bits; + for (auto c = 0ul; c < nVar; ++c) { + const float fv = SU2_TYPE::PassiveValue(f(r, c)); + uint32_t fb; + memcpy(&fb, &fv, sizeof(fb)); + max_abs_bits = std::max(max_abs_bits, fb & 0x7FFFFFFFu); + } + const int e = std::min(127, std::max(-128, static_cast(max_abs_bits >> 23) - 133)); + qs[r] = static_cast(e); + const uint32_t inv_bits = static_cast(127 - e) << 23; + float inv_rscale; + memcpy(&inv_rscale, &inv_bits, sizeof(inv_rscale)); + for (auto c = 0ul; c < nVar; ++c) { + qv[r * nVar + c] = + static_cast(std::max(-128.f, std::min(127.f, roundf(SU2_TYPE::PassiveValue(f(r, c)) * inv_rscale)))); + } + } +} + +/*! + * \brief View of one matrix block, const-correct via the ScalarType template parameter. + * \c CBlockView is read-only; \c CBlockView is mutable + * and exposes \c apply(f) for writing with on-the-fly quantized encoding. + * Evaluates to \c false if the block is absent from the sparsity pattern. + */ +template +struct CBlockView { + using QuantType = std::conditional_t, const int8_t, int8_t>; + + ScalarType* ptr = nullptr; ///< Full-precision block; non-null iff not quantized. + QuantType* qs = nullptr; ///< Per-row binary exponent; non-null iff quantized. + QuantType* qv = nullptr; ///< Quantized values (row-major); non-null iff quantized. + unsigned long nVar = 0; + + /*! \brief False when the block is not present in the sparsity pattern. */ + explicit operator bool() const { return ptr || qs; } + + /*! \brief Return entry (row \p i, col \p j), decoding quantization if necessary. */ + std::remove_const_t operator()(unsigned long i, unsigned long j) const { + using T = std::remove_const_t; + if (ptr) return ptr[i * nVar + j]; + return static_cast(qv[i * nVar + j] * DecodeQuantScale(qs[i])); + } + + /*! + * \brief Write the block from callable \p f(i,j). + * \p Overwrite=true overwrites (or quantizes for Q_LU_SGS off-diagonal blocks); + * \p Overwrite=false accumulates into non-quantized storage only — accumulating into + * quantized storage would require decode-accumulate-encode and is a silent no-op. + * Only enabled for mutable (non-const ScalarType) views. + */ + template > = 0> + void apply(const F& f) const { + if (ptr) { + for (auto i = 0ul; i < nVar; ++i) + for (auto j = 0ul; j < nVar; ++j) { + if constexpr (Overwrite) + ptr[i * nVar + j] = f(i, j); + else + ptr[i * nVar + j] += f(i, j); + } + } else if constexpr (Overwrite) { + if (qs) EncodeQuantBlock(f, qs, qv, nVar); + } + } +}; + /*! * \class CSysMatrix * \ingroup SpLinSys @@ -145,31 +236,50 @@ class CSysMatrix { * the pointers address host or device memory is managed by CSysMatrix. */ struct LDU { - ScalarType* d = nullptr; /*!< \brief Diagonal block values. */ - ScalarType* l = nullptr; /*!< \brief Strictly-lower block values. */ - ScalarType* u = nullptr; /*!< \brief Strictly-upper block values. */ - const unsigned long* row_ptr_l = nullptr; /*!< \brief Row pointers for L (geometry-owned or GPU copy). */ - const unsigned long* col_ind_l = nullptr; /*!< \brief Column indices for L. */ - const unsigned long* row_ptr_u = nullptr; /*!< \brief Row pointers for U. */ - const unsigned long* col_ind_u = nullptr; /*!< \brief Column indices for U. */ - unsigned long nnz_l = 0; /*!< \brief Number of L nonzeros. */ - unsigned long nnz_u = 0; /*!< \brief Number of U nonzeros. */ + ScalarType* d = nullptr; /*!< \brief Diagonal block values. */ + ScalarType* l = nullptr; /*!< \brief Strictly-lower block values. */ + ScalarType* u = nullptr; /*!< \brief Strictly-upper block values. */ + const su2uint* row_ptr_l = nullptr; /*!< \brief Row pointers for L (geometry-owned or GPU copy). */ + const su2uint* col_ind_l = nullptr; /*!< \brief Column indices for L. */ + const su2uint* row_ptr_u = nullptr; /*!< \brief Row pointers for U. */ + const su2uint* col_ind_u = nullptr; /*!< \brief Column indices for U. */ + unsigned long nnz_l = 0; /*!< \brief Number of L nonzeros. */ + unsigned long nnz_u = 0; /*!< \brief Number of U nonzeros. */ }; LDU mat; /*!< \brief Host matrix (values owned via aligned_alloc; pattern from geometry). */ LDU gpu; /*!< \brief Device matrix (all pointers to GPU memory). */ LDU ilu; /*!< \brief ILU factorization, host (values owned; pattern from geometry). */ - bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ - const unsigned long* l_to_u_transp; /*!< \brief L-entry index -> U-entry index of its transpose. */ - const unsigned long* u_to_l_transp; /*!< \brief U-entry index -> L-entry index of its transpose. */ + /*--- Quantized off-diagonal storage (used when quantized_mode == true). ---*/ + using QuantType = int8_t; + + /*! \brief Set by Initialize() when preconditioner == Q_LU_SGS. + * mat.l and mat.u are NOT allocated; off-diagonal blocks live in the + * q_* arrays below. */ +#ifndef CODI_REVERSE_TYPE + bool quantized_mode = false; +#else + static constexpr bool quantized_mode = false; +#endif + QuantType* q_scale_l; /*!< \brief Per-row exponent for L blocks, [nnz_l * nVar]. */ + QuantType* q_blocks_l; /*!< \brief Quantized L block entries, [nnz_l * nVar * nEqn]. */ + QuantType* q_scale_u; /*!< \brief Same as q_scale_l for the upper entries. */ + QuantType* q_blocks_u; /*!< \brief Same as q_blocks_l for the upper entries. */ + QuantType* q_scale_d; /*!< \brief Same as q_scale_l for the diagonal entries, [nPoint * nVar]. + * 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. */ + 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. */ /*! * \brief Lookup table from edges to the L-index in the LDU split. * U-index == edge index by construction (edges are ordered 1:1 with the U pattern). * Therefore, edge_ptr_l == u_to_l_transp, but we keep a separate member for clarity. */ - const unsigned long* edge_ptr_l; + const su2uint* edge_ptr_l; unsigned short ilu_fill_in; /*!< \brief Fill level for the ILU preconditioner. */ @@ -179,11 +289,11 @@ class CSysMatrix { ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/ - mutable vector > + mutable vector> LineletUpper; /*!< \brief Pointers to the upper blocks of the tri-diag system (working memory). */ - mutable vector > + mutable vector> LineletInvDiag; /*!< \brief Inverse of the diagonal blocks of the tri-diag system (working memory). */ - mutable vector > + mutable vector> LineletVector; /*!< \brief Solution and RHS of the tri-diag system (working memory). */ #ifdef USE_MKL @@ -290,7 +400,7 @@ class CSysMatrix { * \param[in,out] matrix - On entry the system matrix, on exit the factorized matrix. * \param[in,out] vec - On entry the rhs, on exit the solution. */ - void Gauss_Elimination(ScalarType* matrix, ScalarType* vec) const; + void GaussElimination(ScalarType* matrix, ScalarType* vec) const; /*! * \brief Invert a small dense matrix. @@ -305,7 +415,7 @@ class CSysMatrix { * \param[in] rhs - Right-hand-side of the linear system. * \return Solution of the linear system (overwritten on rhs). */ - inline void Gauss_Elimination(unsigned long block_i, ScalarType* rhs) const; + inline void GaussElimination(unsigned long block_i, ScalarType* rhs) const; /*! * \brief Inverse diagonal block. @@ -363,6 +473,38 @@ class CSysMatrix { */ void RowProduct(const CSysVector& vec, unsigned long row_i, ScalarType* prod) const; + /*! + * \brief Computes product += A_k * vec using the quantized representation of block k. + * \note Only valid after QuantizeDiagonalBlocks() has been called. + * \param[in] k - Block index in the CSR flat storage. + * \param[in] vec - Input vector (nEqn entries). + * \param[in,out] prod - Accumulation output (nVar entries). + */ + inline void QuantizedMatVecAdd(const QuantType* qs, const QuantType* qv, const ScalarType* vec, + ScalarType* prod) const; + + /*! \brief Quantize one nVar×nVar block (row-major) into the int8 scale+value arrays. + * Called on the hot assembly path (SetBlocks/UpdateBlocks in Q_LU_SGS mode). */ + void QuantizeBlock(const ScalarType* blk, QuantType* qs, QuantType* qv) const; + + /*! \brief Full-row product using quantized L/D/U (Q_LU_SGS SpMV path). */ + inline void QuantizedRowProduct(const CSysVector& vec, unsigned long row_i, ScalarType* prod) const; + + /*! \brief Upper-triangle product using quantized U (Q_LU_SGS backward sweep). */ + inline void QuantizedUpperProduct(const CSysVector& vec, unsigned long row_i, unsigned long col_ub, + ScalarType* prod) const; + + /*! \brief Lower-triangle product using quantized L (Q_LU_SGS forward sweep). */ + inline void QuantizedLowerProduct(const CSysVector& vec, unsigned long row_i, unsigned long col_lb, + ScalarType* prod) const; + + /*! \brief Diagonal product using quantized D (Q_LU_SGS backward sweep). */ + inline void QuantizedDiagonalProduct(const CSysVector& vec, unsigned long row_i, ScalarType* prod) const; + + /*! \brief Gauss elimination on the quantized diagonal block: decodes q_blocks_d into a local + * ScalarType buffer and delegates to the scalar GaussElimination overload. */ + inline void QuantizedGaussElimination(unsigned long block_i, ScalarType* rhs) const; + public: /*! * \brief Constructor of the class. @@ -384,10 +526,20 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. * \param[in] needTranspPtr - If the L/U transpose maps should be built, used for "SetDiagonalAsColumnSum". + * \param[in] grad_mode - Gradient smoothing mode, only used to detect the right preconditioner type. + * \param[in] allow_quant - Quantization is only possible with solvers that "set and forget" the off-diagonal + * blocks of the matrix. Solvers that perform multiple updates would lose too much information, so + * that pattern is not supported with quantization (the code will hit null pointers). It is up to + * the solver to declare whether it will "set and forget". */ void Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry* geometry, const CConfig* config, bool needTranspPtr = false, - bool grad_mode = false); + bool grad_mode = false, bool allow_quant = false); + + /*! + * \brief Compresses off-diagonal blocks into quantized form for use with USE_QUANTIZATION. + */ + void QuantizeDiagonalBlocks(); /*! * \brief Sets to zero all the entries of the sparse matrix. @@ -424,7 +576,7 @@ class CSysMatrix { } /*! - * \brief Get a pointer to the start of block "ij", non-const version + * \brief Get a pointer to the start of block "ij", non-const version. */ FORCEINLINE ScalarType* GetBlock(unsigned long block_i, unsigned long block_j) { const CSysMatrix& const_this = *this; @@ -432,131 +584,98 @@ class CSysMatrix { } /*! - * \brief Gets the value of a particular entry in block "ij". - * \param[in] block_i - Row index. - * \param[in] block_j - Column index. - * \param[in] iVar - Row of the block. - * \param[in] jVar - Column of the block. - * \return Value of the block entry. - */ - FORCEINLINE ScalarType GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, - unsigned short jVar) const { - auto mat_ij = GetBlock(block_i, block_j); - if (!mat_ij) return 0.0; - return mat_ij[iVar * nEqn + jVar]; + * \brief Read-only view of block (block_i, block_j). In Q_LU_SGS mode values are decoded + * on access inside CBlockView::operator()(i,j); no temporary copy is made. + * \return A CBlockView that evaluates to false if the block is absent. + */ + FORCEINLINE CBlockView GetBlockView(unsigned long block_i, unsigned long block_j) const { +#define GET_BLOCK_VIEW_IMPL \ + if (!quantized_mode || block_i == block_j) { \ + return {GetBlock(block_i, block_j), nullptr, nullptr, nVar}; \ + } \ + if (block_j < block_i) { \ + for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) \ + if (mat.col_ind_l[k] == block_j) return {nullptr, &q_scale_l[k * nVar], &q_blocks_l[k * nVar * nVar], nVar}; \ + } else { \ + for (auto k = mat.row_ptr_u[block_i]; k < mat.row_ptr_u[block_i + 1]; ++k) \ + if (mat.col_ind_u[k] == block_j) return {nullptr, &q_scale_u[k * nVar], &q_blocks_u[k * nVar * nVar], nVar}; \ + } \ + return {} + GET_BLOCK_VIEW_IMPL; } /*! - * \brief Set the value of a block (in flat format) in the sparse matrix with scaling. - * \note If the template param Overwrite is false we add to the block (bij += alpha*b). - * \param[in] block_i - Row index. - * \param[in] block_j - Column index. - * \param[in] val_block - Block to set to A(i, j). - * \param[in] alpha - Scale factor. - */ - template ::value> = 0> - inline void SetBlock(unsigned long block_i, unsigned long block_j, const OtherType* val_block, - OtherType alpha = 1.0) { - auto mat_ij = GetBlock(block_i, block_j); - if (!mat_ij) return; - SU2_OMP_SIMD - for (auto iVar = 0ul; iVar < nVar * nEqn; ++iVar) { - mat_ij[iVar] = (Overwrite ? ScalarType(0) : mat_ij[iVar]) + PassiveAssign(alpha * val_block[iVar]); - } - } - - /*! - * \brief Add a scaled block (in flat format) to the sparse matrix (see SetBlock). - * \param[in] block_i - Row index. - * \param[in] block_j - Column index. - * \param[in] val_block - Block to set to A(i, j). - * \param[in] alpha - Scale factor. + * \overload Non const version of GetBlockView. */ - template ::value> = 0> - inline void AddBlock(unsigned long block_i, unsigned long block_j, const OtherType* val_block, - OtherType alpha = 1.0) { - SetBlock(block_i, block_j, val_block, alpha); + FORCEINLINE CBlockView GetBlockView(unsigned long block_i, unsigned long block_j) { + GET_BLOCK_VIEW_IMPL; +#undef GET_BLOCK_VIEW_IMPL } /*! * \brief Set the value of a scaled block in the sparse matrix. - * \note If the template param Overwrite is false we add to the block (bij += alpha*b). + * \note This is an templated overload for C2Dcontainer specialization su2matrix. + * It assumes that MatrixType supports a member type Scalar and access operator(i, j). + * If the template param Overwrite is false we add to the block (bij += alpha*b). * \param[in] block_i - Row index. * \param[in] block_j - Column index. * \param[in] val_block - Block to set to A(i, j). * \param[in] alpha - Scale factor. */ - template - inline void SetBlock(unsigned long block_i, unsigned long block_j, const OtherType* const* val_block, - OtherType alpha = 1.0) { - auto mat_ij = GetBlock(block_i, block_j); - if (!mat_ij) return; - for (auto iVar = 0ul; iVar < nVar; ++iVar) { - for (auto jVar = 0ul; jVar < nEqn; ++jVar) { - *mat_ij = (Overwrite ? ScalarType(0) : *mat_ij) + PassiveAssign(alpha * val_block[iVar][jVar]); - ++mat_ij; - } - } + template + inline void SetBlock(unsigned long block_i, unsigned long block_j, MatrixType& val_block, + std::decay_t alpha = 1.0) { + auto view = GetBlockView(block_i, block_j); + if (!view) return; + view.template apply( + [&](unsigned long i, unsigned long j) { return PassiveAssign(alpha * val_block(i, j)); }); } /*! - * \brief Adds a scaled block to the sparse matrix (see SetBlock). - * \param[in] block_i - Row index. - * \param[in] block_j - Column index. - * \param[in] val_block - Block to add to A(i, j). - * \param[in] alpha - Scale factor. + * \overload val_block is a pointer instead of a matrix type. */ - template - inline void AddBlock(unsigned long block_i, unsigned long block_j, const OtherType* const* val_block, - OtherType alpha = 1.0) { - SetBlock(block_i, block_j, val_block, alpha); + template ::value> = 0> + inline void SetBlock(unsigned long block_i, unsigned long block_j, const OtherType* val_block, + std::decay_t alpha = 1.0) { + auto view = GetBlockView(block_i, block_j); + if (!view) return; + view.template apply( + [&](unsigned long i, unsigned long j) { return PassiveAssign(alpha * val_block[i * nEqn + j]); }); } /*! - * \brief Subtracts the specified block to the sparse matrix (see AddBlock). - * \param[in] block_i - Row index. - * \param[in] block_j - Column index. - * \param[in] val_block - Block to subtract to A(i, j). + * \overload val_block is a double pointer instead of matrix type. */ - template - inline void SubtractBlock(unsigned long block_i, unsigned long block_j, const OtherType* const* val_block) { - AddBlock(block_i, block_j, val_block, OtherType(-1)); + template + inline void SetBlock(unsigned long block_i, unsigned long block_j, const OtherType* const* val_block, + std::decay_t alpha = 1.0) { + auto view = GetBlockView(block_i, block_j); + if (!view) return; + view.template apply( + [&](unsigned long i, unsigned long j) { return PassiveAssign(alpha * val_block[i][j]); }); } /*! - * \brief Set the value of a scaled block in the sparse matrix. - * \note This is an templated overload for C2Dcontainer specialization su2matrix. - * It assumes that MatrixType supports a member type Scalar and access operator[][]. - * If the template param Overwrite is false we add to the block (bij += alpha*b). + * \brief Add a scaled block (in flat format) to the sparse matrix (see SetBlock). * \param[in] block_i - Row index. * \param[in] block_j - Column index. * \param[in] val_block - Block to set to A(i, j). * \param[in] alpha - Scale factor. */ - template - inline void SetBlock(unsigned long block_i, unsigned long block_j, MatrixType& val_block, - typename MatrixType::Scalar alpha = 1.0) { - auto mat_ij = GetBlock(block_i, block_j); - if (!mat_ij) return; - for (auto iVar = 0ul; iVar < nVar; ++iVar) { - for (auto jVar = 0ul; jVar < nEqn; ++jVar) { - *mat_ij = (Overwrite ? ScalarType(0) : *mat_ij) + PassiveAssign(alpha * val_block(iVar, jVar)); - ++mat_ij; - } - } + template + inline void AddBlock(unsigned long block_i, unsigned long block_j, const T& val_block, OtherType alpha = 1.0) { + SetBlock(block_i, block_j, val_block, alpha); } /*! - * \brief Adds a scaled block to the sparse matrix (see SetBlock). + * \brief Subtracts the specified block to the sparse matrix (see AddBlock). * \param[in] block_i - Row index. * \param[in] block_j - Column index. - * \param[in] val_block - Block to add to A(i, j). - * \param[in] alpha - Scale factor. + * \param[in] val_block - Block to subtract to A(i, j). */ - template - inline void AddBlock(unsigned long block_i, unsigned long block_j, MatrixType& val_block, - typename MatrixType::Scalar alpha = 1.0) { - SetBlock(block_i, block_j, val_block, alpha); + template + inline void SubtractBlock(unsigned long block_i, unsigned long block_j, const T& val_block) { + AddBlock(block_i, block_j, val_block, -1); } /*! @@ -589,11 +708,31 @@ class CSysMatrix { template inline void UpdateBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, const MatrixType& block_i, const MatrixType& block_j, OtherType scale = 1) { - ScalarType *bii, *bij, *bji, *bjj; - GetBlocks(iEdge, iPoint, jPoint, bii, bij, bji, bjj); + const auto blkSz = nVar * nEqn; + auto* bii = &mat.d[iPoint * blkSz]; + auto* bjj = &mat.d[jPoint * blkSz]; unsigned long iVar, jVar, offset = 0; + if (quantized_mode) { + assert(OverwriteOffDiag); + /*--- Diagonal: full-precision accumulation. Off-diagonal: quantize on the fly. ---*/ + ScalarType bij_buf[MAXNVAR * MAXNVAR], bji_buf[MAXNVAR * MAXNVAR]; + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++, ++offset) { + bii[offset] += PassiveAssign(block_i[iVar][jVar] * scale); + bjj[offset] -= PassiveAssign(block_j[iVar][jVar] * scale); + bij_buf[offset] = PassiveAssign(block_j[iVar][jVar] * scale); + bji_buf[offset] = -PassiveAssign(block_i[iVar][jVar] * scale); + } + QuantizeBlock(bij_buf, &q_scale_u[iEdge * nVar], &q_blocks_u[iEdge * blkSz]); + const auto k_l = edge_ptr_l[iEdge]; + QuantizeBlock(bji_buf, &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + return; + } + + auto* bij = &mat.u[iEdge * blkSz]; + auto* bji = &mat.l[edge_ptr_l[iEdge] * blkSz]; for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nEqn; jVar++) { bii[offset] += PassiveAssign(block_i[iVar][jVar] * scale); @@ -623,9 +762,9 @@ class CSysMatrix { * \brief SIMD version, does the update for multiple edges and points. * \note Nothing is updated if the mask is 0. */ - template - FORCEINLINE void UpdateBlocks(simd::Array iEdge, simd::Array iPoint, simd::Array jPoint, - const MatTypeSIMD& block_i, const MatTypeSIMD& block_j, simd::Array mask = 1) { + template + FORCEINLINE void SetBlocks(simd::Array iEdge, simd::Array iPoint, simd::Array jPoint, + const MatTypeSIMD& block_i, const MatTypeSIMD& block_j, simd::Array mask = 1) { static_assert(MatTypeSIMD::StaticSize, "This method requires static size blocks."); static_assert(MatTypeSIMD::IsRowMajor, "Block storage is not compatible with matrix."); constexpr size_t blkSz = MatTypeSIMD::StaticSize; @@ -647,20 +786,30 @@ class CSysMatrix { for (size_t k = 0; k < N; ++k) { if (mask[k] == 0) continue; - /*--- Fetch the blocks. ---*/ auto bii = &mat.d[iPoint[k] * blkSz]; auto bjj = &mat.d[jPoint[k] * blkSz]; - auto bij = &mat.u[iEdge[k] * blkSz]; - auto bji = &mat.l[edge_ptr_l[iEdge[k]] * blkSz]; - - /*--- Update, block i was negated during transpose in the - * hope the assignments below become non-temporal stores. ---*/ - SU2_OMP_SIMD - for (size_t i = 0; i < blkSz; ++i) { - bii[i] -= blk_i[k][i]; - bjj[i] -= blk_j[k][i]; - bij[i] = blk_j[k][i]; - bji[i] = blk_i[k][i]; + + if (quantized_mode) { + SU2_OMP_SIMD + for (size_t i = 0; i < blkSz; ++i) { + bii[i] -= blk_i[k][i]; + bjj[i] -= blk_j[k][i]; + } + QuantizeBlock(blk_j[k], &q_scale_u[iEdge[k] * nVar], &q_blocks_u[iEdge[k] * blkSz]); + const auto k_l = edge_ptr_l[iEdge[k]]; + QuantizeBlock(blk_i[k], &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + } else { + auto bij = &mat.u[iEdge[k] * blkSz]; + auto bji = &mat.l[edge_ptr_l[iEdge[k]] * blkSz]; + /*--- Update, block i was negated during transpose in the + * hope the assignments below become non-temporal stores. ---*/ + SU2_OMP_SIMD + for (size_t i = 0; i < blkSz; ++i) { + bii[i] -= blk_i[k][i]; + bjj[i] -= blk_j[k][i]; + bij[i] = blk_j[k][i]; + bji[i] = blk_i[k][i]; + } } } } @@ -679,11 +828,24 @@ class CSysMatrix { inline void SetBlocks(unsigned long iEdge, const MatrixType& block_i, const MatrixType& block_j, OtherType scale = 1) { const auto blkSz = nVar * nEqn; - ScalarType* bij = &mat.u[iEdge * blkSz]; - ScalarType* bji = &mat.l[edge_ptr_l[iEdge] * blkSz]; - unsigned long iVar, jVar, offset = 0; + if (quantized_mode) { + assert(Overwrite); + ScalarType bij_buf[MAXNVAR * MAXNVAR], bji_buf[MAXNVAR * MAXNVAR]; + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++, ++offset) { + bij_buf[offset] = PassiveAssign(block_j[iVar][jVar] * scale); + bji_buf[offset] = -PassiveAssign(block_i[iVar][jVar] * scale); + } + QuantizeBlock(bij_buf, &q_scale_u[iEdge * nVar], &q_blocks_u[iEdge * blkSz]); + const auto k_l = edge_ptr_l[iEdge]; + QuantizeBlock(bji_buf, &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + return; + } + + ScalarType* bij = &mat.u[iEdge * blkSz]; + ScalarType* bji = &mat.l[edge_ptr_l[iEdge] * blkSz]; for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nEqn; jVar++) { bij[offset] = (Overwrite ? ScalarType(0) : bij[offset]) + PassiveAssign(block_j[iVar][jVar] * scale); @@ -738,16 +900,20 @@ class CSysMatrix { for (size_t k = 0; k < N; ++k) { if (mask[k] == 0) continue; - /*--- Fetch the blocks. ---*/ - ScalarType* bij = &mat.u[iEdge[k] * blkSz]; - ScalarType* bji = &mat.l[edge_ptr_l[iEdge[k]] * blkSz]; - - /*--- Update, block i was negated during transpose in the - * hope the assignments below become non-temporal stores. ---*/ - SU2_OMP_SIMD - for (size_t i = 0; i < blkSz; ++i) { - bij[i] = blk_j[k][i]; - bji[i] = blk_i[k][i]; + if (quantized_mode) { + QuantizeBlock(blk_j[k], &q_scale_u[iEdge[k] * nVar], &q_blocks_u[iEdge[k] * blkSz]); + const auto k_l = edge_ptr_l[iEdge[k]]; + QuantizeBlock(blk_i[k], &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + } else { + ScalarType* bij = &mat.u[iEdge[k] * blkSz]; + ScalarType* bji = &mat.l[edge_ptr_l[iEdge[k]] * blkSz]; + /*--- Update, block i was negated during transpose in the + * hope the assignments below become non-temporal stores. ---*/ + SU2_OMP_SIMD + for (size_t i = 0; i < blkSz; ++i) { + bij[i] = blk_j[k][i]; + bji[i] = blk_i[k][i]; + } } } } diff --git a/Common/include/linear_algebra/CSysMatrix.inl b/Common/include/linear_algebra/CSysMatrix.inl index d40d68b3a4f..b5411016af6 100644 --- a/Common/include/linear_algebra/CSysMatrix.inl +++ b/Common/include/linear_algebra/CSysMatrix.inl @@ -137,12 +137,23 @@ FORCEINLINE void CSysMatrix::MatrixMatrixProduct(const ScalarType* m #undef __MATVECPROD_SIGNATURE__ template -FORCEINLINE void CSysMatrix::Gauss_Elimination(unsigned long block_i, ScalarType* rhs) const { +FORCEINLINE void CSysMatrix::GaussElimination(unsigned long block_i, ScalarType* rhs) const { /*--- Copy block, as the algorithm modifies the matrix ---*/ ScalarType block[MAXNVAR * MAXNVAR]; MatrixCopy(&mat.d[block_i * nVar * nVar], block); + GaussElimination(block, rhs); +} - Gauss_Elimination(block, rhs); +template +FORCEINLINE void CSysMatrix::QuantizedGaussElimination(unsigned long block_i, ScalarType* rhs) const { + ScalarType block[MAXNVAR * MAXNVAR]; + const QuantType* __restrict qs = &q_scale_d[block_i * nVar]; + const QuantType* __restrict qv = &q_blocks_d[block_i * nVar * nVar]; + for (auto r = 0ul; r < nVar; ++r) { + const float row_scale = DecodeQuantScale(qs[r]); + for (auto c = 0ul; c < nVar; ++c) block[r * nVar + c] = static_cast(qv[r * nVar + c] * row_scale); + } + GaussElimination(block, rhs); } template @@ -164,18 +175,32 @@ FORCEINLINE const ScalarType* CSysMatrix::InvertDiagonalBlockILUMatr return Uii; } +template +FORCEINLINE void CSysMatrix::QuantizedMatVecAdd(const QuantType* __restrict qs, + const QuantType* __restrict qv, + const ScalarType* __restrict vec, + ScalarType* __restrict prod) const { + for (auto r = 0ul; r < nVar; ++r) { + const float row_scale = DecodeQuantScale(qs[r]); + auto sum = ScalarType(0); + for (auto c = 0ul; c < nVar; ++c) sum += qv[r * nVar + c] * vec[c]; + prod[r] += row_scale * sum; + } +} + template FORCEINLINE void CSysMatrix::RowProduct(const CSysVector& vec, unsigned long row_i, ScalarType* prod) const { for (auto iVar = 0ul; iVar < nVar; iVar++) prod[iVar] = 0.0; - for (auto index = mat.row_ptr_l[row_i]; index < mat.row_ptr_l[row_i + 1]; index++) - MatrixVectorProductAdd(&mat.l[index * nVar * nEqn], &vec[mat.col_ind_l[index] * nEqn], prod); - + for (auto k = mat.row_ptr_l[row_i]; k < mat.row_ptr_l[row_i + 1]; k++) { + MatrixVectorProductAdd(&mat.l[k * nVar * nEqn], &vec[mat.col_ind_l[k] * nEqn], prod); + } MatrixVectorProductAdd(&mat.d[row_i * nVar * nEqn], &vec[row_i * nEqn], prod); - for (auto index = mat.row_ptr_u[row_i]; index < mat.row_ptr_u[row_i + 1]; index++) - MatrixVectorProductAdd(&mat.u[index * nVar * nEqn], &vec[mat.col_ind_u[index] * nEqn], prod); + for (auto k = mat.row_ptr_u[row_i]; k < mat.row_ptr_u[row_i + 1]; k++) { + MatrixVectorProductAdd(&mat.u[k * nVar * nEqn], &vec[mat.col_ind_u[k] * nEqn], prod); + } } template @@ -185,9 +210,10 @@ FORCEINLINE void CSysMatrix::UpperProduct(const CSysVector= nPointDomain) + + if (col_j < col_ub || col_j >= nPointDomain) { MatrixVectorProductAdd(&mat.u[index * nVar * nEqn], &vec[col_j * nEqn], prod); + } } } @@ -198,7 +224,9 @@ FORCEINLINE void CSysMatrix::LowerProduct(const CSysVector= col_lb) MatrixVectorProductAdd(&mat.l[index * nVar * nEqn], &vec[col_j * nEqn], prod); + if (col_j >= col_lb) { + MatrixVectorProductAdd(&mat.l[index * nVar * nEqn], &vec[col_j * nEqn], prod); + } } } @@ -207,3 +235,51 @@ FORCEINLINE void CSysMatrix::DiagonalProduct(const CSysVector +FORCEINLINE void CSysMatrix::QuantizedRowProduct(const CSysVector& vec, unsigned long row_i, + ScalarType* prod) const { + for (auto iVar = 0ul; iVar < nVar; iVar++) prod[iVar] = 0.0; + + for (auto k = mat.row_ptr_l[row_i]; k < mat.row_ptr_l[row_i + 1]; k++) { + QuantizedMatVecAdd(&q_scale_l[k * nVar], &q_blocks_l[k * nVar * nEqn], &vec[mat.col_ind_l[k] * nEqn], prod); + } + QuantizedMatVecAdd(&q_scale_d[row_i * nVar], &q_blocks_d[row_i * nVar * nEqn], &vec[row_i * nEqn], prod); + + for (auto k = mat.row_ptr_u[row_i]; k < mat.row_ptr_u[row_i + 1]; k++) { + QuantizedMatVecAdd(&q_scale_u[k * nVar], &q_blocks_u[k * nVar * nEqn], &vec[mat.col_ind_u[k] * nEqn], prod); + } +} + +template +FORCEINLINE void CSysMatrix::QuantizedUpperProduct(const CSysVector& vec, unsigned long row_i, + unsigned long col_ub, ScalarType* prod) const { + for (auto iVar = 0ul; iVar < nVar; iVar++) prod[iVar] = 0.0; + + for (auto index = mat.row_ptr_u[row_i]; index < mat.row_ptr_u[row_i + 1]; index++) { + auto col_j = mat.col_ind_u[index]; + if (col_j < col_ub || col_j >= nPointDomain) { + QuantizedMatVecAdd(&q_scale_u[index * nVar], &q_blocks_u[index * nVar * nEqn], &vec[col_j * nEqn], prod); + } + } +} + +template +FORCEINLINE void CSysMatrix::QuantizedLowerProduct(const CSysVector& vec, unsigned long row_i, + unsigned long col_lb, ScalarType* prod) const { + for (auto iVar = 0ul; iVar < nVar; iVar++) prod[iVar] = 0.0; + + for (auto index = mat.row_ptr_l[row_i]; index < mat.row_ptr_l[row_i + 1]; index++) { + auto col_j = mat.col_ind_l[index]; + if (col_j >= col_lb) { + QuantizedMatVecAdd(&q_scale_l[index * nVar], &q_blocks_l[index * nVar * nEqn], &vec[col_j * nEqn], prod); + } + } +} + +template +FORCEINLINE void CSysMatrix::QuantizedDiagonalProduct(const CSysVector& vec, + unsigned long row_i, ScalarType* prod) const { + for (auto iVar = 0ul; iVar < nVar; iVar++) prod[iVar] = 0.0; + QuantizedMatVecAdd(&q_scale_d[row_i * nVar], &q_blocks_d[row_i * nVar * nEqn], &vec[row_i * nEqn], prod); +} diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 0cc823bd879..efbad155492 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2517,6 +2517,7 @@ enum ENUM_LINEAR_SOLVER_PREC { LU_SGS, /*!< \brief LU SGS preconditioner. */ LINELET, /*!< \brief Line implicit preconditioner. */ ILU, /*!< \brief ILU(k) preconditioner. */ + Q_LU_SGS, /*!< \brief LU-SGS with quantized (int8) off-diagonal storage; L/U are never allocated as ScalarType. */ PASTIX_ILU=10, /*!< \brief PaStiX ILU(k) preconditioner. */ PASTIX_LU_P, /*!< \brief PaStiX LU as preconditioner. */ PASTIX_LDLT_P, /*!< \brief PaStiX LDLT as preconditioner. */ @@ -2526,6 +2527,7 @@ static const MapType Linear_Solver_Prec_Ma MakePair("LU_SGS", LU_SGS) MakePair("LINELET", LINELET) MakePair("ILU", ILU) + MakePair("Q_LU_SGS", Q_LU_SGS) MakePair("PASTIX_ILU", PASTIX_ILU) MakePair("PASTIX_LU", PASTIX_LU_P) MakePair("PASTIX_LDLT", PASTIX_LDLT_P) diff --git a/Common/include/toolboxes/graph_toolbox.hpp b/Common/include/toolboxes/graph_toolbox.hpp index 170b973ea1f..f1c36969ada 100644 --- a/Common/include/toolboxes/graph_toolbox.hpp +++ b/Common/include/toolboxes/graph_toolbox.hpp @@ -27,6 +27,7 @@ #pragma once +#include "../code_config.hpp" #include "../containers/C2DContainer.hpp" #include "../parallelization/omp_structure.hpp" @@ -333,9 +334,9 @@ class CCompressedSparsePattern { template using CEdgeToNonZeroMap = C2DContainer; -using CCompressedSparsePatternUL = CCompressedSparsePattern; -using CCompressedSparsePatternL = CCompressedSparsePattern; -using CEdgeToNonZeroMapUL = CEdgeToNonZeroMap; +using CCompressedSparsePatternUL = CCompressedSparsePattern; +using CCompressedSparsePatternL = CCompressedSparsePattern; +using CEdgeToNonZeroMapUL = CEdgeToNonZeroMap; /*! * \brief Build a sparse pattern from geometry information, of type FVM or FEM, @@ -661,7 +662,7 @@ T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool includeOuterId /*! * \brief A way to represent one grid color that allows range-for syntax. */ -template +template struct GridColor { static_assert(std::is_integral::value); @@ -679,7 +680,7 @@ struct GridColor { * \brief A way to represent natural coloring {0,1,2,...,size-1} with zero * overhead (behaves like looping with an integer index, after optimization...). */ -template +template struct DummyGridColor { static_assert(std::is_integral::value); diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 2ae45cf340f..33cb5d0914b 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -7395,10 +7395,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { } } switch (Kind_Linear_Solver_Prec) { - case ILU: cout << "Using a ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break; - case LINELET: cout << "Using a linelet preconditioning."<< endl; break; - case LU_SGS: cout << "Using a LU-SGS preconditioning."<< endl; break; - case JACOBI: cout << "Using a Jacobi preconditioning."<< endl; break; + case ILU: cout << "Using ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break; + case LINELET: cout << "Using linelet preconditioning."<< endl; break; + case LU_SGS: cout << "Using LU-SGS preconditioning."<< endl; break; + case Q_LU_SGS: cout << "Using LU-SGS preconditioning with matrix quantization."<< endl; break; + case JACOBI: cout << "Using Jacobi preconditioning."<< endl; break; } break; case SMOOTHER: @@ -7406,6 +7407,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case ILU: cout << "A ILU(" << Linear_Solver_ILU_n << ")"; break; case LINELET: cout << "A Linelet"; break; case LU_SGS: cout << "A LU-SGS"; break; + case Q_LU_SGS: cout << "A quantized LU-SGS"; break; case JACOBI: cout << "A Jacobi"; break; } cout << " method is used for smoothing the linear system." << endl; diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index f6f423761e6..45153456086 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -4128,7 +4128,7 @@ const CGeometry::LDUSparsePattern& CGeometry::GetSparsePattern(ConnectivityType auto& grp = fillLvl == 0 ? (fvm ? finiteVolumePatternFill0 : finiteElementPatternFill0) : (fvm ? finiteVolumePatternFillN : finiteElementPatternFillN); if (grp.empty()) { - grp.csr = buildCSRPattern(*this, type, fillLvl); + grp.csr = buildCSRPattern(*this, type, static_cast(fillLvl)); grp.csr.buildDiagPtr(); grp.l = buildLowerPattern(grp.csr); grp.u = buildUpperPattern(grp.csr); @@ -4136,7 +4136,7 @@ const CGeometry::LDUSparsePattern& CGeometry::GetSparsePattern(ConnectivityType return grp; } -const su2vector& CGeometry::GetLToUTransposeSparsePatternMap(ConnectivityType type) { +const su2vector& CGeometry::GetLToUTransposeSparsePatternMap(ConnectivityType type) { bool fvm = (type == ConnectivityType::FiniteVolume); auto& l_to_u = fvm ? finiteVolumeLToUTranspMap : finiteElementLToUTranspMap; if (l_to_u.empty()) { @@ -4147,7 +4147,7 @@ const su2vector& CGeometry::GetLToUTransposeSparsePatternMap(Conn return l_to_u; } -const su2vector& CGeometry::GetUToLTransposeSparsePatternMap(ConnectivityType type) { +const su2vector& CGeometry::GetUToLTransposeSparsePatternMap(ConnectivityType type) { bool fvm = (type == ConnectivityType::FiniteVolume); auto& u_to_l = fvm ? finiteVolumeUToLTranspMap : finiteElementUToLTranspMap; if (u_to_l.empty()) { diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index f02d908a059..fe72110a411 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -31,6 +31,8 @@ #include "../../include/toolboxes/allocation_toolbox.hpp" #include +#include +#include #include namespace { @@ -83,6 +85,13 @@ CSysMatrix::CSysMatrix() : rank(SU2_MPI::GetRank()), size(SU2_MPI::G ilu.d = nullptr; ilu.u = nullptr; + q_scale_l = nullptr; + q_blocks_l = nullptr; + q_scale_u = nullptr; + q_blocks_u = nullptr; + q_scale_d = nullptr; + q_blocks_d = nullptr; + invM = nullptr; #ifdef USE_MKL @@ -105,6 +114,12 @@ CSysMatrix::~CSysMatrix() { MemoryAllocation::aligned_free(mat.l); MemoryAllocation::aligned_free(mat.u); MemoryAllocation::aligned_free(invM); + MemoryAllocation::aligned_free(q_scale_l); + MemoryAllocation::aligned_free(q_blocks_l); + MemoryAllocation::aligned_free(q_scale_u); + MemoryAllocation::aligned_free(q_blocks_u); + MemoryAllocation::aligned_free(q_scale_d); + MemoryAllocation::aligned_free(q_blocks_d); if (useCuda) { GPUMemoryAllocation::gpu_free(gpu.d); @@ -127,7 +142,7 @@ CSysMatrix::~CSysMatrix() { template void CSysMatrix::Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry* geometry, - const CConfig* config, bool needTranspPtr, bool grad_mode) { + const CConfig* config, bool needTranspPtr, bool grad_mode, bool allow_quant) { SU2_ZONE_SCOPED assert(omp_get_thread_num() == 0 && "Only the master thread is allowed to initialize the matrix."); @@ -161,8 +176,16 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi prec = config->GetKind_Grad_Linear_Solver_Prec(); } + useCuda = config->GetCUDA(); + const bool ilu_needed = (prec == ILU); const bool diag_needed = (prec == JACOBI) || (prec == LINELET); +#ifndef CODI_REVERSE_TYPE + const bool q_lus_needed = allow_quant && !useCuda && (prec == Q_LU_SGS); +#else + /*--- No quantization in adjoint mode for now because TransposeInPlace would get complicated. ---*/ + const bool q_lus_needed = false; +#endif /*--- Basic dimensions. ---*/ nVar = nvar; @@ -175,8 +198,6 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi ptr = MemoryAllocation::aligned_alloc(64, num * sizeof(ScalarType)); }; - useCuda = config->GetCUDA(); - /*--- L/D/U index structures and value arrays. ---*/ { const auto& pat = geometry->GetSparsePattern(type, 0); @@ -188,15 +209,33 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi mat.nnz_u = pat.u.getNumNonZeros(); } allocAndInit(mat.d, nPoint * nVar * nEqn); - allocAndInit(mat.l, mat.nnz_l * nVar * nEqn); - allocAndInit(mat.u, mat.nnz_u * nVar * nEqn); + + if (q_lus_needed) { + /*--- Q_LU_SGS: no full-precision L/U; off-diagonal blocks live in quantized storage. + * L/U are quantized on-the-fly during assembly; diagonal is quantized in Build step. ---*/ +#ifndef CODI_REVERSE_TYPE + quantized_mode = true; +#endif + auto allocQ = [](QuantType*& ptr, unsigned long n) { + ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); + }; + allocQ(q_scale_l, mat.nnz_l * nVar); + allocQ(q_blocks_l, mat.nnz_l * nVar * nEqn); + allocQ(q_scale_u, mat.nnz_u * nVar); + allocQ(q_blocks_u, mat.nnz_u * nVar * nEqn); + allocQ(q_scale_d, nPoint * nVar); + allocQ(q_blocks_d, nPoint * nVar * nEqn); + } else { + allocAndInit(mat.l, mat.nnz_l * nVar * nEqn); + allocAndInit(mat.u, mat.nnz_u * nVar * nEqn); + } if (useCuda) { auto GPUAllocAndInit = [](ScalarType*& ptr, unsigned long num) { ptr = GPUMemoryAllocation::gpu_alloc(num * sizeof(ScalarType)); }; - auto GPUAllocAndCopy = [](const unsigned long*& ptr, const unsigned long* src_ptr, unsigned long num) { - ptr = GPUMemoryAllocation::gpu_alloc_cpy(src_ptr, num * sizeof(unsigned long)); + auto GPUAllocAndCopy = [](const su2uint*& ptr, const su2uint* src_ptr, unsigned long num) { + ptr = GPUMemoryAllocation::gpu_alloc_cpy(src_ptr, num * sizeof(su2uint)); }; GPUAllocAndInit(gpu.d, nPoint * nVar * nEqn); GPUAllocAndInit(gpu.l, mat.nnz_l * nVar * nEqn); @@ -249,7 +288,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi /*--- Set suitable chunk sizes for light static for loops, and heavy dynamic ones, such that threads are approximately evenly loaded. ---*/ - omp_light_size = computeStaticChunkSize((mat.nnz_l + mat.nnz_u + nPoint) * nVar * nEqn, num_threads, OMP_MAX_SIZE_L); + omp_light_size = computeStaticChunkSize(nPoint * nVar * nEqn, num_threads, OMP_MAX_SIZE_L); omp_heavy_size = computeStaticChunkSize(nPointDomain, num_threads, OMP_MAX_SIZE_H); omp_num_parts = config->GetLinear_Solver_Prec_Threads(); @@ -538,21 +577,47 @@ void CSysMatrixComms::Complete(CSysVector& x, CGeometry* geometry, const CCon #endif } +template +void CSysMatrix::QuantizeBlock(const ScalarType* blk, QuantType* qs, QuantType* qv) const { + EncodeQuantBlock([&](unsigned long r, unsigned long c) { return blk[r * nVar + c]; }, qs, qv, nVar); +} + +template +void CSysMatrix::QuantizeDiagonalBlocks() { + SU2_ZONE_SCOPED + + if (quantized_mode) { + /*--- Q_LU_SGS: L/U were quantized during assembly; only the diagonal needs quantization now. ---*/ + SU2_OMP_FOR_DYN(omp_heavy_size) + for (auto i = 0ul; i < nPointDomain; ++i) + QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale_d[i * nVar], &q_blocks_d[i * nVar * nVar]); + END_SU2_OMP_FOR + } +} + template void CSysMatrix::SetValZero() { SU2_ZONE_SCOPED const auto nThreads = static_cast(omp_get_num_threads()); const auto iThread = static_cast(omp_get_thread_num()); - auto zeroChunk = [&](ScalarType* arr, unsigned long n) { + + auto zeroChunk = [&](auto* arr, unsigned long n) { if (n == 0) return; const auto chunk = roundUpDiv(n, nThreads); const auto begin = min(chunk * iThread, n); - const auto mySize = min(chunk, n - begin) * sizeof(ScalarType); + const auto mySize = min(chunk, n - begin) * sizeof(std::remove_pointer_t); if (mySize) memset(&arr[begin], 0, mySize); }; zeroChunk(mat.d, nPoint * nVar * nEqn); - zeroChunk(mat.l, mat.nnz_l * nVar * nEqn); - zeroChunk(mat.u, mat.nnz_u * nVar * nEqn); + if (!quantized_mode) { + zeroChunk(mat.l, mat.nnz_l * nVar * nEqn); + zeroChunk(mat.u, mat.nnz_u * nVar * nEqn); + } else { + zeroChunk(q_scale_l, mat.nnz_l * nVar); + zeroChunk(q_scale_u, mat.nnz_l * nVar); + zeroChunk(q_blocks_l, mat.nnz_l * nVar * nEqn); + zeroChunk(q_blocks_u, mat.nnz_u * nVar * nEqn); + } SU2_OMP_BARRIER } @@ -565,7 +630,7 @@ void CSysMatrix::SetValDiagonalZero() { } template -void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* vec) const { +void CSysMatrix::GaussElimination(ScalarType* matrix, ScalarType* vec) const { #ifdef USE_MKL_LAPACK // With MKL_DIRECT_CALL enabled, this is significantly faster than native code on Intel Architectures. lapack_int ipiv[MAXNVAR]; @@ -585,7 +650,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v 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, "DEBUG Gauss_Elimination"); + RegularizePivot(A(jVar, jVar), jVar, jVar, "DEBUG GaussElimination"); ScalarType weight = A(iVar, jVar) / A(jVar, jVar); @@ -600,7 +665,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v for (auto jVar = iVar + 1; jVar < nVar; jVar++) vec[iVar] -= A(iVar, jVar) * vec[jVar]; /*--- Regularize diagonal if too small ---*/ - RegularizePivot(A(iVar, iVar), iVar, iVar, "DEBUG Gauss_Elimination backsubst"); + RegularizePivot(A(iVar, iVar), iVar, iVar, "DEBUG GaussElimination backsubst"); vec[iVar] /= A(iVar, iVar); } @@ -611,7 +676,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v template void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inverse) const { /*--- This is a generalization of Gaussian elimination for multiple rhs' (the basis vectors). - We could call "Gauss_Elimination" multiple times or fully generalize it for multiple rhs, + We could call "GaussElimination" multiple times or fully generalize it for multiple rhs, the performance of both routines would suffer in both cases without the use of exotic templating. And so it feels reasonable to have some duplication here. ---*/ @@ -673,13 +738,27 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver template void CSysMatrix::DeleteValsRowi(unsigned long block_i, unsigned long row) { SU2_ZONE_SCOPED - for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) - for (auto iVar = 0u; iVar < nVar; iVar++) mat.l[k * nVar * nEqn + row * nEqn + iVar] = 0.0; - auto* d = &mat.d[block_i * nVar * nEqn]; - for (auto iVar = 0u; iVar < nVar; iVar++) d[row * nEqn + iVar] = 0.0; + const auto blkSz = nVar * nEqn; + + auto* d = &mat.d[block_i * blkSz]; + for (auto iVar = 0u; iVar < nEqn; iVar++) d[row * nEqn + iVar] = 0.0; d[row * nEqn + row] = 1.0; - for (auto k = mat.row_ptr_u[block_i]; k < mat.row_ptr_u[block_i + 1]; ++k) - for (auto iVar = 0u; iVar < nVar; iVar++) mat.u[k * nVar * nEqn + row * nEqn + iVar] = 0.0; + + if (quantized_mode) { + for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) { + for (auto iVar = 0u; iVar < nEqn; iVar++) q_blocks_l[k * blkSz + row * nEqn + iVar] = 0; + } + for (auto k = mat.row_ptr_u[block_i]; k < mat.row_ptr_u[block_i + 1]; ++k) { + for (auto iVar = 0u; iVar < nEqn; iVar++) q_blocks_u[k * blkSz + row * nEqn + iVar] = 0; + } + } else { + for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) { + for (auto iVar = 0u; iVar < nEqn; iVar++) mat.l[k * blkSz + row * nEqn + iVar] = 0; + } + for (auto k = mat.row_ptr_u[block_i]; k < mat.row_ptr_u[block_i + 1]; ++k) { + for (auto iVar = 0u; iVar < nEqn; iVar++) mat.u[k * blkSz + row * nEqn + iVar] = 0; + } + } } template @@ -702,11 +781,19 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector& v SU2_OMP_BARRIER - SU2_OMP_FOR_DYN(omp_heavy_size) - for (auto row_i = 0ul; row_i < nPointDomain; row_i++) { - RowProduct(vec, row_i, &prod[row_i * nVar]); + if (quantized_mode) { + SU2_OMP_FOR_DYN(omp_heavy_size) + for (auto row_i = 0ul; row_i < nPointDomain; row_i++) { + QuantizedRowProduct(vec, row_i, &prod[row_i * nVar]); + } + END_SU2_OMP_FOR + } else { + SU2_OMP_FOR_DYN(omp_heavy_size) + for (auto row_i = 0ul; row_i < nPointDomain; row_i++) { + RowProduct(vec, row_i, &prod[row_i * nVar]); + } + END_SU2_OMP_FOR } - END_SU2_OMP_FOR /*--- MPI Parallelization. ---*/ @@ -753,7 +840,7 @@ void CSysMatrix::BuildILUPreconditioner() { if (ilu_fill_in == 0) { /*--- ILU0: Same sparse pattern, copy L and U blocks directly. ---*/ - auto copy = [&](const unsigned long* row_ptr, const ScalarType* mat, ScalarType* ilu) { + auto copy = [&](const su2uint* row_ptr, const ScalarType* mat, ScalarType* ilu) { const unsigned long begin = row_ptr[iPoint] * blockSize; const unsigned long end = row_ptr[iPoint + 1] * blockSize; SU2_OMP_SIMD @@ -764,9 +851,8 @@ void CSysMatrix::BuildILUPreconditioner() { return; } /*--- ILUn: Merge-scan L and U via shared lambda. ---*/ - auto scatterPart = [&](const unsigned long* mat_row_ptr, const unsigned long* mat_col_ind, - const ScalarType* mat_vals, const unsigned long* ilu_row_ptr, - const unsigned long* ilu_col_ind, ScalarType* ilu_vals) { + auto scatterPart = [&](const su2uint* mat_row_ptr, const su2uint* mat_col_ind, const ScalarType* mat_vals, + const su2uint* ilu_row_ptr, const su2uint* ilu_col_ind, ScalarType* ilu_vals) { auto km = mat_row_ptr[iPoint], km_end = mat_row_ptr[iPoint + 1]; for (auto k = ilu_row_ptr[iPoint]; k < ilu_row_ptr[iPoint + 1]; ++k) { const auto jPoint = ilu_col_ind[k]; @@ -976,11 +1062,20 @@ void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector::ComputeLU_SGSPreconditioner(const CSysVector begin;) { - iPoint--; // because of unsigned type - auto idx = iPoint * nVar; - DiagonalProduct(prod, iPoint, dia_prod); // Compute D.x* - UpperProduct(prod, iPoint, row_end, up_prod); // Compute U.x_(n+1) - VectorSubtraction(dia_prod, up_prod, &prod[idx]); // Compute y = D.x*-U.x_(n+1) - Gauss_Elimination(iPoint, &prod[idx]); // Solve D.x* = y + if (quantized_mode) { + for (auto iPoint = row_end; iPoint > begin;) { + iPoint--; + auto idx = iPoint * nVar; + QuantizedDiagonalProduct(prod, iPoint, dia_prod); + QuantizedUpperProduct(prod, iPoint, row_end, up_prod); + VectorSubtraction(dia_prod, up_prod, &prod[idx]); + QuantizedGaussElimination(iPoint, &prod[idx]); + } + } else { + for (auto iPoint = row_end; iPoint > begin;) { + iPoint--; // because of unsigned type + auto idx = iPoint * nVar; + DiagonalProduct(prod, iPoint, dia_prod); // Compute D.x* + UpperProduct(prod, iPoint, row_end, up_prod); // Compute U.x_(n+1) + VectorSubtraction(dia_prod, up_prod, &prod[idx]); // Compute y = D.x*-U.x_(n+1) + GaussElimination(iPoint, &prod[idx]); // Solve D.x* = y + } } } END_SU2_OMP_FOR @@ -1129,7 +1235,7 @@ void CSysMatrix::ComputeLineletPreconditioner(const CSysVector 0; --iElem) { @@ -1244,19 +1350,32 @@ void CSysMatrix::EnforceZeroProjection(unsigned long node_i, const O template void CSysMatrix::SetDiagonalAsColumnSum() { SU2_ZONE_SCOPED + const auto blkSz = nVar * nEqn; SU2_OMP_FOR_DYN(omp_heavy_size) for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { - auto* d_i = &mat.d[iPoint * nVar * nEqn]; - for (auto k = 0ul; k < nVar * nEqn; ++k) d_i[k] = 0.0; - - /*--- For each L entry (iPoint, j): subtract its U-transpose (j, iPoint). ---*/ - for (auto k_l = mat.row_ptr_l[iPoint]; k_l < mat.row_ptr_l[iPoint + 1]; ++k_l) - MatrixSubtraction(d_i, &mat.u[l_to_u_transp[k_l] * nVar * nEqn], d_i); - - /*--- For each U entry (iPoint, j): subtract its L-transpose (j, iPoint). ---*/ - for (auto k_u = mat.row_ptr_u[iPoint]; k_u < mat.row_ptr_u[iPoint + 1]; ++k_u) - MatrixSubtraction(d_i, &mat.l[u_to_l_transp[k_u] * nVar * nEqn], d_i); + auto* d_i = &mat.d[iPoint * blkSz]; + for (auto k = 0ul; k < blkSz; ++k) d_i[k] = 0.0; + + if (!quantized_mode) { + /*--- For each L entry (iPoint, j): subtract its U-transpose (j, iPoint). ---*/ + for (auto k_l = mat.row_ptr_l[iPoint]; k_l < mat.row_ptr_l[iPoint + 1]; ++k_l) + MatrixSubtraction(d_i, &mat.u[l_to_u_transp[k_l] * blkSz], d_i); + + /*--- For each U entry (iPoint, j): subtract its L-transpose (j, iPoint). ---*/ + for (auto k_u = mat.row_ptr_u[iPoint]; k_u < mat.row_ptr_u[iPoint + 1]; ++k_u) + MatrixSubtraction(d_i, &mat.l[u_to_l_transp[k_u] * blkSz], d_i); + } else { + auto subtractTransp = [&](su2uint k_transp, const QuantType* qs, const QuantType* qv) { + const CBlockView view{nullptr, &qs[k_transp * nVar], &qv[k_transp * blkSz], nVar}; + for (auto i = 0ul; i < nVar; ++i) + for (auto j = 0ul; j < nEqn; ++j) d_i[i * nEqn + j] -= view(i, j); + }; + for (auto k_l = mat.row_ptr_l[iPoint]; k_l < mat.row_ptr_l[iPoint + 1]; ++k_l) + subtractTransp(l_to_u_transp[k_l], q_scale_u, q_blocks_u); + for (auto k_u = mat.row_ptr_u[iPoint]; k_u < mat.row_ptr_u[iPoint + 1]; ++k_u) + subtractTransp(u_to_l_transp[k_u], q_scale_l, q_blocks_l); + } } END_SU2_OMP_FOR } diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index aae5d9ce707..a50ed059625 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1569,7 +1569,8 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con if (RequiresTranspose) Jacobian.BuildJacobiPreconditioner(); break; case LU_SGS: - /*--- Nothing to build. ---*/ + case Q_LU_SGS: + /*--- Nothing to build (transpose path not supported for Q_LU_SGS, see CSysMatrix::Initialize). ---*/ break; case PASTIX_ILU: case PASTIX_LU_P: diff --git a/SU2_CFD/include/numerics_simd/util.hpp b/SU2_CFD/include/numerics_simd/util.hpp index 6e84209ee6d..c5fc54c25b3 100644 --- a/SU2_CFD/include/numerics_simd/util.hpp +++ b/SU2_CFD/include/numerics_simd/util.hpp @@ -248,7 +248,7 @@ FORCEINLINE void updateLinearSystem(Int iEdge, vector.UpdateBlocks(iPoint, jPoint, flux, updateMask); if(implicit) { auto wasActive = AD::BeginPassive(); - matrix.UpdateBlocks(iEdge, iPoint, jPoint, jac_i, jac_j, updateMask); + matrix.SetBlocks(iEdge, iPoint, jPoint, jac_i, jac_j, updateMask); AD::EndPassive(wasActive); } } diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 93529597cdc..b2dd6441318 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -327,20 +327,58 @@ class CFVMFlowSolverBase : public CSolver { /*! * \brief Compute the viscous contribution for a particular edge. - * \note The convective residual methods include a call to this for each edge, - * this allows convective and viscous loops to be "fused". + * \note The convective residual methods include a call to this for each edge, this allows convective and + * viscous loops to be "fused". Only the residual is applied here, the Jacobians are returned so that + * the caller can update the system matrix in a single operation together with the convective part + * (a requirement of quantized matrix storage). * \param[in] iEdge - Edge for which the flux and Jacobians are to be computed. * \param[in] geometry - Geometrical definition of the problem. * \param[in] solver_container - Container vector with all the solutions. * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. + * \return The viscous Jacobians (null for inviscid solvers). */ - inline virtual void Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config) { } - void Viscous_Residual_impl(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config); + inline virtual CNumerics::ResidualType<> Viscous_Residual(unsigned long iEdge, CGeometry *geometry, + CSolver **solver_container, CNumerics *numerics, + CConfig *config) { + return CNumerics::ResidualType<>(nullptr, nullptr, nullptr); + } + CNumerics::ResidualType<> Viscous_Residual_impl(unsigned long iEdge, CGeometry *geometry, + CSolver **solver_container, CNumerics *numerics, CConfig *config); using CSolver::Viscous_Residual; /*--- Silence warning ---*/ + /*! + * \brief Update the Jacobian for one edge with the fused convective and viscous contributions. + * \note Both contributions must be applied at once because in quantized mode the + * off-diagonal blocks of the matrix can only be overwritten, not accumulated. + * \param[in] iEdge - Edge index for the off-diagonal blocks. + * \param[in] iPoint, jPoint - Points connected by the edge (diagonal blocks). + * \param[in] conv - Convective residual/Jacobians (added to i, subtracted from j). + * \param[in] visc - Viscous residual/Jacobians (subtracted from i, added to j), may hold null Jacobians. + */ + inline void UpdateJacobian(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, + const CNumerics::ResidualType<>& conv, const CNumerics::ResidualType<>& visc) { + /*--- Lazy element-wise difference, presented with the [i][j] access the matrix expects. ---*/ + struct CJacobianDifference { + const su2double* const* conv; + const su2double* const* visc; + struct Row { + const su2double *c, *v; + su2double operator[](unsigned long j) const { return c[j] - v[j]; } + }; + Row operator[](unsigned long i) const { return {conv[i], visc[i]}; } + }; + if (visc.jacobian_i != nullptr) { + const CJacobianDifference jac_i{conv.jacobian_i, visc.jacobian_i}; + const CJacobianDifference jac_j{conv.jacobian_j, visc.jacobian_j}; + if (ReducerStrategy) Jacobian.SetBlocks(iEdge, jac_i, jac_j); + else Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, jac_i, jac_j); + } else { + if (ReducerStrategy) Jacobian.SetBlocks(iEdge, conv.jacobian_i, conv.jacobian_j); + else Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, conv.jacobian_i, conv.jacobian_j); + } + } + /*! * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over a nonlinear * iteration for stability. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 28ca74d8149..080704d698d 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -317,7 +317,7 @@ void CFVMFlowSolverBase::HybridParallelInitialization(const CConfig& confi if (!coloring.empty()) { /*--- If the reducer strategy is used we are not constrained by group * size as we have no other edge loops in the Euler/NS solvers. ---*/ - auto groupSize = ReducerStrategy ? 1ul : geometry.GetEdgeColorGroupSize(); + auto groupSize = static_cast(ReducerStrategy ? 1ul : geometry.GetEdgeColorGroupSize()); auto nColor = coloring.getOuterSize(); EdgeColoring.reserve(nColor); @@ -445,11 +445,11 @@ void CFVMFlowSolverBase::SetPrimitive_Limiter(CGeometry* geometry, const C } template -void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config) { +CNumerics::ResidualType<> CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeometry *geometry, + CSolver **solver_container, + CNumerics *numerics, CConfig *config) { SU2_ZONE_SCOPED - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); const bool backscatter = config->GetSBSParam().StochasticBackscatter; const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || @@ -518,17 +518,14 @@ void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeome if (ReducerStrategy) { EdgeFluxes.SubtractBlock(iEdge, residual); - if (implicit) - Jacobian.UpdateBlocksSub(iEdge, residual.jacobian_i, residual.jacobian_j); } else { LinSysRes.SubtractBlock(iPoint, residual); LinSysRes.AddBlock(jPoint, residual); - - if (implicit) - Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } + /*--- The Jacobians are applied by the caller, fused with the convective contribution. ---*/ + return residual; } template @@ -1296,13 +1293,14 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve auto ModifyJacobian = [&](const unsigned long jPoint) { su2double jac[MAXNVAR * MAXNVAR], newJac[MAXNVAR * MAXNVAR]; - auto* block = Jacobian.GetBlock(iPoint, jPoint); - for (auto iVar = 0u; iVar < nVar * nVar; iVar++) jac[iVar] = block[iVar]; + const auto view = Jacobian.GetBlockView(iPoint, jPoint); + if (!view) return; + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto jVar = 0u; jVar < nVar; jVar++) jac[iVar * nVar + jVar] = view(iVar, jVar); CBlasStructure().gemm(nVar, nVar, nVar, mat, jac, newJac, config); - for (auto iVar = 0u; iVar < nVar * nVar; iVar++) - block[iVar] = SU2_TYPE::GetValue(newJac[iVar]); + Jacobian.SetBlock(iPoint, jPoint, newJac); }; ModifyJacobian(iPoint); for (size_t iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { diff --git a/SU2_CFD/include/solvers/CIncNSSolver.hpp b/SU2_CFD/include/solvers/CIncNSSolver.hpp index dc2d80e86e6..712855d0a20 100644 --- a/SU2_CFD/include/solvers/CIncNSSolver.hpp +++ b/SU2_CFD/include/solvers/CIncNSSolver.hpp @@ -61,9 +61,10 @@ class CIncNSSolver final : public CIncEulerSolver { * \param[in] solver_container - Container vector with all the solutions. * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. + * \return The viscous Jacobians, to be applied by the caller together with the convective part. */ - void Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config) override; + CNumerics::ResidualType<> Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, + CNumerics *numerics, CConfig *config) override; /*! * \brief Computes the wall shear stress (Tau_Wall) on the surface using a wall function. diff --git a/SU2_CFD/include/solvers/CNSSolver.hpp b/SU2_CFD/include/solvers/CNSSolver.hpp index 4ac0ee4ae7d..405ff6de9a7 100644 --- a/SU2_CFD/include/solvers/CNSSolver.hpp +++ b/SU2_CFD/include/solvers/CNSSolver.hpp @@ -111,9 +111,10 @@ class CNSSolver final : public CEulerSolver { * \param[in] solver_container - Container vector with all the solutions. * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. + * \return The viscous Jacobians, to be applied by the caller together with the convective part. */ - void Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config) override; + CNumerics::ResidualType<> Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, + CNumerics *numerics, CConfig *config) override; /*! * \brief Computes the wall shear stress (Tau_Wall) on the surface using a wall function. diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 7d6871f8ca4..26ced4a77a8 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -63,7 +63,7 @@ CScalarSolver::CScalarSolver(CGeometry* geometry, CConfig* config, if (ReducerStrategy && (coloring.getOuterSize() > 1)) geometry->SetNaturalEdgeColoring(); if (!coloring.empty()) { - auto groupSize = ReducerStrategy ? 1ul : geometry->GetEdgeColorGroupSize(); + auto groupSize = static_cast(ReducerStrategy ? 1ul : geometry->GetEdgeColorGroupSize()); auto nColor = coloring.getOuterSize(); EdgeColoring.reserve(nColor); @@ -105,7 +105,7 @@ void CScalarSolver::CommonPreprocessing(CGeometry *geometry, const if (!ReducerStrategy && !Output) { LinSysRes.SetValZero(); if (implicit) { - Jacobian.SetValZero(); + Jacobian.SetValDiagonalZero(); } else { SU2_OMP_BARRIER } @@ -291,7 +291,7 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** } else { LinSysRes.AddBlock(iPoint, residual); LinSysRes.SubtractBlock(jPoint, residual); - if (implicit) Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + if (implicit) Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } /*--- Apply convective flux correction to negate the effects of flow divergence in case of incompressible flow. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 3fc3c2a7bc9..942e2e25877 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3492,7 +3492,7 @@ class CSolver { * \param[in] rhs - Right hand side. * \param[in] nVar - Number of variables. */ - void Gauss_Elimination(su2double** A, + void GaussElimination(su2double** A, su2double* rhs, unsigned short nVar); diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index 361df6d4c6b..88ee3214a63 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -2094,7 +2094,7 @@ void CAdjEulerSolver::Smooth_Sensitivity(CGeometry *geometry, CSolver **solver_c A[iVertex][iVertex+1] = 0.0; A[iVertex][iVertex-1] = 0.0; - Gauss_Elimination(A, b, (unsigned short)nVertex); + GaussElimination(A, b, (unsigned short)nVertex); /*--- Set the new value of the sensitiviy ---*/ diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index e89f022920c..d02385633d2 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -164,7 +164,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, false, true); } else { if (rank == MASTER_NODE) @@ -1983,7 +1983,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Compute the residual ---*/ - auto residual = numerics->ComputeResidual(config); + auto conv_residual = numerics->ComputeResidual(config); /*--- Set the final value of the Roe dissipation coefficient ---*/ @@ -1995,23 +1995,18 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Update residual value ---*/ if (ReducerStrategy) { - EdgeFluxes.SetBlock(iEdge, residual); - if (implicit) - Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + EdgeFluxes.SetBlock(iEdge, conv_residual); + } else { + LinSysRes.AddBlock(iPoint, conv_residual); + LinSysRes.SubtractBlock(jPoint, conv_residual); } - else { - LinSysRes.AddBlock(iPoint, residual); - LinSysRes.SubtractBlock(jPoint, residual); - /*--- Set implicit computation ---*/ - if (implicit) - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); - } + /*--- Viscous contribution, returns its Jacobians so that the matrix is updated once. ---*/ - /*--- Viscous contribution. ---*/ + const auto visc_residual = Viscous_Residual( + iEdge, geometry, solver_container, numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); - Viscous_Residual(iEdge, geometry, solver_container, - numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); + if (implicit) UpdateJacobian(iEdge, iPoint, jPoint, conv_residual, visc_residual); } END_SU2_OMP_FOR } // end color loop diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 31cb41c52d9..27ad9c8fdff 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -258,7 +258,7 @@ void CFEASolver::HybridParallelInitialization(CGeometry* geometry) { if (!coloring.empty()) { /*--- We are not constrained by the color group size when using locks. ---*/ - auto groupSize = LockStrategy? 1ul : geometry->GetElementColorGroupSize(); + auto groupSize = static_cast(LockStrategy ? 1ul : geometry->GetElementColorGroupSize()); auto nColor = coloring.getOuterSize(); ElemColoring.reserve(nColor); diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index 52b591770ca..2a9e4790760 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -180,13 +180,11 @@ void CHeatSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, /*--- Need to clear EdgeFluxes and Jacobian when only the viscous part is called for solid heat transfer, * for the weakly coupled energy equation the convection part does this by setting instead of incrementing. ---*/ - if (!Output && !flow && ReducerStrategy) { - EdgeFluxes.SetValZero(); - if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) { - Jacobian.SetValZero(); - } else { - SU2_OMP_BARRIER - } + if (!Output && !flow) { + if (ReducerStrategy) EdgeFluxes.SetValZero(); + if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) Jacobian.SetValZero(); + + SU2_OMP_BARRIER } } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 2c4c9b05e1c..4951551851f 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -155,7 +155,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, false, true); } else { if (rank == MASTER_NODE) @@ -1013,7 +1013,7 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ if(!ReducerStrategy && !Output) { LinSysRes.SetValZero(); - if (implicit) Jacobian.SetValZero(); + if (implicit) Jacobian.SetValDiagonalZero(); else {SU2_OMP_BARRIER} // because of "nowait" in LinSysRes } } @@ -1199,30 +1199,25 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co /*--- Compute residuals, and Jacobians ---*/ - auto residual = numerics->ComputeResidual(config); + auto conv_residual = numerics->ComputeResidual(config); - if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; + if (bounded_scalar) EdgeMassFluxes[iEdge] = conv_residual[0]; /*--- Update residual value ---*/ if (ReducerStrategy) { - EdgeFluxes.SetBlock(iEdge, residual); - if (implicit) - Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + EdgeFluxes.SetBlock(iEdge, conv_residual); + } else { + LinSysRes.AddBlock(iPoint, conv_residual); + LinSysRes.SubtractBlock(jPoint, conv_residual); } - else { - LinSysRes.AddBlock(iPoint, residual); - LinSysRes.SubtractBlock(jPoint, residual); - /*--- Set implicit computation ---*/ - if (implicit) - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); - } + /*--- Viscous contribution, returns its Jacobians so that the matrix is updated once. ---*/ - /*--- Viscous contribution. ---*/ + const auto visc_residual = Viscous_Residual( + iEdge, geometry, solver_container, numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); - Viscous_Residual(iEdge, geometry, solver_container, - numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); + if (implicit) UpdateJacobian(iEdge, iPoint, jPoint, conv_residual, visc_residual); } END_SU2_OMP_FOR } // end color loop @@ -1377,30 +1372,26 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont /*--- Compute the residual ---*/ - auto residual = numerics->ComputeResidual(config); + auto conv_residual = numerics->ComputeResidual(config); - if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; + if (bounded_scalar) EdgeMassFluxes[iEdge] = conv_residual[0]; /*--- Update residual value ---*/ if (ReducerStrategy) { - EdgeFluxes.SetBlock(iEdge, residual); - if (implicit) - Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + EdgeFluxes.SetBlock(iEdge, conv_residual); } else { - LinSysRes.AddBlock(iPoint, residual); - LinSysRes.SubtractBlock(jPoint, residual); - - /*--- Set implicit computation ---*/ - if (implicit) - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + LinSysRes.AddBlock(iPoint, conv_residual); + LinSysRes.SubtractBlock(jPoint, conv_residual); } - /*--- Viscous contribution. ---*/ + /*--- Viscous contribution, returns its Jacobians so that the matrix is updated once. ---*/ + + const auto visc_residual = Viscous_Residual( + iEdge, geometry, solver_container, numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); - Viscous_Residual(iEdge, geometry, solver_container, - numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); + if (implicit) UpdateJacobian(iEdge, iPoint, jPoint, conv_residual, visc_residual); } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 2d6dd26786d..2532215cf8d 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -275,8 +275,9 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, SU2_OMP_SAFE_GLOBAL_ACCESS(GetStreamwise_Periodic_Properties(geometry, config, iMesh);) } -void CIncNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config) { +CNumerics::ResidualType<> CIncNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, + CSolver **solver_container, CNumerics *numerics, + CConfig *config) { const bool energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); /*--- Contribution to heat flux due to enthalpy diffusion for multicomponent and reacting flows ---*/ @@ -286,7 +287,7 @@ void CIncNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CS Compute_Enthalpy_Diffusion(iEdge, geometry, solver_container, numerics, n_species, implicit); } - Viscous_Residual_impl(iEdge, geometry, solver_container, numerics, config); + return Viscous_Residual_impl(iEdge, geometry, solver_container, numerics, config); } void CIncNSSolver::Compute_Enthalpy_Diffusion(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 37de9a83207..93bc7037e1f 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -282,7 +282,7 @@ void CNEMOEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver if(!ReducerStrategy && !Output) { LinSysRes.SetValZero(); - if (implicit) Jacobian.SetValZero(); + if (implicit) Jacobian.SetValDiagonalZero(); } } @@ -460,7 +460,7 @@ void CNEMOEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_c LinSysRes.AddBlock(iPoint, residual); LinSysRes.SubtractBlock(jPoint, residual); if (implicit) { - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } } } @@ -625,7 +625,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con LinSysRes.AddBlock(iPoint, residual); LinSysRes.SubtractBlock(jPoint, residual); if (implicit) { - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } } } @@ -1894,7 +1894,7 @@ void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, // /*--- Jacobian contribution for implicit integration ---*/ // if (implicit) -// Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); +// Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); // } } } diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index b3578a0d554..bae5591d946 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -798,7 +798,7 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, Jacobian_i[iVar][jVar] += Jacobian_j[iVar][kVar]*dVdU[kVar][jVar]*Area; /*--- Apply to the linear system ---*/ - Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); + Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); } } else { diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 367d7db04bb..643cb9a4f6c 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -178,9 +178,10 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, cons return nonPhysicalPoints; } -void CNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, - CNumerics *numerics, CConfig *config) { - Viscous_Residual_impl(iEdge, geometry, solver_container, numerics, config); +CNumerics::ResidualType<> CNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, + CSolver **solver_container, CNumerics *numerics, + CConfig *config) { + return Viscous_Residual_impl(iEdge, geometry, solver_container, numerics, config); } void CNSSolver::Buffet_Monitoring(const CGeometry *geometry, const CConfig *config) { diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 72685eb9bc5..e014c794a52 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -27,6 +27,9 @@ #include "../../include/solvers/CSolver.hpp" + +#include + #include "../../include/gradients/computeGradientsGreenGauss.hpp" #include "../../include/gradients/computeGradientsLeastSquares.hpp" #include "../../include/limiters/computeLimiters.hpp" @@ -542,9 +545,11 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, if (implicit_periodic) { + const auto block = Jacobian.GetBlockView(iPoint, iPoint); + for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { - jacBlock[iVar][jVar] = Jacobian.GetBlock(iPoint, iPoint, iVar, jVar); + jacBlock[iVar][jVar] = block(iVar, jVar); } } @@ -553,21 +558,15 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, if (rotate_periodic) { for (iVar = 0; iVar < nVar; iVar++) { if (nDim == 2) { - jacBlock[1][iVar] = (rotMatrix2D[0][0]*Jacobian.GetBlock(iPoint, iPoint, 1, iVar) + - rotMatrix2D[0][1]*Jacobian.GetBlock(iPoint, iPoint, 2, iVar)); - jacBlock[2][iVar] = (rotMatrix2D[1][0]*Jacobian.GetBlock(iPoint, iPoint, 1, iVar) + - rotMatrix2D[1][1]*Jacobian.GetBlock(iPoint, iPoint, 2, iVar)); + jacBlock[1][iVar] = rotMatrix2D[0][0]*block(1, iVar) + rotMatrix2D[0][1]*block(2, iVar); + jacBlock[2][iVar] = rotMatrix2D[1][0]*block(1, iVar) + rotMatrix2D[1][1]*block(2, iVar); } else { - - jacBlock[1][iVar] = (rotMatrix3D[0][0]*Jacobian.GetBlock(iPoint, iPoint, 1, iVar) + - rotMatrix3D[0][1]*Jacobian.GetBlock(iPoint, iPoint, 2, iVar) + - rotMatrix3D[0][2]*Jacobian.GetBlock(iPoint, iPoint, 3, iVar)); - jacBlock[2][iVar] = (rotMatrix3D[1][0]*Jacobian.GetBlock(iPoint, iPoint, 1, iVar) + - rotMatrix3D[1][1]*Jacobian.GetBlock(iPoint, iPoint, 2, iVar) + - rotMatrix3D[1][2]*Jacobian.GetBlock(iPoint, iPoint, 3, iVar)); - jacBlock[3][iVar] = (rotMatrix3D[2][0]*Jacobian.GetBlock(iPoint, iPoint, 1, iVar) + - rotMatrix3D[2][1]*Jacobian.GetBlock(iPoint, iPoint, 2, iVar) + - rotMatrix3D[2][2]*Jacobian.GetBlock(iPoint, iPoint, 3, iVar)); + jacBlock[1][iVar] = rotMatrix3D[0][0]*block(1, iVar) + rotMatrix3D[0][1]*block(2, iVar) + + rotMatrix3D[0][2]*block(3, iVar); + jacBlock[2][iVar] = rotMatrix3D[1][0]*block(1, iVar) + rotMatrix3D[1][1]*block(2, iVar) + + rotMatrix3D[1][2]*block(3, iVar); + jacBlock[3][iVar] = rotMatrix3D[2][0]*block(1, iVar) + rotMatrix3D[2][1]*block(2, iVar) + + rotMatrix3D[2][2]*block(3, iVar); } } } @@ -1997,7 +1996,7 @@ void CSolver::SetResidual_RMS(const CGeometry *geometry, const CConfig *config) /*--- Set the L2 Norm residual in all the processors. ---*/ - vector rbuf_res(nVar); + vector rbuf_res(nVar * nDim); unsigned long Global_nPointDomain = 0; if (config->GetComm_Level() == COMM_FULL) { @@ -2028,21 +2027,22 @@ void CSolver::SetResidual_RMS(const CGeometry *geometry, const CConfig *config) /*--- Set the Maximum residual in all the processors. ---*/ if (config->GetComm_Level() == COMM_FULL) { - - const unsigned long nProcessor = size; - - su2activematrix rbuf_residual(nProcessor,nVar); - su2matrix rbuf_point(nProcessor,nVar); - su2activematrix rbuf_coord(nProcessor*nVar, nDim); - - SU2_MPI::Allgather(Residual_Max.data(), nVar, MPI_DOUBLE, rbuf_residual.data(), nVar, MPI_DOUBLE, SU2_MPI::GetComm()); - SU2_MPI::Allgather(Point_Max.data(), nVar, MPI_UNSIGNED_LONG, rbuf_point.data(), nVar, MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); - SU2_MPI::Allgather(Point_Max_Coord.data(), nVar*nDim, MPI_DOUBLE, rbuf_coord.data(), nVar*nDim, MPI_DOUBLE, SU2_MPI::GetComm()); - + SU2_MPI::Allreduce(Residual_Max.data(), rbuf_res.data(), nVar, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); for (unsigned short iVar = 0; iVar < nVar; iVar++) { - for (auto iProcessor = 0ul; iProcessor < nProcessor; iProcessor++) { - AddRes_Max(iVar, rbuf_residual(iProcessor,iVar), rbuf_point(iProcessor,iVar), rbuf_coord[iProcessor*nVar+iVar]); + if (Residual_Max[iVar] < rbuf_res[iVar]) { + Point_Max[iVar] = 0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Point_Max_Coord(iVar, iDim) = std::numeric_limits::lowest(); + } } + Residual_Max[iVar] = rbuf_res[iVar]; + } + vector rbuf_point(nVar); + SU2_MPI::Allreduce(Point_Max.data(), rbuf_point.data(), nVar, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(Point_Max_Coord.data(), rbuf_res.data(), nVar*nDim, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + Point_Max = std::move(rbuf_point); + for (unsigned short iVar = 0; iVar < nVar * nDim; iVar++) { + Point_Max_Coord.data()[iVar] = rbuf_res[iVar]; } } @@ -2331,7 +2331,7 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) { *geometry, *config, 0, nVar, umusclKappa, solution, gradient, solMin, solMax, limiter); } -void CSolver::Gauss_Elimination(su2double** A, su2double* rhs, unsigned short nVar) { +void CSolver::GaussElimination(su2double** A, su2double* rhs, unsigned short nVar) { SU2_ZONE_SCOPED short iVar, jVar, kVar; diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index e32cd1b071f..e108fb978dc 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -103,7 +103,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 100 - flatplate.test_vals = [-6.543288, -1.065157, 0.001196, 0.029390, 2.361500, -2.332100, 0.000000, 0.000000] + flatplate.test_vals = [-6.543299, -1.065161, 0.001196, 0.029390, 2.361500, -2.332100, 0.000000, 0.000000] test_list.append(flatplate) # Laminar cylinder (steady) @@ -119,7 +119,7 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.476618, -1.014896, 0.594402, 70.856965, 0.000000] + cylinder_lowmach.test_vals = [-6.476618, -1.014896, 0.594402, 70.856953, 0.000000] cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0] test_list.append(cylinder_lowmach) @@ -128,7 +128,7 @@ def main(): poiseuille.cfg_dir = "navierstokes/poiseuille" poiseuille.cfg_file = "lam_poiseuille.cfg" poiseuille.test_iter = 10 - poiseuille.test_vals = [-5.046182, 0.652932, 0.008488, 13.734488, 0.000000] + poiseuille.test_vals = [-7.595535, -6.993332, 0.007670, 13.591507, 0] test_list.append(poiseuille) # 2D Poiseuille flow (inlet profile file) @@ -136,7 +136,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.005131, -7.582325, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.005122, -7.582423, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262530, -0.000000, 2.089953] test_list.append(poiseuille_profile) @@ -145,7 +145,7 @@ def main(): periodic2d.cfg_dir = "navierstokes/periodic2D" periodic2d.cfg_file = "config.cfg" periodic2d.test_iter = 1400 - periodic2d.test_vals = [-10.817608, -8.363542, -8.287458, -5.334102, -1.088411, -2945.200000] + periodic2d.test_vals = [-10.817608, -8.363540, -8.287459, -5.334101, -1.088412, -2945.200000] test_list.append(periodic2d) ########################## @@ -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.038011, -16.332088, 1.080346, 0.018385, 20.000000, -2.873515, 0.000000, -14.250270, 0.000000] + turb_naca0012_sa.test_vals = [-12.038027, -16.332088, 1.080346, 0.018385, 20.000000, -2.873302, 0.000000, -14.250270, 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.093908, -15.250756, -5.906323, 1.070413, 0.015775, -2.855199, 0.000000] + turb_naca0012_sst.test_vals = [-12.093908, -15.250825, -5.906324, 1.070413, 0.015775, -2.856106, 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.080880, -14.837176, -5.732907, 1.000893, 0.019109, -2.119818] + turb_naca0012_sst_sust.test_vals = [-12.080858, -14.837175, -5.732907, 1.000893, 0.019109, -2.120045] 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) @@ -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.083072, -5.374686, -8.880093, -4.073514, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-11.083070, -5.374688, -8.880090, -4.073523, 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) @@ -367,7 +367,7 @@ def main(): inc_lam_cylinder.cfg_dir = "incomp_navierstokes/cylinder" inc_lam_cylinder.cfg_file = "incomp_cylinder.cfg" inc_lam_cylinder.test_iter = 10 - inc_lam_cylinder.test_vals = [-4.161196, -3.573053, 0.025533, 4.944647] + inc_lam_cylinder.test_vals = [-4.161196, -3.573053, 0.025533, 4.944646] test_list.append(inc_lam_cylinder) # Buoyancy-driven cavity @@ -404,7 +404,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.758114, -10.974548, -0.000004, -0.028637, 5.000000, -4.080548, 2.000000, -4.490129] + inc_turb_naca0012.test_vals = [-4.758114, -10.974548, -0.000004, -0.028637, 5.000000, -4.080155, 2.000000, -4.491724] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST @@ -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.106209, -16.303012, -16.484904, -15.006575, -17.858050, -14.024869, 5.609100] + inc_weakly_coupled.test_vals = [-18.106234, -16.302995, -16.484896, -15.006575, -17.858050, -14.024855, 5.609100] test_list.append(inc_weakly_coupled) ###################################### @@ -570,7 +570,7 @@ def main(): transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg" transonic_stator_restart.test_iter = 20 - transonic_stator_restart.test_vals = [-4.367141, -2.487739, -2.079071, 1.728176, -1.464952, 3.225028, -471620.000000, 94.839000, -0.051126] + transonic_stator_restart.test_vals = [-4.367141, -2.487739, -2.079071, 1.728176, -1.464951, 3.225028, -471620.000000, 94.839000, -0.051125] transonic_stator_restart.test_vals_aarch64 = [-4.442510, -2.561369, -2.165778, 1.652750, -1.355494, 3.172712, -471620.000000, 94.843000, -0.043825] test_list.append(transonic_stator_restart) @@ -592,7 +592,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 5 - uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624456] + uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624446] uniform_flow.unsteady = True uniform_flow.multizone = True test_list.append(uniform_flow) @@ -663,7 +663,7 @@ def main(): slinc_steady.cfg_dir = "sliding_interface/incompressible_steady" slinc_steady.cfg_file = "config.cfg" slinc_steady.test_iter = 19 - slinc_steady.test_vals = [19.000000, -1.144102, -1.424987] + slinc_steady.test_vals = [19.000000, -1.144102, -1.424986] slinc_steady.test_vals_aarch64 = [19.000000, -1.154874, -1.378120] slinc_steady.multizone = True test_list.append(slinc_steady) @@ -746,7 +746,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.631546, 0.000000, 0.000000] + mms_fvm_inc_ns.test_vals = [-7.414945, -7.631547, 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 a6c8c03b01b..4db36093031 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -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.665653, -6.239112, 0.000000] + discadj_incomp_cylinder.test_vals = [20.000000, -1.665652, -6.239091, 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.023526, 0.000000] + discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.845989, -1.023525, 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.775278, -3.089107, -7.143663, 0.000000, -0.896760] + discadj_incomp_turb_NACA0012_sst.test_vals = [-3.775275, -3.089116, -7.143654, 0.000000, -0.896764] 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.183380] + discadj_cylinder.test_vals = [2.183381] discadj_cylinder.unsteady = True discadj_cylinder.enabled_with_tsan = False test_list.append(discadj_cylinder) @@ -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.496079, 1.441818, 0.000000] + pywrapper_CFD_AD_MeshDisp.test_vals = [30.000000, -2.496158, 1.441842, 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/incomp_navierstokes/sphere/sphere.cfg b/TestCases/incomp_navierstokes/sphere/sphere.cfg index 0150ea5c0ba..2f9eef68ff8 100644 --- a/TestCases/incomp_navierstokes/sphere/sphere.cfg +++ b/TestCases/incomp_navierstokes/sphere/sphere.cfg @@ -51,7 +51,7 @@ MARKER_MONITORING= ( wall ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= GREEN_GAUSS -CFL_NUMBER= 100.0 +CFL_NUMBER= 25.0 CFL_ADAPT= NO ITER=50000 VENKAT_LIMITER_COEFF= 0.01 @@ -59,7 +59,7 @@ VENKAT_LIMITER_COEFF= 0.01 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % LINEAR_SOLVER= FGMRES -LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_PREC= Q_LU_SGS LINEAR_SOLVER_ERROR= 1E-01 LINEAR_SOLVER_ITER= 3 @@ -100,5 +100,5 @@ VOLUME_FILENAME= flow SURFACE_FILENAME= surface_flow OUTPUT_WRT_FREQ= 100 WRT_VOLUME_OVERWRITE= YES -SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_VELOCITY-Z, LIFT, DRAG) +SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_VELOCITY-Z, LIFT, DRAG, LINSOL_RESIDUAL) VOLUME_OUTPUT= ( SOLUTION,PRIMITIVE,RESIDUAL,MULTIGRID, RANK ) diff --git a/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg b/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg index fbdf927f053..553e71b0ce7 100644 --- a/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg +++ b/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg @@ -35,7 +35,7 @@ REF_ORIGIN_MOMENT_Y = 0.00 REF_ORIGIN_MOMENT_Z = 0.00 REF_LENGTH= 1.0 REF_AREA= 1.0 -REF_DIMENSIONALIZATION= DIMENSIONAL +REF_DIMENSIONALIZATION= FREESTREAM_VEL_EQ_MACH % ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% % @@ -77,7 +77,7 @@ ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % LINEAR_SOLVER= FGMRES -LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_PREC= Q_LU_SGS LINEAR_SOLVER_ERROR= 1E-4 LINEAR_SOLVER_ITER= 5 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index a9825672229..042b8f1853f 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -369,7 +369,7 @@ def main(): poiseuille.cfg_dir = "navierstokes/poiseuille" poiseuille.cfg_file = "lam_poiseuille.cfg" poiseuille.test_iter = 10 - poiseuille.test_vals = [0.648196, 0.000199, 13.639173, 0.000000] + poiseuille.test_vals = [-6.994786, 0.000197, 13.596249, 0] poiseuille.tol = 0.001 test_list.append(poiseuille) @@ -378,7 +378,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.004284, -7.577596, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.004277, -7.577681, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.007498, -7.226926, -0.000000, 2.089953] poiseuille_profile.tol = [0.001, 0.001, 1e-5, 1e-5, 1e-5] test_list.append(poiseuille_profile) @@ -492,7 +492,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.037535, -16.376949, 1.080346, 0.018385, 20.000000, -1.564088, 20.000000, -4.180928, 0.000000] + turb_naca0012_sa.test_vals = [-12.037552, -16.376949, 1.080346, 0.018385, 20.000000, -1.564090, 20.000000, -4.180926, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.037489, -16.376949, 1.080346, 0.018385, 20.000000, -1.564143, 20.000000, -4.180945, 0.000000] turb_naca0012_sa.timeout = 3200 test_list.append(turb_naca0012_sa) @@ -502,7 +502,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.094759, -15.251093, -5.906365, 1.070413, 0.015775, -2.376032, 0.000000] + turb_naca0012_sst.test_vals = [-12.094758, -15.251092, -5.906365, 1.070413, 0.015775, -2.375975, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000] turb_naca0012_sst.timeout = 3200 test_list.append(turb_naca0012_sst) @@ -512,7 +512,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.082056, -14.837177, -5.733436, 1.000893, 0.019109, -2.240933] + turb_naca0012_sst_sust.test_vals = [-12.082110, -14.837177, -5.733436, 1.000893, 0.019109, -2.240999] turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073964, -14.836726, -5.732390, 1.000050, 0.019144, -2.229074] turb_naca0012_sst_sust.timeout = 3200 test_list.append(turb_naca0012_sst_sust) @@ -581,7 +581,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.056236, -5.334117, -8.842317, -4.067917, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-11.056238, -5.334118, -8.842319, -4.067907, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) @@ -646,7 +646,7 @@ def main(): inc_lam_sphere.cfg_dir = "incomp_navierstokes/sphere" inc_lam_sphere.cfg_file = "sphere.cfg" inc_lam_sphere.test_iter = 5 - inc_lam_sphere.test_vals = [-8.190948, -8.992588, 0.121003, 25.782691] + inc_lam_sphere.test_vals = [-7.600533, -8.244915, -8.361301, -9.325293, 0.121003, 25.782687, -1.881890] test_list.append(inc_lam_sphere) # Buoyancy-driven cavity @@ -1191,7 +1191,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 5 - uniform_flow.test_vals = [5.000000, 0.000000, -0.195001, -10.624449] + uniform_flow.test_vals = [5.000000, 0.000000, -0.195001, -10.624453] uniform_flow.unsteady = True uniform_flow.multizone = True test_list.append(uniform_flow) @@ -1363,7 +1363,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.330741, -4.152826, 0, 75] + dyn_fsi.test_vals = [-4.330741, -4.152826, 0.000000, 76.000000] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) @@ -1453,7 +1453,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-12.094759, -15.251093, -5.906365, 1.070413, 0.015775, -2.376032, 0.000000] + pywrapper_turb_naca0012_sst.test_vals = [-12.094758, -15.251092, -5.906365, 1.070413, 0.015775, -2.375975, 0.000000] pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000] pywrapper_turb_naca0012_sst.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f") pywrapper_turb_naca0012_sst.timeout = 3200 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 2b8283f49a4..c3a07802f6b 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -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 = [-2.717803, -3.199669, -2.499149] + da_sp_pinArray_cht_2d_dp_hf.test_vals = [-11.620815, -6.488915, -13.316503] da_sp_pinArray_cht_2d_dp_hf.multizone = True test_list.append(da_sp_pinArray_cht_2d_dp_hf) 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 10f18610446..814ca3d83ff 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.95, 0.00 +19, -6.30, 48.94, 0.00 20, -5.67, 44.59, 0.00 21, -5.01, 39.84, 0.00 22, -4.30, 34.65, 0.00 @@ -95,23 +95,23 @@ 93, 50.49, 53.24, 0.00 94, 63.17, 53.66, 0.00 95, 73.70, 48.01, 0.00 -96, 73.02, 33.47, 0.00 -97, 67.53, 18.28, 0.00 +96, 73.03, 33.47, 0.00 +97, 67.54, 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.59, -25.33, 0.00 -102, 62.72, -28.75, 0.00 +101, 93.60, -25.33, 0.00 +102, 62.71, -28.74, 0.00 103, 27.14, -17.68, 0.00 -104, 23.11, -19.64, 0.00 +104, 23.12, -19.64, 0.00 105, 5.43, -5.72, 0.00 -106, -0.12, 0.15, 0.00 +106, -0.11, 0.15, 0.00 107, -12.46, 18.52, 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.23, 0.00 +112, -39.31, 108.24, 0.00 113, -41.08, 125.40, 0.00 114, -43.50, 146.70, 0.00 115, -44.37, 164.95, 0.00 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 88a3cfe2b92..7265c2b39bc 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -190,7 +190,7 @@ def main(): poiseuille.cfg_dir = "navierstokes/poiseuille" poiseuille.cfg_file = "lam_poiseuille.cfg" poiseuille.test_iter = 10 - poiseuille.test_vals = [-5.050753, 0.648333, 0.012273, 13.643141, 0.000000] + poiseuille.test_vals = [-7.597002, -6.994814, 0.012154, 13.596296, 0] test_list.append(poiseuille) # 2D Poiseuille flow (inlet profile file) @@ -198,7 +198,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.003743, -7.573444, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.003752, -7.573460, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns test_list.append(poiseuille_profile) @@ -283,7 +283,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.037309, -16.384159, 1.080346, 0.018385, 20.000000, -3.456846, 20.000000, -4.641251, 0.000000] + turb_naca0012_sa.test_vals = [-12.037295, -16.384159, 1.080346, 0.018385, 20.000000, -3.457520, 20.000000, -4.641248, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.037297, -16.384158, 1.080346, 0.018385, 20.000000, -3.455886, 20.000000, -4.641247, 0.000000] turb_naca0012_sa.timeout = 3200 test_list.append(turb_naca0012_sa) @@ -293,7 +293,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.094445, -15.251083, -5.906366, 1.070413, 0.015775, -3.178548, 0.000000] + turb_naca0012_sst.test_vals = [-12.094336, -15.251083, -5.906366, 1.070413, 0.015775, -3.178695, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000] turb_naca0012_sst.timeout = 3200 test_list.append(turb_naca0012_sst) @@ -312,7 +312,7 @@ def main(): turb_naca0012_sst_sust_restart.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust_restart.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust_restart.test_iter = 10 - turb_naca0012_sst_sust_restart.test_vals = [-12.080455, -14.837169, -5.733461, 1.000893, 0.019109, -2.634140] + turb_naca0012_sst_sust_restart.test_vals = [-12.080498, -14.837169, -5.733461, 1.000893, 0.019109, -2.634000] turb_naca0012_sst_sust_restart.test_vals_aarch64 = [-12.074189, -14.836725, -5.732398, 1.000050, 0.019144, -3.315560] turb_naca0012_sst_sust_restart.timeout = 3200 test_list.append(turb_naca0012_sst_sust_restart) @@ -344,7 +344,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.054279, -5.328901, -8.835585, -4.056810, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-11.054282, -5.328899, -8.835603, -4.056836, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) @@ -946,7 +946,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 2 - uniform_flow.test_vals = [2.000000, 0.000000, -0.230639, -13.253272] + uniform_flow.test_vals = [2.000000, 0.000000, -0.230639, -13.252043] uniform_flow.test_vals_aarch64 = [2.000000, 0.000000, -0.230641, -13.249000] uniform_flow.tol = 0.000001 uniform_flow.unsteady = True @@ -1089,7 +1089,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, -5.077002, -5.379449, -9.247803, -9.319973, -9.184965, 608.35, -0.012973, 0, 30] fsi_cht.multizone = True test_list.append(fsi_cht) @@ -1558,7 +1558,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-12.094445, -15.251083, -5.906366, 1.070413, 0.015775, -3.178548, 0.000000] + pywrapper_turb_naca0012_sst.test_vals = [-12.094336, -15.251083, -5.906366, 1.070413, 0.015775, -3.178695, 0.000000] pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000] pywrapper_turb_naca0012_sst.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_turb_naca0012_sst.timeout = 3200 diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 4f73d7fb1fa..4e039eea09b 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -237,7 +237,7 @@ def main(): tutorial_trans_flatplate_T3A.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A" tutorial_trans_flatplate_T3A.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3A.test_iter = 20 - tutorial_trans_flatplate_T3A.test_vals = [-5.790137, -2.054834, -3.894659, -0.255074, -1.747098, 5.119341, -3.493237, 0.393262] + tutorial_trans_flatplate_T3A.test_vals = [-5.790137, -2.054834, -3.894659, -0.255074, -1.747087, 5.119341, -3.493237, 0.393262] tutorial_trans_flatplate_T3A.test_vals_aarch64 = [-5.808996, -2.070606, -3.969765, -0.277943, -1.953289, 1.708472, -3.514943, 0.357411] tutorial_trans_flatplate_T3A.no_restart = True test_list.append(tutorial_trans_flatplate_T3A) @@ -247,7 +247,7 @@ def main(): tutorial_trans_flatplate_T3Am.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A-" tutorial_trans_flatplate_T3Am.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3Am.test_iter = 20 - tutorial_trans_flatplate_T3Am.test_vals = [-5.587332, -1.700868, -3.093872, -0.102783, -3.750523, 3.287643, -2.394575, 1.119623] + tutorial_trans_flatplate_T3Am.test_vals = [-5.587389, -1.700868, -3.093936, -0.102834, -3.750523, 3.287643, -2.394575, 1.119623] tutorial_trans_flatplate_T3Am.test_vals_aarch64 = [-5.540938, -1.681627, -2.878831, -0.058224, -3.695533, 3.413628, -2.385345, 1.103633] tutorial_trans_flatplate_T3Am.no_restart = True test_list.append(tutorial_trans_flatplate_T3Am) diff --git a/TestCases/vandv.py b/TestCases/vandv.py index 01bace8011e..2adaa2e6c88 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -63,7 +63,7 @@ def main(): flatplate_sst1994m.cfg_dir = "vandv/rans/flatplate" flatplate_sst1994m.cfg_file = "turb_flatplate_sst.cfg" flatplate_sst1994m.test_iter = 5 - flatplate_sst1994m.test_vals = [-13.040081, -10.136488, -10.941822, -7.977797, -10.323859, -4.732970, 0.002801] + flatplate_sst1994m.test_vals = [-13.041178, -10.137005, -10.939481, -7.988319, -10.323869, -4.732381, 0.002801] flatplate_sst1994m.test_vals_aarch64 = [-13.021715, -9.534786, -10.401912, -7.501836, -9.750800, -4.850665, 0.002807] test_list.append(flatplate_sst1994m) @@ -72,7 +72,7 @@ def main(): bump_sst1994m.cfg_dir = "vandv/rans/bump_in_channel" bump_sst1994m.cfg_file = "turb_bump_sst.cfg" bump_sst1994m.test_iter = 5 - bump_sst1994m.test_vals = [-11.928608, -10.096068, -9.512742, -6.445912, -11.774007, -6.978081, 0.004931] + bump_sst1994m.test_vals = [-11.928336, -10.095735, -9.512439, -6.445676, -11.773921, -6.995865, 0.004931] bump_sst1994m.test_vals_aarch64 = [-13.042689, -10.812982, -10.604523, -7.655547, -10.816257, -5.308083, 0.004911] test_list.append(bump_sst1994m) @@ -99,7 +99,7 @@ def main(): dsma661_sa.cfg_dir = "vandv/rans/dsma661" dsma661_sa.cfg_file = "dsma661_sa_config.cfg" dsma661_sa.test_iter = 5 - dsma661_sa.test_vals = [-11.240967, -8.243826, -8.958188, -5.895885, -10.737669, 0.155687, 0.024232] + dsma661_sa.test_vals = [-11.237268, -8.242347, -9.007582, -5.887860, -10.737674, 0.155687, 0.024232] dsma661_sa.test_vals_aarch64 = [-11.293183, -8.241775, -9.083761, -6.011398, -10.737680, 0.155687, 0.024232] test_list.append(dsma661_sa) @@ -108,7 +108,7 @@ def main(): dsma661_sst.cfg_dir = "vandv/rans/dsma661" dsma661_sst.cfg_file = "dsma661_sst_config.cfg" dsma661_sst.test_iter = 5 - dsma661_sst.test_vals = [-11.023511, -8.156964, -9.060365, -5.934828, -10.651368, -7.898758, 0.155882, 0.023344] + dsma661_sst.test_vals = [-11.022853, -8.157766, -9.043824, -5.929729, -10.651421, -7.904970, 0.155882, 0.023344] dsma661_sst.test_vals_aarch64 = [-10.977195, -8.403731, -8.747068, -5.808899, -10.522786, -7.369851, 0.155875, 0.023353] test_list.append(dsma661_sst) diff --git a/UnitTests/Common/linear_algebra/quantization_tests.cpp b/UnitTests/Common/linear_algebra/quantization_tests.cpp new file mode 100644 index 00000000000..d074a6f16e2 --- /dev/null +++ b/UnitTests/Common/linear_algebra/quantization_tests.cpp @@ -0,0 +1,113 @@ +/*! + * \file quantization_tests.cpp + * \brief Unit tests for the int8 row-scaled block quantization used by Q_LU_SGS. + * \author 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 . + */ + +#include "catch.hpp" +#include "../../../Common/include/linear_algebra/CSysMatrix.hpp" + +/*--- Row-major access into a flat quantized block, matching CBlockView's decode. ---*/ +static su2double Decode(const int8_t* qs, const int8_t* qv, unsigned long nVar, unsigned long r, unsigned long c) { + return static_cast(qv[r * nVar + c]) * DecodeQuantScale(qs[r]); +} + +TEST_CASE("Quantization round-trip is stable for well-behaved values", "[LinearAlgebra]") { + /*--- EncodeQuantBlock always encodes a full nVar x nVar block (one scale per row), + * so every row is set to the same values here and only row 0 is inspected. Values + * are moderate: none close enough to the row's own max to saturate, and none land + * exactly on the int8 min (-128), so the derived per-row scale does not shift + * between successive encode/decode passes. ---*/ + constexpr unsigned long nVar = 4; + const su2double row[nVar] = {3.0, -2.0, 1.5, 0.75}; + auto f = [&](unsigned long, unsigned long c) { return row[c]; }; + + int8_t qs1[nVar], qv1[nVar * nVar]; + EncodeQuantBlock(f, qs1, qv1, nVar); + + auto decoded = [&](unsigned long, unsigned long c) { return Decode(qs1, qv1, nVar, 0, c); }; + + int8_t qs2[nVar], qv2[nVar * nVar]; + EncodeQuantBlock(decoded, qs2, qv2, nVar); + + CHECK(qs2[0] == qs1[0]); + for (unsigned long c = 0; c < nVar; ++c) CHECK(qv2[c] == qv1[c]); + + /*--- A second decode/encode cycle from the now-stable representation must + * reproduce the exact same codes again. ---*/ + auto decoded2 = [&](unsigned long, unsigned long c) { return Decode(qs2, qv2, nVar, 0, c); }; + int8_t qs3[nVar], qv3[nVar * nVar]; + EncodeQuantBlock(decoded2, qs3, qv3, nVar); + + CHECK(qs3[0] == qs2[0]); + for (unsigned long c = 0; c < nVar; ++c) CHECK(qv3[c] == qv2[c]); +} + +TEST_CASE("Quantization saturates and truncates with a known, bounded error", "[LinearAlgebra]") { + /*--- One row engineered so that, at the scale it forces (2^0 = 1 here, since the + * largest magnitude in the row is in [64, 128)): + * col 0: 127.9 -> rounds to 128, clamped to 127 (int8 max), error ~0.9 (close to + * the theoretical worst case: clamping can only ever push a value that + * rounds to +128 down to +127, an error that approaches but never reaches + * 1 full scale unit). + * col 1: -127.9 -> rounds to -128 (int8 min), no clamping needed since -128 is a + * representable int8 value; error ~0.1 (near-exact). Clamping can only + * ever trigger on the positive side because int8 is asymmetric + * ([-128, 127]) while the encoding is symmetric around the row's max abs. + * col 2: 0.2 -> rounds to 0: small values are truncated away entirely. + * col 3: 60.0 -> reconstructs exactly, since the row's scale (2^0 = 1) is an + * exact power of two and 60 fits in int8. ---*/ + constexpr unsigned long nVar = 4; + const su2double row[nVar] = {127.9, -127.9, 0.2, 60.0}; + auto f = [&](unsigned long, unsigned long c) { return row[c]; }; + + int8_t qs[nVar], qv[nVar * nVar]; + EncodeQuantBlock(f, qs, qv, nVar); + + CHECK(qs[0] == 0); + CHECK(static_cast(qv[0]) == 127); // Saturated at the int8 maximum. + CHECK(static_cast(qv[1]) == -128); // Hits the int8 minimum exactly. + CHECK(static_cast(qv[2]) == 0); // Truncated to zero. + CHECK(static_cast(qv[3]) == 60); // Exact, mid-range value. + + const su2double scale = DecodeQuantScale(qs[0]); + CHECK(scale == Approx(1.0)); + + const su2double error0 = std::abs(row[0] - Decode(qs, qv, nVar, 0, 0)); + const su2double error1 = std::abs(row[1] - Decode(qs, qv, nVar, 0, 1)); + const su2double error2 = std::abs(row[2] - Decode(qs, qv, nVar, 0, 2)); + const su2double error3 = std::abs(row[3] - Decode(qs, qv, nVar, 0, 3)); + + /*--- The quantization error for any entry is strictly bounded by one scale unit: + * 0.5 from rounding to the nearest representable level, plus (only for values that + * saturate on the positive side) less than another 0.5 from clamping 128 down to 127. ---*/ + CHECK(error0 < scale); + CHECK(error1 < scale); + CHECK(error2 < scale); + CHECK(error3 == Approx(0.0).margin(1e-12)); + + CHECK(error0 == Approx(0.9).margin(1e-5)); + CHECK(error1 == Approx(0.1).margin(1e-5)); + CHECK(error2 == Approx(0.2).margin(1e-5)); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 76eecd1dc57..f8c22511b5f 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -17,7 +17,8 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/windowing.cpp', - 'Common/toolboxes/random_toolbox_tests.cpp']) + 'Common/toolboxes/random_toolbox_tests.cpp', + 'Common/linear_algebra/quantization_tests.cpp']) # Reverse-mode (algorithmic differentiation) tests: su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp', diff --git a/config_template.cfg b/config_template.cfg index 998360550b3..e6f0083fe25 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1648,7 +1648,11 @@ DISCADJ_LIN_SOLVER= FGMRES % Use CUDA GPU Acceleration for FGMRES Linear Solver Only ENABLE_CUDA=NO % -% Preconditioner of the Krylov linear solver or type of smoother (ILU, LU_SGS, LINELET, JACOBI) +% Preconditioner of the Krylov linear solver or type of smoother (ILU, LU_SGS, Q_LU_SGS, LINELET, JACOBI) +% Q_LU_SGS triggers the use of quantization to reduce the size of the sparse matrix (for compressible flow +% the matrix becomes 3x smaller relative to mixed-precision mode). This is only used by the compressible +% and incompressible solvers, others fallback silently to LU_SGS. A suitable nondimensionalization mode +% MUST be used otherwise the solver is very likely to diverge. LINEAR_SOLVER_PREC= ILU % % Same for discrete adjoint (JACOBI or ILU), replaces LINEAR_SOLVER_PREC in SU2_*_AD codes.