diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 3544f028267..934f88feb22 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -639,7 +639,8 @@ class CConfig { su2double Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ su2double Deform_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ su2double Linear_Solver_Smoother_Relaxation; /*!< \brief Relaxation factor for iterative linear smoothers. */ - unsigned long Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the implicit formulation. */ + unsigned long Linear_Solver_Iter; + unsigned long Linear_Solver_Prec_Freeze; /*!< \brief Reuse the finest-grid preconditioner for this many solves. */ /*!< \brief Max iterations of the linear solver for the implicit formulation. */ unsigned long Deform_Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the implicit formulation. */ unsigned long Linear_Solver_Restart_Frequency; /*!< \brief Restart frequency of the linear solver for the implicit formulation. */ unsigned long Linear_Solver_Restart_Deflation; /*!< \brief Number of vectors used for deflated restarts. */ @@ -1098,6 +1099,7 @@ class CConfig { su2double ParMETIS_tolerance; /*!< \brief Load balancing tolerance for ParMETIS. */ long ParMETIS_pointWgt; /*!< \brief Load balancing weight given to points. */ long ParMETIS_edgeWgt; /*!< \brief Load balancing weight given to edges. */ + su2double ParMETIS_anisoWgt; /*!< \brief Strength of the anisotropy-aware ParMETIS edge weights. 0 disables them. */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ bool DiscreteAdjoint, /*!< \brief AD-based discrete adjoint mode. */ DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ @@ -4380,6 +4382,12 @@ class CConfig { */ unsigned long GetLinear_Solver_Iter(void) const { return Linear_Solver_Iter; } + /*! + * \brief Number of consecutive linear solves that reuse one finest-grid preconditioner. + * \return Freeze period, 1 meaning rebuild on every solve. + */ + unsigned long GetLinear_Solver_Prec_Freeze(void) const { return Linear_Solver_Prec_Freeze; } + /*! * \brief Get max number of iterations of the linear solver for the implicit formulation. * \return Max number of iterations of the linear solver for the implicit formulation. @@ -10166,6 +10174,11 @@ class CConfig { */ long GetParMETIS_EdgeWeight() const { return ParMETIS_edgeWgt; } + /*! + * \brief Get the strength of the anisotropy-aware ParMETIS edge weights (0 disables them). + */ + passivedouble GetParMETIS_AnisoWeight() const { return SU2_TYPE::GetValue(ParMETIS_anisoWgt); } + /*! * \brief Find the marker index (if any) that is part of a given interface pair. * \param[in] iInterface - Number of the interface pair being tested, starting at 0. diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index 1f9bc851b92..4191b0cf11c 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/Common/include/linear_algebra/CSysSolve.hpp @@ -136,6 +136,13 @@ class CSysSolve { /*!< \brief Inner solver for nested preconditioning. */ std::unique_ptr> inner_solver; + /*--- Preconditioner freezing on coarse multigrid levels. The factorization lives in the + * CSysMatrix (not in the short-lived CPreconditioner object built in Solve), so simply + * skipping Build() reuses the previous one. This instance belongs to one solver on one + * grid level, so the counter is naturally per-level. See MG_COARSE_PREC_FREEZE. ---*/ + mutable unsigned long precSolveCount = 0; /*!< \brief Linear solves done by this instance. */ + mutable bool buildPrecThisSolve = true; /*!< \brief Decision for the current solve, shared by all threads. */ + /*! * \brief sign transfer function * \param[in] x - value having sign prescribed diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index d4d0e7f0fa0..bcc766f08da 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. */ @@ -1126,6 +1126,16 @@ struct CMGOptions { 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). */ + bool MG_Implicit_Lines_Isotropic{false}; /*!< \brief Use isotropic (vs anisotropic) agglomeration along implicit lines. */ + unsigned long MG_Implicit_Lines_Max_Group{0}; /*!< \brief Max number of parallel implicit lines merged into one coarse + CV tangential to the wall. 0 = dimension-appropriate default + (2 in 2D, 4 in 3D). See CMultiGridGeometry::AgglomerateImplicitLines. */ + su2double MG_Implicit_Lines_Min_AR{2.0}; /*!< \brief Smallest local cell aspect ratio for which a node still counts as + part of a stretched layer. Ends a line where the mesh stops being + stretched along it, and decides which boundaries carry a layer normal + to them. See CMultiGridGeometry::AgglomerateImplicitLines. */ + unsigned long MG_Coarse_Prec_Freeze{1}; /*!< \brief On MG levels > 0, reuse the linear-solver preconditioner for this many consecutive solves. 1 = rebuild every solve. */ + su2double MG_Correction_Limit{0.0}; /*!< \brief Max relative change of any solution component from one prolongated FAS correction. 0 = no limit. */ 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 c21116c6f3e..23273086cb4 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2065,12 +2065,36 @@ void CConfig::SetConfig_Options() { 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_COARSE_PREC_FREEZE\n DESCRIPTION: On multigrid levels above MESH_0, reuse the linear-solver preconditioner + * (e.g. the ILU factorization) for this many consecutive linear solves instead of rebuilding it every time. + * 1 reproduces the previous behaviour exactly. DEFAULT: 1 \ingroup Config*/ + addUnsignedLongOption("MG_COARSE_PREC_FREEZE", MGOptions.MG_Coarse_Prec_Freeze, 1); + /*!\brief LINEAR_SOLVER_PREC_FREEZE\n DESCRIPTION: Reuse the linear-solver preconditioner (e.g. the ILU factorization) + * for this many consecutive solves on the finest grid, instead of rebuilding it every time. Applies to MESH_0 and to + * single-grid runs; MG_COARSE_PREC_FREEZE covers the coarse levels. The fine-grid preconditioner drives the outer + * nonlinear convergence, so raise this one with more care. 1 rebuilds every solve. DEFAULT: 1 \ingroup Config*/ + addUnsignedLongOption("LINEAR_SOLVER_PREC_FREEZE", Linear_Solver_Prec_Freeze, 1); + /*!\brief MG_CORRECTION_LIMIT\n DESCRIPTION: Largest relative change any solution component may undergo from a single + * prolongated multigrid correction, e.g. 0.1 caps it at 10%. The whole correction vector at a point is scaled by one + * factor so its direction is preserved. 0 disables the limiter (previous behaviour). DEFAULT: 0 \ingroup Config*/ + addDoubleOption("MG_CORRECTION_LIMIT", MGOptions.MG_Correction_Limit, 0.0); + /*!\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: 50 \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_IMPLICIT_LINES_ISOTROPIC\n DESCRIPTION: Use isotropic agglomeration along implicit lines (4 cells per coarse CV) instead of anisotropic (2 cells per coarse CV). DEFAULT: NO \ingroup Config*/ + addBoolOption("MG_IMPLICIT_LINES_ISOTROPIC", MGOptions.MG_Implicit_Lines_Isotropic, false); + /*!\brief MG_IMPLICIT_LINES_MAX_GROUP\n DESCRIPTION: Maximum number of parallel implicit lines merged tangential to + * the wall into one coarse CV (2D: always 2; 3D: e.g. 4 for a wall quad/hex corner, 3 for a triangular prism apex). + * 0 uses the dimension-appropriate default (2 in 2D, 4 in 3D). DEFAULT: 0 \ingroup Config*/ + addUnsignedLongOption("MG_IMPLICIT_LINES_MAX_GROUP", MGOptions.MG_Implicit_Lines_Max_Group, 0); + /*!\brief MG_IMPLICIT_LINES_MIN_AR\n DESCRIPTION: Smallest local cell aspect ratio for which a node still counts as part + * of a stretched layer, measured from the ratio of dual-grid edge weights. Ends an implicit line where the mesh stops + * being stretched along it, instead of letting it run to the far field, and decides which boundaries carry a layer + * normal to them and may therefore seed lines. 1.0 disables both tests. DEFAULT: 2.0 \ingroup Config*/ + addDoubleOption("MG_IMPLICIT_LINES_MIN_AR", MGOptions.MG_Implicit_Lines_Min_AR, 2.0); /*!\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); @@ -3053,6 +3077,14 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: ParMETIS load balancing weight for edges (equiv. to neighbors) */ addLongOption("PARMETIS_EDGE_WEIGHT", ParMETIS_edgeWgt, 1); + /* DESCRIPTION: Strength of the anisotropy-aware ParMETIS edge weights. ParMETIS is otherwise given no edge weights at + * all, so every edge is equally cheap to cut and partition boundaries slice straight through the stretched cells of a + * boundary layer, splitting the wall-normal columns that implicit-line agglomeration and line-implicit smoothing rely + * on. Weighting an edge by the inverse of its length makes the short wall-normal edges expensive to cut and pushes the + * cuts into the tangential direction instead. On a mesh without stretching all edges are of similar length, the + * weights come out uniform, and the partitioning is the same as with no weights at all. 0 disables the weights. + * DEFAULT: 0 */ + addDoubleOption("PARMETIS_ANISO_WEIGHT", ParMETIS_anisoWgt, 0.0); /*--- options that are used in the Hybrid RANS/LES Simulations ---*/ /*!\par CONFIG_CATEGORY:Hybrid_RANSLES Options\ingroup Config*/ diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 360da1aaaa1..981ac966c66 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -25,6 +25,7 @@ * License along with SU2. If not, see . */ +#include #include #include "../../include/geometry/CGeometry.hpp" @@ -4356,6 +4357,19 @@ const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config) unsigned long maxNPoints = 0, sumNPoints = 0; + /*--- Why each line stopped growing. A line that ends because the mesh has become isotropic is a + * line that has done its job; one that ends because no neighbour was well enough aligned means + * the walk lost the wall-normal direction and the line is short despite the mesh still being + * stretched. The two call for opposite responses on coarse grids, so they are counted apart. ---*/ + unsigned long nStopIsotropic = 0, nStopNoNeighbour = 0, nStopCap = 0; + /*--- "No neighbour" has two very different causes: every candidate was already taken by another + * line (a competition/ordering problem), or candidates were free but none lay within 45 deg of + * the current direction (a geometry problem). Only the second says the mesh lost its + * wall-normal structure. For the latter, also accumulate the best alignment on offer, which + * says whether the 45 deg threshold is merely too tight or the direction is truly lost. ---*/ + unsigned long nStopAllTaken = 0, nStopMisaligned = 0; + su2double sumBestCos = 0.0; + if (nLinelet != 0) { /*--- Define the basic linelets, starting from each vertex, preventing duplication of points. ---*/ @@ -4375,11 +4389,23 @@ const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config) } li.linelets.resize(nLinelet); - /*--- Create the linelet structure. ---*/ + /*--- Grow the lines breadth first: each pass advances every still-growing line by one point. + * Growing them depth first - running one line to completion before starting the next - lets an + * early line exhaust its own column, turn sideways (the 45 deg test permits it) and consume the + * points its neighbours needed, starving them into one- and two-point stubs. That is harmless on + * the fine grid, where the boundary layer is deeper than MAX_LINELET_POINTS so no line ever runs + * out of vertical room, but on agglomerated grids the layer is shallower than the cap and the + * starvation is severe. Advancing in lockstep makes the lines compete on equal terms for the + * layer they are all entitled to. ---*/ - nLinelet = 0; - for (auto& linelet : li.linelets) { - while (linelet.size() < CLineletInfo::MAX_LINELET_POINTS) { + std::vector growing(nLinelet, 1); + + for (unsigned long step = 1; step < CLineletInfo::MAX_LINELET_POINTS; ++step) { + bool anyGrew = false; + + for (auto iLine = 0ul; iLine < nLinelet; ++iLine) { + if (!growing[iLine]) continue; + auto& linelet = li.linelets[iLine]; const auto iPoint = linelet.back(); /*--- Compute the value of the max and min weights to detect if this region is isotropic. ---*/ @@ -4398,26 +4424,35 @@ const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config) } /*--- Isotropic, stop this linelet. ---*/ - if (min_weight / max_weight > CLineletInfo::ALPHA_ISOTROPIC()) break; + if (min_weight / max_weight > CLineletInfo::ALPHA_ISOTROPIC()) { + growing[iLine] = 0; + ++nStopIsotropic; + continue; + } /*--- Otherwise, add the closest valid neighbor. ---*/ su2double min_dist2 = std::numeric_limits::max(); auto next_Point = iPoint; const auto* iCoord = nodes->GetCoord(iPoint); + unsigned long nFreeCandidates = 0; + su2double bestCos = -1.0; for (const auto jPoint : nodes->GetPoints(iPoint)) { if (li.lineletIdx[jPoint] == CLineletInfo::NO_LINELET && nodes->GetDomain(jPoint)) { + ++nFreeCandidates; const auto* jCoord = nodes->GetCoord(jPoint); const su2double d2 = GeometryToolbox::SquaredDistance(nDim, iCoord, jCoord); su2double cosTheta = 1; + su2double dij[3] = {0.0}; + GeometryToolbox::Distance(nDim, jCoord, iCoord, dij); if (linelet.size() > 1) { const auto* kCoord = nodes->GetCoord(linelet[linelet.size() - 2]); - su2double dij[3] = {0.0}, dki[3] = {0.0}; + su2double dki[3] = {0.0}; GeometryToolbox::Distance(nDim, iCoord, kCoord, dki); - GeometryToolbox::Distance(nDim, jCoord, iCoord, dij); cosTheta = GeometryToolbox::DotProduct(3, dki, dij) / sqrt(d2 * GeometryToolbox::SquaredNorm(nDim, dki)); } + bestCos = max(bestCos, cosTheta); if (d2 < min_dist2 && cosTheta > 0.7071) { next_Point = jPoint; min_dist2 = d2; @@ -4426,28 +4461,61 @@ const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config) } /*--- Did not find a suitable point. ---*/ - if (next_Point == iPoint) break; + if (next_Point == iPoint) { + growing[iLine] = 0; + ++nStopNoNeighbour; + if (nFreeCandidates == 0) { + ++nStopAllTaken; + } else { + ++nStopMisaligned; + sumBestCos += bestCos; + } + continue; + } linelet.push_back(next_Point); - li.lineletIdx[next_Point] = nLinelet; + li.lineletIdx[next_Point] = iLine; + anyGrew = true; } - ++nLinelet; - maxNPoints = max(maxNPoints, linelet.size()); - sumNPoints += linelet.size(); + if (!anyGrew) break; + } + + /*--- A line that never stopped advancing ran into the length cap. ---*/ + for (auto iLine = 0ul; iLine < nLinelet; ++iLine) { + if (growing[iLine]) ++nStopCap; + maxNPoints = max(maxNPoints, li.linelets[iLine].size()); + sumNPoints += li.linelets[iLine].size(); } } /*--- Average linelet size over all ranks. ---*/ - unsigned long globalNPoints, globalNLineLets; + unsigned long globalNPoints, globalNLineLets, globalMaxNPoints; + unsigned long stopCounts[5] = {nStopIsotropic, nStopNoNeighbour, nStopCap, nStopAllTaken, nStopMisaligned}; + unsigned long globalStop[5] = {}; SU2_MPI::Allreduce(&sumNPoints, &globalNPoints, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); SU2_MPI::Allreduce(&nLinelet, &globalNLineLets, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); - - if (rank == MASTER_NODE) { - std::cout << "Computed linelet structure, " + SU2_MPI::Allreduce(&maxNPoints, &globalMaxNPoints, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(stopCounts, globalStop, 5, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + su2double globalSumBestCos = 0.0; + SU2_MPI::Allreduce(&sumBestCos, &globalSumBestCos, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + if (rank == MASTER_NODE && globalNLineLets > 0) { + const auto pct = [&](unsigned long n) { return 100.0 * passivedouble(n) / globalNLineLets; }; + std::cout << "Computed linelet structure on MG level " << MGLevel << ", " << static_cast(passivedouble(globalNPoints) / globalNLineLets) - << " points in each line (average)." << std::endl; + << " points in each line (average), " << globalMaxNPoints << " longest, " << globalNLineLets + << " lines.\n" + << " Line ends because: " << pct(globalStop[0]) << "% mesh became isotropic, " << pct(globalStop[1]) + << "% no aligned neighbour, " << pct(globalStop[2]) << "% hit the " << CLineletInfo::MAX_LINELET_POINTS + << "-point cap.\n" + << " of the 'no aligned neighbour': " << pct(globalStop[3]) + << "% all candidates already claimed by another line, " << pct(globalStop[4]) + << "% candidates free but misaligned"; + if (globalStop[4] > 0) + std::cout << " (best cos on offer " << globalSumBestCos / globalStop[4] << ", need > 0.7071)"; + std::cout << "." << std::endl; } /*--- Color the linelets for OpenMP parallelization and visualization. ---*/ diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 684da742b13..c056215e17c 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -88,11 +88,28 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } + /*--- STEP 0: agglomerate the stretched layer above viscous walls along implicit lines, wall CV + * included. This runs before the general boundary agglomeration below so that the wall control + * volume and the layers stacked on top of it share one footprint; letting the general scheme + * claim the wall first would fix a footprint chosen without any knowledge of the lines, and the + * stack above it could then only be misaligned with its own base. Everything it claims is + * already marked agglomerated, so the boundary and interior passes below simply skip it. ---*/ + if (config->GetMGOptions().MG_Implicit_Lines) { + AgglomerateImplicitLines(Index_CoarseCV, fine_grid, config, MGQueue_InnerCV); + } + /*--- 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. Carrying one does not put a point on a boundary, it only + * records that the point is mirrored on another rank. A point whose only markers are + * SEND_RECEIVE is an interior point, and is left to the domain pass (STEP 2) which is + * where a serial run would agglomerate it too. Points that do sit on a physical boundary + * are still reached here through their physical marker. ---*/ + 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 +135,18 @@ 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. SEND_RECEIVE markers are deliberately not counted: including them would make + an ordinary wall point look like a ridge, and a wall/symmetry ridge look like a corner + (which the counter > 2 rule below then refuses to agglomerate at all), so a point would be + classified differently depending only on where the partition happens to cut. ---*/ 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 +185,21 @@ 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 feature where two surface markers meet). + Always allow it to attempt agglomeration here; SetBoundAgglomeration() enforces + the actual ridge-ridge rule downstream: it may only pair with a neighboring ridge + point that carries the identical physical marker pair. A mismatched marker pair + usually indicates a genuine sharp corner in the geometry and is correctly left + un-merged (falls through to the singleton leftover loop). ---*/ + 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 @@ -281,6 +301,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,11 +336,6 @@ 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. ---*/ auto iteration = 0ul; @@ -433,8 +453,30 @@ 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. ---*/ + + 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 ((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); @@ -651,28 +693,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); @@ -766,18 +811,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 +835,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 +849,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. ---*/ @@ -1283,229 +1332,552 @@ su2double CMultiGridGeometry::ComputeLocalCurvature(const CGeometry* fine_grid, 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 su2double ANGLE_THRESHOLD_DEG = 20.0; /*!< Stop a line if the 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); + const bool ISOTROPIC = config->GetMGOptions().MG_Implicit_Lines_Isotropic; const unsigned long nPointFine = fine_grid->GetnPoint(); + const unsigned long starting_Index_CoarseCV = Index_CoarseCV; + const bool DEBUG_OUTPUT = (rank == MASTER_NODE); + + /*--- How many parallel lines one coarse CV may span tangential to the wall. In 2D a wall "face" is + * an edge with 2 end nodes, in 3D a quadrilateral with 4 corner nodes, which is the number of + * lines that must be bundled to coarsen by 2 in every wall-tangential direction. ---*/ + unsigned long max_group = config->GetMGOptions().MG_Implicit_Lines_Max_Group; + if (max_group == 0) max_group = (nDim == 2) ? 2 : 4; + + /*--- Smallest local cell aspect ratio that still counts as a stretched layer. ---*/ + const su2double MIN_AR = config->GetMGOptions().MG_Implicit_Lines_Min_AR; + const bool USE_AR = (MIN_AR > 1.0); + + /*--- Strength of the coupling across the dual face between a node and one of its neighbours. For a + * cell of streamwise size dx and wall-normal size dy this is 1/dy across the wall-normal face and + * 1/dx across the tangential one, so the ratio of the largest weight at a node to the smallest is + * the local cell aspect ratio. That makes the aspect ratio available from the dual grid alone, + * which SetControlVolume builds on every multigrid level, whereas CGeometry::Aspect_Ratio exists + * only on MESH_0. The same quantity already decides where LINELET preconditioner lines stop, in + * CGeometry::GetLineletInfo. ---*/ + auto edgeWeight = [&](unsigned long iPoint, unsigned short 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)); + return 0.5 * area * (1.0 / fine_grid->nodes->GetVolume(iPoint) + 1.0 / fine_grid->nodes->GetVolume(jPoint)); + }; + + /*--- Weakest coupling at a node, i.e. the denominator of the local aspect ratio. ---*/ + auto minEdgeWeight = [&](unsigned long iPoint) { + su2double wmin = std::numeric_limits::max(); + for (auto iNeigh = 0u; iNeigh < fine_grid->nodes->GetnPoint(iPoint); iNeigh++) + wmin = std::min(wmin, edgeWeight(iPoint, iNeigh)); + return wmin; + }; + + /*--- Aspect ratio of the mesh at iPoint measured along the edge to jPoint, and the neighbour the + * stiffest edge leads to. Taking the weight of one specific edge over the weakest edge at the + * node, rather than the largest over the smallest, keeps the measure directional: a mesh graded + * in the streamwise direction reads as stretched to the undirected form even far from any wall, + * which is why GetLineletInfo's min/max test cannot be used to decide where a line ends. ---*/ + auto aspectRatioAlong = [&](unsigned long iPoint, unsigned long jPoint) { + const su2double wmin = minEdgeWeight(iPoint); + if (wmin <= 0.0) return su2double(1.0); + for (auto iNeigh = 0u; iNeigh < fine_grid->nodes->GetnPoint(iPoint); iNeigh++) + if (fine_grid->nodes->GetPoint(iPoint, iNeigh) == jPoint) return edgeWeight(iPoint, iNeigh) / wmin; + return su2double(1.0); + }; + + /*================================================================================================== + * PHASE A - build the implicit lines. + * + * Lines are grown one step at a time across ALL lines simultaneously rather than one line to + * completion at a time, and a node is claimed globally the moment any line takes it. Growing them + * one-at-a-time lets an early line run the full depth of the layer and consume nodes that a later, + * neighbouring line needed, so that later line terminates after a step or two; the lines then have + * wildly different lengths and cannot be bundled into columns of uniform depth. Advancing in + * lockstep makes all lines compete for each layer on equal terms, which on an extruded prismatic + * layer reproduces the mesh's own structure: every line reaches the same depth. + *================================================================================================*/ + vector> lines; /*!< lines[i] = [wall_node, interior_1, interior_2, ...] */ + vector dir; /*!< Current marching direction of each line, nDim per line. */ + vector claimed(nPointFine, 0); /*!< Node already belongs to some line. */ + + /*--- Nodes that sit on a boundary with a boundary condition on it. A line must not grow into one, + * because those nodes belong to the boundary agglomeration and a stack that absorbed one would + * straddle two boundaries. CPoint's Boundary flag cannot answer this on its own: it is set for + * every marker a node belongs to, SEND_RECEIVE included, so on a partitioned mesh it is also + * true for the ordinary interior nodes of the send fringe. Testing it directly would stop every + * line that reaches the fringe one layer short of the partition, leaving the top of those + * columns to isotropic agglomeration purely because of where the mesh was cut. ---*/ + vector onPhysicalBoundary(nPointFine, 0); + for (auto iPoint = 0ul; iPoint < nPointFine; ++iPoint) { + if (!fine_grid->nodes->GetBoundary(iPoint)) continue; + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + if (fine_grid->nodes->GetVertex(iPoint, iMarker) != -1) { + onPhysicalBoundary[iPoint] = 1; + break; + } + } + } - /*--- 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; - + /*--- Seed a line at iPoint growing away from the boundary along unitNormal. ---*/ + auto seedLine = [&](unsigned long iPoint, const su2double* unitNormal) { + lines.push_back({iPoint}); + claimed[iPoint] = 1; + for (unsigned short d = 0; d < nDim; ++d) dir.push_back(unitNormal[d]); + }; + + /*--- Unit normal of the boundary at a vertex, false if the marker does not reach iPoint. ---*/ + auto vertexNormal = [&](unsigned long iPoint, unsigned short iMarker, su2double* unitNormal) { + const long ChildVertex = fine_grid->nodes->GetVertex(iPoint, iMarker); + if (ChildVertex == -1) return false; + fine_grid->vertex[iMarker][ChildVertex]->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; + }; + + /*--- Viscous walls always carry a stretched layer, so they seed unconditionally. Running them + * first also settles the nodes where a wall meets another boundary: the wall claims them, and + * the line there is a wall line. ---*/ 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; 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 (claimed[iPoint]) continue; /*--- A node on two wall markers must seed only one line. ---*/ - /*--- 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 (!vertexNormal(iPoint, iMarker, Normal)) continue; + seedLine(iPoint, Normal); + } + } + + /*================================================================================================== + * Boundaries other than viscous walls that nevertheless carry a stretched layer normal to + * themselves. A symmetry plane laid in the same surface as a wall, such as the one ahead of a + * flat plate's leading edge or the flat floor sections either side of a bump, is meshed with the + * very same normal spacing as the wall it continues. Seeding only from walls leaves the mesh + * above it to isotropic agglomeration, so the coarse grid changes character across the line where + * the two meet even though the fine grid does not, and that shows up as a residual there. + * + * A node qualifies when its stiffest edge is both stretched and points along the boundary normal, + * which is what distinguishes a layer growing off this boundary from one merely passing by: on + * the side planes of a bump the mesh is just as stretched, but in the wall-normal direction that + * runs ALONG the plane, and those nodes belong to the wall's own lines. + * + * The decision is then taken per marker rather than per node. Seeding individual qualifying nodes + * on a marker that mostly does not qualify scatters isolated lines across a face whose neighbours + * seed nothing, and those become one-line bundles, i.e. coarse CVs one fine CV wide that do not + * coarsen tangentially at all. Measured on a flat plate and a 3D bump the two populations are far + * apart - boundaries with a layer normal to them qualify at 100%, while side planes, inlets, + * outlets and far fields come in at 13% and below - so any threshold near a half separates them. + *================================================================================================*/ + if (USE_AR) { + const auto nMarkerFine = fine_grid->GetnMarker(); + + /*--- Counts are kept per marker of the configuration file, not per local marker. Ranks do not + * agree on either the number of markers or their order, because the SEND_RECEIVE markers of a + * partition are appended to its own list, so the same index means a different boundary on + * another rank and a reduction over it would add unrelated boundaries together. The + * configuration file list is the same everywhere. ---*/ + const auto nMarkerCfg = config->GetnMarker_CfgFile(); + vector nValid(nMarkerCfg, 0), nQualified(nMarkerCfg, 0); + + /*--- True if the mesh at iPoint is stretched along the boundary normal, i.e. this boundary has a + * layer growing off it in the same way a viscous wall does. ---*/ + auto hasLayerNormalTo = [&](unsigned long iPoint, const su2double* unitNormal) { + su2double wmax = 0.0, wmin = std::numeric_limits::max(); + unsigned long jStiffest = ULONG_MAX; + for (auto iNeigh = 0u; iNeigh < fine_grid->nodes->GetnPoint(iPoint); iNeigh++) { + const su2double w = edgeWeight(iPoint, iNeigh); + if (w > wmax) { + wmax = w; + jStiffest = fine_grid->nodes->GetPoint(iPoint, iNeigh); } + wmin = std::min(wmin, w); + } + if (jStiffest == ULONG_MAX || wmin <= 0.0) return false; + if (wmax / wmin < 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; + }; + + /*--- Markers that may be tested at all. Periodic boundaries are left out: the two halves are the + * same physical location under a transform and have their own matching, which a line running + * into one would interfere with. ---*/ + auto canSeed = [&](unsigned short iMarker) { + const auto bc = config->GetMarker_All_KindBC(iMarker); + if (bc == SEND_RECEIVE || bc == PERIODIC_BOUNDARY) return false; + return (bc != HEAT_FLUX && bc != ISOTHERMAL && bc != CHT_WALL_INTERFACE && bc != SMOLUCHOWSKI_MAXWELL); + }; + + /*--- Position of a local marker in the configuration file list. Only meaningful for the markers + * canSeed accepts: a SEND_RECEIVE marker is named per partition and is not in that list. ---*/ + 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 (!vertexNormal(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 has to be taken on the whole + * of it or two ranks could disagree about the same boundary. ---*/ + 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); + } - L.push_back(best_neighbor); + for (auto iMarker = 0u; iMarker < nMarkerFine; iMarker++) { + if (!canSeed(iMarker)) continue; + const auto iCfg = cfgOfMarker[iMarker]; + if (nValid[iCfg] == 0) continue; + if (2 * nQualified[iCfg] < nValid[iCfg]) continue; /*--- Fewer than half, not a layer. ---*/ + + 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 (claimed[iPoint]) continue; + + su2double Normal[MAXNDIM] = {0.0}; + if (!vertexNormal(iPoint, iMarker, Normal)) continue; + /*--- The marker carries a layer, but this node still has to be in it. ---*/ + if (!hasLayerNormalTo(iPoint, Normal)) continue; + seedLine(iPoint, Normal); + } + } + } - /*--- 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; + if (lines.empty()) return; + + vector growing(lines.size(), 1); + for (bool any_grew = true; any_grew;) { + any_grew = false; + for (unsigned long li = 0; li < lines.size(); ++li) { + if (!growing[li]) continue; + if (lines[li].size() >= MAX_LINE_LENGTH) { + growing[li] = 0; + continue; + } + const auto current = lines[li].back(); + su2double best_dot = -2.0; + unsigned long best_neighbor = ULONG_MAX; + + for (auto jPoint : fine_grid->nodes->GetPoints(current)) { + /*--- Halo nodes stay out: their parent is dictated by the rank that owns them and arrives + * through the MPI relay, so a line claiming one here would fight that assignment. A line + * therefore still ends at the partition itself, but now only there, instead of one layer + * earlier at the fringe of owned nodes. ---*/ + if (!fine_grid->nodes->GetDomain(jPoint)) continue; + if (onPhysicalBoundary[jPoint]) continue; + if (fine_grid->nodes->GetAgglomerate(jPoint)) continue; + if (claimed[jPoint]) continue; + + 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; + + const su2double dot = GeometryToolbox::DotProduct(nDim, vec, &dir[li * nDim]); + if (dot > best_dot) { + best_dot = dot; + best_neighbor = jPoint; + } + } + + if (best_neighbor == ULONG_MAX || best_dot < cos_threshold) { + growing[li] = 0; + continue; + } - current = best_neighbor; + /*--- End the line where the mesh stops being stretched along it. Without this the only limits + * are the direction cone and MAX_LINE_LENGTH, so a line leaves the boundary layer and keeps + * going into the far field, stacking coarse CVs along a direction the fine grid does not + * single out. Ordinary agglomeration handles that region better. ---*/ + if (USE_AR && aspectRatioAlong(current, best_neighbor) < MIN_AR) { + growing[li] = 0; + continue; } - /*--- Accept only lines with at least 2 interior nodes (length >= 3 including wall) ---*/ - if (L.size() >= 3) { - lines.push_back(std::move(L)); + su2double step[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(best_neighbor), fine_grid->nodes->GetCoord(current), + step); + const su2double slen = GeometryToolbox::Norm(nDim, step); + if (slen <= 0.0) { + growing[li] = 0; + continue; } + for (unsigned short d = 0; d < nDim; ++d) dir[li * nDim + d] = step[d] / slen; + + lines[li].push_back(best_neighbor); + claimed[best_neighbor] = 1; + any_grew = true; } } + /*--- A line needs at least one interior node to contribute anything. Drop the rest and release + * their wall seed so ordinary boundary agglomeration can treat it normally. ---*/ + { + vector> kept; + kept.reserve(lines.size()); + for (auto& L : lines) { + if (L.size() >= 2) + kept.push_back(std::move(L)); + else + claimed[L[0]] = 0; + } + lines = std::move(kept); + } if (lines.empty()) return; - if (rank == MASTER_NODE) { - cout << "Implicit line agglomeration: detected " << lines.size() << " lines." << endl; + /*================================================================================================== + * PHASE B - partition the lines into bundles. + * + * Every line must end up in exactly one bundle, and a bundle must be a compact patch on the wall: + * in 3D the four lines rising from the corners of one wall quadrilateral, in 2D the two lines from + * the ends of one wall edge. Selecting, for each line independently, some set of neighbours to + * merge with does not do this - the relation is not symmetric, so line 1 claiming {2,3,4} does not + * stop line 2 from claiming {1,3,5}, and the bundles overlap and fight over nodes. + * + * Building the partition by repeated pairwise matching avoids that by construction. One matching + * round pairs adjacent lines into 2-bundles (the wall edge); a second round pairs adjacent + * 2-bundles into 4-bundles (the wall quadrilateral). Each round is a matching, so membership is + * mutually exclusive at every stage and the result is a true partition. It also needs nothing but + * point-to-point connectivity, so it works identically on every multigrid level - boundary face + * connectivity does not exist on agglomerated grids, so a literal "same quadrilateral" test would + * only ever work for the first coarsening. + * + * Lines may only be bundled when their wall nodes carry the same set of physical markers, so that + * a bundle never straddles a boundary-condition change (the same rule ordinary agglomeration uses: + * ridges merge only with ridges, valleys only with valleys). + *================================================================================================*/ + const unsigned long nLines = lines.size(); + + /*--- Marker signature of each line's wall node. ---*/ + vector> sig(nLines); + for (unsigned long li = 0; li < nLines; ++li) { + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) continue; + if (fine_grid->nodes->GetVertex(lines[li][0], iMarker) != -1) sig[li].push_back(iMarker); + } + sort(sig[li].begin(), sig[li].end()); } - /*--- 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; + /*--- Line adjacency, inherited from the wall nodes' mesh connectivity. ---*/ + unordered_map lineOfWallNode; + lineOfWallNode.reserve(nLines); + for (unsigned long li = 0; li < nLines; ++li) lineOfWallNode[lines[li][0]] = li; - 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); + vector> adj(nLines); + for (unsigned long li = 0; li < nLines; ++li) { + for (auto jPoint : fine_grid->nodes->GetPoints(lines[li][0])) { + const auto it = lineOfWallNode.find(jPoint); + if (it != lineOfWallNode.end() && it->second != li) adj[li].push_back(it->second); } + } - vector line_processed(lines.size(), 0); - - /*--- A) Cross-line merges: parents with multiple lines ---*/ - for (auto& [parent, line_ids] : parent_to_lines) { - if (line_ids.size() < 2) continue; - - 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; - - 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; - - const auto a = L1[idx1], b = L1[idx2]; - const auto c = L2[idx1], d = L2[idx2]; - - /*--- 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)) - 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; - - /*--- 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; - } + /*--- Round 1: match adjacent lines into pairs. ---*/ + vector> bundles; + vector bundleOf(nLines, -1); + bundles.reserve(nLines); + for (unsigned long li = 0; li < nLines; ++li) { + if (bundleOf[li] >= 0) continue; + const long b = static_cast(bundles.size()); + bundles.push_back({li}); + bundleOf[li] = b; + if (max_group < 2) continue; + for (auto lj : adj[li]) { + if (bundleOf[lj] >= 0 || sig[lj] != sig[li]) continue; + bundles[b].push_back(lj); + bundleOf[lj] = b; + break; + } + } - /*--- 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); + /*--- Further rounds: merge adjacent bundles while they still fit. In 3D this turns pairs into + * quadrilaterals; in 2D max_group is 2 so it only absorbs leftover singletons. ---*/ + for (bool changed = true; changed;) { + changed = false; + + vector> badj(bundles.size()); + for (unsigned long li = 0; li < nLines; ++li) + for (auto lj : adj[li]) + if (bundleOf[li] != bundleOf[lj]) badj[bundleOf[li]].push_back(bundleOf[lj]); + for (auto& v : badj) { + sort(v.begin(), v.end()); + v.erase(unique(v.begin(), v.end()), v.end()); + } - 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; + vector consumed(bundles.size(), 0); + vector> merged; + merged.reserve(bundles.size()); + for (unsigned long b = 0; b < bundles.size(); ++b) { + if (consumed[b]) continue; + consumed[b] = 1; + auto group = bundles[b]; + for (auto h : badj[b]) { + if (consumed[h]) continue; + if (group.size() + bundles[h].size() > max_group) continue; + if (sig[bundles[h].front()] != sig[group.front()]) continue; + consumed[h] = 1; + group.insert(group.end(), bundles[h].begin(), bundles[h].end()); + changed = true; + break; } + merged.push_back(std::move(group)); } - /*--- 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); + bundles = std::move(merged); + for (unsigned long b = 0; b < bundles.size(); ++b) + for (auto li : bundles[b]) bundleOf[li] = static_cast(b); + } - Index_CoarseCV++; - any_work = true; + /*================================================================================================== + * PHASE C - extrude each bundle into a stack of coarse control volumes. + * + * The bundle's wall nodes become one coarse CV, and each successive layer of the bundle's lines + * becomes the next, so the coarse grid inherits the layer structure of the fine grid and a line + * relaxation remains meaningful on it. Because Phase A made the lines node-disjoint and Phase B + * made the bundles a partition, no two bundles can ever contend for the same node, so a stack is + * never interrupted part way up. + * + * The lines in a bundle need not be equally long. A stack therefore keeps rising for as long as + * enough of its lines have nodes left, narrowing where the shorter ones end, instead of stopping + * where the shortest one does and abandoning everything the taller ones still had. + * + * The multigrid queue is deliberately not updated here: the sync loop that follows the boundary + * agglomeration removes every point already marked agglomerated, so removing them a second time + * from this function would be an error. + *================================================================================================*/ + const unsigned long nBlock = ISOTROPIC ? 2 : 1; + + map bundle_size_histogram; + unsigned long nStacks = 0, nTruncated = 0; + + for (const auto& members : bundles) { + bundle_size_histogram[members.size()]++; + + /*--- The wall CV. Claiming it here, before ordinary boundary agglomeration runs, is what keeps + * the whole stack aligned: the layer above a wall CV has exactly the same footprint. ---*/ + bool valid = true; + for (auto li : members) + if (!GeometricalCheck(lines[li][0], fine_grid, config)) valid = false; + if (!valid) continue; + + for (unsigned long c = 0; c < members.size(); ++c) { + const auto p = lines[members[c]][0]; + fine_grid->nodes->SetParent_CV(p, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, c, p); } + nodes->SetnChildren_CV(Index_CoarseCV, static_cast(members.size())); + Index_CoarseCV++; + nStacks++; + + /*--- The interior layers, in lockstep. A line that runs out simply stops contributing and the + * others carry on without it, so a bundle is no longer cut down to its shortest member: the + * stack narrows as it rises instead of ending. The CVs still form one connected column, which + * is what a line relaxation on the coarse grid needs. + * + * Two situations end the stack rather than narrowing it further. Dropping below two lines + * would extrude a column one line wide, thinner than anything the domain pass would build + * there, and a set of lines that is no longer connected on the wall would put two separated + * columns into a single CV. In both cases the nodes above are better left to ordinary domain + * agglomeration. A bundle that only ever had one line is exempt from the first rule: it is + * one line wide by construction, and stopping early would gain nothing. ---*/ + + /*--- Are the still-growing lines one connected patch on the wall? Takes positions into members, + * which is at most max_group long, so the quadratic scan over Phase B's adjacency is cheap. ---*/ + auto isConnected = [&adj, &members](const vector& act) { + if (act.size() <= 1) return true; + vector seen(act.size(), 0); + vector stack{0}; + seen[0] = 1; + unsigned long nSeen = 1; + while (!stack.empty()) { + const auto cur = stack.back(); + stack.pop_back(); + const auto& neighbors = adj[members[act[cur]]]; + for (unsigned long k = 0; k < act.size(); ++k) { + if (seen[k]) continue; + if (find(neighbors.begin(), neighbors.end(), members[act[k]]) != neighbors.end()) { + seen[k] = 1; + nSeen++; + stack.push_back(k); + } + } + } + return nSeen == act.size(); + }; - pair_idx++; - if (!any_work) break; + const unsigned long minActive = std::min(2, members.size()); - /*--- 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; + vector placed(members.size(), 1); /*!< First node of each line not yet in a CV. */ + vector active, group; + + for (unsigned long first = 1;; first += nBlock) { + /*--- The lines that still have a whole block left at this height. ---*/ + active.clear(); + for (unsigned long m = 0; m < members.size(); ++m) + if (first + nBlock <= lines[members[m]].size()) active.push_back(m); + + if (active.size() < minActive) break; + if (!isConnected(active)) break; + + group.clear(); + group.reserve(active.size() * nBlock); + for (auto m : active) + for (unsigned long b = 0; b < nBlock; ++b) group.push_back(lines[members[m]][first + b]); + + valid = true; + for (auto p : group) + if (!GeometricalCheck(p, fine_grid, config)) valid = false; + if (!valid) break; + + for (unsigned long c = 0; c < group.size(); ++c) { + fine_grid->nodes->SetParent_CV(group[c], Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, c, group[c]); } + nodes->SetnChildren_CV(Index_CoarseCV, static_cast(group.size())); + Index_CoarseCV++; + for (auto m : active) placed[m] = first + nBlock; } - if (!any_more) break; + + /*--- Whatever each line still carries above the last CV it contributed to. ---*/ + for (unsigned long m = 0; m < members.size(); ++m) nTruncated += lines[members[m]].size() - placed[m]; + } + + if (DEBUG_OUTPUT) { + unsigned long nLineNodes = 0; + for (const auto& L : lines) nLineNodes += L.size(); + cout << " Implicit lines: " << nLines << " lines, " << nStacks << " stacks, bundle sizes "; + for (const auto& h : bundle_size_histogram) cout << h.first << "x" << h.second << " "; + cout << "\n Coarse CVs from lines: " << (Index_CoarseCV - starting_Index_CoarseCV) << " covering " + << (nLineNodes - nTruncated) << "/" << nLineNodes << " line nodes"; + if (nTruncated > 0) cout << " (" << nTruncated << " left to domain agglomeration)"; + cout << endl; } } diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index f7d2d36f198..38bf0ebea8a 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -7379,7 +7379,7 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { /*--- Some recommended defaults for the various ParMETIS options. ---*/ - idx_t wgtflag = 2; + idx_t wgtflag = 2; /*--- Weights on the vertices only, raised to 3 below if edge weights are built. ---*/ idx_t numflag = 0; idx_t ncon = 1; real_t ubvec = 1.0 + config->GetParMETIS_Tolerance(); @@ -7412,6 +7412,149 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { vwgt[iPoint] = wp + we * (xadj[iPoint + 1] - xadj[iPoint]); } + /*--- Cost of cutting each edge of the graph. + * + * Without these ParMETIS is given no edge weights at all and every edge is equally cheap to + * cut, so nothing stops a partition boundary from running straight through the stretched cells + * of a boundary layer and splitting the wall-normal columns that implicit-line agglomeration + * and line-implicit smoothing are built on. + * + * In a stretched cell the wall-normal spacing is the small one, so the short edges are exactly + * the ones that should stay inside a partition. Weighting an edge by the inverse of its length + * therefore makes cutting across the layer expensive and leaves the long tangential edges as + * the cheap place to cut. Length is used rather than the face area over volume ratio that + * measures the same thing elsewhere, because the dual grid does not exist yet at this point of + * the setup: this runs before SetControlVolume, and only the coordinates are available. + * + * On a mesh with no stretching every edge is of similar length, so the weights come out + * uniform and minimizing their sum is the same problem as minimizing the number of cut edges. + * Such a mesh is therefore partitioned exactly as it is with no weights at all. ---*/ + + vector adjwgt; + const su2double anisoWgt = config->GetParMETIS_AnisoWeight(); + + if (anisoWgt > 0.0) { + const auto firstIdx = pointPartitioner.GetFirstIndexOnRank(rank); + const auto lastIdx = pointPartitioner.GetLastIndexOnRank(rank); + auto isLocal = [&](unsigned long g) { return g >= firstIdx && g < lastIdx; }; + + /*--- The graph is split linearly and its entries are global indices, so an edge near a linear + * partition boundary has one end that is not stored here. Those are few, of the order of a + * percent of the entries, and each is asked for from the rank the linear partitioner says + * owns it. The same request lists are used twice, once for coordinates and once for the + * longest edge at the point, which only its owner can work out. ---*/ + vector> wanted(size); + for (auto gPoint : adjacency) + if (!isLocal(gPoint)) wanted[pointPartitioner.GetRankContainingIndex(gPoint)].push_back(gPoint); + + vector nSend(size, 0), nRecv(size, 0), sDisp(size + 1, 0), rDisp(size + 1, 0); + for (int r = 0; r < size; ++r) { + auto& w = wanted[r]; + sort(w.begin(), w.end()); + w.erase(unique(w.begin(), w.end()), w.end()); + nSend[r] = static_cast(w.size()); + } + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, comm); + for (int r = 0; r < size; ++r) { + sDisp[r + 1] = sDisp[r] + nSend[r]; + rDisp[r + 1] = rDisp[r] + nRecv[r]; + } + + vector sendIdx(sDisp[size]), recvIdx(rDisp[size]); + for (int r = 0; r < size; ++r) copy(wanted[r].begin(), wanted[r].end(), sendIdx.begin() + sDisp[r]); + SU2_MPI::Alltoallv(sendIdx.data(), nSend.data(), sDisp.data(), MPI_UNSIGNED_LONG, recvIdx.data(), nRecv.data(), + rDisp.data(), MPI_UNSIGNED_LONG, comm); + + /*--- Round one, the coordinates of the points that were asked for. ---*/ + map> remoteCoord; + { + vector sendBuf(static_cast(rDisp[size]) * nDim), + recvBuf(static_cast(sDisp[size]) * nDim); + for (size_t i = 0; i < recvIdx.size(); ++i) + for (unsigned short iDim = 0; iDim < nDim; ++iDim) + sendBuf[i * nDim + iDim] = nodes->GetCoord(recvIdx[i] - firstIdx, iDim); + + vector nS(size), nR(size), sD(size), rD(size); + for (int r = 0; r < size; ++r) { + nS[r] = nRecv[r] * nDim; + nR[r] = nSend[r] * nDim; + sD[r] = rDisp[r] * nDim; + rD[r] = sDisp[r] * nDim; + } + SU2_MPI::Alltoallv(sendBuf.data(), nS.data(), sD.data(), MPI_DOUBLE, recvBuf.data(), nR.data(), rD.data(), + MPI_DOUBLE, comm); + for (size_t i = 0; i < sendIdx.size(); ++i) { + array c = {0.0, 0.0, 0.0}; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) c[iDim] = recvBuf[i * nDim + iDim]; + remoteCoord[sendIdx[i]] = c; + } + } + + /*--- Length of every edge of the local part of the graph, and the longest edge at each point + * this rank owns. ---*/ + vector edgeLen(adjacency.size(), 0.0); + vector maxLen(nPoint, 0.0); + + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; ++k) { + const auto gPoint = adjacency[k]; + array tmp = {0.0, 0.0, 0.0}; + const su2double* coord_j = nullptr; + if (isLocal(gPoint)) { + coord_j = nodes->GetCoord(gPoint - firstIdx); + } else { + const auto it = remoteCoord.find(gPoint); + if (it == remoteCoord.end()) continue; + tmp = it->second; + coord_j = tmp.data(); + } + edgeLen[k] = GeometryToolbox::Distance(nDim, nodes->GetCoord(iPoint), coord_j); + maxLen[iPoint] = max(maxLen[iPoint], edgeLen[k]); + } + } + + /*--- Round two, the longest edge at each of the remote points. ---*/ + map remoteMaxLen; + { + vector sendBuf(rDisp[size]), recvBuf(sDisp[size]); + for (size_t i = 0; i < recvIdx.size(); ++i) sendBuf[i] = maxLen[recvIdx[i] - firstIdx]; + SU2_MPI::Alltoallv(sendBuf.data(), nRecv.data(), rDisp.data(), MPI_DOUBLE, recvBuf.data(), nSend.data(), + sDisp.data(), MPI_DOUBLE, comm); + for (size_t i = 0; i < sendIdx.size(); ++i) remoteMaxLen[sendIdx[i]] = recvBuf[i]; + } + + /*--- An edge is expensive to cut when it is much shorter than the other edges meeting it, which + * is the definition of the local cell aspect ratio and is exactly the situation inside a + * boundary layer, where the short edges are the wall-normal ones. Comparing an edge only + * against its own neighbourhood, rather than against a global length, is what keeps the + * measure a ratio: scaling the whole mesh, or refining one region of it isotropically, leaves + * every weight unchanged. Weighting by absolute length instead would make any small cell + * expensive to cut and would steer the partitioner away from refined regions that are not + * stretched at all. Averaging the two ends keeps the weight of an edge symmetric, which + * ParMETIS requires. ---*/ + const idx_t MAX_EDGE_WEIGHT = 1000; + adjwgt.resize(adjacency.size(), 1); + + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; ++k) { + if (edgeLen[k] <= 0.0) continue; + const auto gPoint = adjacency[k]; + su2double maxLen_j = 0.0; + if (isLocal(gPoint)) { + maxLen_j = maxLen[gPoint - firstIdx]; + } else { + const auto it = remoteMaxLen.find(gPoint); + if (it == remoteMaxLen.end()) continue; + maxLen_j = it->second; + } + const su2double ratio = 0.5 * (maxLen[iPoint] + maxLen_j) / edgeLen[k]; + const su2double w = 1.0 + anisoWgt * (ratio - 1.0); + adjwgt[k] = static_cast(min(max(w, 1.0), MAX_EDGE_WEIGHT)); + } + } + wgtflag = 3; /*--- Weights on both the vertices and the edges. ---*/ + } + /*--- Create some structures that ParMETIS needs to output the partitioning. ---*/ idx_t edgecut; @@ -7420,9 +7563,9 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { /*--- Calling ParMETIS ---*/ if (rank == MASTER_NODE) cout << "Calling ParMETIS..."; - auto err = - ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjacency.data(), vwgt.data(), nullptr, &wgtflag, &numflag, - &ncon, &nparts, tpwgts.data(), &ubvec, options, &edgecut, part.data(), &comm); + auto err = ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjacency.data(), vwgt.data(), + adjwgt.empty() ? nullptr : adjwgt.data(), &wgtflag, &numflag, &ncon, &nparts, + tpwgts.data(), &ubvec, options, &edgecut, part.data(), &comm); if (err != METIS_OK) SU2_MPI::Error("Partitioning failed.", CURRENT_FUNCTION); if (rank == MASTER_NODE) { cout << " graph partitioning complete (" << edgecut << " edge cuts)." << endl; diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 3c9176d8716..aa6ca89515b 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1435,6 +1435,36 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con const bool nested = SetupInnerSolver(KindSolver, config); + /*--- Decide whether to rebuild the preconditioner for this solve. + * + * On coarse multigrid levels the Jacobian changes little between smoothing sweeps, yet the + * factorization is rebuilt from scratch on every one of them. Profiling a V-cycle shows the + * ILU build is the single most expensive zone in the cycle, so MG_COARSE_PREC_FREEZE lets a + * factorization be reused for several consecutive solves. The factorization itself lives in + * the CSysMatrix, which outlives the CPreconditioner object created below, so not calling + * Build() is all that is needed to reuse it. + * + * Restricted to the standard solver mode (mesh deformation and gradient smoothing are left + * alone). Coarse levels and the finest grid have separate periods because the finest-grid + * preconditioner drives the outer nonlinear convergence and so carries more risk. The first + * solve on this instance always builds (count 0), which matters because the factorization is + * otherwise uninitialized. + * + * The decision is taken by one thread and read by all of them, because Build() is internally + * OpenMP-parallel and every thread must make the same choice. ---*/ + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + unsigned long freeze = 1; + if (lin_sol_mode == LINEAR_SOLVER_MODE::STANDARD && geometry != nullptr) { + freeze = (geometry->GetMGLevel() != MESH_0) ? config->GetMGOptions().MG_Coarse_Prec_Freeze + : config->GetLinear_Solver_Prec_Freeze(); + freeze = std::max(1, freeze); + } + buildPrecThisSolve = (precSolveCount % freeze == 0); + precSolveCount++; + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + /*--- Stop the recording for the linear solver ---*/ bool TapeActive = NO; @@ -1472,7 +1502,7 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con const auto kindPrec = static_cast(KindPrecond); auto* normal_prec = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); - normal_prec->Build(); + if (buildPrecThisSolve) normal_prec->Build(); CPreconditioner* nested_prec = nullptr; if (nested) { diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index eee8c7235fd..8b607fb5c01 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -321,6 +321,13 @@ class CMultiGridIntegration final : public CIntegration { static constexpr int MAX_MG_LEVELS = 10; + /*--- Upper bound on nVar for the small per-point scratch arrays used by the restriction and + * prolongation loops, so they can live on the stack instead of being heap-allocated on every + * call. Must be >= the largest MAXNVAR of any variable class integrated by this class, + * currently CNEMOEulerVariable::MAXNVAR = 25. ---*/ + 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. */ passivedouble mg_initial_smooth_rms = 0.0; /*!< \brief Initial RMS residual before current smoothing phase (FAS). */ diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index fecbf0491b7..ac0a8114a32 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -63,6 +63,33 @@ 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, so the same loop serves both the FAS + * correction and the Full-MG solution handoff. + * + * The loop covers all coarse points, halos included. Halo coarse CVs own the fine halo points as + * children (CMultiGridGeometry sets Children_CV for received CVs), so injecting from them is what + * fills the fine-grid halo entries of the prolongated field. Restricting the loop to domain points + * leaves those entries at whatever the last solver update left there (zero, for LinSysRes), which + * is wrong for any operator that reads the prolongated field at neighbours across a partition + * boundary - the Jacobi smoother in SmoothProlongated_Correction does exactly that. The caller + * must therefore have synchronized the coarse-grid field being read before calling this. + \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 +516,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 +540,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(); @@ -653,6 +685,19 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetKind_TimeIntScheme(EULER_IMPLICIT);) } + /*--- NOTE: the coarse-grid residual computed just above is evaluated at the restricted + * solution, i.e. at exactly the state the first pre-smoothing sweep of the recursive call + * below re-evaluates it at, so it looks like that sweep could reuse LinSysRes (and, if the + * Jacobian were assembled here, the Jacobian too) and skip its own Preprocessing and + * Space_Integration. It cannot, as things stand: Space_Integration is not a pure producer of + * LinSysRes/Jacobian. BC_Sym_Plane (which serves both SYMMETRY_PLANE and EULER_WALL) also + * projects Res_TruncError and Solution_Old onto the wall tangent plane, and in the current + * ordering that projection is what makes the FAS forcing term written by SetForcing_Term + * below, and the Solution_Old written by Set_OldSolution, wall-consistent before their first + * use. Reusing the residual moves both projections to the wrong side of the writes. + * Factoring those side effects out of Space_Integration would make the reuse safe and save + * one full residual evaluation per coarse level per cycle. ---*/ + /*--- Recursive call to MultiGrid_Cycle (this routine). ---*/ /*--- Execute multigrid cycles sequentially to ensure deterministic recursion order ---*/ /*--- This prevents accumulation of floating-point variations across recursive calls ---*/ @@ -877,14 +922,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 +947,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 +974,22 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS } } - /*--- MPI the set solution old ---*/ + /*--- MPI the set solution old. Required: the loop above only writes domain points, and + * ProlongateField below injects from every coarse point including halos in order to fill the + * fine-grid halo entries of the correction. ---*/ 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 (held in Solution_Old) onto the fine + * grid and store it in LinSysRes, which SetProlongated_Correction then damps + * and adds to the fine-grid solution. ---*/ + + 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); + }); } @@ -957,6 +1003,10 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ const unsigned short nVar = solver->GetnVar(); + /*--- Seeded over all points, halos included: the restore loop below reads Residual_Old at the + * vertices of the physical markers, and on a partitioned mesh some of those are halo points + * owned by another rank. ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); iPoint++) { const auto* Residual_Old = solver->LinSysRes.GetBlock(iPoint); @@ -969,10 +1019,13 @@ 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 (sum the residuals of direct neighbors). + * Halo points are deliberately not smoothed here: their own neighbor stencil is incomplete + * on this rank, so the average would be meaningless, and the halo exchange at the end of + * each sweep overwrites them with the value their owner computed anyway. ---*/ - 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 +1038,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 +1053,23 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } END_SU2_OMP_FOR - /*--- Restore original residuals (without average) at boundary points. ---*/ + /*--- Restore original residuals (without average) at physical boundary points. + * + * SEND_RECEIVE is excluded: carrying such a marker does not put a point on a boundary, it + * only records that the point is mirrored on another rank. Restoring those points froze the + * correction on the whole send fringe, which is exactly the ring of domain points that have + * a halo neighbour, so the smoothing this function applied depended on where the mesh + * happened to be partitioned rather than on the geometry alone. + * + * Note this removes one source of rank-dependence, not all of them: the coarse grids are + * agglomerated per rank, so the multigrid operator itself still differs between partition + * counts and a run on 1 and on N ranks is not expected to match bit for bit. ---*/ 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 +1081,19 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } } + /*--- Refresh the halo entries of the correction with the values their owner ranks just + * computed. The next sweep averages LinSysRes over the neighbours of every domain point, + * and across a partition boundary those neighbours are halo points, so this has to run + * once per sweep rather than once at the end. It comes after the restore so that a halo + * point sitting on a physical boundary mirrors its owner's restored value. + * + * The barrier is required: the restore loop above only carries an implicit barrier for the + * markers that pass the test, so if the last marker is skipped there is none. ---*/ + + SU2_OMP_BARRIER + CSysMatrixComms::Initiate(solver->LinSysRes, geometry, config); + CSysMatrixComms::Complete(solver->LinSysRes, geometry, config); + } /*--- Record final correction norm for debugging output. ---*/ @@ -1035,19 +1112,48 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet /*--- Use the adaptive damping factor uniformly across all prolongation levels. ---*/ const su2double factor = config->GetDamp_Correc_Prolong(); + /*--- Optional cap on how much one coarse-grid correction may move the solution at a point. + * Without it the only guard below is the NaN check, and a correction can drive a cell + * non-physical in a single application - measured on the turbulent flat plate, a W-cycle + * correction cuts wall-adjacent density by 25%, from which the energy equation never + * recovers. The cap is relative and per point, and the whole correction vector is scaled by + * one factor so its direction is preserved (scaling components independently would rotate + * the correction and break the coupling between the equations). + * + * Components that are negligible against the largest one at that point (transverse momentum + * in a freestream cell, say) carry no meaningful relative bound and are skipped, otherwise + * they would veto every correction. ---*/ + const su2double limit = config->GetMGOptions().MG_Correction_Limit; + const bool limiting = (limit > 0.0); + SU2_OMP_FOR_STAT(roundUpDiv(geo_fine->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Fine = 0ul; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) { auto* Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); auto* Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine); + + /*--- Prevent a fine grid divergence due to a coarse grid divergence ---*/ for (auto iVar = 0u; iVar < nVar; iVar++) { - /*--- Prevent a fine grid divergence due to a coarse grid divergence ---*/ if (Residual_Fine[iVar] != Residual_Fine[iVar]) Residual_Fine[iVar] = 0.0; + } - su2double correction = factor * Residual_Fine[iVar]; - - Solution_Fine[iVar] += correction; + su2double omega = 1.0; + if (limiting) { + su2double ref = 0.0; + for (auto iVar = 0u; iVar < nVar; iVar++) + ref = max(ref, fabs(Solution_Fine[iVar])); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double scale = fabs(Solution_Fine[iVar]); + if (scale < 1e-6 * ref) continue; + const su2double correction = fabs(factor * Residual_Fine[iVar]); + if (correction > limit * scale) + omega = min(omega, limit * scale / correction); + } } + + for (auto iVar = 0u; iVar < nVar; iVar++) + Solution_Fine[iVar] += omega * factor * Residual_Fine[iVar]; } END_SU2_OMP_FOR @@ -1171,26 +1277,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 +1392,15 @@ 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; + /*--- Row-major scratch plus the row pointers SetGradient expects. ---*/ + 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 +1415,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..0f051fdc2b6 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]; + if (RunTime_EqSystem == RUNTIME_TURB_SYS) { + /*--- CFL scaling of turbulence during the warmup phase if FMG. ---*/ + 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/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/config_template.cfg b/config_template.cfg index 535ae564510..2fdad98b989 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1740,6 +1740,11 @@ MG_IMPLICIT_LINES= NO % Increase to extend the line deeper into the boundary layer (default 20). MG_IMPLICIT_LINES_MAX_LENGTH= 20 % +% Use isotropic (vs anisotropic) agglomeration for implicit lines (NO, YES) +% Anisotropic (NO): Pair cells normal to wall (2 cells per coarse CV, reduces mesh ~2x) +% Isotropic (YES): Pair cells in all directions (4 cells per coarse CV, reduces mesh ~4x) +MG_IMPLICIT_LINES_ISOTROPIC= NO +% % 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