diff --git a/Common/include/geometry/CMultiGridGeometry.hpp b/Common/include/geometry/CMultiGridGeometry.hpp index ba635b40f8c..e1ba5a80fde 100644 --- a/Common/include/geometry/CMultiGridGeometry.hpp +++ b/Common/include/geometry/CMultiGridGeometry.hpp @@ -48,7 +48,16 @@ class CMultiGridGeometry final : public CGeometry { * \return TRUE or FALSE depending if the control volume can be agglomerated. */ bool SetBoundAgglomeration(unsigned long CVPoint, vector marker_seed, const CGeometry* fine_grid, - const CConfig* config) const; + const CConfig* config, const vector& mixedBC) const; + + /*! + * \brief Find nodes where two boundary conditions of different type meet. These are never + * agglomerated, since a coarse CV holding one would average both conditions. + * \param[in] fine_grid - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \return One flag per fine grid point, set where that point must stay on its own. + */ + vector FindMixedBoundaryNodes(const CGeometry* fine_grid, const CConfig* config) const; /*! * \brief Determine if a Point can be agglomerated using geometrical criteria. @@ -78,16 +87,71 @@ class CMultiGridGeometry final : public CGeometry { su2double ComputeLocalCurvature(const CGeometry* fine_grid, unsigned long iPoint, unsigned short iMarker) const; /*! - * \brief Agglomerate high-aspect-ratio interior cells along implicit lines from wall vertices. + * \brief Pave the domain with advancing fronts extruded from the boundary patches. * \param[in,out] Index_CoarseCV - Current coarse CV index, incremented as new coarse CVs are created. * \param[in] fine_grid - Fine grid geometry. * \param[in] config - Configuration. - * \param[in,out] MGQueue_InnerCV - Queue for domain agglomeration; processed points are removed. + * \param[in] iMesh - Multigrid level being built, used to label the summary. + * \return Summary of the paving, empty except on the master rank. + */ + string AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const CGeometry* fine_grid, const CConfig* config, + unsigned short iMesh); + + /*! + * \brief Weakest and strongest dual-grid coupling at each node, and the neighbour across the + * strongest edge. Their ratio is the local cell aspect ratio, available on every MG level. + */ + struct CNodeStiffness { + vector wMin, wMax; /*!< \brief Weakest and strongest edge coupling at each node. */ + vector jStiffest; /*!< \brief Neighbour across the strongest edge. */ + + /*!< \brief Local aspect ratio at a node, 1.0 where it could not be measured. */ + su2double AspectRatio(unsigned long iPoint) const { + return (wMin[iPoint] > 0.0) ? wMax[iPoint] / wMin[iPoint] : su2double(1.0); + } + }; + + /*! + * \brief Measure the dual-grid coupling at every node of a grid. + * \param[in] fine_grid - Grid to measure. + * \return Weakest/strongest coupling per node. + */ + CNodeStiffness ComputeNodeStiffness(const CGeometry* fine_grid) const; + + /*! + * \brief Boundary nodes that seed a front, with the direction each starts marching in. + */ + struct CFrontSeeds { + vector node; /*!< \brief Seed node on the boundary. */ + vector> normal; /*!< \brief Unit normal there, pointing into the domain. */ + }; + + /*! + * \brief Collect the boundary nodes that seed an advancing front: those on a viscous wall, or on a + * boundary carrying a stretched layer normal to itself. + * \param[in] fine_grid - Fine grid geometry. + * \param[in] config - Definition of the particular problem. + * \param[in] stiff - Node coupling from ComputeNodeStiffness. + * \return Seed nodes and their inward boundary normals. */ - void AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const CGeometry* fine_grid, const CConfig* config, - CMultiGridQueue& MGQueue_InnerCV); + CFrontSeeds SeedFrontNodes(const CGeometry* fine_grid, const CConfig* config, const CNodeStiffness& stiff) const; + + /*! + * \brief Partition the seed nodes into compact surface patches by repeated pairwise matching. Each + * patch is the footprint of one front and fixes the shape of the stack above it. + * \param[in] seeds - Seed nodes from SeedFrontNodes. + * \param[in] fine_grid - Fine grid geometry. + * \param[in] config - Definition of the particular problem. + * \param[in] mixedBC - Nodes that must stay on their own, from FindMixedBoundaryNodes. + * \return One vector of indices into seeds.node per patch. + */ + vector> BuildFrontPatches(const CFrontSeeds& seeds, const CGeometry* fine_grid, + const CConfig* config, const vector& mixedBC) const; public: + /*!< \brief Paving summary for this level. */ + string pavingReport; + /*--- This is to suppress Woverloaded-virtual, omitting it has no negative impact. ---*/ using CGeometry::SetBoundControlVolume; using CGeometry::SetControlVolume; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 3d5d57c1767..bb71aea8aea 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1116,7 +1116,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne struct CMGOptions { su2double MG_Smooth_Res_Threshold{0.0}; /*!< \brief RMS reduction threshold for MG smoothing early exit. */ su2double MG_Smooth_Coeff{0.0}; /*!< \brief Jacobi smoother coefficient for coarse-grid correction. */ - unsigned long MG_Min_MeshSize{0}; /*!< \brief Minimum CVs on coarsest MG level. */ + unsigned long MG_Min_MeshSize{0}; /*!< \brief Minimum CVs on coarsest MG level, per MPI rank. */ std::vector MG_PreSmooth; /*!< \brief Multigrid pre-smoothing iterations per level. */ std::vector MG_PostSmooth; /*!< \brief Multigrid post-smoothing iterations per level. */ std::vector MG_CorrecSmooth; /*!< \brief Multigrid Jacobi correction-smoothing per level. */ @@ -1125,7 +1125,6 @@ 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. */ - unsigned long MG_Implicit_Lines_MaxLength{20}; /*!< \brief Maximum nodes on a wall-normal implicit line (including wall seed). */ 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 diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 4a2f292bd61..29161a74338 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2061,16 +2061,16 @@ void CConfig::SetConfig_Options() { addDoubleOption("MG_SMOOTH_RES_THRESHOLD", MGOptions.MG_Smooth_Res_Threshold, 0.9); /*!\brief MG_SMOOTH_OUTPUT\n DESCRIPTION: Print compact per-cycle smoothing iteration summary. DEFAULT: NO \ingroup Config*/ addBoolOption("MG_SMOOTH_OUTPUT", MGOptions.MG_Smooth_Output, false); - /*!\brief MG_SMOOTH_STAGNATION_TOL\n DESCRIPTION: Stop smoothing if current_rms >= previous_rms * this value. Values < 1.0 enable early exit on stagnation, 1.0 only exits on defect growth. DEFAULT: 0.99 \ingroup Config*/ + /*!\brief MG_SMOOTH_STAGNATION_TOL\n DESCRIPTION: Stop smoothing if current_rms >= previous_rms * this value. + * Values < 1.0 enable early exit on stagnation, 1.0 only exits on defect growth. DEFAULT: 0.99 \ingroup Config*/ addDoubleOption("MG_SMOOTH_STAGNATION_TOL", MGOptions.MG_Smooth_StagnationTol, 0.99); /*!\brief MG_SMOOTH_COEFF\n DESCRIPTION: Smoothing coefficient for the correction prolongation Jacobi smoother. DEFAULT: 1.25 \ingroup Config*/ addDoubleOption("MG_SMOOTH_COEFF", MGOptions.MG_Smooth_Coeff, 1.25); - /*!\brief MG_MIN_MESHSIZE\n DESCRIPTION: Minimum number of CVs on the coarsest multigrid level. Levels that would produce fewer CVs are not created. DEFAULT: 50 \ingroup Config*/ + /*!\brief MG_MIN_MESHSIZE\n DESCRIPTION: Minimum number of CVs on the coarsest multigrid level, checked per MPI rank + * (i.e. on the smallest partition). Levels that would produce fewer CVs on any rank are not created. DEFAULT: 500 \ingroup Config*/ addUnsignedLongOption("MG_MIN_MESHSIZE", MGOptions.MG_Min_MeshSize, 500); /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Enable agglomeration along implicit lines from wall seeds. DEFAULT: NO \ingroup Config*/ addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false); - /*!\brief MG_IMPLICIT_LINES_MAX_LENGTH\n DESCRIPTION: Maximum number of nodes on a wall-normal implicit agglomeration line (including the wall seed node). DEFAULT: 20 \ingroup Config*/ - addUnsignedLongOption("MG_IMPLICIT_LINES_MAX_LENGTH", MGOptions.MG_Implicit_Lines_MaxLength, 20); /*!\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); @@ -2085,7 +2085,8 @@ void CConfig::SetConfig_Options() { /*!\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*/ addUnsignedLongOption("MG_STARTUP_STAGNATION_ITER", MGOptions.MG_Startup_Stagnation_Iter, 5); - /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ + /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). + * If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p); /*!\par CONFIG_CATEGORY: Spatial Discretization \ingroup Config*/ diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 684da742b13..bd1877de9a4 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -88,11 +88,35 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } + /*--- STEP 0: pave the domain with advancing fronts rising from the boundaries, wall CV included. + * It runs before the boundary agglomeration so that the wall CV and the layers above it share one + * footprint. Everything it claims is marked agglomerated, so the passes below skip it. The coarse + * CVs it creates occupy [firstLineCV, endLineCV), which is how the repair passes identify them. ---*/ + const auto firstLineCV = Index_CoarseCV; + if (config->GetMGOptions().MG_Implicit_Lines) { + pavingReport = AgglomerateImplicitLines(Index_CoarseCV, fine_grid, config, iMesh); + } + const auto endLineCV = Index_CoarseCV; + + /*--- Points carrying a physical boundary condition. This does not include SEND_RECEIVE. ---*/ + vector onPhysBoundary(fine_grid->GetnPoint(), 0); + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) + onPhysBoundary[fine_grid->vertex[iMarker][iVertex]->GetNode()] = 1; + } + /*--- Nodes where two different boundary conditions meet. ---*/ + const auto mixedBC = FindMixedBoundaryNodes(fine_grid, config); + /*--- STEP 1: The first step is the boundary agglomeration. ---*/ for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { /*--- Skip periodic boundaries: do not agglomerate on periodic markers. ---*/ if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) continue; + /*--- Skip SEND_RECEIVE markers: they record that a point is mirrored on another rank, not that it + * lies on a boundary. Such points are left to the domain pass, as in a serial run. ---*/ + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); @@ -118,13 +142,15 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un unsigned short copy_marker[3] = {}; marker_seed.push_back(iMarker); - /*--- For a particular point in the fine grid we save all the markers - that are in that point ---*/ + /*--- For a particular point in the fine grid we save all the physical markers that are in + that point. ---*/ for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker(); jMarker++) { - const string Marker_Tag = config->GetMarker_All_TagBound(iMarker); + if (config->GetMarker_All_KindBC(jMarker) == SEND_RECEIVE) continue; if (fine_grid->nodes->GetVertex(iPoint, jMarker) != -1) { - copy_marker[counter] = jMarker; + /*--- Count every physical marker, the counter > 2 test needs the true count, but only + store the first few, which is all the matching rules ever look at. ---*/ + if (counter < 3) copy_marker[counter] = jMarker; counter++; if (jMarker != iMarker) { @@ -163,23 +189,16 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un euler_wall_agglomerated[marker_seed[0]]++; } } - - /*--- Note that if the (single) marker is a SEND_RECEIVE, then the node is actually an interior point. - In that case it can only be agglomerated with another interior point. ---*/ - if (config->GetMarker_All_KindBC(marker_seed[0]) == SEND_RECEIVE) { - agglomerate_seed = true; - } } - /*--- Note that in 2D, this is a corner and we do not agglomerate unless one of them is SEND_RECEIVE. ---*/ - /*--- In 3D, we agglomerate if the 2 markers are the same. ---*/ + /*--- Two physical markers meet here. ---*/ if (counter == 2) { - if (nDim == 2) { - agglomerate_seed = ((config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) || - (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE)); - } - /*--- agglomerate if both markers are the same. ---*/ - if (nDim == 3) agglomerate_seed = (copy_marker[0] == copy_marker[1]); + /*--- In 2D that is a genuine corner in the geometry, which is never agglomerated. A wall + point merely split by a partition interface no longer reaches this branch: it counts + one physical marker and is handled above as the valley point it is. ---*/ + if (nDim == 2) agglomerate_seed = false; + /*--- In 3D, this is a ridge point (an edge where two surface markers meet). ---*/ + if (nDim == 3) agglomerate_seed = true; /*--- Euler walls: check curvature-based agglomeration criterion for both markers ---*/ // only in 3d because in 2d it's a corner @@ -208,6 +227,11 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un if (counter > 2) agglomerate_seed = false; + /*--- ...and so is a node where two markers of DIFFERENT type meet, whatever the count and + * whatever the dimension. A coarse CV holding such a node would average two conditions the + * fine grid applies separately. ---*/ + if (mixedBC[iPoint]) agglomerate_seed = false; + /*--- If the seed (parent) can be agglomerated, we try to agglomerate connected childs to the parent ---*/ /*--- Note that in 2D we allow a maximum of 4 nodes to be agglomerated ---*/ @@ -217,7 +241,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un for (auto CVPoint : fine_grid->nodes->GetPoints(iPoint)) { /*--- The new point can be agglomerated ---*/ - if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) { + if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config, mixedBC)) { /*--- We set the value of the parent ---*/ fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); @@ -231,12 +255,12 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- In 3D, we agglomerate exactly 2 nodes if the nodes are on the surface edge. ---*/ if ((nDim == 3) && (counter == 2)) break; /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes (counter==1 in 3D). ---*/ - if (nChildren == maxAgglomSize) break; + if (nChildren >= maxAgglomSize) break; } } - /*--- Only take into account indirect neighbors for 3D faces, not 2D. ---*/ - if (nDim == 3) { + /*--- Indirect neighbors only for 3D faces. ---*/ + if ((nDim == 3) && (nChildren < maxAgglomSize)) { Suitable_Indirect_Neighbors.clear(); if (fine_grid->nodes->GetAgglomerate_Indirect(iPoint)) @@ -247,7 +271,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un for (auto CVPoint : Suitable_Indirect_Neighbors) { /*--- The new point can be agglomerated ---*/ - if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) { + if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config, mixedBC)) { /*--- We set the value of the parent ---*/ fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); @@ -263,7 +287,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); nChildren++; /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes. ---*/ - if (nChildren == maxAgglomSize) break; + if (nChildren >= maxAgglomSize) break; } } } @@ -281,6 +305,11 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un i.e. make one coarse CV with a single child. ---*/ for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + /*--- As in STEP 1, a SEND_RECEIVE marker does not make a point a boundary point. Turning the + leftovers of those markers into single-child coarse CVs here would strand every interior point + along a partition interface before the domain pass below ever gets to see it. ---*/ + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); @@ -311,12 +340,19 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } - /*--- Agglomerate high-aspect-ratio interior nodes along implicit lines from walls. ---*/ - if (config->GetMGOptions().MG_Implicit_Lines) { - AgglomerateImplicitLines(Index_CoarseCV, fine_grid, config, MGQueue_InnerCV); - } + /*--- STEP 2: Agglomerate the domain points. A seed grows one node at a time, taking the candidate + * that shares the most edges with the current members. ---*/ + vector inCV(fine_grid->GetnPoint(), 0); + vector isCandidate(fine_grid->GetnPoint(), 0); + vector members, candidates; + members.reserve(maxAgglomSize); - /*--- STEP 2: Agglomerate the domain points. ---*/ + /*--- A local frame at the seed: up to nDim of its incident edges, as mutually orthogonal as possible, + * each with its own length. An offset is resolved onto these and each component divided by that + * direction's spacing. ---*/ + vector> frameDir, edgeDir; + vector frameLen, edgeLen; + vector edgeUsed; auto iteration = 0ul; while (!MGQueue_InnerCV.EmptyQueue() && (iteration < fine_grid->GetnPoint())) { @@ -328,79 +364,138 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un if ((!fine_grid->nodes->GetAgglomerate(iPoint)) && (fine_grid->nodes->GetDomain(iPoint)) && (GeometricalCheck(iPoint, fine_grid, config))) { - unsigned short nChildren = 1; - - /*--- We set an index for the parent control volume ---*/ - - fine_grid->nodes->SetParent_CV(iPoint, Index_CoarseCV); - - /*--- We add the seed point (child) to the parent control volume ---*/ - - nodes->SetChildren_CV(Index_CoarseCV, 0, iPoint); - - /*--- Update the queue with the seed point (remove the seed and - increase the priority of its neighbors) ---*/ - - MGQueue_InnerCV.Update(iPoint, fine_grid); - - /*--- Now we do a sweep over all the nodes that surround the seed point ---*/ - - for (auto CVPoint : fine_grid->nodes->GetPoints(iPoint)) { - /*--- Determine if the CVPoint can be agglomerated ---*/ + members.clear(); + candidates.clear(); + unsigned short nChildren = 0; + + /*--- Take a node into the CV and let the frontier grow with it. ---*/ + auto addMember = [&](unsigned long CVPoint) { + fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); + nChildren++; - if ((!fine_grid->nodes->GetAgglomerate(CVPoint)) && (fine_grid->nodes->GetDomain(CVPoint)) && - (GeometricalCheck(CVPoint, fine_grid, config))) { - /*--- We set the value of the parent ---*/ + if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint)) nodes->SetAgglomerate_Indirect(Index_CoarseCV, true); - fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); + /*--- Remove it from the queue and raise the priority of its neighbours. ---*/ + MGQueue_InnerCV.Update(CVPoint, fine_grid); - /*--- We set the value of the child ---*/ + members.push_back(CVPoint); + inCV[CVPoint] = 1; - nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); - nChildren++; + for (auto jPoint : fine_grid->nodes->GetPoints(CVPoint)) { + if (inCV[jPoint] || isCandidate[jPoint]) continue; + if (fine_grid->nodes->GetAgglomerate(jPoint) || !fine_grid->nodes->GetDomain(jPoint)) continue; + if (!GeometricalCheck(jPoint, fine_grid, config)) continue; + isCandidate[jPoint] = 1; + candidates.push_back(jPoint); + } + }; - /*--- Update the queue with the new control volume (remove the CV and - increase the priority of its neighbors) ---*/ + addMember(iPoint); - MGQueue_InnerCV.Update(CVPoint, fine_grid); + edgeDir.clear(); + edgeLen.clear(); + su2double seedShortest = std::numeric_limits::max(); + for (auto jPoint : fine_grid->nodes->GetPoints(iPoint)) { + su2double e[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(jPoint), fine_grid->nodes->GetCoord(iPoint), e); + const su2double len = GeometryToolbox::Norm(nDim, e); + if (len <= 0.0) continue; + std::array u{}; + for (unsigned short d = 0; d < nDim; ++d) u[d] = e[d] / len; + edgeDir.push_back(u); + edgeLen.push_back(len); + seedShortest = std::min(seedShortest, len); + } + /*--- A seed with no usable edge cannot be measured against anything, leave the scale at one. ---*/ + if (seedShortest == std::numeric_limits::max()) seedShortest = 1.0; + + /*--- The shortest edge goes in first, so that the direction the mesh is stretched along is the + * one measured against its own small spacing; the rest are taken in order of how orthogonal + * they are to what is already in the frame. ---*/ + frameDir.clear(); + frameLen.clear(); + edgeUsed.assign(edgeDir.size(), 0); + while (frameDir.size() < nDim) { + long pick = -1; + su2double bestScore = -1.0; + for (size_t i = 0; i < edgeDir.size(); ++i) { + if (edgeUsed[i]) continue; + su2double score; + if (frameDir.empty()) { + score = 1.0 / edgeLen[i]; + } else { + score = 2.0; + for (const auto& f : frameDir) { + su2double a = 0.0; + for (unsigned short d = 0; d < nDim; ++d) a += f[d] * edgeDir[i][d]; + score = std::min(score, 1.0 - fabs(a)); + } + } + if (score > bestScore) { + bestScore = score; + pick = static_cast(i); + } } - if (nChildren == maxAgglomSize) break; + if (pick < 0) break; + edgeUsed[pick] = 1; + frameDir.push_back(edgeDir[pick]); + frameLen.push_back(edgeLen[pick]); } + /*--- Without a full frame the offset cannot be resolved, fall back to one isotropic scale. ---*/ + const bool haveFrame = (frameDir.size() == nDim); + + while (nChildren < maxAgglomSize) { + /*--- Centroid of what the CV holds so far, used only to break ties. ---*/ + su2double centroid[MAXNDIM] = {0.0}; + for (auto jPoint : members) { + const auto* coord = fine_grid->nodes->GetCoord(jPoint); + for (auto iDim = 0u; iDim < nDim; iDim++) centroid[iDim] += coord[iDim] / su2double(members.size()); + } - /*--- Identify the indirect neighbors ---*/ - - Suitable_Indirect_Neighbors.clear(); - if (fine_grid->nodes->GetAgglomerate_Indirect(iPoint)) - SetSuitableNeighbors(Suitable_Indirect_Neighbors, iPoint, Index_CoarseCV, fine_grid); - - /*--- Now we do a sweep over all the indirect nodes that can be added ---*/ - - for (auto CVPoint : Suitable_Indirect_Neighbors) { - // if we have reached the maximum, get out. - if (nChildren == maxAgglomSize) break; - /*--- The new point can be agglomerated ---*/ - - if ((!fine_grid->nodes->GetAgglomerate(CVPoint)) && (fine_grid->nodes->GetDomain(CVPoint))) { - /*--- We set the value of the parent ---*/ - - fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); - - /*--- We set the indirect agglomeration information ---*/ - - if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint)) nodes->SetAgglomerate_Indirect(Index_CoarseCV, true); - - /*--- We set the value of the child ---*/ - - nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); - nChildren++; - - /*--- Update the queue with the new control volume (remove the CV and - increase the priority of the neighbors) ---*/ + unsigned long best = std::numeric_limits::max(); + unsigned short best_shared = 0; + su2double best_dist = std::numeric_limits::max(); + + for (auto CVPoint : candidates) { + if (inCV[CVPoint]) continue; + + unsigned short shared = 0; + for (auto jPoint : fine_grid->nodes->GetPoints(CVPoint)) shared += inCV[jPoint]; + + /*--- Distance to the centroid in cells rather than in metres: the offset is scaled by the + * seed edge pointing most nearly along it. ---*/ + su2double off[MAXNDIM] = {0.0}; + for (unsigned short d = 0; d < nDim; ++d) off[d] = fine_grid->nodes->GetCoord(CVPoint)[d] - centroid[d]; + const su2double offLen = GeometryToolbox::Norm(nDim, off); + + su2double dist = 0.0; + if (haveFrame) { + for (size_t k = 0; k < frameDir.size(); ++k) { + su2double p = 0.0; + for (unsigned short d = 0; d < nDim; ++d) p += off[d] * frameDir[k][d]; + const su2double q = p / frameLen[k]; + dist += q * q; + } + } else { + const su2double r = offLen / seedShortest; + dist = r * r; + } - MGQueue_InnerCV.Update(CVPoint, fine_grid); + if ((shared > best_shared) || ((shared == best_shared) && (dist < best_dist))) { + best = CVPoint; + best_shared = shared; + best_dist = dist; + } } + + if (best == std::numeric_limits::max()) break; + addMember(best); } + for (auto jPoint : members) inCV[jPoint] = 0; + for (auto jPoint : candidates) isCandidate[jPoint] = 0; + /*--- Update the number of control of childrens ---*/ nodes->SetnChildren_CV(Index_CoarseCV, nChildren); @@ -433,11 +528,90 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un SetPoint_Connectivity(fine_grid); + /*--- The connectivity just built only knows about coarse CVs of this rank: the halo CVs do not + exist until the MPI relay below runs, and the relay cannot run earlier because it broadcasts the + parent indices that the merge here is still free to change. So a CV touching a partition boundary + may look isolated while actually having neighbors on the other rank, and merging it would be + wrong. Mark those CVs and leave them alone; a genuinely isolated CV in the interior is unaffected. + Note this deliberately keeps the conservative outcome the sentinel used to produce by accident, + but only for the CVs that really do border another rank rather than for every CV near one. ---*/ + + /*--- Coarse CVs to leave exactly as the agglomeration made them: those holding a node where two + * different boundary conditions meet. Both repair passes remove one-child CVs, and a deliberately + * isolated junction is one, so they must be protected as a merge TARGET as well as a source. ---*/ + vector mustStayAlone(nPointDomain, false); + /*--- ...and which coarse CVs hold a boundary node at all. A boundary node is never agglomerated with + * an interior one, which the paving already respects. The repair passes below would undo it from + * the other end, since a boundary CV's neighbours include the interior CV sitting on top of it. ---*/ + vector cvOnBoundary(nPointDomain, false); + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + const auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + if (mixedBC[iFinePoint]) mustStayAlone[iCoarsePoint] = true; + if (onPhysBoundary[iFinePoint]) cvOnBoundary[iCoarsePoint] = true; + } + + /*--- Which physical boundaries each coarse CV sits on, one bit per marker. cvOnBoundary records only + * THAT a CV touches a boundary, so a one-child CV on boundary A could be merged into a neighbour on + * boundary B; SetVertex would then give the merged CV marker A while its body sits on B. Comparing + * marker sets is strictly stronger and keeps a merge inside one boundary. ---*/ + vector cvMarkerMask(nPointDomain, 0); + { + vector bitOfMarker(fine_grid->GetnMarker(), -1); + int nBits = 0; + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + if (nBits < 64) bitOfMarker[iMarker] = nBits; + nBits++; + } + if (nBits <= 64) { + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + const auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (bitOfMarker[iMarker] < 0) continue; + if (fine_grid->nodes->GetVertex(iFinePoint, iMarker) >= 0) + cvMarkerMask[iCoarsePoint] |= 1ULL << bitOfMarker[iMarker]; + } + } + } else { + /*--- More markers than bits: fall back to the boolean test. ---*/ + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) + cvMarkerMask[iCoarsePoint] = cvOnBoundary[iCoarsePoint] ? 1ULL : 0ULL; + } + } + + /*--- A boundary CV built by the paving is the base of a stack and must keep its footprint. Merging + * one into a neighbour leaves the column headless and widens the base, which on the next level + * stops it seeding a front and lets an interior front swallow it mid-stack. The repair passes + * below are for the interior singleton left where a line narrows, not for these. ---*/ + auto isStackBase = [&](unsigned long iCoarsePoint) { + return cvOnBoundary[iCoarsePoint] && (iCoarsePoint >= firstLineCV) && (iCoarsePoint < endLineCV); + }; + + vector touchesPartition(nPointDomain, false); + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + const auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + for (auto iFinePoint_Neighbor : fine_grid->nodes->GetPoints(iFinePoint)) { + if (fine_grid->nodes->GetParent_CV(iFinePoint_Neighbor) == std::numeric_limits::max()) { + touchesPartition[iCoarsePoint] = true; + break; + } + } + if (touchesPartition[iCoarsePoint]) break; + } + } + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { - if (nodes->GetnPoint(iCoarsePoint) == 1) { + if (mustStayAlone[iCoarsePoint]) continue; + if ((nodes->GetnPoint(iCoarsePoint) == 1) && !touchesPartition[iCoarsePoint]) { /*--- Find the neighbor of the isolated point. This neighbor is the right control volume ---*/ const auto iCoarsePoint_Complete = nodes->GetPoint(iCoarsePoint, 0); + if (mustStayAlone[iCoarsePoint_Complete]) continue; + if (isStackBase(iCoarsePoint) || isStackBase(iCoarsePoint_Complete)) continue; + if (cvMarkerMask[iCoarsePoint] != cvMarkerMask[iCoarsePoint_Complete]) continue; /*--- Check if merging would exceed the maximum agglomeration size ---*/ auto nChildren_Target = nodes->GetnChildren_CV(iCoarsePoint_Complete); @@ -451,6 +625,9 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un for (auto jCoarsePoint : nodes->GetPoints(iCoarsePoint_Complete)) { if (nChildrenToRedistribute == 0) break; + if (mustStayAlone[jCoarsePoint]) continue; + if (isStackBase(jCoarsePoint)) continue; + if (cvMarkerMask[jCoarsePoint] != cvMarkerMask[iCoarsePoint_Complete]) continue; auto nChildren_Neighbor = nodes->GetnChildren_CV(jCoarsePoint); if (nChildren_Neighbor < maxAgglomSize) { @@ -501,6 +678,114 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } + /*--- The pass above only rescues a coarse CV that has a single coarse neighbor, i.e. one enclosed + entirely within another. It does nothing for a CV that was left with a single fine-grid child but + sits between several coarse neighbors, which happens routinely with the implicit-line stacks above: + a bundle whose members reach the top of the boundary layer at slightly different heights (graded + spacing, local curvature) narrows one line at a time (see AgglomerateImplicitLines, Phase C), and + the fine node just past where a line dropped out is left to fend for itself. It usually has no + unclaimed neighbor left to agglomerate with by the time the ordinary domain pass (STEP 2) reaches + it, so it becomes a coarse CV of its own, surrounded on multiple sides by CVs it cannot join under + the single-neighbor rule above. Left alone, this is exactly the "hole" that turns what should be a + compact block into an L shape with an isolated singleton sitting in the missing corner. Merge such + a CV into whichever neighbor currently has the fewest children (and still has room under + maxAgglomSize): the smallest neighbor is the one most likely to be the under-filled block the + singleton belongs with, e.g. the 3-node L that this merge completes into a 4-node square, and + folding into it keeps the coarse grid from accumulating disproportionately large CVs. + + Unlike the pass above this one deliberately does not skip CVs that border another rank. That guard + exists because "has exactly one coarse neighbour" is measured on connectivity this rank cannot yet + see in full, so it misfires on a CV whose other neighbours merely live across the partition. The + trigger here is a child count, which is complete locally: a CV owning one fine point owns one fine + point no matter how the mesh was cut. The merge is local too - every candidate neighbour comes from + the connectivity built above and is therefore an owned CV, every child is an owned point, and the + MPI relay below broadcasts the result afterwards, so the ranks stay consistent. Requiring at least + two local neighbours below also rules out the one case that could not be repaired locally, a CV + whose neighbours are all halo: a halo CV mirrors another rank's decision and must not be added to + here. Keeping the guard cost 455 of 473 unrepaired singletons on a four-rank 3D bump. ---*/ + + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { + if (nodes->GetnChildren_CV(iCoarsePoint) != 1) continue; + if (mustStayAlone[iCoarsePoint]) continue; + if (isStackBase(iCoarsePoint)) continue; + if (nodes->GetnPoint(iCoarsePoint) <= 1) continue; /*--- Already handled above, or truly islanded. ---*/ + + /*--- The smallest neighbour is the one most likely to be an under-filled block, and taking the + smallest also means a neighbour with room is always preferred over a full one: anything below + maxAgglomSize necessarily has fewer children than anything at it. When every neighbour is full + the smallest is taken anyway, one child over the limit. maxAgglomSize is a quality target rather + than a capacity - Children_CV grows on demand - and a CV of nine is a far smaller blemish on the + coarse grid than the CV of one it removes. ---*/ + unsigned long best_neighbor = std::numeric_limits::max(); + unsigned short best_nChildren = 0; + for (auto jCoarsePoint : nodes->GetPoints(iCoarsePoint)) { + if (mustStayAlone[jCoarsePoint]) continue; + if (isStackBase(jCoarsePoint)) continue; + if (cvMarkerMask[jCoarsePoint] != cvMarkerMask[iCoarsePoint]) continue; + const auto nChildren_Neighbor = nodes->GetnChildren_CV(jCoarsePoint); + /*--- Skip neighbors already emptied by an earlier merge in this same pass. ---*/ + if (nChildren_Neighbor == 0) continue; + if ((best_neighbor == std::numeric_limits::max()) || (nChildren_Neighbor < best_nChildren)) { + best_nChildren = nChildren_Neighbor; + best_neighbor = jCoarsePoint; + } + } + if (best_neighbor == std::numeric_limits::max()) continue; /*--- Every neighbor was emptied. ---*/ + + const auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, 0); + nodes->SetChildren_CV(best_neighbor, best_nChildren, iFinePoint); + nodes->SetnChildren_CV(best_neighbor, best_nChildren + 1); + fine_grid->nodes->SetParent_CV(iFinePoint, best_neighbor); + nodes->SetnChildren_CV(iCoarsePoint, 0); + } + + /*--- Compact the coarse numbering. Both repair passes empty a control volume by moving its children + elsewhere, but nPointDomain was fixed before them, so an emptied index survives as a control volume + that owns no fine points at all. It ends up with EPS volume, (0,0,0) coordinates and, because + neighbours are derived through children, no neighbours either. Nothing is then adjacent to it on the + next level, so it cannot be agglomerated there and becomes a one-child CV again, and so on down the + hierarchy: the repair would otherwise trade a bad CV on this level for a degenerate one on every + level above. Squeezing the empty slots out keeps the CV count honest and stops that propagation. + + Renumbering is safe here because only three things refer to coarse indices at this point - the + children lists, the indirect-agglomeration flags, and the fine grid's parent indices - and the MPI + relay below has not run yet, so no other rank has seen these numbers. Only parents of *owned* fine + points are remapped; halo points are assigned by the relay afterwards. Every surviving parent keeps + a valid new index because a CV that still has a child is never removed. ---*/ + + { + constexpr auto NO_INDEX = std::numeric_limits::max(); + vector newIndex(nPointDomain, NO_INDEX); + unsigned long nKept = 0; + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) + if (nodes->GetnChildren_CV(iCoarsePoint) > 0) newIndex[iCoarsePoint] = nKept++; + + if (nKept < nPointDomain) { + /*--- Move the survivors down. newIndex[i] <= i, so a forward sweep never lands on a CV that has + not been moved yet. ---*/ + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { + const auto iNew = newIndex[iCoarsePoint]; + if ((iNew == NO_INDEX) || (iNew == iCoarsePoint)) continue; + const auto nChildren = nodes->GetnChildren_CV(iCoarsePoint); + for (auto iChildren = 0u; iChildren < nChildren; iChildren++) + nodes->SetChildren_CV(iNew, iChildren, nodes->GetChildren_CV(iCoarsePoint, iChildren)); + nodes->SetnChildren_CV(iNew, nChildren); + nodes->SetAgglomerate_Indirect(iNew, nodes->GetAgglomerate_Indirect(iCoarsePoint)); + } + for (auto iCoarsePoint = nKept; iCoarsePoint < nPointDomain; iCoarsePoint++) + nodes->SetnChildren_CV(iCoarsePoint, 0); + + for (auto iFinePoint = 0ul; iFinePoint < fine_grid->GetnPointDomain(); iFinePoint++) { + const auto iParent = fine_grid->nodes->GetParent_CV(iFinePoint); + if (iParent != NO_INDEX) fine_grid->nodes->SetParent_CV(iFinePoint, newIndex[iParent]); + } + + nPointDomain = nKept; + nPoint = nKept; + Index_CoarseCV = nKept; + } + } + /*--- Reset the neighbor information. ---*/ nodes->ResetPoints(); @@ -651,28 +936,31 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- Console output with the summary of the agglomeration ---*/ unsigned long nPointFine = fine_grid->GetnPointDomain(); - unsigned long Global_nPointCoarse, Global_nPointFine; + unsigned long Global_nPointCoarse, Global_nPointFine, Min_nPointCoarse; SU2_MPI::Allreduce(&nPointDomain, &Global_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); SU2_MPI::Allreduce(&nPointFine, &Global_nPointFine, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&nPointDomain, &Min_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_MIN, SU2_MPI::GetComm()); SetGlobal_nPointDomain(Global_nPointCoarse); if (iMesh != MESH_0) { - /*--- Note: CFL at the coarse levels have a large impact on convergence, - this should be rewritten to use adaptive CFL. ---*/ - const su2double Coeff = 1.5; - const su2double CFL = config->GetCFL(iMesh - 1) / Coeff; - config->SetCFL(iMesh, CFL); + /*--- Initialize coarse-level CFL from config. MG_CFL_SCALING will + apply per-level reductions during the multigrid cycle. ---*/ + config->SetCFL(iMesh, config->GetCFL(MESH_0)); } const su2double ratio = su2double(Global_nPointFine) / su2double(Global_nPointCoarse); - if (Global_nPointCoarse < config->GetMGOptions().MG_Min_MeshSize) { + /*--- Stop coarsening once the smallest per-rank partition falls below the minimum, + not just the summed total, since each rank runs its own MG hierarchy locally + and a partition that agglomerates down to too few CVs degenerates the operator + on that rank even if other ranks still have plenty of points. ---*/ + if (Min_nPointCoarse < config->GetMGOptions().MG_Min_MeshSize) { if (rank == MASTER_NODE) - cout << "MG level " << iMesh << " has only " << Global_nPointCoarse - << " CVs (< MG_MIN_MESHSIZE=" << config->GetMGOptions().MG_Min_MeshSize << "). Reducing MG levels to " - << iMesh - 1 << "." << endl; + cout << "MG level " << iMesh << " has only " << Min_nPointCoarse + << " CVs on the smallest partition (< MG_MIN_MESHSIZE=" << config->GetMGOptions().MG_Min_MeshSize + << "). Reducing MG levels to " << iMesh - 1 << "." << endl; config->SetMGLevels(iMesh - 1); } else if (rank == MASTER_NODE) { PrintingToolbox::CTablePrinter MGTable(&std::cout); @@ -751,8 +1039,35 @@ bool CMultiGridGeometry::GeometricalCheck(unsigned long iPoint, const CGeometry* return (Volume); } +vector CMultiGridGeometry::FindMixedBoundaryNodes(const CGeometry* fine_grid, const CConfig* config) const { + vector mixed(fine_grid->GetnPoint(), 0); + + /*--- The first physical condition seen at each node, -1 until one is. A second one of a different + * KIND is what makes the node mixed. A second marker of the same kind - two wall patches meeting + * - is not, and neither is SEND_RECEIVE, which records where the partition runs and says nothing + * about the boundary condition. ---*/ + vector firstBC(fine_grid->GetnPoint(), -1); + + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + const auto bc = static_cast(config->GetMarker_All_KindBC(iMarker)); + for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { + const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); + if (firstBC[iPoint] < 0) + firstBC[iPoint] = bc; + else if (firstBC[iPoint] != bc) + mixed[iPoint] = 1; + } + } + return mixed; +} + bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, vector marker_seed, - const CGeometry* fine_grid, const CConfig* config) const { + const CGeometry* fine_grid, const CConfig* config, + const vector& mixedBC) const { + /*--- A node where two boundary conditions of different type meet is never merged with anything. ---*/ + if (mixedBC[CVPoint]) return false; + bool agglomerate_CV = false; /*--- Basic condition, the point has not been previously agglomerated, it belongs to the domain, @@ -766,18 +1081,22 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, vectornodes->GetBoundary(CVPoint)) { - /*--- Identify the markers of the vertex that we want to agglomerate ---*/ - - // count number of markers on the agglomeration candidate - for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker() && counter < 3; jMarker++) { + /*--- Identify the physical markers of the vertex that we want to agglomerate. SEND_RECEIVE + markers are skipped for the same reason as on the seed side: they say nothing about the + boundary condition the candidate carries, only that it is mirrored on another rank. A + candidate whose markers are all SEND_RECEIVE therefore ends up with counter == 0 and is + rejected below, which is the answer a serial run gives for the interior point it really is. ---*/ + + for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker(); jMarker++) { + if (config->GetMarker_All_KindBC(jMarker) == SEND_RECEIVE) continue; if (fine_grid->nodes->GetVertex(CVPoint, jMarker) != -1) { - copy_marker[counter] = jMarker; + if (counter < 3) copy_marker[counter] = jMarker; counter++; } } - /*--- The basic condition is that the agglomerated vertex must have the same physical marker, - but eventually a send-receive condition ---*/ + /*--- The basic condition is that the agglomerated vertex must have the same physical marker + as the seed. ---*/ /*--- Only one marker in the vertex that is going to be agglomerated ---*/ @@ -786,17 +1105,11 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, vectorGetMarker_All_KindBC(other_marker) == SEND_RECEIVE) { - agglomerate_CV = true; - } - } } /*--- If there are two markers in the vertex that is going to be aglomerated ---*/ @@ -806,7 +1119,7 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, vectornodes->GetPoints(iFinePoint)) { const auto iParent = fine_grid->nodes->GetParent_CV(iFinePoint_Neighbor); + /*--- Skip neighbors whose parent is not known yet. The first call to this function happens + during construction, before the MPI relay has assigned parents to the fine grid's halo + points, so those still hold the sentinel. Letting it through would add one fake neighbor to + every coarse CV along a partition boundary, which both corrupts nNeighbor and hides the CV + from the isolated-CV repair. The driver calls this again once the relay has run. ---*/ + if (iParent == std::numeric_limits::max()) continue; /*--- If it is not the target coarse point, it is a coarse neighbor. ---*/ if (iParent != iCoarsePoint) { /*--- Avoid duplicates. ---*/ @@ -1280,232 +1599,903 @@ su2double CMultiGridGeometry::ComputeLocalCurvature(const CGeometry* fine_grid, return max_angle; } -void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const CGeometry* fine_grid, - const CConfig* config, CMultiGridQueue& MGQueue_InnerCV) { - /*--- Parameters ---*/ - const su2double ANGLE_THRESHOLD_DEG = 20.0; /*!< Stop line if direction deviates more than this. */ - const unsigned long MAX_LINE_LENGTH = config->GetMGOptions().MG_Implicit_Lines_MaxLength; - const su2double cos_threshold = cos(ANGLE_THRESHOLD_DEG * PI_NUMBER / 180.0); +CMultiGridGeometry::CNodeStiffness CMultiGridGeometry::ComputeNodeStiffness(const CGeometry* fine_grid) const { + /*--- Coupling across the dual face between a node and a neighbour. For a cell of size dx by dy this + * is 1/dy across the wall-normal face and 1/dx across the tangential one, so the ratio of largest + * to smallest weight at a node is the local aspect ratio. Unlike CGeometry::Aspect_Ratio this is + * available on every multigrid level, since SetControlVolume builds the dual grid on each. ---*/ + const auto nPointFine = fine_grid->GetnPoint(); + + CNodeStiffness stiff; + stiff.wMin.assign(nPointFine, 0.0); + stiff.wMax.assign(nPointFine, 0.0); + stiff.jStiffest.assign(nPointFine, std::numeric_limits::max()); + + for (auto iPoint = 0ul; iPoint < nPointFine; ++iPoint) { + su2double wmin = std::numeric_limits::max(), wmax = 0.0; + auto jStiffest = std::numeric_limits::max(); + + for (auto iNeigh = 0u; iNeigh < fine_grid->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = fine_grid->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = fine_grid->nodes->GetEdge(iPoint, iNeigh); + const su2double area = GeometryToolbox::Norm(nDim, fine_grid->edges->GetNormal(iEdge)); + const su2double w = + 0.5 * area * (1.0 / fine_grid->nodes->GetVolume(iPoint) + 1.0 / fine_grid->nodes->GetVolume(jPoint)); + if (w > wmax) { + wmax = w; + jStiffest = jPoint; + } + wmin = std::min(wmin, w); + } + + /*--- A node with no neighbours keeps the zeroed defaults, so AspectRatio reads 1. ---*/ + if (jStiffest != std::numeric_limits::max()) { + stiff.wMin[iPoint] = wmin; + stiff.wMax[iPoint] = wmax; + stiff.jStiffest[iPoint] = jStiffest; + } + } + return stiff; +} - const unsigned long nPointFine = fine_grid->GetnPoint(); +namespace { + +/*--- Unit normal of a boundary at a vertex, false if the marker does not reach iPoint. Boundary + * normals point INTO the domain, so this doubles as a starting direction to march in. ---*/ +bool VertexUnitNormal(const CGeometry* grid, unsigned short nDim, unsigned long iPoint, unsigned short iMarker, + su2double* unitNormal) { + const long iVertex = grid->nodes->GetVertex(iPoint, iMarker); + if (iVertex == -1) return false; + grid->vertex[iMarker][iVertex]->GetNormal(unitNormal); + const su2double nrm = GeometryToolbox::Norm(nDim, unitNormal); + if (nrm <= 0.0) return false; + for (unsigned short d = 0; d < nDim; ++d) unitNormal[d] /= nrm; + return true; +} - /*--- Collect implicit lines starting at viscous (no-slip) wall vertices only. - * Seeding from non-wall boundaries (farfield, inlet, outlet, symmetry) would - * claim interior BL cells before the wall lines can reach them, leaving - * wall-seeded lines with length < 3 (discarded). Restricting to viscous walls - * ensures the boundary-layer cells are agglomerated wall-first. - * Each line: [wall_node, interior_1, interior_2, ...]. - * The wall node (index 0) is already agglomerated by boundary agglomeration; - * only interior nodes (index >= 1) are paired into coarse CVs. ---*/ - vector> lines; +/*--- Are a and b neighbours in the fine grid? ---*/ +bool IsAdjacent(const CGeometry* grid, unsigned long a, unsigned long b) { + const auto& pts = grid->nodes->GetPoints(a); + return std::find(pts.begin(), pts.end(), b) != pts.end(); +} - for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { - /*--- Only seed lines from viscous (no-slip) wall markers. - * Non-wall boundaries (farfield, inlet, outlet, symmetry) must NOT seed - * lines because they would prematurely claim boundary-layer interior nodes. ---*/ - const auto bc = config->GetMarker_All_KindBC(iMarker); - if (bc != HEAT_FLUX && bc != ISOTHERMAL && bc != CHT_WALL_INTERFACE && bc != SMOLUCHOWSKI_MAXWELL) continue; +/*--- Is a footprint one connected patch? A set that falls into pieces is not the extrusion of + * anything, so an arriving footprint and a piece left by a split both have to pass this. ---*/ +bool IsConnectedLayer(const CGeometry* grid, const vector& layer) { + if (layer.size() < 2) return true; + vector seen(layer.size(), 0); + vector stk{0}; + seen[0] = 1; + size_t nSeen = 1; + while (!stk.empty()) { + const auto cur = stk.back(); + stk.pop_back(); + for (size_t k = 0; k < layer.size(); ++k) { + if (seen[k] || !IsAdjacent(grid, layer[cur], layer[k])) continue; + seen[k] = 1; + nSeen++; + stk.push_back(k); + } + } + return nSeen == layer.size(); +} + +/*--- Is the new layer topologically identical to the old? They are index-aligned, so phi maps + * old[k] to new[k], and the layer is valid when phi is an isomorphism of the induced subgraphs. ---*/ +bool LayerIsIsomorphic(const CGeometry* grid, const vector& oldL, const vector& newL) { + const auto n = oldL.size(); + if (newL.size() != n) return false; + + for (size_t k = 0; k < n; ++k) { + unsigned nOld = 0, nNew = 0; + for (size_t l = 0; l < n; ++l) { + nOld += IsAdjacent(grid, newL[k], oldL[l]); + nNew += IsAdjacent(grid, oldL[k], newL[l]); + } + /*--- Exactly one partner each way, and it has to be the one phi names. ---*/ + if ((nOld != 1) || (nNew != 1)) return false; + if (!IsAdjacent(grid, oldL[k], newL[k])) return false; + } + + for (size_t k = 0; k < n; ++k) + for (size_t l = k + 1; l < n; ++l) + if (IsAdjacent(grid, oldL[k], oldL[l]) != IsAdjacent(grid, newL[k], newL[l])) return false; + + return true; +} + +/*--- A boundary node ends a front only if the step runs INTO it, i.e. roughly along its normal. One + * running ALONG a boundary is a legitimate interior node of the stack. ---*/ +bool EntersBoundary(const CGeometry* grid, const CConfig* config, unsigned short nDim, unsigned long jPoint, + const su2double* stepDir, su2double cosBoundary) { + for (unsigned short iMarker = 0; iMarker < grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + su2double n[3] = {0.0}; /*--- nDim is at most 3. ---*/ + if (!VertexUnitNormal(grid, nDim, jPoint, iMarker, n)) continue; + if (fabs(GeometryToolbox::DotProduct(nDim, n, stepDir)) >= cosBoundary) return true; + } + return false; +} +/*--- Fine layers the next coarse CV of this front holds: two, so the stack coarsens by the same + * factor along the marching direction as the footprint does across it. One only if a second + * layer would exceed the agglomeration size limit. ---*/ +unsigned long BlockFor(short int maxAgglomSize, const vector& layer) { + return (layer.size() * 2 > static_cast(maxAgglomSize)) ? 1 : 2; +} + +/*--- Rank-independent name for a set of nodes: the smallest global index in it, +1 so that 0 is + * free to mean "nothing". A footprint is claimed by one front and no other, so it is unique. ---*/ +unsigned long TagOfSet(const CGeometry* grid, const vector& set) { + unsigned long t = std::numeric_limits::max(); + for (auto p : set) t = std::min(t, grid->nodes->GetGlobalIndex(p)); + return t + 1; +} + +} // namespace + +CMultiGridGeometry::CFrontSeeds CMultiGridGeometry::SeedFrontNodes(const CGeometry* fine_grid, const CConfig* config, + const CNodeStiffness& stiff) const { + /*--- Fraction of a marker's nodes that must sit in a layer before the whole marker may seed. ---*/ + constexpr passivedouble QUALIFIED_FRACTION = 0.5; + constexpr passivedouble ANGLE_THRESHOLD_DEG = 30.0; + /*--- Smallest local cell aspect ratio for which a node still counts as part of a stretched layer. + * Not exposed as an option: the two populations it separates are decades apart - a boundary + * carrying a layer normal to itself measures in the hundreds or more, one in isotropic mesh + * measures about one - so every value between them selects the same set of markers. ---*/ + constexpr passivedouble MIN_AR = 2.0; + const su2double cos_threshold = cos(ANGLE_THRESHOLD_DEG * PI_NUMBER / 180.0); + + const auto nMarkerFine = fine_grid->GetnMarker(); + constexpr auto NO_POINT = std::numeric_limits::max(); + + CFrontSeeds seeds; + vector taken(fine_grid->GetnPoint(), 0); + + auto isWall = [&](unsigned short bc) { + return (bc == HEAT_FLUX) || (bc == ISOTHERMAL) || (bc == CHT_WALL_INTERFACE) || (bc == SMOLUCHOWSKI_MAXWELL); + }; + + /*--- True if the mesh at iPoint is stretched along the boundary normal, i.e. this boundary has a + * layer growing off it the way a viscous wall does. On the side planes of a bump the mesh is just + * as stretched, but in a direction running ALONG the plane, and those nodes belong to the wall. ---*/ + auto hasLayerNormalTo = [&](unsigned long iPoint, const su2double* unitNormal) { + const auto jStiffest = stiff.jStiffest[iPoint]; + if (jStiffest == NO_POINT) return false; + if (stiff.AspectRatio(iPoint) < MIN_AR) return false; + + su2double vec[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(jStiffest), fine_grid->nodes->GetCoord(iPoint), vec); + const su2double len = GeometryToolbox::Norm(nDim, vec); + if (len <= 0.0) return false; + for (unsigned short d = 0; d < nDim; ++d) vec[d] /= len; + return fabs(GeometryToolbox::DotProduct(nDim, vec, unitNormal)) >= cos_threshold; + }; + + auto seedMarker = [&](unsigned short iMarker, bool requireLayer) { for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); + if (!fine_grid->nodes->GetDomain(iPoint)) continue; + if (fine_grid->nodes->GetAgglomerate(iPoint)) continue; + if (taken[iPoint]) continue; /*--- A node on two markers must seed only one front. ---*/ - /*--- Get vertex normal to seed the line direction ---*/ - const long ChildVertex = fine_grid->nodes->GetVertex(iPoint, iMarker); - if (ChildVertex == -1) continue; su2double Normal[MAXNDIM] = {0.0}; - fine_grid->vertex[iMarker][ChildVertex]->GetNormal(Normal); - - /*--- Normalize the direction ---*/ - su2double prev_dir[MAXNDIM] = {0.0}; - su2double norm_prev = 0.0; - for (unsigned short d = 0; d < nDim; ++d) { - prev_dir[d] = Normal[d]; - norm_prev += Normal[d] * Normal[d]; - } - if (norm_prev <= 0.0) continue; - norm_prev = sqrt(norm_prev); - for (unsigned short d = 0; d < nDim; ++d) prev_dir[d] /= norm_prev; - - /*--- Build the implicit line by following the best-aligned interior neighbor ---*/ - vector L; - L.push_back(iPoint); - auto current = iPoint; - - while (L.size() < MAX_LINE_LENGTH) { - su2double best_dot = -2.0; - unsigned long best_neighbor = ULONG_MAX; - - for (auto jPoint : fine_grid->nodes->GetPoints(current)) { - if (jPoint == current) continue; - if (!fine_grid->nodes->GetDomain(jPoint)) continue; - if (fine_grid->nodes->GetBoundary(jPoint)) continue; - if (fine_grid->nodes->GetAgglomerate(jPoint)) continue; - - /*--- Compute normalized direction to candidate ---*/ - su2double vec[MAXNDIM] = {0.0}; - GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(jPoint), fine_grid->nodes->GetCoord(current), vec); - const su2double len = GeometryToolbox::Norm(nDim, vec); - if (len <= 0.0) continue; - for (unsigned short d = 0; d < nDim; ++d) vec[d] /= len; - - /*--- Alignment with previous direction ---*/ - const su2double dot = GeometryToolbox::DotProduct(nDim, vec, prev_dir); - if (dot > best_dot) { - best_dot = dot; - best_neighbor = jPoint; - } - } + if (!VertexUnitNormal(fine_grid, nDim, iPoint, iMarker, Normal)) continue; + if (requireLayer && !hasLayerNormalTo(iPoint, Normal)) continue; + + std::array n0{}; + for (unsigned short d = 0; d < nDim; ++d) n0[d] = Normal[d]; + seeds.node.push_back(iPoint); + seeds.normal.push_back(n0); + taken[iPoint] = 1; + } + }; + + /*--- Viscous walls always carry a stretched layer, so they seed unconditionally. Running them first + * also settles nodes where a wall meets another boundary: the wall claims them. ---*/ + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) + if (isWall(config->GetMarker_All_KindBC(iMarker))) seedMarker(iMarker, false); + + /*--- Non-wall boundaries that still carry a layer normal to themselves, such as a symmetry plane in + * the same surface as a wall. The verdict is per marker, not per node: seeding isolated qualifying + * nodes scatters one-node patches that do not coarsen tangentially at all. ---*/ + /*--- Counted per configuration-file marker, not per local marker: ranks agree on neither the + * number nor the order of local markers, because each partition appends its own SEND_RECEIVE + * markers, so the same index means a different boundary elsewhere. ---*/ + const auto nMarkerCfg = config->GetnMarker_CfgFile(); + vector nValid(nMarkerCfg, 0), nQualified(nMarkerCfg, 0); + + auto canSeed = [&](unsigned short iMarker) { + const auto bc = config->GetMarker_All_KindBC(iMarker); + /*--- Periodic boundaries are left out: the two halves are the same physical location under a + * transform and have their own matching, which a front running into one would disturb. ---*/ + return (bc != SEND_RECEIVE) && (bc != PERIODIC_BOUNDARY) && !isWall(bc); + }; + + vector cfgOfMarker(nMarkerFine, 0); + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) + if (canSeed(iMarker)) + cfgOfMarker[iMarker] = config->GetMarker_CfgFile_TagBound(config->GetMarker_All_TagBound(iMarker)); + + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) { + if (!canSeed(iMarker)) continue; + for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { + const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); + if (!fine_grid->nodes->GetDomain(iPoint)) continue; + su2double Normal[MAXNDIM] = {0.0}; + if (!VertexUnitNormal(fine_grid, nDim, iPoint, iMarker, Normal)) continue; + nValid[cfgOfMarker[iMarker]]++; + if (hasLayerNormalTo(iPoint, Normal)) nQualified[cfgOfMarker[iMarker]]++; + } + } - if (best_neighbor == ULONG_MAX || best_dot < cos_threshold) break; + /*--- A marker is generally split over several ranks, so the verdict must be taken on all of it. + * Every rank reaches these collectives, including one that owns no boundary at all. ---*/ + if (nMarkerCfg > 0) { + vector tmp(nMarkerCfg); + SU2_MPI::Allreduce(nValid.data(), tmp.data(), nMarkerCfg, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + nValid.swap(tmp); + SU2_MPI::Allreduce(nQualified.data(), tmp.data(), nMarkerCfg, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + nQualified.swap(tmp); + } + + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) { + if (!canSeed(iMarker)) continue; + const auto iCfg = cfgOfMarker[iMarker]; + if (nValid[iCfg] == 0) continue; + if (su2double(nQualified[iCfg]) < QUALIFIED_FRACTION * su2double(nValid[iCfg])) continue; + seedMarker(iMarker, true); + } - L.push_back(best_neighbor); + return seeds; +} - /*--- Update direction for next step ---*/ - GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(best_neighbor), fine_grid->nodes->GetCoord(current), - prev_dir); - const su2double len = GeometryToolbox::Norm(nDim, prev_dir); - if (len <= 0.0) break; - for (unsigned short d = 0; d < nDim; ++d) prev_dir[d] /= len; +vector> CMultiGridGeometry::BuildFrontPatches(const CFrontSeeds& seeds, + const CGeometry* fine_grid, const CConfig* config, + const vector& mixedBC) const { + /*--- Repeated pairwise matching, one round per doubling, partitions the seeds into compact patches: + * a boundary edge in 2D, a boundary quadrilateral in 3D. Matching keeps membership mutually + * exclusive, which choosing neighbours per seed would not. It needs only point connectivity, so it + * works on every multigrid level. This patch alone fixes the footprint of the stack above it. ---*/ + const auto nSeeds = seeds.node.size(); + const unsigned long max_group = (nDim == 2) ? 2 : 4; + + /*--- Physical markers each seed lies on, ascending. Seeds may only be matched when these agree, so + * a patch never straddles a change of boundary condition. ---*/ + const auto nMarkerFine = fine_grid->GetnMarker(); + vector> sig(nSeeds); + for (unsigned long si = 0; si < nSeeds; ++si) + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) + if ((config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) && + (fine_grid->nodes->GetVertex(seeds.node[si], iMarker) != -1)) + sig[si].push_back(iMarker); + + /*--- Seed-to-seed adjacency, inherited from the boundary nodes' mesh connectivity. ---*/ + vector seedOfNode(fine_grid->GetnPoint(), -1); + for (unsigned long si = 0; si < nSeeds; ++si) seedOfNode[seeds.node[si]] = static_cast(si); + + vector> adj(nSeeds); + for (unsigned long si = 0; si < nSeeds; ++si) + for (auto jPoint : fine_grid->nodes->GetPoints(seeds.node[si])) { + const auto sj = seedOfNode[jPoint]; + if ((sj >= 0) && (static_cast(sj) != si)) adj[si].push_back(static_cast(sj)); + } - current = best_neighbor; - } + /*--- Global point index of each seed, used everywhere below as the deterministic sort key. Local + * indices depend on the partitioning, so ordering by them would make the coarse grid depend on + * the rank count. ---*/ + vector sgkey(nSeeds); + for (unsigned long si = 0; si < nSeeds; ++si) sgkey[si] = fine_grid->nodes->GetGlobalIndex(seeds.node[si]); - /*--- Accept only lines with at least 2 interior nodes (length >= 3 including wall) ---*/ - if (L.size() >= 3) { - lines.push_back(std::move(L)); + vector> groups; + vector groupOf(nSeeds); + groups.reserve(nSeeds); + + /*--- Every seed starts as its own group; the rounds below merge them. ---*/ + for (unsigned long si = 0; si < nSeeds; ++si) { + groupOf[si] = si; + groups.push_back({si}); + } + + const unsigned nRounds = (max_group <= 1) ? 0 : ((max_group <= 2) ? 1 : 2); + + /*--- One admissible merge of two groups, weighted by how many seed-to-seed adjacencies they share. + * A group lying ALONGSIDE this one touches it along its whole length and shares two, whereas one + * continuing in the same direction touches at an end and shares one. So weight 2 is the square + * and weight 1 is the strip, and the strip extrudes into a coarse CV elongated in one + * boundary-tangential direction. ---*/ + struct CMerge { + unsigned long g, h; /*!< \brief The two groups, g < h. */ + unsigned long weight; /*!< \brief Shared adjacencies: 2 makes a square, 1 makes a strip. */ + unsigned long keyG, keyH; /*!< \brief Their global-index keys, the deterministic tie-break. */ + }; + + vector merges; + vector touched, nShared(nSeeds, 0); + vector gkey; + vector consumed; + + for (unsigned round = 0; round < nRounds; ++round) { + const auto nGroups = groups.size(); + + /*--- Sort key of each group: the smallest global point index it holds. ---*/ + gkey.assign(nGroups, std::numeric_limits::max()); + for (unsigned long g = 0; g < nGroups; ++g) + for (auto si : groups[g]) gkey[g] = std::min(gkey[g], sgkey[si]); + + /*--- Every merge this round could make. Counted once per unordered pair: adjacency is symmetric, + * so the count from g's side equals the count from h's, and taking only h > g avoids both. ---*/ + merges.clear(); + for (unsigned long g = 0; g < nGroups; ++g) { + /*--- A node where two different boundary conditions meet stays a patch of its own, so the front + * rising from it is one node wide. Both sides of a merge are tested: skipping only the mixed + * group would still let an ordinary group reach out and take it. ---*/ + if (mixedBC[seeds.node[groups[g].front()]]) continue; + touched.clear(); + for (auto si : groups[g]) + for (auto sj : adj[si]) { + const auto h = groupOf[sj]; + if (h <= g) continue; + if (mixedBC[seeds.node[groups[h].front()]]) continue; + if (groups[g].size() + groups[h].size() > max_group) continue; + if (sig[groups[h].front()] != sig[groups[g].front()]) continue; + if (nShared[h]++ == 0) touched.push_back(h); + } + for (auto h : touched) { + merges.push_back({g, h, nShared[h], gkey[g], gkey[h]}); + nShared[h] = 0; } } + + /*--- Best merges first over all groups at once. Sweeping groups in index order lets one with no + * square partner take a weight-1 merge and consume a group a later one needed, and the failures + * cascade. Ordering globally makes every square before the first strip is considered. ---*/ + std::sort(merges.begin(), merges.end(), [](const CMerge& a, const CMerge& b) { + if (a.weight != b.weight) return a.weight > b.weight; + if (a.keyG != b.keyG) return a.keyG < b.keyG; + return a.keyH < b.keyH; + }); + + consumed.assign(nGroups, 0); + vector> merged; + merged.reserve(nGroups); + + for (const auto& m : merges) { + if (consumed[m.g] || consumed[m.h]) continue; + consumed[m.g] = consumed[m.h] = 1; + auto group = groups[m.g]; + group.insert(group.end(), groups[m.h].begin(), groups[m.h].end()); + merged.push_back(std::move(group)); + } + /*--- Whatever found no partner passes through unchanged. ---*/ + for (unsigned long g = 0; g < nGroups; ++g) + if (!consumed[g]) merged.push_back(std::move(groups[g])); + + groups = std::move(merged); + for (unsigned long g = 0; g < groups.size(); ++g) + for (auto si : groups[g]) groupOf[si] = g; } - if (lines.empty()) return; + return groups; +} - if (rank == MASTER_NODE) { - cout << "Implicit line agglomeration: detected " << lines.size() << " lines." << endl; +string CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const CGeometry* fine_grid, + const CConfig* config, unsigned short iMesh) { + /*--- Paving by advancing fronts. Each boundary patch rises into the domain keeping its footprint, so + * a coarse CV never spans two patches. A front stops at a boundary or when the next layer is not + * isomorphic to the current one. Each coarse CV is the footprint two layers deep. ---*/ + const auto starting_Index_CoarseCV = Index_CoarseCV; + const auto nPointFine = fine_grid->GetnPoint(); + const auto nMarkerFine = fine_grid->GetnMarker(); + constexpr auto NO_POINT = std::numeric_limits::max(); + const short int maxAgglomSize = (nDim == 2) ? 4 : 8; + + /*--- How nearly parallel a step must be to a boundary's normal to count as running INTO that + * boundary rather than along it. This is not a limit on where a front may go - it is only how + * "the front has reached a boundary" is recognised. ---*/ + constexpr passivedouble BOUNDARY_ALIGN_DEG = 30.0; + const su2double cos_boundary = cos(BOUNDARY_ALIGN_DEG * PI_NUMBER / 180.0); + /*--- Weight of the new step direction when the front's marching direction is updated. The direction + * only ever RANKS candidates, it never rejects one, so this is a preference and not a limit. ---*/ + constexpr passivedouble DIR_BLEND = 0.5; + + const auto stiff = ComputeNodeStiffness(fine_grid); + + /*--- PHASE 1. SeedFrontNodes must be reached by every rank, including one that owns no boundary: + * it takes a collective to agree on which markers carry a layer. ---*/ + const auto seeds = SeedFrontNodes(fine_grid, config, stiff); + const auto mixedBC = FindMixedBoundaryNodes(fine_grid, config); + const auto patches = BuildFrontPatches(seeds, fine_grid, config, mixedBC); + + /*--- Nodes on a boundary carrying a boundary condition; a front must not grow into one or the stack + * would straddle two boundaries. CPoint's Boundary flag is also set by SEND_RECEIVE, so it would + * stop every front one layer short of a partition. ---*/ + vector onPhysicalBoundary(nPointFine, 0); + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) + onPhysicalBoundary[fine_grid->vertex[iMarker][iVertex]->GetNode()] = 1; } - /*--- Advancing-front greedy pairing with cross-line merging. - * For each pair stage k, process interior positions (1+2k, 1+2k+1). - * When two lines share the same wall-node parent CV, merge their pairs - * into a single 4-child coarse CV. Otherwise create 2-child coarse CVs. ---*/ - vector reserved(nPointFine, 0); - unsigned pair_idx = 0; - - while (true) { - bool any_work = false; - - /*--- Build map: wall parent CV -> list of line indices ---*/ - unordered_map> parent_to_lines; - parent_to_lines.reserve(lines.size()); - for (unsigned long li = 0; li < lines.size(); ++li) { - const auto& L = lines[li]; - if (L.empty()) continue; - const auto idx2 = 1 + 2 * pair_idx + 1; - if (L.size() <= idx2) continue; // no pair at this stage - const auto pW = fine_grid->nodes->GetParent_CV(L[0]); - parent_to_lines[pW].push_back(li); + /*================================================================================================== + * PHASE 2 - advance every front, one layer per round. + *================================================================================================*/ + + /*--- One front node's chosen successor, before contention over it is resolved. ---*/ + struct CStep { + unsigned long node; /*!< Candidate successor. */ + unsigned long from; /*!< The front node that proposed it. */ + unsigned long key; /*!< Global index of "from", the final deterministic tie-break. */ + su2double score; /*!< Alignment of the step with the front's marching direction. */ + su2double dist; /*!< Length of the step. */ + su2double dir[MAXNDIM]; /*!< Unit step direction, reused to update the front's direction. */ + }; + + /*--- One advancing front. Fronts are not a fixed set: one handed over from a neighbouring rank is + * appended while the rounds are running, so the loops below are bounded by fronts.size(). ---*/ + struct CFront { + vector nodes; /*!< \brief Current footprint. */ + vector pending; /*!< \brief Nodes buffered for the coarse CV being built. */ + vector handTo; /*!< \brief Halo nodes the stack should continue onto. */ + vector prop; /*!< \brief This round's proposed successors. */ + std::array dir{}; /*!< \brief Marching direction. */ + unsigned long tag = 0; /*!< \brief Rank-independent name, see TagOfSet. */ + unsigned long handTag = 0; /*!< \brief Name the handed-over piece travels under. After a + * split this is NOT tag: the pieces are separate stacks. */ + unsigned long depth = 0; /*!< \brief Layers laid. */ + unsigned long nBlock = 0; /*!< \brief Fine layers the next coarse CV holds. */ + unsigned long pendingLayers = 0; /*!< \brief Layers currently buffered. */ + char alive = 1; + char failed = 0; + char keepLocal = 0; /*!< \brief Handed over only part of its footprint, so it marches on here. */ + }; + vector fronts; + + auto addFront = [&](const vector& layer, const std::array& dir, + unsigned long frontTag, unsigned long block) { + CFront F; + F.nodes = F.pending = layer; + F.dir = dir; + F.tag = frontTag; + F.nBlock = block; + F.pendingLayers = 1; + fronts.push_back(std::move(F)); + return fronts.size() - 1; + }; + + vector claimed(nPointFine, 0); + + /*--- Best free neighbour of n to step onto, ranked by alignment with dir. Local and halo candidates + * are ranked separately: a halo node cannot be claimed here, but it is where the stack would go + * next, so it is what gets handed over. ---*/ + struct CCandidate { + unsigned long node = std::numeric_limits::max(); + unsigned long halo = std::numeric_limits::max(); + su2double dot = -2.0, len = 0.0, dir[MAXNDIM] = {0.0}; + }; + + auto bestSuccessor = [&](unsigned long n, const su2double* marchDir) { + CCandidate c; + su2double haloDot = -2.0; + + for (auto jPoint : fine_grid->nodes->GetPoints(n)) { + su2double vec[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(jPoint), fine_grid->nodes->GetCoord(n), vec); + const su2double len = GeometryToolbox::Norm(nDim, vec); + if (len <= 0.0) continue; + for (unsigned short d = 0; d < nDim; ++d) vec[d] /= len; + + /*--- The direction only RANKS candidates; a front is never stopped for turning. ---*/ + const su2double dot = GeometryToolbox::DotProduct(nDim, vec, marchDir); + const bool admissible = + !(onPhysicalBoundary[jPoint] && EntersBoundary(fine_grid, config, nDim, jPoint, vec, cos_boundary)) && + GeometricalCheck(jPoint, fine_grid, config); + + /*--- Halo parents are assigned by the owning rank through the MPI relay, so claiming one here + * would fight that assignment. Admissibility is still checked, as the owner would. ---*/ + if (!fine_grid->nodes->GetDomain(jPoint)) { + if (dot > haloDot && admissible) { + haloDot = dot; + c.halo = jPoint; + } + continue; + } + if (fine_grid->nodes->GetAgglomerate(jPoint) || claimed[jPoint] || !admissible) continue; + + if (dot > c.dot) { + c.dot = dot; + c.node = jPoint; + c.len = len; + for (unsigned short d = 0; d < nDim; ++d) c.dir[d] = vec[d]; + } + } + return c; + }; + + /*--- Bid table: only an index per mesh point, the bids themselves in a compact vector. A whole + * CStep per point would be hundreds of megabytes that is empty almost everywhere. ---*/ + constexpr unsigned NOBID = std::numeric_limits::max(); + vector bidIdx(nPointFine, NOBID); + vector bids; + vector bidOwner; + + /*--- Scratch for the layer under construction, hoisted so a front does not allocate per layer. ---*/ + vector newLayer; + + /*--- Summed over all ranks for the one-line report at the end. ---*/ + enum { P_STACKS, P_LAYERS, P_COVERED, P_COUNT }; + unsigned long ct[P_COUNT] = {0}; + + /*--- One footprint node arriving from a neighbouring rank, to be regrouped by tag. ---*/ + struct CInherited { + unsigned long tag; /*!< \brief The front it belongs to, the same name on both ranks. */ + unsigned long node; /*!< \brief Local index of the node, owned by this rank. */ + su2double dir[MAXNDIM]; /*!< \brief The marching direction the stack arrives with. */ + }; + vector inherited; + + /*--- Where each halo node sits in this rank's SEND_RECEIVE receive lists, so a handover can be + * packed against the right vertex without searching. A halo node appears in exactly one. ---*/ + vector haloMarker(nPointFine, -1); + vector haloVertex(nPointFine, 0); + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + if (!((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0))) + continue; + const auto MarkerR = iMarker + 1; + for (auto iVertex = 0ul; iVertex < fine_grid->nVertex[MarkerR]; iVertex++) { + const auto p = fine_grid->vertex[MarkerR][iVertex]->GetNode(); + haloMarker[p] = static_cast(MarkerR); + haloVertex[p] = iVertex; } + } - vector line_processed(lines.size(), 0); + auto markFail = [&](unsigned long f) { fronts[f].failed = 1; }; + + /*--- How many fine layers the next coarse CV of this front holds: two, so the stack coarsens by the + * same factor along the marching direction as the footprint does across it. One only if a second + * layer would exceed the agglomeration size limit. ---*/ + /*--- Turn everything buffered for this front into one coarse control volume. ---*/ + auto emit = [&](unsigned long f) { + if (fronts[f].pending.empty()) return; + for (unsigned long c = 0; c < fronts[f].pending.size(); ++c) { + const auto p = fronts[f].pending[c]; + fine_grid->nodes->SetParent_CV(p, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, c, p); + if (fine_grid->nodes->GetAgglomerate_Indirect(p)) nodes->SetAgglomerate_Indirect(Index_CoarseCV, true); + } + nodes->SetnChildren_CV(Index_CoarseCV, static_cast(fronts[f].pending.size())); + Index_CoarseCV++; + ct[P_COVERED] += fronts[f].pending.size(); + + fronts[f].pending.clear(); + fronts[f].pendingLayers = 0; + fronts[f].nBlock = BlockFor(maxAgglomSize, fronts[f].nodes); + }; + + /*--- The boundary layer of every front. Claiming it before ordinary boundary agglomeration runs is + * what keeps the stack aligned: every layer above has exactly the same footprint. ---*/ + for (const auto& patch : patches) { + /*--- A one-node patch may seed a front and marches as a stack one node wide. Leaving it to + * ordinary agglomeration is worse: its whole column then goes unpaved, and such seeds sit where + * markers meet, cutting full-height stripes through the paved region. ---*/ + bool valid = !patch.empty(); + for (auto si : patch) { + const auto p = seeds.node[si]; + if (!GeometricalCheck(p, fine_grid, config) || fine_grid->nodes->GetAgglomerate(p)) valid = false; + } + if (!valid) continue; + + vector layer0; + std::array n0{}; + unsigned long frontTag = std::numeric_limits::max(); + for (auto si : patch) { + layer0.push_back(seeds.node[si]); + frontTag = std::min(frontTag, fine_grid->nodes->GetGlobalIndex(seeds.node[si])); + for (unsigned short d = 0; d < nDim; ++d) n0[d] += seeds.normal[si][d]; + } + const su2double nrm = GeometryToolbox::Norm(nDim, n0.data()); + /*--- A patch whose members' normals cancel has no direction to march in. ---*/ + if (nrm <= 0.0) continue; + for (unsigned short d = 0; d < nDim; ++d) n0[d] /= nrm; + + /*--- The boundary layer is a coarse CV on its own: a boundary node is never merged with an + * interior one, so only the first advance is a single layer. This also isolates a junction + * node without a rule of its own, since such a node is a patch of one. ---*/ + const auto f = addFront(layer0, n0, frontTag + 1, 1); + for (auto p : layer0) { + claimed[p] = 1; + } + ct[P_STACKS]++; + ct[P_LAYERS]++; + emit(f); + } - /*--- A) Cross-line merges: parents with multiple lines ---*/ - for (auto& [parent, line_ids] : parent_to_lines) { - if (line_ids.size() < 2) continue; + for (unsigned long layer = 1;; ++layer) { + /*--- Every rank runs the same number of rounds: each ends in a collective handover exchange, so + * one dropping out early would hang the others. ---*/ + int aliveLocal = 0; + for (unsigned long f = 0; f < fronts.size(); ++f) aliveLocal |= fronts[f].alive; + int aliveGlobal = 0; + SU2_MPI::Allreduce(&aliveLocal, &aliveGlobal, 1, MPI_INT, MPI_MAX, SU2_MPI::GetComm()); + if (aliveGlobal == 0) break; + + for (auto& F : fronts) F.failed = F.keepLocal = F.handTag = 0; + for (const auto& b : bids) bidIdx[b.node] = NOBID; + bids.clear(); + bidOwner.clear(); + + /*--- (a) Every alive front proposes a successor for each of its nodes. A front that cannot fill a + * whole layer proposes NOTHING: it is retiring this round anyway, and letting its partial bids + * stand would let a dying front displace a healthy one out of nodes it can still use. ---*/ + for (unsigned long f = 0; f < fronts.size(); ++f) { + if (!fronts[f].alive) continue; + fronts[f].prop.clear(); + fronts[f].handTo.clear(); + + for (auto n : fronts[f].nodes) { + const auto c = bestSuccessor(n, fronts[f].dir.data()); + + /*--- Nothing free here but the stack continues across the interface; the split test below + * decides whether the whole layer goes over. ---*/ + if ((c.node == NO_POINT) && (c.halo != NO_POINT)) { + fronts[f].handTo.push_back(c.halo); + continue; + } + /*--- No successor at all: a boundary, a partition, another front, or unusable mesh. ---*/ + if (c.node == NO_POINT) { + markFail(f); + break; + } - for (size_t k = 0; k + 1 < line_ids.size(); k += 2) { - const auto li1 = line_ids[k]; - const auto li2 = line_ids[k + 1]; - if (line_processed[li1] || line_processed[li2]) continue; + CStep s{c.node, n, fine_grid->nodes->GetGlobalIndex(n), c.dot, c.len, {}}; + for (unsigned short d = 0; d < nDim; ++d) s.dir[d] = c.dir[d]; + fronts[f].prop.push_back(s); + } - const auto& L1 = lines[li1]; - const auto& L2 = lines[li2]; - const auto idx1 = 1 + 2 * pair_idx; - const auto idx2 = idx1 + 1; - if (L1.size() <= idx2 || L2.size() <= idx2) continue; + /*--- An interface can cut a footprint. If all of it crosses, the stack is handed over intact + * and this front is finished. If only part crosses, the footprint is SPLIT: the piece whose + * successors are local marches on here, the rest is handed across, and both are renamed. ---*/ + if (fronts[f].failed) { + fronts[f].prop.clear(); + fronts[f].handTo.clear(); + } else if (!fronts[f].handTo.empty()) { + /*--- fronts[f].prop is built in the order of fronts[f].nodes, so this is the piece that stays, in the same + * order, and phi still runs index for index between it and the layer it proposes. ---*/ + vector narrow; + for (const auto& s : fronts[f].prop) narrow.push_back(s.from); + + /*--- A cut can leave the local piece in two disconnected halves - a square footprint cut + * diagonally does exactly that - and that is not a layer. Drop it and hand over the rest; + * the stack still survives on the far side instead of ending here. ---*/ + if (!narrow.empty() && !IsConnectedLayer(fine_grid, narrow)) { + narrow.clear(); + } - const auto a = L1[idx1], b = L1[idx2]; - const auto c = L2[idx1], d = L2[idx2]; + fronts[f].handTag = TagOfSet(fine_grid, fronts[f].handTo); + + if (narrow.empty()) { + fronts[f].prop.clear(); + } else { + /*--- Close the coarse CV that is open on the WIDE footprint before narrowing, so that no CV + * ever ends up holding two layers of different shape. ---*/ + fronts[f].nodes = narrow; + emit(f); + fronts[f].nBlock = BlockFor(maxAgglomSize, fronts[f].nodes); + fronts[f].tag = TagOfSet(fine_grid, fronts[f].nodes); + fronts[f].keepLocal = 1; + } + } + } - /*--- Skip if any node is already claimed ---*/ - if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b) || - fine_grid->nodes->GetAgglomerate(c) || fine_grid->nodes->GetAgglomerate(d)) + /*--- (b) Contention resolved from bids that were all collected before any was granted, so the + * outcome is a pure function of the proposals and does not depend on the order the fronts are + * visited in. That is what makes the coarse grid reproducible. ---*/ + auto better = [](const CStep& a, const CStep& b) { + if (a.score != b.score) return a.score > b.score; + if (a.dist != b.dist) return a.dist < b.dist; + return a.key < b.key; + }; + + for (unsigned long f = 0; f < fronts.size(); ++f) { + if (!fronts[f].alive || fronts[f].failed) continue; + for (const auto& s : fronts[f].prop) { + /*--- A front that has lost a bid is retiring and must not displace a healthy one with the + * rest of its layer. What it placed before losing stays, which is conservative. ---*/ + if (fronts[f].failed) break; + + if (bidIdx[s.node] == NOBID) { + bidIdx[s.node] = static_cast(bids.size()); + bids.push_back(s); + bidOwner.push_back(f); continue; - if (reserved[a] || reserved[b] || reserved[c] || reserved[d]) continue; + } - /*--- Geometrical quality check ---*/ - if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config) || - !GeometricalCheck(c, fine_grid, config) || !GeometricalCheck(d, fine_grid, config)) - continue; + const auto k = bidIdx[s.node]; + const auto g = bidOwner[k]; + /*--- Two nodes of the SAME front reaching for one successor is a pinch: the layer would come + * out narrower than the front, which all-or-nothing does not allow. ---*/ + if (better(s, bids[k])) { + markFail(g); + bids[k] = s; + bidOwner[k] = f; + } else { + markFail(f); + } - /*--- Guard against duplicate indices ---*/ - if (a == b || a == c || a == d || b == c || b == d || c == d) { - for (auto other_li : line_ids) line_processed[other_li] = 1; - continue; + /*--- A head-on meeting stops BOTH fronts, or the winner overshoots into the other's + * territory. A glancing contact is not a meeting and only costs the loser. ---*/ + if ((g != f) && (GeometryToolbox::DotProduct(nDim, fronts[f].dir.data(), fronts[g].dir.data()) < 0.0)) { + markFail(f); + markFail(g); } + } + } - /*--- Create 4-child coarse CV ---*/ - fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 0, a); - fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 1, b); - fine_grid->nodes->SetParent_CV(c, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 2, c); - fine_grid->nodes->SetParent_CV(d, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 3, d); - nodes->SetnChildren_CV(Index_CoarseCV, 4); - - reserved[a] = reserved[b] = reserved[c] = reserved[d] = 1; - MGQueue_InnerCV.RemoveCV(a); - MGQueue_InnerCV.RemoveCV(b); - MGQueue_InnerCV.RemoveCV(c); - MGQueue_InnerCV.RemoveCV(d); + /*--- (c) All-or-nothing acceptance: a front takes the whole layer or none of it and retires. A + * bid only becomes a claim here, so a retiring front never has to give anything back and the + * nodes it was reaching for stay available to ordinary agglomeration. ---*/ + for (unsigned long f = 0; f < fronts.size(); ++f) { + if (!fronts[f].alive) continue; + + newLayer.clear(); + if (!fronts[f].failed) { + /*--- Built in proposal order, so newLayer[k] is the successor of nodes[k] and the two carry phi. ---*/ + for (const auto& s : fronts[f].prop) { + const auto k = bidIdx[s.node]; + if ((k != NOBID) && (bidOwner[k] == f)) newLayer.push_back(s.node); + } + /*--- Every bid of a front that was not marked failed must have been granted. ---*/ + if ((newLayer.size() != fronts[f].nodes.size()) || !LayerIsIsomorphic(fine_grid, fronts[f].nodes, newLayer)) + markFail(f); + } - Index_CoarseCV++; - line_processed[li1] = line_processed[li2] = 1; - for (auto other_li : line_ids) - if (other_li != li1 && other_li != li2) line_processed[other_li] = 1; - any_work = true; + if (fronts[f].failed) { + /*--- Nothing to give back: a bid only becomes a claim on acceptance. ---*/ + fronts[f].alive = 0; + emit(f); + continue; + } + + /*--- Turn the marching direction towards the mean of the steps just taken. ---*/ + su2double mean[MAXNDIM] = {0.0}; + for (const auto& s : fronts[f].prop) + for (unsigned short d = 0; d < nDim; ++d) mean[d] += s.dir[d]; + const su2double meanNrm = GeometryToolbox::Norm(nDim, mean); + if (meanNrm > 0.0) { + su2double blended[MAXNDIM] = {0.0}; + for (unsigned short d = 0; d < nDim; ++d) + blended[d] = (1.0 - DIR_BLEND) * fronts[f].dir[d] + DIR_BLEND * mean[d] / meanNrm; + const su2double bNrm = GeometryToolbox::Norm(nDim, blended); + if (bNrm > 0.0) + for (unsigned short d = 0; d < nDim; ++d) fronts[f].dir[d] = blended[d] / bNrm; } + + for (auto p : newLayer) { + claimed[p] = 1; + } + fronts[f].nodes = std::move(newLayer); + fronts[f].depth++; + ct[P_LAYERS]++; + + fronts[f].pending.insert(fronts[f].pending.end(), fronts[f].nodes.begin(), fronts[f].nodes.end()); + fronts[f].pendingLayers++; + if (fronts[f].pendingLayers >= fronts[f].nBlock) emit(f); } - /*--- B) Single-line 2-child merges for remaining lines ---*/ - for (unsigned long li = 0; li < lines.size(); ++li) { - if (line_processed[li]) continue; - const auto& L = lines[li]; - const auto idx1 = 1 + 2 * pair_idx; - const auto idx2 = idx1 + 1; - if (L.size() <= idx2) continue; - - const auto a = L[idx1], b = L[idx2]; - if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b)) continue; - if (reserved[a] || reserved[b]) continue; - if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config)) continue; - - /*--- Create 2-child coarse CV ---*/ - fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 0, a); - fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 1, b); - nodes->SetnChildren_CV(Index_CoarseCV, 2); - - reserved[a] = reserved[b] = 1; - MGQueue_InnerCV.RemoveCV(a); - MGQueue_InnerCV.RemoveCV(b); + /*--- (d) Hand stacks across partition interfaces. A front that runs into the halo cannot go on + * here, so the footprint is sent to the owning rank, which picks the stack up and carries on. + * What crosses is the footprint, not a coarse CV, so both halves keep the same shape. It travels + * the reverse of the usual halo direction: packed against the RECEIVE marker and sent to the + * rank that marker receives from. ---*/ + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + if (!((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0))) + continue; - Index_CoarseCV++; - any_work = true; + const auto MarkerS = iMarker, MarkerR = iMarker + 1; + const auto send_to = config->GetMarker_All_SendRecv(MarkerS) - 1; + const auto receive_from = abs(config->GetMarker_All_SendRecv(MarkerR)) - 1; + const auto nVertexS = fine_grid->nVertex[MarkerS]; + const auto nVertexR = fine_grid->nVertex[MarkerR]; + + /*--- Packed against the halo vertices, i.e. what this rank wants the neighbour to continue. + * Tag and direction go separately: the AD MPI wrapper has no byte type to send a struct. ---*/ + vector tagOut(nVertexR, 0), tagIn(nVertexS, 0); + vector dirOut(nVertexR * nDim, 0.0), dirIn(nVertexS * nDim, 0.0); + + for (auto& F : fronts) + for (auto p : F.handTo) { + if (haloMarker[p] != static_cast(MarkerR)) continue; + const auto v = haloVertex[p]; + /*--- Two fronts reaching for one node: the lower tag takes it, the same way on both ranks. ---*/ + if ((tagOut[v] != 0) && (tagOut[v] <= F.handTag)) continue; + tagOut[v] = F.handTag; + for (unsigned short d = 0; d < nDim; ++d) dirOut[v * nDim + d] = F.dir[d]; + } + + SU2_MPI::Sendrecv(tagOut.data(), nVertexR, MPI_UNSIGNED_LONG, receive_from, 2, tagIn.data(), nVertexS, + MPI_UNSIGNED_LONG, send_to, 2, SU2_MPI::GetComm(), MPI_STATUS_IGNORE); + SU2_MPI::Sendrecv(dirOut.data(), nVertexR * nDim, MPI_DOUBLE, receive_from, 3, dirIn.data(), nVertexS * nDim, + MPI_DOUBLE, send_to, 3, SU2_MPI::GetComm(), MPI_STATUS_IGNORE); + + for (auto iVertex = 0ul; iVertex < nVertexS; iVertex++) { + if (tagIn[iVertex] == 0) continue; + inherited.push_back({tagIn[iVertex], fine_grid->vertex[MarkerS][iVertex]->GetNode(), {}}); + for (unsigned short d = 0; d < nDim; ++d) inherited.back().dir[d] = dirIn[iVertex * nDim + d]; + } } - pair_idx++; - if (!any_work) break; + /*--- A front that handed its WHOLE footprint over is finished here; the neighbour owns the rest + * of the stack. One that handed over only a piece keeps marching on what was left of it. ---*/ + for (unsigned long f = 0; f < fronts.size(); ++f) { + if (fronts[f].handTo.empty()) continue; + fronts[f].handTo.clear(); + if (fronts[f].keepLocal) continue; + fronts[f].alive = 0; + emit(f); + } - /*--- Check if any line still has pairs at the next stage ---*/ - bool any_more = false; - for (const auto& L : lines) { - if (L.size() > 1 + 2 * pair_idx + 1) { - any_more = true; - break; + /*--- Adopt what the neighbours sent, tags ascending so arrival order cannot change the outcome. + * A footprint whose nodes are not all free is dropped and the stack simply ends. ---*/ + std::sort(inherited.begin(), inherited.end(), + [](const CInherited& a, const CInherited& b) { return a.tag < b.tag; }); + + for (size_t i = 0; i < inherited.size();) { + size_t j = i; + while ((j < inherited.size()) && (inherited[j].tag == inherited[i].tag)) ++j; + + vector layer0; + bool ok = true; + for (size_t k = i; k < j; ++k) { + const auto p = inherited[k].node; + if (claimed[p] || fine_grid->nodes->GetAgglomerate(p) || !GeometricalCheck(p, fine_grid, config)) ok = false; + layer0.push_back(p); } + /*--- The footprint has to arrive whole and connected, the same test any other layer passes. ---*/ + if (ok && !IsConnectedLayer(fine_grid, layer0)) ok = false; + + if (ok) { + std::array d0{}; + for (unsigned short d = 0; d < nDim; ++d) d0[d] = inherited[i].dir[d]; + /*--- An inherited layer is an interior one, so it is NOT subject to the single-layer rule the + * boundary layer gets: it opens an ordinary two-deep coarse CV and waits for its partner. ---*/ + const auto nf = addFront(layer0, d0, inherited[i].tag, BlockFor(maxAgglomSize, layer0)); + for (auto p : layer0) { + claimed[p] = 1; + } + ct[P_LAYERS]++; + if (fronts[nf].pendingLayers >= fronts[nf].nBlock) emit(nf); + } + i = j; } - if (!any_more) break; + inherited.clear(); + } + + /*--- Nothing should be left buffered, but a front retired outside the loop would strand its nodes + * with a parent index that was never assigned. ---*/ + for (unsigned long f = 0; f < fronts.size(); ++f) emit(f); + + /*--- A rank with no fronts leaves dmin at its sentinel so it stays out of the MPI_MIN. ---*/ + unsigned long dmin = std::numeric_limits::max(), dmax = 0; + for (unsigned long f = 0; f < fronts.size(); ++f) { + if (fronts[f].nodes.empty()) continue; + dmin = std::min(dmin, fronts[f].depth); + dmax = std::max(dmax, fronts[f].depth); } + + /*--- Every rank must reach these. ---*/ + unsigned long tot[P_COUNT] = {0}, pairTot[2] = {0}, depthMin = 0, depthMax = 0; + unsigned long pair[2] = {Index_CoarseCV - starting_Index_CoarseCV, seeds.node.size()}; + SU2_MPI::Allreduce(ct, tot, P_COUNT, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(pair, pairTot, 2, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&dmin, &depthMin, 1, MPI_UNSIGNED_LONG, MPI_MIN, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&dmax, &depthMax, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm()); + if (depthMin == std::numeric_limits::max()) depthMin = 0; + + if (rank != MASTER_NODE) return {}; + + stringstream out; + out << " MG level " << iMesh << " paving: " << tot[P_STACKS] << " fronts from " << pairTot[1] << " seeds, " + << pairTot[0] << " CVs covering " << tot[P_COVERED] << " nodes in " << tot[P_LAYERS] << " layers, depth " + << depthMin << " to " << depthMax << "\n"; + return out.str(); } diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index eee8c7235fd..b4d552d5dcc 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -321,12 +321,15 @@ class CMultiGridIntegration final : public CIntegration { static constexpr int MAX_MG_LEVELS = 10; + static constexpr unsigned short MAXNVAR = 25; + static constexpr unsigned short MAXNDIM = 3; + /*--- Early-exit smoothing state (shared across OMP threads via master write + barrier). ---*/ - bool mg_early_exit_flag = false; /*!< \brief Shared flag for early exit across OMP threads. */ + bool mg_early_exit_flag = false; /*!< \brief Shared flag for early exit across OMP threads. */ passivedouble mg_initial_smooth_rms = 0.0; /*!< \brief Initial RMS residual before current smoothing phase (FAS). */ passivedouble mg_prev_smooth_rms = 0.0; /*!< \brief RMS residual from previous smoothing step; used for stagnation detection. */ - passivedouble mg_fine_rms_ema = 0.0; /*!< \brief EMA of fine-grid pre-smooth RMS across cycles; cross-cycle trend signal. */ - passivedouble last_crossCycleRatio = 1.0; /*!< \brief crossCycleRatio from the most recent cycle; stored for display only. */ + passivedouble mg_fine_rms_ema = 0.0; /*!< \brief EMA of fine-grid pre-smooth RMS across cycles; cross-cycle trend signal. */ + passivedouble last_crossCycleRatio = 1.0; /*!< \brief crossCycleRatio from the most recent cycle; stored for display only. */ /*--- Actual iteration counts per MG level, filled each cycle for the compact output summary. ---*/ unsigned short lastPreSmoothIters[MAX_MG_LEVELS+1] = {}; diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 5904ed9febb..fa16c685437 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -809,11 +809,14 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { /*--- Loop over all the new grid ---*/ + string pavingReports; + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { /*--- Create main agglomeration structure ---*/ - geometry[iMGlevel] = new CMultiGridGeometry(geometry[iMGlevel-1], config, iMGlevel); + auto* coarse_grid = new CMultiGridGeometry(geometry[iMGlevel-1], config, iMGlevel); + geometry[iMGlevel] = coarse_grid; /*--- Protect against the situation that we were not able to complete the agglomeration for this level, i.e., there weren't enough points. @@ -824,6 +827,7 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { geometry[iMGlevel] = nullptr; break; } + pavingReports += coarse_grid->pavingReport; /*--- Compute points surrounding points. ---*/ @@ -850,6 +854,9 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { } + /*--- Held back so they do not interleave with the multigrid level table. ---*/ + if (rank == MASTER_NODE) cout << pavingReports; + if (config->GetWrt_MultiGrid()) geometry[MESH_0]->ColorMGLevels(config->GetnMGLevels(), geometry); /*--- For unsteady simulations, initialize the grid volumes diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index fecbf0491b7..5c85dabb4fb 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -63,6 +63,24 @@ inline passivedouble ComputeLinSysResRMS(const CSolver* solver) { return sqrt(result); } +/*!\cond PRIVATE + * Prolongate a coarse-grid field onto the fine grid via constant injection: every fine + * child gets its parent's value. \c getCoarse returns the coarse-grid block of a point + * and \c setFine writes it to a fine-grid point. + \endcond */ +template +void ProlongateField(CGeometry* geo_coarse, GetCoarse getCoarse, SetFine setFine) { + + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + setFine(Point_Fine, getCoarse(Point_Coarse)); + } + } + END_SU2_OMP_FOR +} + } // anonymous namespace void CMultiGridIntegration::adaptDampingFactors(CConfig* config, passivedouble crossCycleRatio) { @@ -489,8 +507,8 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, END_SU2_OMP_SAFE_GLOBAL_ACCESS } - /*--- Print compact smoothing summary when MG_SMOOTH_OUTPUT= YES. ---*/ - if (mgOptsZone.MG_Smooth_Output) { + /*--- Print compact smoothing summary when MG_SMOOTH_OUTPUT= YES and MGLEVEL > 0. ---*/ + if ((mgOptsZone.MG_Smooth_Output) && (nMGLevels > 0)) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS if (SU2_MPI::GetRank() == MASTER_NODE) { @@ -513,8 +531,13 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, return ss.str(); }; + const string eqName = (RunTime_EqSystem == RUNTIME_FLOW_SYS) ? "Flow" : + (RunTime_EqSystem == RUNTIME_TURB_SYS) ? "Turb" : + (RunTime_EqSystem == RUNTIME_SPECIES_SYS) ? "Species" : + (RunTime_EqSystem == RUNTIME_TRANS_SYS) ? "Trans" : "Other"; + PrintingToolbox::CTablePrinter table(&std::cout); - table.AddColumn("Smoother", 13); + table.AddColumn("Smoother [" + eqName + "]", 13 + 7); for (unsigned short i = 0; i <= nMGLevels; ++i) table.AddColumn("Level " + std::to_string(i), 38); table.PrintHeader(); @@ -877,14 +900,13 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS SU2_ZONE_SCOPED const unsigned short nVar = sol_coarse->GetnVar(); - su2activevector Solution(nVar); SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { - su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); + su2double Solution[MAXNVAR] = {0.0}; - Solution = su2double(0); + su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); /*--- Accumulate children contributions with stable ordering ---*/ /*--- Process all children in sequential order to ensure deterministic FP summation ---*/ @@ -903,8 +925,7 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] += Solution_Coarse[iVar]; - for (auto iVar = 0u; iVar < nVar; iVar++) - sol_coarse->GetNodes()->SetSolution_Old(Point_Coarse, Solution.data()); + sol_coarse->GetNodes()->SetSolution_Old(Point_Coarse, Solution); } END_SU2_OMP_FOR @@ -931,19 +952,19 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS } } - /*--- MPI the set solution old ---*/ + /*--- MPI the set solution old. ---*/ sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD); sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD); - SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { - for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - sol_fine->LinSysRes.SetBlock(Point_Fine, sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse)); - } - } - END_SU2_OMP_FOR + /*--- Interpolate the coarse-grid correction onto the fine + * grid and store in LinSysRes. ---*/ + + ProlongateField(geo_coarse, + [&](unsigned long iPoint) { return sol_coarse->GetNodes()->GetSolution_Old(iPoint); }, + [&](unsigned long Point_Fine, const su2double* value) { + sol_fine->LinSysRes.SetBlock(Point_Fine, value); + }); } @@ -969,10 +990,10 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ for (auto iSmooth = 0u; iSmooth < val_nSmooth; iSmooth++) { - /*--- Loop over all mesh points (sum the residuals of direct neighbors). ---*/ + /*--- Loop over the domain points, exclude halo points ---*/ - SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) - for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { + SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPointDomain(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geometry->GetnPointDomain(); ++iPoint) { solver->GetNodes()->SetResidualSumZero(iPoint); @@ -985,10 +1006,10 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } END_SU2_OMP_FOR - /*--- Loop over all mesh points (update residuals with the neighbor averages). ---*/ + /*--- Loop over the domain points (update residuals with the neighbor averages). ---*/ - SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) - for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { + SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPointDomain(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geometry->GetnPointDomain(); ++iPoint) { su2double factor = 1.0/(1.0+val_smooth_coeff*su2double(geometry->nodes->GetnPoint(iPoint))); @@ -1000,12 +1021,13 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } END_SU2_OMP_FOR - /*--- Restore original residuals (without average) at boundary points. ---*/ + /*--- Restore original residuals at physical boundary points. ---*/ for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && (config->GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { + (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) && + (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE)) { SU2_OMP_FOR_STAT(32) for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { @@ -1017,6 +1039,13 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } } + /*--- Refresh the halo entries of the correction with the values their owner ranks just + * computed. ---*/ + + SU2_OMP_BARRIER + CSysMatrixComms::Initiate(solver->LinSysRes, geometry, config); + CSysMatrixComms::Complete(solver->LinSysRes, geometry, config); + } /*--- Record final correction norm for debugging output. ---*/ @@ -1171,26 +1200,23 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar CGeometry *geo_coarse, CConfig *config, unsigned short iMesh) { SU2_ZONE_SCOPED - const su2double *Residual_Fine; - const unsigned short nVar = sol_coarse->GetnVar(); const su2double factor = config->GetDamp_Res_Restric(); - su2activevector Residual(nVar); - SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { sol_coarse->GetNodes()->SetRes_TruncErrorZero(Point_Coarse); - Residual = su2double(0); + su2double RestrictedDefect[MAXNVAR] = {0.0}; + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); + const su2double* Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); for (auto iVar = 0u; iVar < nVar; iVar++) - Residual[iVar] += factor * Residual_Fine[iVar]; + RestrictedDefect[iVar] += factor * Residual_Fine[iVar]; } - sol_coarse->GetNodes()->AddRes_TruncError(Point_Coarse, Residual.data()); + sol_coarse->GetNodes()->AddRes_TruncError(Point_Coarse, RestrictedDefect); } END_SU2_OMP_FOR @@ -1289,17 +1315,14 @@ void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSyst const unsigned short nDim = geo_coarse->GetnDim(); const unsigned short nVar = sol_coarse->GetnVar(); - auto **Gradient = new su2double* [nVar]; - for (auto iVar = 0u; iVar < nVar; iVar++) - Gradient[iVar] = new su2double [nDim]; - SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { - su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); - for (auto iVar = 0u; iVar < nVar; iVar++) - for (auto iDim = 0u; iDim < nDim; iDim++) - Gradient[iVar][iDim] = 0.0; + su2double GradientData[MAXNVAR][MAXNDIM] = {{0.0}}; + su2double* Gradient[MAXNVAR]; + for (auto iVar = 0u; iVar < nVar; iVar++) Gradient[iVar] = GradientData[iVar]; + + su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { unsigned long Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); @@ -1314,10 +1337,6 @@ void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSyst } END_SU2_OMP_FOR - for (auto iVar = 0u; iVar < nVar; iVar++) - delete [] Gradient[iVar]; - delete [] Gradient; - } void CMultiGridIntegration::NonDimensional_Parameters(CGeometry **geometry, CSolver ***solver_container, diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index fab97842048..29c71b2ed36 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -49,6 +49,20 @@ void CSingleGridIntegration::SingleGrid_Iteration(CGeometry ****geometry, CSolve CGeometry* geometry_fine = geometry[iZone][iInst][FinestMesh]; CSolver** solvers_fine = solver_container[iZone][iInst][FinestMesh]; + /*--- CFL scaling of turbulence during the warmpup phase if FMG. ---*/ + if ((RunTime_EqSystem == RUNTIME_TURB_SYS) && (FinestMesh != MESH_0)) { + const su2double turbReduction = SU2_TYPE::GetValue(config[iZone]->GetCFLRedCoeff_Turb()); + const su2double turbCFL = SU2_TYPE::GetValue(config[iZone]->GetCFL(FinestMesh)) * turbReduction; + auto* turbSolver = solvers_fine[Solver_Position]; + + SU2_OMP_SAFE_GLOBAL_ACCESS(turbSolver->SetCFL_Local_Stats(turbCFL);) + SU2_OMP_FOR_STAT(roundUpDiv(geometry_fine->GetnPoint(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geometry_fine->GetnPoint(); ++iPoint) { + turbSolver->GetNodes()->SetLocalCFL(iPoint, turbCFL); + } + END_SU2_OMP_FOR + } + /*--- Preprocessing ---*/ solvers_fine[Solver_Position]->Preprocessing(geometry_fine, solvers_fine, config[iZone], diff --git a/TestCases/euler/CRM/inv_CRM_JST.cfg b/TestCases/euler/CRM/inv_CRM_JST.cfg index 63c749021f2..60080a0254b 100644 --- a/TestCases/euler/CRM/inv_CRM_JST.cfg +++ b/TestCases/euler/CRM/inv_CRM_JST.cfg @@ -43,11 +43,11 @@ MARKER_MONITORING= ( fuselage , Wing , HTP ) % ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------% % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 5.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -EXT_ITER= 99999 +ITER= 1000 LINEAR_SOLVER= FGMRES LINEAR_SOLVER_ERROR= 1E-1 LINEAR_SOLVER_ITER= 3 @@ -55,12 +55,16 @@ LINEAR_SOLVER_ITER= 3 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % MGLEVEL= 3 +MG_MIN_MESHSIZE= 100 MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_SMOOTH_OUTPUT= YES +MG_SMOOTH_EARLY_EXIT= YES +MG_PRE_SMOOTH= ( 5, 5, 5, 5 ) +MG_POST_SMOOTH= ( 5, 5, 5, 5 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 +MG_CFL_SCALING= 0.5, 0.5, 0.5 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/euler/channel/inv_channel.cfg b/TestCases/euler/channel/inv_channel.cfg index 5161a04413d..8145521b147 100644 --- a/TestCases/euler/channel/inv_channel.cfg +++ b/TestCases/euler/channel/inv_channel.cfg @@ -47,11 +47,11 @@ MARKER_MONITORING= ( upper_wall, lower_wall ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= GREEN_GAUSS -CFL_NUMBER= 6.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -EXT_ITER= 999999 +ITER= 1000 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % @@ -63,12 +63,15 @@ LINEAR_SOLVER_ITER= 3 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % MGLEVEL= 3 -MGCYCLE= V_CYCLE +MGCYCLE= W_CYCLE +MG_SMOOTH_OUTPUT= NO +MG_SMOOTH_EARLY_EXIT= YES MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) MG_POST_SMOOTH= ( 4, 4, 4, 4 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) MG_DAMP_RESTRICTION= 0.5 MG_DAMP_PROLONGATION= 0.5 +MG_CFL_SCALING= 0.5, 0.5, 0.5 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/euler/channel/inv_channel_RK.cfg b/TestCases/euler/channel/inv_channel_RK.cfg index 6d44b6ad5e0..f0a487a55f5 100644 --- a/TestCases/euler/channel/inv_channel_RK.cfg +++ b/TestCases/euler/channel/inv_channel_RK.cfg @@ -44,11 +44,11 @@ MARKER_MONITORING= ( upper_wall, lower_wall ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 1.0 +CFL_NUMBER= 1 CFL_ADAPT= NO -CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +CFL_ADAPT_PARAM= ( 0.9, 1.1, 1.0, 100.0, 0.1 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -ITER= 110 +ITER= 1000 LINEAR_SOLVER= BCGSTAB LINEAR_SOLVER_ERROR= 1E-1 LINEAR_SOLVER_ITER= 10 @@ -57,11 +57,11 @@ LINEAR_SOLVER_ITER= 10 % MGLEVEL= 3 MGCYCLE= W_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) -MG_CORRECTION_SMOOTH= ( 4, 4, 4, 4 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2 ) +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % @@ -94,4 +94,4 @@ GRAD_OBJFUNC_FILENAME= of_grad SURFACE_FILENAME= surface_flow SURFACE_ADJ_FILENAME= surface_adjoint OUTPUT_WRT_FREQ= 250 -SCREEN_OUTPUT = (INNER_ITER, RMS_DENSITY, RMS_ENERGY, LIFT, DRAG) +SCREEN_OUTPUT = (INNER_ITER, RMS_DENSITY, RMS_ENERGY, LIFT, DRAG, AVG_CFL) diff --git a/TestCases/euler/oneram6/inv_ONERAM6.cfg b/TestCases/euler/oneram6/inv_ONERAM6.cfg index 6517b957871..59c19e14a53 100644 --- a/TestCases/euler/oneram6/inv_ONERAM6.cfg +++ b/TestCases/euler/oneram6/inv_ONERAM6.cfg @@ -49,15 +49,15 @@ MARKER_DESIGNING = ( WING ) % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES OBJECTIVE_FUNCTION= DRAG -CFL_NUMBER= 5.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) ITER= 99999 LINEAR_SOLVER= FGMRES LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-6 -LINEAR_SOLVER_ITER= 2 +LINEAR_SOLVER_ERROR= 1E-1 +LINEAR_SOLVER_ITER= 5 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------% % @@ -92,7 +92,7 @@ TIME_DISCRE_ADJFLOW= EULER_IMPLICIT % --------------------------- CONVERGENCE PARAMETERS --------------------------& % CONV_RESIDUAL_MINVAL= -12 -CONV_STARTITER= 25 +CONV_STARTITER= 10 CONV_CAUCHY_ELEMS= 100 CONV_CAUCHY_EPS= 1E-10 diff --git a/TestCases/euler/wedge/inv_wedge_HLLC.cfg b/TestCases/euler/wedge/inv_wedge_HLLC.cfg index 64fb1967366..4dce47518b8 100644 --- a/TestCases/euler/wedge/inv_wedge_HLLC.cfg +++ b/TestCases/euler/wedge/inv_wedge_HLLC.cfg @@ -55,11 +55,11 @@ LINEAR_SOLVER_ITER= 3 % MGLEVEL= 3 MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/fixed_cl/naca0012/inv_NACA0012.cfg b/TestCases/fixed_cl/naca0012/inv_NACA0012.cfg index b9c7cc06667..ade966ede21 100644 --- a/TestCases/fixed_cl/naca0012/inv_NACA0012.cfg +++ b/TestCases/fixed_cl/naca0012/inv_NACA0012.cfg @@ -69,11 +69,11 @@ LINEAR_SOLVER_ITER= 3 % MGLEVEL= 3 MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/navierstokes/cylinder/cylinder_lowmach.cfg b/TestCases/navierstokes/cylinder/cylinder_lowmach.cfg index 1c175581e3a..402b8c9158b 100644 --- a/TestCases/navierstokes/cylinder/cylinder_lowmach.cfg +++ b/TestCases/navierstokes/cylinder/cylinder_lowmach.cfg @@ -42,7 +42,7 @@ MARKER_MONITORING= ( cylinder ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 4.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) @@ -58,12 +58,12 @@ LINEAR_SOLVER_ITER= 3 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % MGLEVEL= 3 -MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MGCYCLE= W_CYCLE +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/navierstokes/cylinder/lam_cylinder.cfg b/TestCases/navierstokes/cylinder/lam_cylinder.cfg index 70f8692f64b..8d696643fbf 100644 --- a/TestCases/navierstokes/cylinder/lam_cylinder.cfg +++ b/TestCases/navierstokes/cylinder/lam_cylinder.cfg @@ -42,7 +42,7 @@ MARKER_MONITORING= ( cylinder ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 500.0 +CFL_NUMBER= 1000.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) @@ -59,11 +59,11 @@ LINEAR_SOLVER_ITER= 3 % MGLEVEL= 3 MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/navierstokes/flatplate/lam_flatplate.cfg b/TestCases/navierstokes/flatplate/lam_flatplate.cfg index f213b10a0a7..0be7e3573cf 100644 --- a/TestCases/navierstokes/flatplate/lam_flatplate.cfg +++ b/TestCases/navierstokes/flatplate/lam_flatplate.cfg @@ -60,11 +60,11 @@ LINEAR_SOLVER_ITER= 3 % MGLEVEL= 3 MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % @@ -92,3 +92,4 @@ RESTART_ADJ_FILENAME= restart_adj GRAD_OBJFUNC_FILENAME= of_grad OUTPUT_WRT_FREQ= 250 SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_ENERGY, DRAG, SURFACE_MASSFLOW, TOTAL_HEATFLUX, TOTAL_HEATFLUX_ON_SURFACE) +WRT_PERFORMANCE= YES diff --git a/config_template.cfg b/config_template.cfg index 535ae564510..e5fdd436aba 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1736,10 +1736,6 @@ MG_MIN_MESHSIZE= 500 % Enable agglomeration along implicit lines seeded from viscous walls (NO, YES) MG_IMPLICIT_LINES= NO % -% Maximum nodes on a wall-normal implicit agglomeration line, including the wall seed. -% Increase to extend the line deeper into the boundary layer (default 20). -MG_IMPLICIT_LINES_MAX_LENGTH= 20 -% % Number of iterations spent on each mesh during the Full Multigrid (FMG) startup phase. After % this many iterations the solution is prolongated to the next finer mesh. It is also the length % of the CFL ramp every level is brought up over, the finest grid included, which is what keeps