diff --git a/Common/include/containers/container_decorators.hpp b/Common/include/containers/container_decorators.hpp index 331e87b9a0f4..82b89814960a 100644 --- a/Common/include/containers/container_decorators.hpp +++ b/Common/include/containers/container_decorators.hpp @@ -126,7 +126,8 @@ class C3DContainerDecorator { private: Storage m_storage; - Index m_innerSz; + /*--- One, not zero, so rows() on a container that was never resized does not divide by it. ---*/ + Index m_innerSz = 1; public: C3DContainerDecorator() = default; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bcd2e9c06328..506d4cf83678 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -80,6 +80,7 @@ const unsigned int MAX_NUMBER_FFD = 15; /*!< \brief Maximum number of FFDB enum: unsigned int{MAX_SOLS = 14}; /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */ const unsigned int MAX_TERMS = 7; /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */ const unsigned int MAX_ZONES = 3; /*!< \brief Maximum number of zones. */ +const unsigned short MAX_MGLEVELS = 10; /*!< \brief Maximum number of coarse multigrid levels, which bounds the per-level arrays of the multigrid integration. */ const unsigned int MAX_FE_KINDS = 7; /*!< \brief Maximum number of Finite Elements. */ const unsigned int NO_RK_ITER = 0; /*!< \brief No Runge-Kutta iteration. */ @@ -1234,6 +1235,8 @@ struct CMGOptions { bool MG_Smooth_Output{false}; /*!< \brief Output compact per-cycle smoothing summary. */ su2double MG_Smooth_StagnationTol{0.0}; /*!< \brief Stagnation early exit: stop if current_rms >= prev_rms * tol. 0 = disabled. */ bool MG_Implicit_Lines{false}; /*!< \brief Enable implicit-lines agglomeration from walls. */ + bool MG_Linear_Prolongation{false}; /*!< \brief Prolong the correction with a limited least-squares + gradient instead of piecewise-constant injection. */ unsigned long MG_Startup_Iter{100}; /*!< \brief Iterations per mesh during FMG startup, and the length of each level's CFL ramp. 0 = no iteration budget. */ su2double MG_Startup_Convergence{-2.0}; /*!< \brief FMG: orders of magnitude (log10) that CONV_FIELD must drop on the active level before promoting to the next finer one. Negative is a @@ -2875,6 +2878,7 @@ enum class MPI_QUANTITIES { COORDINATES , /*!< \brief Vertex coordinates communication. */ COORDINATES_OLD , /*!< \brief Old vertex coordinates communication. */ MAX_LENGTH , /*!< \brief Maximum length communication. */ + WALL_DISTANCE , /*!< \brief Wall distance and roughness of the nearest wall communication. */ GRID_VELOCITY , /*!< \brief Grid velocity communication. */ SOLUTION_EDDY , /*!< \brief Turbulent solution plus eddy viscosity communication. */ STOCH_SOURCE_LANG , /*!< \brief Stochastic source term for Langevin equations communication. */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index fd4207b0c5fe..39d9e23b8f18 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2104,6 +2104,10 @@ void CConfig::SetConfig_Options() { /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Pave the coarse grid with advancing fronts raised from boundaries * that carry a stretched layer normal to themselves. DEFAULT: NO \ingroup Config*/ addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false); + /*!\brief MG_LINEAR_PROLONGATION\n DESCRIPTION: Prolong the multigrid correction with a limited least-squares + * gradient over the coarse control volume instead of injecting the parent value into every child. Coarse CVs on + * walls and symmetry planes keep the constant operator. DEFAULT: NO \ingroup Config*/ + addBoolOption("MG_LINEAR_PROLONGATION", MGOptions.MG_Linear_Prolongation, false); /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Max number of iterations spent on each mesh during the Full * Multigrid (FMG) startup phase. DEFAULT: 100 \ingroup Config*/ addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100); @@ -2112,8 +2116,9 @@ void CConfig::SetConfig_Options() { * DEFAULT: -2 \ingroup Config*/ addDoubleOption("MG_STARTUP_CONVERGENCE", MGOptions.MG_Startup_Convergence, -2.0); /*!\brief MG_STARTUP_STAGNATION\n DESCRIPTION: Full-MG promotion on stagnation. If the active level's residual ratio - * between successive iterations exceeds this value for MG_STARTUP_STAGNATION_ITER consecutive iterations, promote to - * the next finer level without waiting out MG_STARTUP_ITER. 0 disables it. DEFAULT: 0.99 \ingroup Config*/ + * between successive iterations stays between this value and its inverse for MG_STARTUP_STAGNATION_ITER consecutive + * iterations, promote to the next finer level without waiting out MG_STARTUP_ITER. 0 disables it. + * DEFAULT: 0.99 \ingroup Config*/ addDoubleOption("MG_STARTUP_STAGNATION", MGOptions.MG_Startup_Stagnation, 0.99); /*!\brief MG_STARTUP_STAGNATION_ITER\n DESCRIPTION: Consecutive stalled iterations required before Full-MG promotes * on stagnation. 0 disables it, as MG_STARTUP_STAGNATION= 0 does. DEFAULT: 5 \ingroup Config*/ @@ -4979,6 +4984,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Kappa_2nd_AdjFlow = jst_adj_coeff[0]; Kappa_4th_AdjFlow = jst_adj_coeff[1]; + /*--- The multigrid integration carries per-level arrays of this size. ---*/ + + if (nMGLevels > MAX_MGLEVELS) { + SU2_MPI::Error("MGLEVEL is larger than the supported maximum of " + std::to_string(MAX_MGLEVELS) + ".", + CURRENT_FUNCTION); + } + /*--- Fill MG smooth vectors to size nMGLevels+1. Use parsed values (truncating or extending by repeat) or defaults if not set. ---*/ diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 360da1aaaa1e..a9adc2cd8198 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -627,6 +627,10 @@ void CGeometry::GetCommCountAndType(const CConfig* config, MPI_QUANTITIES commTy COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE::DOUBLE; break; + case MPI_QUANTITIES::WALL_DISTANCE: + COUNT_PER_POINT = 2; + MPI_TYPE = COMM_TYPE::DOUBLE; + break; case MPI_QUANTITIES::NEIGHBORS: COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT; @@ -718,6 +722,10 @@ void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, MPI_QU case MPI_QUANTITIES::MAX_LENGTH: bufDSend[buf_offset] = nodes->GetMaxLength(iPoint); break; + case MPI_QUANTITIES::WALL_DISTANCE: + bufDSend[buf_offset] = nodes->GetWall_Distance(iPoint); + bufDSend[buf_offset + 1] = nodes->GetRoughnessHeight(iPoint); + break; case MPI_QUANTITIES::NEIGHBORS: bufSSend[buf_offset] = geometry->nodes->GetnNeighbor(iPoint); break; @@ -811,6 +819,10 @@ void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, MPI_QU case MPI_QUANTITIES::MAX_LENGTH: nodes->SetMaxLength(iPoint, bufDRecv[buf_offset]); break; + case MPI_QUANTITIES::WALL_DISTANCE: + nodes->SetWall_Distance(iPoint, bufDRecv[buf_offset]); + nodes->SetRoughnessHeight(iPoint, bufDRecv[buf_offset + 1]); + break; case MPI_QUANTITIES::NEIGHBORS: nodes->SetnNeighbor(iPoint, bufSRecv[buf_offset]); break; @@ -2741,6 +2753,16 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE); geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE); geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]); + + /*--- SetCoord centred a halo agglomerate on the partial child list this rank holds. Take + the owner's coordinate, otherwise coarse stencils depend on the partitioning. Deformation + runs this inside a parallel region, so only one thread may reach the exchange. ---*/ + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + geometry_container[iMesh]->InitiateComms(geometry_container[iMesh], config, MPI_QUANTITIES::COORDINATES); + geometry_container[iMesh]->CompleteComms(geometry_container[iMesh], config, MPI_QUANTITIES::COORDINATES); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS } /*--- Compute the global surface areas for all markers. ---*/ @@ -4538,6 +4560,47 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); return 2 * Vol / GeometryToolbox::Norm(3, Normal); } + +/*--- Volume-averages the wall distance and roughness of the children onto a coarse grid, then sets the + * nearest-neighbor distance of its viscous wall vertices. ---*/ +void RestrictWallDistance(const CGeometry* geo_fine, CGeometry* geo_coarse, const CConfig* config) { + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geo_coarse->GetnPointDomain(); iPoint++) { + su2double dist = 0.0, roughness = 0.0, vol = 0.0; + for (auto iChild = 0u; iChild < geo_coarse->nodes->GetnChildren_CV(iPoint); iChild++) { + const auto jPoint = geo_coarse->nodes->GetChildren_CV(iPoint, iChild); + const su2double volChild = geo_fine->nodes->GetVolume(jPoint); + dist += geo_fine->nodes->GetWall_Distance(jPoint) * volChild; + roughness += geo_fine->nodes->GetRoughnessHeight(jPoint) * volChild; + vol += volChild; + } + geo_coarse->nodes->SetWall_Distance(iPoint, (vol > 0.0) ? su2double(dist / vol) : su2double(0.0)); + geo_coarse->nodes->SetRoughnessHeight(iPoint, (vol > 0.0) ? su2double(roughness / vol) : su2double(0.0)); + } + END_SU2_OMP_FOR + + /*--- A halo agglomerate only holds the children on this rank, so take the owner's values. ---*/ + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + geo_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::WALL_DISTANCE); + geo_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::WALL_DISTANCE); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) { + const auto viscous = config->GetViscous_Wall(iMarker); + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + const auto iPoint = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + const su2double dist = (viscous && geo_coarse->nodes->GetDomain(iPoint)) + ? NearestNeighborDistance(geo_coarse, config, iPoint) + : geo_coarse->nodes->GetWall_Distance(iPoint); + geo_coarse->vertex[iMarker][iVertex]->SetNearestNeighborDistance(dist); + } + END_SU2_OMP_FOR + } +} } // namespace void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container, @@ -4638,6 +4701,13 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo } END_SU2_OMP_FOR } + + /*--- The Full-MG startup solves the turbulence model on coarse grids. ---*/ + + for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { + RestrictWallDistance(geometry_container[iZone][iInst][iMesh - 1], geometry_container[iZone][iInst][iMesh], + config); + } } } } diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index d3bca729bf5b..86b8c272871f 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1619,7 +1619,7 @@ CMultiGridGeometry::CFrontSeeds CMultiGridGeometry::SeedFrontNodes(const CGeomet /*--- Paving order: a no-slip wall first, a slip wall next, everything else last. ---*/ auto tierOfBC = [](unsigned short bc) -> char { if ((bc == HEAT_FLUX) || (bc == ISOTHERMAL) || (bc == CHT_WALL_INTERFACE) || (bc == SMOLUCHOWSKI_MAXWELL)) return 0; - return (bc == EULER_WALL) ? 1 : 2; + return ((bc == EULER_WALL) || (bc == SYMMETRY_PLANE)) ? 1 : 2; }; for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { @@ -2322,6 +2322,63 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co byLayer[iColumn][layerOf[iPoint]].push_back(iPoint); } + /*--- A boundary run with an odd number of seeds, and every node where two boundaries meet, + * leaves one column a single node wide. Two such columns side by side are emitted together, + * so the domain gets square cells instead of a file of narrow ones. ---*/ + vector partnerOf(nColumn, -1); + { + auto isNarrow = [&](unsigned long iColumn) { + return isSeeded[iColumn] && (byLayer[iColumn].size() > 1) && (byLayer[iColumn][0].size() == 1); + }; + + /*--- Both columns must stand shoulder to shoulder over every layer they share, otherwise the + * cells the merge makes are not compact. ---*/ + auto aligned = [&](unsigned long iColumn, unsigned long jColumn) { + const auto nCommon = std::min(byLayer[iColumn].size(), byLayer[jColumn].size()); + for (auto iLayer = 1ul; iLayer < nCommon; ++iLayer) { + const auto iPoint = byLayer[iColumn][iLayer].front(); + const auto jPoint = byLayer[jColumn][iLayer].front(); + bool touch = false; + for (auto kPoint : fine_grid->nodes->GetPoints(iPoint)) touch = touch || (kPoint == jPoint); + if (!touch) return false; + } + return nCommon > 1; + }; + + /*--- Taken in global index order, so the pairing does not depend on the partitioning. ---*/ + vector narrow; + for (auto iColumn = 0ul; iColumn < nColumn; ++iColumn) + if (isNarrow(iColumn)) narrow.push_back(iColumn); + std::sort(narrow.begin(), narrow.end(), [&](unsigned long a, unsigned long b) { + return fine_grid->nodes->GetGlobalIndex(byLayer[a][0].front()) < + fine_grid->nodes->GetGlobalIndex(byLayer[b][0].front()); + }); + + for (auto iColumn : narrow) { + if (partnerOf[iColumn] >= 0) continue; + auto best = NO_COLUMN; + auto bestKey = std::numeric_limits::max(); + + for (auto jPoint : fine_grid->nodes->GetPoints(byLayer[iColumn][1].front())) { + const auto jColumn = columnOf[jPoint]; + if ((jColumn == NO_COLUMN) || (jColumn == iColumn)) continue; + if ((partnerOf[jColumn] >= 0) || !isNarrow(jColumn)) continue; + if (tierOf[jColumn] != tierOf[iColumn]) continue; + if (!aligned(iColumn, jColumn)) continue; + + const auto key = fine_grid->nodes->GetGlobalIndex(byLayer[jColumn][0].front()); + if (key < bestKey) { + bestKey = key; + best = jColumn; + } + } + + if (best == NO_COLUMN) continue; + partnerOf[iColumn] = static_cast(best); + partnerOf[best] = static_cast(iColumn); + } + } + auto emitGroup = [&](const vector& group) { nodes->SetChildren_CV(Index_CoarseCV, group); for (auto iPoint : group) { @@ -2334,6 +2391,7 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co auto minDepth = std::numeric_limits::max(), maxDepth = 0ul; vector group; + vector> paired; /*--- Hand a set of nodes out as coarse CVs, each connected and within the size limit. ---*/ vector inSet(nPointFine, 0); @@ -2379,7 +2437,8 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co auto iLayer = 0ul; if (isSeeded[iColumn]) { - /*--- The boundary row is a coarse CV of its own, which fixes the footprint above it. ---*/ + /*--- The boundary row is a coarse CV of its own, which fixes the footprint above it. A + * paired column keeps its own row, which may hold a boundary condition of its own. ---*/ const auto baseCV = emitGroup(layers[0]); ct[P_CVS]++; ct[P_COVERED] += layers[0].size(); @@ -2387,19 +2446,35 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co iLayer = 1; } + /*--- Of a pair, the column of lower index emits what stands above both boundary rows. ---*/ + if ((partnerOf[iColumn] >= 0) && (partnerOf[iColumn] < static_cast(iColumn))) continue; + + const auto* emitted = &layers; + if (partnerOf[iColumn] > static_cast(iColumn)) { + const auto& other = byLayer[partnerOf[iColumn]]; + paired = layers; + for (auto k = 1ul; k < other.size(); ++k) { + if (k < paired.size()) + paired[k].insert(paired[k].end(), other[k].begin(), other[k].end()); + else + paired.push_back(other[k]); + } + emitted = &paired; + } + /*--- Above it, consecutive layers are blocked so the coarse cell coarsens by the same ratio * along the column as the patch does across it. ---*/ - while (iLayer < layers.size()) { + while (iLayer < emitted->size()) { group.clear(); - const auto block = BlockFor(maxAgglomSize, layers[iLayer].size()); - for (auto k = 0ul; (k < block) && (iLayer < layers.size()); ++k) { - if (!group.empty() && (group.size() + layers[iLayer].size() > static_cast(maxAgglomSize))) break; - group.insert(group.end(), layers[iLayer].begin(), layers[iLayer].end()); + const auto block = BlockFor(maxAgglomSize, (*emitted)[iLayer].size()); + for (auto k = 0ul; (k < block) && (iLayer < emitted->size()); ++k) { + if (!group.empty() && (group.size() + (*emitted)[iLayer].size() > static_cast(maxAgglomSize))) break; + group.insert(group.end(), (*emitted)[iLayer].begin(), (*emitted)[iLayer].end()); iLayer++; } /*--- A single layer wider than the limit still has to go somewhere. ---*/ if (group.empty()) { - group = layers[iLayer]; + group = (*emitted)[iLayer]; iLayer++; } emitConnected(group); diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index d73b4962128f..155f98523d38 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -26,6 +26,7 @@ */ #include "CIntegration.hpp" +#include "../../../Common/include/containers/container_decorators.hpp" /*! * \class CMultiGridIntegration @@ -190,7 +191,20 @@ class CMultiGridIntegration final : public CIntegration { * \param[in] config - Definition of the particular problem. */ void GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh); + + /*! + * \brief Least-squares gradient of the coarse-grid correction, limited so no child value + * leaves the range of the coarse stencil. + * \param[in] sol_coarse - Solver holding the correction in Solution_Old. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the coarse mesh. + */ + void ComputeProlongationGradient(CSolver *sol_coarse, CGeometry *geo_coarse, CGeometry *geo_fine, + const CConfig *config, unsigned short iMesh); /*! * \brief Do an implicit smoothing of the prolongated correction. @@ -319,7 +333,8 @@ class CMultiGridIntegration final : public CIntegration { passivedouble lastRMS[2], char& exitReason, passivedouble& worstStepRatio, unsigned short& worstStep); - static constexpr int MAX_MG_LEVELS = 10; + /*--- CConfig rejects a larger MGLEVEL, so the per-level arrays below always fit. ---*/ + static constexpr int MAX_MG_LEVELS = MAX_MGLEVELS; /*--- Bounds the fixed-size stack buffers in the restriction and prolongation kernels, * independently of the CSysMatrix limit of the same name. ---*/ @@ -367,6 +382,10 @@ class CMultiGridIntegration final : public CIntegration { enum class MGStartupPromote { NONE, BUDGET, CONVERGENCE, STAGNATION }; MGStartupPromote mg_startup_promote_reason = MGStartupPromote::NONE; + /*! \brief Limited least-squares gradient of the correction, indexed by coarse level. + * Allocated on first use, only when MG_LINEAR_PROLONGATION is on. */ + vector prolongGradient; + vector mg_startup_conv_start; /*!< \brief Field values when the active level became active. */ vector mg_startup_conv_prev; /*!< \brief Field values on the previous iteration. */ unsigned long mg_startup_stall_count = 0; /*!< \brief Consecutive iterations without useful reduction. */ diff --git a/SU2_CFD/include/variables/CAdjEulerVariable.hpp b/SU2_CFD/include/variables/CAdjEulerVariable.hpp index 1faf911816aa..cb039f291358 100644 --- a/SU2_CFD/include/variables/CAdjEulerVariable.hpp +++ b/SU2_CFD/include/variables/CAdjEulerVariable.hpp @@ -102,6 +102,11 @@ class CAdjEulerVariable : public CVariable { Solution_Old(iPoint,iDim+1) = val_velocity[iDim]*Solution(iPoint,0); } + /*! + * \brief Index of the adjoint momentum in the solution vector. + */ + inline short GetVelocityIndex() const final { return 1; } + /*! * \brief Set the momentum part of the truncation error to zero. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index abd63a3afd27..ebd817fa1339 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -298,6 +298,11 @@ class CEulerVariable : public CFlowVariable { Solution_Old(iPoint,iDim+1) = val_velocity[iDim]*Solution(iPoint,0); } + /*! + * \brief Index of the momentum in the conservative solution vector. + */ + inline short GetVelocityIndex() const final { return 1; } + /*! * \brief Set the momentum part of the truncation error to zero. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index aaa6c437b665..7fc1860a826a 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -213,6 +213,11 @@ class CIncEulerVariable : public CFlowVariable { Solution_Old(iPoint,iDim+1) = val_velocity[iDim]; } + /*! + * \brief Index of the velocity in the solution vector. + */ + inline short GetVelocityIndex() const final { return 1; } + /*! * \brief Set the momentum part of the truncation error to zero. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index e0d4199144c4..37229a595945 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -298,6 +298,12 @@ class CNEMOEulerVariable : public CFlowVariable { } } + /*! + * \brief Index of the momentum in the conservative solution vector, which the primitive + * index Velocity() does not give. + */ + inline short GetVelocityIndex() const final { return static_cast(nSpecies); } + /*! * \brief A virtual member. * \return Value of the vibrational-electronic temperature. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 739a1af75fa5..13224df85125 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1607,6 +1607,12 @@ class CVariable { */ inline virtual void SetVelocity_Old(unsigned long iPoint, const su2double *val_velocity) {} + /*! + * \brief Index of the first velocity or momentum component in the solution vector. + * \return The index, or -1 when the solution carries no velocity. + */ + inline virtual short GetVelocityIndex() const { return -1; } + /*! * \brief A virtual member. * \param[in] laminarViscosity diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 29fb325fa233..dbe6458b9aa2 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -881,10 +881,6 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { geometry[iMGlevel]->SetCoord(geometry[iMGlevel-1]); - /*--- Find closest, most normal, neighbor to a surface point ---*/ - - geometry[iMGlevel]->FindNormal_Neighbor(config); - /*--- Store our multigrid index. ---*/ geometry[iMGlevel]->SetMGLevel(iMGlevel); @@ -934,6 +930,18 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + /*--- SetCoord centred a halo agglomerate on the partial child list this rank holds. Take + the owner's coordinate, otherwise coarse stencils depend on the partitioning. ---*/ + + if (iMGlevel > MESH_0) { + geometry[iMGlevel]->InitiateComms(geometry[iMGlevel], config, MPI_QUANTITIES::COORDINATES); + geometry[iMGlevel]->CompleteComms(geometry[iMGlevel], config, MPI_QUANTITIES::COORDINATES); + + /*--- Find closest, most normal, neighbor to a surface point ---*/ + + geometry[iMGlevel]->FindNormal_Neighbor(config); + } + /*--- Compute the max length. ---*/ if (!fea) { diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 607a0915d3c8..b6908a389b92 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -28,6 +28,7 @@ #include "../../include/integration/CMultiGridIntegration.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include @@ -50,6 +51,53 @@ static su2double applyGlobalTrend(su2double factor, passivedouble crossCycleRati return max(su2double{CLAMP_MIN}, min(su2double{CLAMP_MAX}, factor)); } +/*!\cond PRIVATE + * Inverts the symmetric least-squares matrix. False means the stencil is too degenerate + * to define a gradient. + \endcond */ +static bool invertLeastSquaresMatrix(unsigned short nDim, const su2double A[3][3], su2double inv[3][3]) { + + /*--- Compare against the trace so the test does not depend on the mesh units. ---*/ + + su2double trace = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) trace += A[iDim][iDim]; + if (trace <= 0.0) return false; + + const su2double scale = trace / su2double(nDim); + + /*--- det/scale^nDim is the conditioning of the fit. At 1e-8 the inverse still carries about + * nine digits; below it the gradient is noise and the control volume stays constant. ---*/ + constexpr passivedouble REL_TOL = 1e-8; + + if (nDim == 2) { + const su2double det = A[0][0]*A[1][1] - A[0][1]*A[0][1]; + if (det <= REL_TOL * scale * scale) return false; + + inv[0][0] = A[1][1]/det; + inv[0][1] = -A[0][1]/det; + inv[1][0] = inv[0][1]; + inv[1][1] = A[0][0]/det; + return true; + } + + const su2double a = A[0][0], b = A[0][1], c = A[0][2]; + const su2double d = A[1][1], e = A[1][2], f = A[2][2]; + + const su2double det = a*(d*f - e*e) - b*(b*f - e*c) + c*(b*e - d*c); + if (det <= REL_TOL * scale * scale * scale) return false; + + inv[0][0] = (d*f - e*e)/det; + inv[0][1] = (c*e - b*f)/det; + inv[0][2] = (b*e - c*d)/det; + inv[1][1] = (a*f - c*c)/det; + inv[1][2] = (b*c - a*e)/det; + inv[2][2] = (a*d - b*b)/det; + inv[1][0] = inv[0][1]; + inv[2][0] = inv[0][2]; + inv[2][1] = inv[1][2]; + return true; +} + inline passivedouble ComputeLinSysResRMS(const CSolver* solver) { passivedouble result = 0; for (unsigned short iVar = 0; iVar < solver->GetnVar(); ++iVar) { @@ -72,7 +120,9 @@ void CMultiGridIntegration::adaptDampingFactors(CConfig* config, passivedouble c applyGlobalTrend(config->GetDamp_Correc_Prolong(), crossCycleRatio)); } -CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { } +CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { + prolongGradient.resize(MAX_MG_LEVELS + 1); +} void CMultiGridIntegration::MonitorFullMG_Startup(const vector >& convFields, const CConfig *config) { @@ -101,14 +151,15 @@ void CMultiGridIntegration::MonitorFullMG_Startup(const vector 0.0 && !mg_startup_conv_prev.empty()) { + const passivedouble band = fabs(log10(stall_tol)); bool stalled = true; for (auto iField = 0ul; iField < nFields; iField++) - stalled = stalled && (convFields[iField].second - mg_startup_conv_prev[iField] >= log10(stall_tol)); + stalled = stalled && (fabs(convFields[iField].second - mg_startup_conv_prev[iField]) <= band); if (stalled) mg_startup_stall_count++; @@ -451,6 +502,15 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, geometry[iZone][iInst][FinestMesh], config[iZone], true); } + /*--- The scalar solvers run on the active startup level and need its vorticity and strain rate. ---*/ + + if (fmg_warmup) { + solver_container[iZone][iInst][FinestMesh][Solver_Position]->Preprocessing(geometry[iZone][iInst][FinestMesh], + solver_container[iZone][iInst][FinestMesh], + config[iZone], FinestMesh, NO_RK_ITER, + RunTime_EqSystem, true); + } + /*--- Computes primitive variables and gradients in the finest mesh (useful for the next solver (turbulence) and output ---*/ solver_container[iZone][iInst][MESH_0][Solver_Position]->Preprocessing(geometry[iZone][iInst][MESH_0], @@ -670,7 +730,8 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, /*--- Compute prolongated solution, and smooth the correction $u^(new)_k = u_k + Smooth(I^k_(k+1)(u_(k+1)-I^(k+1)_k u_k))$ ---*/ - GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, + iMesh+1); const auto& mgOpts = config->GetMGOptions(); SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, mgOpts.MG_CorrecSmooth[iMesh], mgOpts.MG_Smooth_Coeff, config, iMesh); @@ -873,7 +934,8 @@ void CMultiGridIntegration::PostSmoothing(unsigned short RunTime_EqSystem, void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh) { SU2_ZONE_SCOPED const unsigned short nVar = sol_coarse->GetnVar(); @@ -941,13 +1003,203 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS /*--- Interpolate the coarse-grid correction onto the fine * grid and store in LinSysRes. ---*/ - /*--- Halos too: the correction smoother reads them before its first exchange. ---*/ - SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) - for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { + if (!config->GetMGOptions().MG_Linear_Prolongation) { + + /*--- Halos too: the correction smoother reads them before its first exchange. ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { + const auto* Correction = sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + sol_fine->LinSysRes.SetBlock(Point_Fine, Correction); + } + } + END_SU2_OMP_FOR + + return; + } + + /*--- The coarse coordinate is the volume-weighted centroid of the children, so the + * reconstruction below averages back to the coarse correction exactly. ---*/ + + ComputeProlongationGradient(sol_coarse, geo_coarse, geo_fine, config, iMesh); + + const unsigned short nDim = geo_coarse->GetnDim(); + const auto& gradient = prolongGradient[iMesh]; + + /*--- A halo agglomerate holds a partial child list, so only owners write, and the + * exchange below fills the fine halos the correction smoother reads. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + const auto* Correction = sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse); + const auto* Coord_Coarse = geo_coarse->nodes->GetCoord(Point_Coarse); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - sol_fine->LinSysRes.SetBlock(Point_Fine, Correction); + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_fine->nodes->GetCoord(Point_Fine), Coord_Coarse, Distance); + + su2double Solution[MAXNVAR] = {0.0}; + for (auto iVar = 0u; iVar < nVar; iVar++) { + Solution[iVar] = Correction[iVar]; + for (auto iDim = 0u; iDim < nDim; iDim++) + Solution[iVar] += gradient(Point_Coarse, iVar, iDim) * Distance[iDim]; + } + sol_fine->LinSysRes.SetBlock(Point_Fine, Solution); + } + } + END_SU2_OMP_FOR + + SU2_OMP_BARRIER + CSysMatrixComms::Initiate(sol_fine->LinSysRes, geo_fine, config); + CSysMatrixComms::Complete(sol_fine->LinSysRes, geo_fine, config); + +} + +void CMultiGridIntegration::ComputeProlongationGradient(CSolver *sol_coarse, CGeometry *geo_coarse, + CGeometry *geo_fine, const CConfig *config, + unsigned short iMesh) { + SU2_ZONE_SCOPED + + const unsigned short nVar = sol_coarse->GetnVar(); + const unsigned short nDim = geo_coarse->GetnDim(); + const auto nPointDomain = geo_coarse->GetnPointDomain(); + auto* nodes = sol_coarse->GetNodes(); + + /*--- Every thread of this rank takes this branch together, and the loops below would all be + * empty anyway. It also keeps the sizing test off a container that was never resized. ---*/ + if (nPointDomain == 0) return; + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + auto& grad = prolongGradient[iMesh]; + if ((grad.length() != nPointDomain) || (grad.rows() != nVar) || (grad.cols() != nDim)) + grad.resize(nPointDomain, nVar, nDim, 0.0); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + auto& gradient = prolongGradient[iMesh]; + + /*--- Unweighted fit: inverse-distance weights would let the short wall-normal edges of a + * stretched agglomerate dominate it. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain, omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) = 0.0; + + const auto* Coord_i = geo_coarse->nodes->GetCoord(iPoint); + const auto* Correction_i = nodes->GetSolution_Old(iPoint); + + su2double Amat[MAXNDIM][MAXNDIM] = {{0.0}}; + su2double rhs[MAXNVAR][MAXNDIM] = {{0.0}}; + + for (auto jPoint : geo_coarse->nodes->GetPoints(iPoint)) { + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_coarse->nodes->GetCoord(jPoint), Coord_i, Distance); + + for (auto iDim = 0u; iDim < nDim; iDim++) + for (auto jDim = 0u; jDim < nDim; jDim++) Amat[iDim][jDim] += Distance[iDim]*Distance[jDim]; + + const auto* Correction_j = nodes->GetSolution_Old(jPoint); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double delta = Correction_j[iVar] - Correction_i[iVar]; + for (auto iDim = 0u; iDim < nDim; iDim++) rhs[iVar][iDim] += Distance[iDim]*delta; + } + } + + su2double Ainv[MAXNDIM][MAXNDIM] = {{0.0}}; + if (!invertLeastSquaresMatrix(nDim, Amat, Ainv)) continue; + + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) { + su2double value = 0.0; + for (auto jDim = 0u; jDim < nDim; jDim++) value += Ainv[iDim][jDim]*rhs[iVar][jDim]; + gradient(iPoint, iVar, iDim) = value; + } + } + END_SU2_OMP_FOR + + /*--- The children of a wall agglomerate are all wall nodes, so their velocity correction is + * the constrained one the caller already imposed and only its gradient has to go. The + * remaining variables carry no wall condition and reconstruct along the surface. ---*/ + + const short iVel = nodes->GetVelocityIndex(); + + if (iVel >= 0) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + + const auto kindBC = config->GetMarker_All_KindBC(iMarker); + if (!config->GetViscous_Wall(iMarker) && (kindBC != EULER_WALL) && (kindBC != SYMMETRY_PLANE)) continue; + + SU2_OMP_FOR_STAT(32) + for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + const auto iPoint = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + if (iPoint >= nPointDomain) continue; + for (auto iVar = iVel; iVar < min(iVel + nDim, nVar); iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) = 0.0; + } + END_SU2_OMP_FOR + } + } + + /*--- Barth-Jespersen limiter over the children. One factor per variable scales the whole + * gradient, which is what keeps the children averaging back to the coarse value. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain, omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + + const auto* Coord_i = geo_coarse->nodes->GetCoord(iPoint); + const auto* Correction_i = nodes->GetSolution_Old(iPoint); + + su2double Correction_Min[MAXNVAR], Correction_Max[MAXNVAR]; + for (auto iVar = 0u; iVar < nVar; iVar++) { + Correction_Min[iVar] = Correction_i[iVar]; + Correction_Max[iVar] = Correction_i[iVar]; + } + + for (auto jPoint : geo_coarse->nodes->GetPoints(iPoint)) { + const auto* Correction_j = nodes->GetSolution_Old(jPoint); + for (auto iVar = 0u; iVar < nVar; iVar++) { + Correction_Min[iVar] = min(Correction_Min[iVar], Correction_j[iVar]); + Correction_Max[iVar] = max(Correction_Max[iVar], Correction_j[iVar]); + } + } + + su2double limiter[MAXNVAR]; + for (auto iVar = 0u; iVar < nVar; iVar++) limiter[iVar] = 1.0; + + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(iPoint); iChildren++) { + + const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(iPoint, iChildren); + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_fine->nodes->GetCoord(Point_Fine), Coord_i, Distance); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + + su2double delta = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) delta += gradient(iPoint, iVar, iDim)*Distance[iDim]; + + /*--- Scaled with the correction so the cut-off tracks its magnitude. ---*/ + const su2double tol = EPS*max(fabs(Correction_Max[iVar]), fabs(Correction_Min[iVar])) + EPS*EPS; + + if (delta > tol) + limiter[iVar] = min(limiter[iVar], (Correction_Max[iVar] - Correction_i[iVar])/delta); + else if (delta < -tol) + limiter[iVar] = min(limiter[iVar], (Correction_Min[iVar] - Correction_i[iVar])/delta); + } + } + + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double factor = max(su2double(0.0), min(su2double(1.0), limiter[iVar])); + for (auto iDim = 0u; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) *= factor; } } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index 29c71b2ed360..3fb340cda924 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -87,6 +87,13 @@ void CSingleGridIntegration::SingleGrid_Iteration(CGeometry ****geometry, CSolve Time_Integration(geometry_fine, solvers_fine, config[iZone], NO_RK_ITER, RunTime_EqSystem); + /*--- Coarse-level residuals are per-rank sums of squares unless reduced, which the smoothing + * early exit already does. ---*/ + + if ((FinestMesh != MESH_0) && !config[iZone]->GetMGOptions().MG_Smooth_EarlyExit) { + solvers_fine[Solver_Position]->SetResidual_RMS(geometry_fine, config[iZone], true); + } + /*--- Postprocessing ---*/ solvers_fine[Solver_Position]->Postprocessing(geometry_fine, solvers_fine, config[iZone], FinestMesh); diff --git a/config_template.cfg b/config_template.cfg index 88533d33c689..d26c2a90453b 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1766,6 +1766,12 @@ MG_DAMP_RESTRICTION= 0.5 % Damping factor for the correction prolongation MG_DAMP_PROLONGATION= 0.5 % +% Prolong the correction with a limited least-squares gradient over the coarse control +% volume instead of injecting the parent value into every child. The reconstruction is +% exact for a linear correction and the children still average back to the coarse value. +% Control volumes on walls and symmetry planes keep the constant operator (NO, YES) +MG_LINEAR_PROLONGATION= NO +% % Enable early exit from multigrid smoothing when the residual RMS drops % below MG_SMOOTH_RES_THRESHOLD * initial_rms (NO, YES) MG_SMOOTH_EARLY_EXIT= YES @@ -1816,8 +1822,9 @@ MG_STARTUP_ITER= 100 MG_STARTUP_CONVERGENCE= -2 % % Full-MG promotion on stagnation: if the active level's residual ratio between -% successive iterations exceeds this for MG_STARTUP_STAGNATION_ITER consecutive -% iterations, promote to the next finer level early. 0 disables (default 0.99). +% successive iterations stays between this value and its inverse for +% MG_STARTUP_STAGNATION_ITER consecutive iterations, promote to the next finer level +% early. A growing residual does not count as stagnation. 0 disables (default 0.99). MG_STARTUP_STAGNATION= 0.99 % % Consecutive stalled iterations required before Full-MG promotes on stagnation.