diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 3284f952736..c40dc37f3d6 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -9733,6 +9733,14 @@ class CConfig { */ bool GetIntegrated_HeatFlux() const { return Integrated_HeatFlux; } + /*! + * \brief Whether the factorized least-squares gradient metric terms are cached and reused + * across evaluations. This is the default behavior, except for periodic boundaries + * (their metric and RHS accumulations are fused in one exchange) and the discrete + * adjoint (the coordinate dependence of the metrics must remain on the tape). + */ + bool GetLSQMetricCaching() const { return (nMarker_PerBound == 0) && !DiscreteAdjoint; } + /*! * \brief Get Compute Average. * \return YES if start computing averages diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 97d3f845537..3dae2cbb235 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -217,6 +217,12 @@ class CGeometry { unsigned long edgeColorGroupSize{1}; /*!< \brief Size of the edge groups within each color. */ unsigned long elemColorGroupSize{1}; /*!< \brief Size of the element groups within each color. */ + /*--- Cached least-squares gradient metric terms (see computeGradientsLeastSquares.hpp). ---*/ + + su2activematrix LSQMetricCache[2]; /*!< \brief Cached LSQ metrics S = inv(A), upper triangle stored + row-wise, [0] unweighted, [1] inverse-distance weighted. */ + bool LSQMetricCacheValid[2] = {false, false}; /*!< \brief Validity of the cached LSQ metrics per weighting. */ + ColMajorMatrix CoarseGridColor_; /*!< \brief Coarse grid levels, colorized. */ public: @@ -1919,6 +1925,28 @@ class CGeometry { */ void SetNaturalEdgeColoring(); + /*! + * \brief Get the cached least-squares metric terms (S = inv(A), upper triangle row-wise). + * \param[in] weighted - False for unweighted, true for inverse-distance weighting. + */ + inline su2activematrix& GetLSQMetricCache(bool weighted) { return LSQMetricCache[weighted]; } + inline const su2activematrix& GetLSQMetricCache(bool weighted) const { return LSQMetricCache[weighted]; } + + /*! + * \brief Check whether the cached least-squares metric terms are valid for a weighting. + */ + inline bool LSQMetricCacheIsValid(bool weighted) const { return LSQMetricCacheValid[weighted]; } + + /*! + * \brief Declare the cached least-squares metric terms valid for a weighting. + */ + inline void SetLSQMetricCacheValid(bool weighted) { LSQMetricCacheValid[weighted] = true; } + + /*! + * \brief Invalidate the cached least-squares metric terms (e.g. if node coordinates change). + */ + inline void InvalidateLSQMetricCache() { LSQMetricCacheValid[0] = LSQMetricCacheValid[1] = false; } + /*! * \brief Get the group size used in edge coloring. * \return Group size. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c2633f84ec4..5d04ef364b6 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3886,6 +3886,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i MUSCL_Flow = false; } } + if (MUSCL_AdjFlow && (Kind_ConvNumScheme_AdjFlow == SPACE_CENTERED)) { if (OptionIsSet("MUSCL_ADJFLOW")) { SU2_MPI::Error("Centered schemes do not use MUSCL reconstruction (use MUSCL_ADJFLOW= NO).", CURRENT_FUNCTION); diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 684da742b13..2c0c64272dc 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1000,6 +1000,9 @@ void CMultiGridGeometry::MatchActuator_Disk(const CConfig* config) { void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned short action) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + /*--- Coarse coordinates change with the fine grid, any cached LSQ metrics are stale. ---*/ + if (action != ALLOCATE) InvalidateLSQMetricCache(); + /*--- Compute the area of the coarse volume ---*/ for (auto iCoarsePoint = 0u; iCoarsePoint < nPoint; iCoarsePoint++) { nodes->SetVolume(iCoarsePoint, 0.0); diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 199c8cec03f..03945668dc7 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -6759,6 +6759,9 @@ void CPhysicalGeometry::SetControlVolume(CConfig* config, unsigned short action) BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { /*--- The following is difficult to parallelize with threads. ---*/ + /*--- Node coordinates changed, any cached LSQ gradient metrics are stale. ---*/ + if (action != ALLOCATE) InvalidateLSQMetricCache(); + su2double my_DomainVolume = 0.0; for (auto iElem = 0ul; iElem < nElem; iElem++) { const auto nNodes = elem[iElem]->GetnNodes(); diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp index 296f433ff40..44e3ee95bc0 100644 --- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -32,6 +32,14 @@ namespace detail { +/*! + * \brief Flattened index of entry (iDim,jDim), iDim <= jDim, of an upper triangular + * matrix stored row-wise (the layout of the cached LSQ metric terms). + */ +FORCEINLINE constexpr size_t lsqCacheIdx(size_t nDim, size_t iDim, size_t jDim) { + return iDim * nDim - (iDim * (iDim - 1)) / 2 + (jDim - iDim); +} + /*! * \brief Prepare Smatrix for 2D. * \ingroup FvmAlgos @@ -66,6 +74,45 @@ FORCEINLINE void computeSmatrix(su2double r11, su2double r12, su2double r13, Smatrix[2][2] = (z33*z33)/detR2; } +/*! + * \brief Factorize the accumulated normal matrix A of the least-squares problem + * (Cholesky) and form S = inv(A), the entries r* are the unique entries of A. + * A (nearly) singular matrix results in S = 0, i.e. a zero gradient. + * \ingroup FvmAlgos + */ +template +FORCEINLINE void invertNormalMatrix(su2double r11, su2double r12, su2double r13, + su2double r22, su2double r23_a, su2double r23_b, + su2double r33, su2double Smatrix[][nDim]) +{ + const auto eps = pow(std::numeric_limits::epsilon(),2); + + r11 = sqrt(max(r11, eps)); + r12 /= r11; + r22 = sqrt(max(r22 - r12*r12, eps)); + + su2double r23 = 0.0; + if (nDim == 3) { + r13 /= r11; + r23 = r23_a/r22 - r23_b*r12/(r11*r22); + r33 = sqrt(max(r33 - r23*r23 - r13*r13, eps)); + } + else { + r13 = 0.0; + r33 = 1.0; + } + + /*--- Compute determinant ---*/ + + const su2double detR2 = pow(r11*r22*r33, 2); + + /*--- S matrix := inv(R)*traspose(inv(R)), detect singular matrix ---*/ + + if (detR2 > eps) { + computeSmatrix(r11, r12, r13, r22, r23, r33, detR2, Smatrix); + } +} + /*! * \brief Solve the least-squares problem for one point. * \ingroup FvmAlgos @@ -79,9 +126,7 @@ FORCEINLINE void solveLeastSquares(size_t iPoint, const RMatrixType& Rmatrix, GradientType& gradient) { - const auto eps = pow(std::numeric_limits::epsilon(),2); - - /*--- Entries of upper triangular matrix R. ---*/ + /*--- Entries of the normal matrix A. ---*/ if (periodic) { AD::StartPreacc(); @@ -90,14 +135,10 @@ FORCEINLINE void solveLeastSquares(size_t iPoint, AD::SetPreaccIn(Rmatrix(iPoint,1,1)); } - su2double r11 = Rmatrix(iPoint,0,0); - su2double r12 = Rmatrix(iPoint,0,1); - su2double r22 = Rmatrix(iPoint,1,1); - su2double r13 = 0.0, r23 = 0.0, r33 = 1.0; - - r11 = sqrt(max(r11, eps)); - r12 /= r11; - r22 = sqrt(max(r22 - r12*r12, eps)); + const su2double r11 = Rmatrix(iPoint,0,0); + const su2double r12 = Rmatrix(iPoint,0,1); + const su2double r22 = Rmatrix(iPoint,1,1); + su2double r13 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0; if (nDim == 3) { if (periodic) { @@ -109,27 +150,13 @@ FORCEINLINE void solveLeastSquares(size_t iPoint, r13 = Rmatrix(iPoint,0,2); r33 = Rmatrix(iPoint,2,2); - const auto r23_a = Rmatrix(iPoint,1,2); - const auto r23_b = Rmatrix(iPoint,2,1); - - r13 /= r11; - r23 = r23_a/r22 - r23_b*r12/(r11*r22); - r33 = sqrt(max(r33 - r23*r23 - r13*r13, eps)); + r23_a = Rmatrix(iPoint,1,2); + r23_b = Rmatrix(iPoint,2,1); } - /*--- Compute determinant ---*/ - - const su2double detR2 = pow(r11*r22*r33, 2); - - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - su2double Smatrix[nDim][nDim] = {{0.0}}; - /*--- Detect singular matrix ---*/ - - if (detR2 > eps) { - computeSmatrix(r11, r12, r13, r22, r23, r33, detR2, Smatrix); - } + invertNormalMatrix(r11, r12, r13, r22, r23_a, r23_b, r33, Smatrix); if (periodic) { /*--- Stop preacc here as gradient is in/out. ---*/ @@ -162,6 +189,233 @@ FORCEINLINE void solveLeastSquares(size_t iPoint, } } +/*! + * \brief Assemble, factorize, and store the least-squares gradient metric terms + * (S = inv(A), upper triangle row-wise) of a grid in the geometry cache. + * \ingroup FvmAlgos + * \note The metric terms depend only on the node coordinates and the weighting. They are + * computed during the geometry preprocessing (see CDriver::InitializeGeometry) and, + * on moving/deforming grids, recomputed on the first gradient evaluation after the + * dual grid update invalidates them (CGeometry::SetControlVolume). The function is + * safe to call from inside or outside an OpenMP parallel region, and returns + * immediately if the cache is already valid. + */ +template +void computeLSQMetrics(CGeometry& geometry, bool weighted) +{ + if (geometry.LSQMetricCacheIsValid(weighted)) return; + + const size_t nPointDomain = geometry.GetnPointDomain(); + +#ifdef HAVE_OMP + constexpr size_t OMP_MAX_CHUNK = 512; + + const size_t chunkSize = computeStaticChunkSize(nPointDomain, + omp_get_max_threads(), OMP_MAX_CHUNK); +#endif + + auto& metricCache = geometry.GetLSQMetricCache(weighted); + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + if (metricCache.size() == 0) metricCache.resize(nPointDomain, nDim*(nDim+1)/2); + } END_SU2_OMP_SAFE_GLOBAL_ACCESS + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) { + const auto coord_i = geometry.nodes->GetCoord(iPoint); + + /*--- Accumulate the unique entries of the normal matrix A. ---*/ + + su2double r11 = 0.0, r12 = 0.0, r13 = 0.0, r22 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0; + + for (auto jPoint : geometry.nodes->GetPoints(iPoint)) { + su2double dist_ij[nDim] = {0.0}; + GeometryToolbox::Distance(nDim, geometry.nodes->GetCoord(jPoint), coord_i, dist_ij); + + su2double weight = 1.0; + if (weighted) { + const su2double dist2 = GeometryToolbox::SquaredNorm(nDim, dist_ij); + if (dist2 <= 0.0) continue; + weight = 1.0 / dist2; + } + + r11 += dist_ij[0]*dist_ij[0]*weight; + r12 += dist_ij[0]*dist_ij[1]*weight; + r22 += dist_ij[1]*dist_ij[1]*weight; + + if (nDim == 3) { + r13 += dist_ij[0]*dist_ij[2]*weight; + r23_a += dist_ij[1]*dist_ij[2]*weight; + r23_b += dist_ij[0]*dist_ij[2]*weight; + r33 += dist_ij[2]*dist_ij[2]*weight; + } + } + + su2double Smatrix[nDim][nDim] = {{0.0}}; + + invertNormalMatrix(r11, r12, r13, r22, r23_a, r23_b, r33, Smatrix); + + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = iDim; jDim < nDim; ++jDim) + metricCache(iPoint, lsqCacheIdx(nDim, iDim, jDim)) = Smatrix[iDim][jDim]; + } + END_SU2_OMP_FOR + + /*--- Declare the cache valid and make sure the edge coloring is available before the + * first cached evaluation, building it here avoids a race on its lazy construction. ---*/ + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + geometry.GetEdgeColoring(); + geometry.SetLSQMetricCacheValid(weighted); + } END_SU2_OMP_SAFE_GLOBAL_ACCESS +} + +/*! + * \brief Fast least-squares gradient evaluation reusing the cached metric terms. + * \ingroup FvmAlgos + * \note Requires valid metric terms for this weighting in the geometry cache (see + * computeLSQMetrics, shared by all solvers, one slot per weighting). Only the + * right-hand side b = sum_k w*dist*(u_k - u_i) is accumulated (in an edge loop, + * i.e. each edge is visited once since its contribution is identical for both + * end points), followed by the product S*b per point. Not compatible with periodic + * boundaries. + */ +template +void computeGradientsLeastSquaresCached(CSolver* solver, + MPI_QUANTITIES kindMpiComm, + CGeometry& geometry, + const CConfig& config, + bool weighted, + const FieldType& field, + const size_t varBegin, + const size_t varEnd, + const int idxVel, + GradientType& gradient) +{ + const auto& metricCache = geometry.GetLSQMetricCache(weighted); + const size_t nPoint = geometry.GetnPoint(); + const size_t nPointDomain = geometry.GetnPointDomain(); + +#ifdef HAVE_OMP + constexpr size_t OMP_MAX_CHUNK = 512; + + const size_t chunkSize = computeStaticChunkSize(nPointDomain, + omp_get_max_threads(), OMP_MAX_CHUNK); +#endif + + /*--- Clear the right-hand-side accumulators, including halo points, which + * receive edge contributions (discarded when halos are communicated). ---*/ + + SU2_OMP_FOR_STAT(2048) + for (size_t iPoint = 0; iPoint < nPoint; ++iPoint) + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) = 0.0; + END_SU2_OMP_FOR + + /*--- Accumulate the RHS in a loop over the edges: the contribution of edge {i,j} is + * w*dist_ij*(u_j - u_i) for BOTH end points. A race-free edge coloring is required + * with multiple threads, the "natural" coloring (single color, used with the + * reducer strategy) forces the fallback to a thread-safe loop over nodes. ---*/ + + const auto& coloring = geometry.GetEdgeColoring(); + + const bool safeColoring = (omp_get_max_threads() == 1) || (coloring.getOuterSize() > 1); + + if (safeColoring) { + const size_t groupSize = geometry.GetEdgeColorGroupSize(); + + for (auto iColor = 0ul; iColor < coloring.getOuterSize(); ++iColor) { + const auto* edgeIndices = coloring.innerIdx(iColor); + const auto nEdgesColor = coloring.getNumNonZeros(iColor); + + SU2_OMP_FOR_DYN(nextMultiple(size_t(32), groupSize)) + for (auto k = 0ul; k < nEdgesColor; ++k) { + const auto iEdge = edgeIndices[k]; + const auto iPoint = geometry.edges->GetNode(iEdge, 0); + const auto jPoint = geometry.edges->GetNode(iEdge, 1); + + su2double dist_ij[nDim] = {0.0}; + GeometryToolbox::Distance(nDim, geometry.nodes->GetCoord(jPoint), + geometry.nodes->GetCoord(iPoint), dist_ij); + + su2double weight = 1.0; + if (weighted) { + const su2double dist2 = GeometryToolbox::SquaredNorm(nDim, dist_ij); + if (dist2 <= 0.0) continue; + weight = 1.0 / dist2; + } + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) { + const su2double delta_ij = weight * (field(jPoint,iVar) - field(iPoint,iVar)); + + for (size_t iDim = 0; iDim < nDim; ++iDim) { + const su2double contrib = dist_ij[iDim] * delta_ij; + gradient(iPoint, iVar, iDim) += contrib; + gradient(jPoint, iVar, iDim) += contrib; + } + } + } + END_SU2_OMP_FOR + } + } + else { + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) { + const auto coord_i = geometry.nodes->GetCoord(iPoint); + + for (auto jPoint : geometry.nodes->GetPoints(iPoint)) { + su2double dist_ij[nDim] = {0.0}; + GeometryToolbox::Distance(nDim, geometry.nodes->GetCoord(jPoint), coord_i, dist_ij); + + su2double weight = 1.0; + if (weighted) { + const su2double dist2 = GeometryToolbox::SquaredNorm(nDim, dist_ij); + if (dist2 <= 0.0) continue; + weight = 1.0 / dist2; + } + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) { + const su2double delta_ij = weight * (field(jPoint,iVar) - field(iPoint,iVar)); + + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) += dist_ij[iDim] * delta_ij; + } + } + } + END_SU2_OMP_FOR + } + + /*--- Multiply the RHS by the cached S matrix. ---*/ + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) { + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) { + su2double Cvector[nDim] = {0.0}; + + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = 0; jDim < nDim; ++jDim) + Cvector[iDim] += metricCache(iPoint, lsqCacheIdx(nDim, min(iDim,jDim), max(iDim,jDim))) * + gradient(iPoint, iVar, jDim); + + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) = Cvector[iDim]; + } + } + END_SU2_OMP_FOR + + /*--- Compute the corrections for symmetry planes and Euler walls. ---*/ + + correctGradientsSymmetry(geometry, config, varBegin, varEnd, idxVel, gradient); + + /*--- Obtain the gradients at halo points from the MPI ranks that own them. ---*/ + + if (solver != nullptr) { + solver->InitiateComms(&geometry, &config, kindMpiComm); + solver->CompleteComms(&geometry, &config, kindMpiComm); + } +} + /*! * \brief Compute the gradient of a field using inverse-distance-weighted or * unweighted Least-Squares approximation. @@ -192,10 +446,22 @@ void computeGradientsLeastSquares(CSolver* solver, const size_t varEnd, const int idxVel, GradientType& gradient, - RMatrixType& Rmatrix) + RMatrixType& Rmatrix, + bool useCaching) { const bool periodic = (solver != nullptr) && (config.GetnMarker_Periodic() > 0); + /*--- Use the cached metric terms, rebuilding them if the coordinates changed since the + * geometry preprocessing (or if this combination was not covered by it). The caching + * does not support the mid-computation periodic accumulations. ---*/ + + if (useCaching && !periodic) { + computeLSQMetrics(geometry, weighted); + computeGradientsLeastSquaresCached(solver, kindMpiComm, geometry, config, weighted, + field, varBegin, varEnd, idxVel, gradient); + return; + } + const size_t nPointDomain = geometry.GetnPointDomain(); #ifdef HAVE_OMP @@ -348,15 +614,35 @@ void computeGradientsLeastSquares(CSolver* solver, const size_t varEnd, const int idxVel, GradientType& gradient, - RMatrixType& Rmatrix) { + RMatrixType& Rmatrix, + bool useCaching = false) { switch (geometry.GetnDim()) { case 2: detail::computeGradientsLeastSquares<2>(solver, kindMpiComm, kindPeriodicComm, geometry, config, - weighted, field, varBegin, varEnd, idxVel, gradient, Rmatrix); + weighted, field, varBegin, varEnd, idxVel, gradient, Rmatrix, useCaching); break; case 3: detail::computeGradientsLeastSquares<3>(solver, kindMpiComm, kindPeriodicComm, geometry, config, - weighted, field, varBegin, varEnd, idxVel, gradient, Rmatrix); + weighted, field, varBegin, varEnd, idxVel, gradient, Rmatrix, useCaching); + break; + default: + SU2_MPI::Error("Too many dimensions to compute gradients.", CURRENT_FUNCTION); + break; + } +} + +/*! + * \brief Compute (if not already valid) the cached least-squares gradient metric terms of + * a grid for one type of weighting, see detail::computeLSQMetrics. + * \ingroup FvmAlgos + */ +inline void computeLSQGradientMetrics(CGeometry& geometry, bool weighted) { + switch (geometry.GetnDim()) { + case 2: + detail::computeLSQMetrics<2>(geometry, weighted); + break; + case 3: + detail::computeLSQMetrics<3>(geometry, weighted); break; default: SU2_MPI::Error("Too many dimensions to compute gradients.", CURRENT_FUNCTION); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 5b1c49b75d9..41e4c08ddb9 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -425,7 +425,8 @@ void CFVMFlowSolverBase::SetPrimitive_Gradient_LS(CGeometry* geometry, con const auto comm = reconstruction? MPI_QUANTITIES::PRIMITIVE_GRAD_REC : MPI_QUANTITIES::PRIMITIVE_GRADIENT; computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, - primitives, 0, nPrimVarGrad, prim_idx.Velocity(), gradient, rmatrix); + primitives, 0, nPrimVarGrad, prim_idx.Velocity(), gradient, rmatrix, + config->GetLSQMetricCaching()); } template diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 63d48590098..cc120cf3f10 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -32,6 +32,8 @@ #include "../../../Common/include/geometry/CPhysicalGeometry.hpp" #include "../../../Common/include/geometry/CMultiGridGeometry.hpp" +#include "../../include/gradients/computeGradientsLeastSquares.hpp" + #include "../../include/solvers/CSolverFactory.hpp" #include "../../include/solvers/CFEM_DG_EulerSolver.hpp" @@ -670,6 +672,33 @@ void CDriver::InitializeGeometry(CConfig* config, CGeometry **&geometry, bool du } } + /*--- Precompute the least-squares gradient metric terms (S = inv(A)) required by the + * numerical settings, they depend only on the grid coordinates and the weighting. + * Combinations not covered here (e.g. auxiliary variable gradients) are computed on + * first use, and on moving/deforming grids the terms are recomputed after each mesh + * update. ---*/ + + if (!dummy && !fem_solver && config->GetLSQMetricCaching()) { + + const auto kindGrad = config->GetKind_Gradient_Method(); + const bool lsqGrad = (kindGrad == LEAST_SQUARES) || (kindGrad == WEIGHTED_LEAST_SQUARES); + + const auto kindGradRecon = config->GetKind_Gradient_Method_Recon(); + const bool lsqGradRecon = config->GetReconstructionGradientRequired() && + ((kindGradRecon == LEAST_SQUARES) || (kindGradRecon == WEIGHTED_LEAST_SQUARES)); + + if (lsqGrad) computeLSQGradientMetrics(*geometry[MESH_0], kindGrad == WEIGHTED_LEAST_SQUARES); + if (lsqGradRecon) computeLSQGradientMetrics(*geometry[MESH_0], kindGradRecon == WEIGHTED_LEAST_SQUARES); + + /*--- The coarse multigrid levels only compute gradients for the viscous fluxes. ---*/ + + if (lsqGrad && config->GetViscous()) { + for (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { + computeLSQGradientMetrics(*geometry[iMesh], kindGrad == WEIGHTED_LEAST_SQUARES); + } + } + } + /*--- If activated by the compile directive, perform a partition analysis. ---*/ #if PARTITION if (!dummy){ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e014c794a52..9bfbd89a7ea 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -2175,7 +2175,8 @@ void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) auto& rmatrix = base_nodes->GetRmatrix(); computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, - weighted, solution, 0, base_nodes->GetnAuxVar(), -1, gradient, rmatrix); + weighted, solution, 0, base_nodes->GetnAuxVar(), -1, gradient, rmatrix, + config->GetLSQMetricCaching()); } void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, const CConfig *config, short idxVel, bool reconstruction) { @@ -2209,7 +2210,8 @@ void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, const CConfig *config auto& gradient = reconstruction? base_nodes->GetGradient_Reconstruction() : base_nodes->GetGradient(); const auto comm = reconstruction? MPI_QUANTITIES::SOLUTION_GRAD_REC : MPI_QUANTITIES::SOLUTION_GRADIENT; - computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, solution, 0, nVar, idxVel, gradient, rmatrix); + computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, solution, 0, nVar, idxVel, + gradient, rmatrix, config->GetLSQMetricCaching()); } void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) {