diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index ff68a7a024b9..992a0d6e30ad 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -122,6 +122,8 @@ jobs: with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} + - name: Lower ASLR entropy for TSan + run: sudo sysctl vm.mmap_rnd_bits=28 - name: Build uses: docker://ghcr.io/su2code/su2/build-su2-tsan:261011-1107 with: @@ -324,6 +326,8 @@ jobs: done chmod a+x $BIN_FOLDER/* ls -lahR $BIN_FOLDER + - name: Lower ASLR entropy for TSan + run: sudo sysctl vm.mmap_rnd_bits=28 - name: Run Tests in Container uses: docker://ghcr.io/su2code/su2/test-su2-tsan:261011-1107 env: diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index e85f783f1de6..5b3dd507b34e 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -601,7 +601,8 @@ class CConfig { MUSCL_Heat, /*!< \brief MUSCL scheme for the (fvm) heat equation.*/ MUSCL_AdjFlow, /*!< \brief MUSCL scheme for the adj flow equations.*/ MUSCL_AdjTurb, /*!< \brief MUSCL scheme for the adj turbulence equations.*/ - MUSCL_Species; /*!< \brief MUSCL scheme for the species equations.*/ + MUSCL_Species, /*!< \brief MUSCL scheme for the species equations.*/ + MUSCL_FirstOrderWallEnds; /*!< \brief First-order reconstruction at no-slip wall ends and sharp wall edges.*/ su2double MUSCL_Kappa, /*!< \brief Blending coefficient for U-MUSCL scheme (for the runtime eq. system). */ MUSCL_Kappa_Flow, /*!< \brief Blending coefficient for U-MUSCL scheme for the flow equations.*/ MUSCL_Kappa_Turb, /*!< \brief Blending coefficient for U-MUSCL scheme for the turbulence equations.*/ @@ -1117,6 +1118,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. */ + bool ParMETIS_anisotropyWgt; /*!< \brief Weight graph edges by cell anisotropy and keep anisotropic lines on one rank. */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ bool DiscreteAdjoint, /*!< \brief AD-based discrete adjoint mode. */ DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ @@ -1177,6 +1179,8 @@ class CConfig { unsigned short *MG_PreSmooth_p{nullptr}, *MG_PostSmooth_p{nullptr}, *MG_CorrecSmooth_p{nullptr}; unsigned short nMG_CflScaling_p{0}; su2double *MG_CflScaling_p{nullptr}; + unsigned short nMG_LinSolIter_p{0}; + unsigned short *MG_LinSolIter_p{nullptr}; ENUM_STREAMWISE_PERIODIC Kind_Streamwise_Periodic; /*!< \brief Kind of Streamwise periodic flow (pressure drop or massflow) */ bool Streamwise_Periodic_Temperature; /*!< \brief Use real periodicity for Energy equation or otherwise outlet source term. */ @@ -1752,6 +1756,10 @@ class CConfig { */ bool GetCFL_Adapt(void) const { return CFL_Adapt; } + /*! + * \brief Whether the adaptive CFL follows the nonlinear residual trend instead of the linear solve. + */ + /*! * \brief Get the outlier mitigation parameters. */ @@ -4859,6 +4867,11 @@ class CConfig { */ bool GetMUSCL_Flow(void) const { return MUSCL_Flow; } + /*! + * \brief Get if edges at no-slip wall ends and sharp wall edges are reconstructed to first order. + */ + bool GetMUSCL_FirstOrderWallEnds(void) const { return MUSCL_FirstOrderWallEnds; } + /*! * \brief Get if the upwind scheme used MUSCL or not. * \return MUSCL scheme. @@ -10334,6 +10347,11 @@ class CConfig { */ long GetParMETIS_EdgeWeight() const { return ParMETIS_edgeWgt; } + /*! + * \brief Whether ParMETIS weights graph edges by cell anisotropy and partitions with coordinates. + */ + bool GetParMETIS_AnisotropyWeight() const { return ParMETIS_anisotropyWgt; } + /*! * \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/containers/container_decorators.hpp b/Common/include/containers/container_decorators.hpp index 331e87b9a0f4..82b89814960a 100644 --- a/Common/include/containers/container_decorators.hpp +++ b/Common/include/containers/container_decorators.hpp @@ -126,7 +126,8 @@ class C3DContainerDecorator { private: Storage m_storage; - Index m_innerSz; + /*--- One, not zero, so rows() on a container that was never resized does not divide by it. ---*/ + Index m_innerSz = 1; public: C3DContainerDecorator() = default; diff --git a/Common/include/geometry/CMultiGridGeometry.hpp b/Common/include/geometry/CMultiGridGeometry.hpp index 3658dc0ea2f5..6e38ef2fea7a 100644 --- a/Common/include/geometry/CMultiGridGeometry.hpp +++ b/Common/include/geometry/CMultiGridGeometry.hpp @@ -110,7 +110,8 @@ class CMultiGridGeometry final : public CGeometry { vector> normal; /*!< \brief Unit normal there, pointing into the domain. */ vector strength; /*!< \brief Local anisotropy, the ordering key. */ vector tier; /*!< \brief Paving order: 0 viscous wall, 1 Euler - wall, 2 anything else. */ + wall, 2 anything else; one more once the + mixed nodes are put first. */ unsigned long nRefusedCurvature = 0; /*!< \brief Euler wall nodes the curvature limit kept out. */ }; diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index 52d8ea6da08e..e69b1a04a9b0 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -337,6 +337,42 @@ class CPhysicalGeometry final : public CGeometry { */ void PrepareAdjacency(const CConfig* config); +#if defined(HAVE_MPI) && defined(HAVE_PARMETIS) + /*! + * \brief ParMETIS edge weights that make cutting across the short edges of anisotropic cells expensive. + * \param[in] maxWeight - Cap on the weight of any edge. + * \return One weight per entry of the adjacency array. + */ + vector ComputeAnisotropyEdgeWeights(long maxWeight) const; + + /*! + * \brief Reads, for any global points of the linear partition, a value each rank holds for its own points. + * \param[in] wanted - Global indices to read. + * \param[in] localValue - Value of every local point. + * \return The value of each wanted point. + */ + vector FetchFromOwners(const vector& wanted, + const vector& localValue) const; + + /*! + * \brief Joins points into lines along the strong anisotropic edges and names each line by its lowest global index. + * \param[in] adjwgt - Anisotropy weight of every entry of the adjacency array. + * \param[out] lineEdge - Whether each entry of the adjacency array is a line edge. + * \return Global index of the root of the line of every local point, the point itself when it is on no line. + */ + vector LinkAnisotropicLines(const vector& adjwgt, vector& lineEdge) const; + + /*! + * \brief Moves each line to the partition that holds most of its points. + * \param[in] root - Line roots from LinkAnisotropicLines. + * \param[in,out] part - Partition of every local point. + * \param[out] longestLine - Number of points of the longest line. + * \return Number of points that changed partition, summed over all ranks. + */ + unsigned long KeepLinesTogether(const vector& root, vector& part, + unsigned long& longestLine) const; +#endif + /*! * \brief Find repeated nodes between two elements to identify the common face. * \param[in] first_elem - Identification of the first element. diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index d0aa6124d96a..663fdeeddf92 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -34,6 +34,7 @@ #include "CPastixWrapper.hpp" #include "../toolboxes/graph_toolbox.hpp" +#include #include #include #include @@ -363,6 +364,19 @@ class CSysMatrix { * into ilu_level_ptr / d_precond_level_idx below for the GPU triangular solves. */ CCompressedSparsePatternUL levels_ilu; + /*!< \brief Per-thread progress counter, padded to a cache line. */ + struct alignas(64) LevelProgress { + std::atomic steps{0}; + LevelProgress() = default; + LevelProgress(const LevelProgress& other) : steps(other.steps.load()) {} + }; + + /*!< \brief Steps each thread must see completed by every other thread, per (step, thread, other thread), + * for the forward and backward level scheduled sweeps. Empty when barriers separate the levels. */ + std::vector level_waits_fwd, level_waits_bwd; + mutable std::vector level_progress; + unsigned long level_num_threads = 0; + /*!< \brief Coloring of the (domain-only) ILU dependency graph, used only by the GPU iterative * factorization (see IluFactorColorKernel and ilu_color_ptr / d_ilu_color_idx below). The * host/OMP path and the GPU triangular solves use levels_ilu instead. */ @@ -666,6 +680,25 @@ class CSysMatrix { void ComputeJacobiPreconditionerGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Compute the wait tables of the point-to-point synchronized level scheduled sweeps. + */ + void BuildLevelWaits(); + + /*! + * \brief Apply a function to the rows of all levels in increasing or decreasing order, waiting only for + * the threads owning the rows depended on. Must be called by all threads of the team. + */ + template + void LevelSweep(bool forward, const F& row_func) const; + + /*! + * \brief Whether the calling thread team matches the one the wait tables were built for. + */ + inline bool UseLevelWaits() const { + return !level_waits_fwd.empty() && static_cast(omp_get_num_threads()) == level_num_threads; + } + /*! * \brief Build the ILU preconditioner on the device, from the device copy of the matrix. * \note Requires the device matrix to be up to date, see HtDTransfer. diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index 98ab059804ce..ca0f9db8e74d 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/Common/include/linear_algebra/CSysSolve.hpp @@ -455,9 +455,10 @@ class CSysSolve { * \param[in,out] LinSysSol - Linear system solution * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] maxIter - Iteration cap replacing LINEAR_SOLVER_ITER in standard mode, 0 keeps the config value. */ unsigned long Solve(MatrixType& Jacobian, const CSysVector& LinSysRes, CSysVector& LinSysSol, - CGeometry* geometry, const CConfig* config); + CGeometry* geometry, const CConfig* config, unsigned long maxIter = 0); /*! * \brief Solve the adjoint linear system using a Krylov subspace method diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1af43985c126..4fbfd376b1ad 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -80,6 +80,7 @@ const unsigned int MAX_NUMBER_FFD = 15; /*!< \brief Maximum number of FFDB enum: unsigned int{MAX_SOLS = 14}; /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */ const unsigned int MAX_TERMS = 7; /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */ const unsigned int MAX_ZONES = 3; /*!< \brief Maximum number of zones. */ +const unsigned short MAX_MGLEVELS = 10; /*!< \brief Maximum number of coarse multigrid levels, which bounds the per-level arrays of the multigrid integration. */ const unsigned int MAX_FE_KINDS = 7; /*!< \brief Maximum number of Finite Elements. */ const unsigned int NO_RK_ITER = 0; /*!< \brief No Runge-Kutta iteration. */ @@ -1246,10 +1247,13 @@ struct CMGOptions { std::vector MG_PostSmooth; /*!< \brief Multigrid post-smoothing iterations per level. */ std::vector MG_CorrecSmooth; /*!< \brief Multigrid Jacobi correction-smoothing per level. */ std::vector MG_CflScaling; /*!< \brief Per-level CFL scaling factors relative to the previous (finer) level. Entry [i] scales level i+1 from level i. Size = nMGLevels. */ + std::vector MG_LinSolIter; /*!< \brief Linear solver iteration cap per MG level (entry 0 = finest). */ bool MG_Smooth_EarlyExit{false}; /*!< \brief Enable early exit for MG smoothing iterations. */ bool MG_Smooth_Output{false}; /*!< \brief Output compact per-cycle smoothing summary. */ su2double MG_Smooth_StagnationTol{0.0}; /*!< \brief Stagnation early exit: stop if current_rms >= prev_rms * tol. 0 = disabled. */ bool MG_Implicit_Lines{false}; /*!< \brief Enable implicit-lines agglomeration from walls. */ + bool MG_Linear_Prolongation{false}; /*!< \brief Prolong the correction with a limited least-squares + gradient instead of piecewise-constant injection. */ unsigned long MG_Startup_Iter{100}; /*!< \brief Iterations per mesh during FMG startup, and the length of each level's CFL ramp. 0 = no iteration budget. */ su2double MG_Startup_Convergence{-2.0}; /*!< \brief FMG: orders of magnitude (log10) that CONV_FIELD must drop on the active level before promoting to the next finer one. Negative is a @@ -2940,6 +2944,7 @@ enum class MPI_QUANTITIES { COORDINATES , /*!< \brief Vertex coordinates communication. */ COORDINATES_OLD , /*!< \brief Old vertex coordinates communication. */ MAX_LENGTH , /*!< \brief Maximum length communication. */ + WALL_DISTANCE , /*!< \brief Wall distance and roughness of the nearest wall communication. */ GRID_VELOCITY , /*!< \brief Grid velocity communication. */ SOLUTION_EDDY , /*!< \brief Turbulent solution plus eddy viscosity communication. */ STOCH_SOURCE_LANG , /*!< \brief Stochastic source term for Langevin equations communication. */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 517b8a9579bc..13269b272eb8 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1930,6 +1930,7 @@ void CConfig::SetConfig_Options() { addULongArrayOption("OUTLIER_MITIGATION_PARAM", 4, true, outlierMitigationParam); /* DESCRIPTION: Activate The adaptive CFL number. */ addBoolOption("CFL_ADAPT", CFL_Adapt, false); + /* DESCRIPTION: Adapt the CFL on the finest-grid nonlinear residuals instead of the linear solve, with a cut on divergence. */ /* !\brief CFL_ADAPT_PARAM * DESCRIPTION: Parameters of the adaptive CFL number (factor down, factor up, CFL limit (min and max)[, acceptable linear residual][, starting iteration]). * Parameters in square brackets are optional, parameter "starting iteration" only valid with parameter "acceptable linear residual". @@ -2129,6 +2130,10 @@ void CConfig::SetConfig_Options() { /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Pave the coarse grid with advancing fronts raised from boundaries * that carry a stretched layer normal to themselves. DEFAULT: NO \ingroup Config*/ addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false); + /*!\brief MG_LINEAR_PROLONGATION\n DESCRIPTION: Prolong the multigrid correction with a limited least-squares + * gradient over the coarse control volume instead of injecting the parent value into every child. Coarse CVs on + * walls and symmetry planes keep the constant operator. DEFAULT: NO \ingroup Config*/ + addBoolOption("MG_LINEAR_PROLONGATION", MGOptions.MG_Linear_Prolongation, false); /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Max number of iterations spent on each mesh during the Full * Multigrid (FMG) startup phase. DEFAULT: 100 \ingroup Config*/ addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100); @@ -2137,8 +2142,9 @@ void CConfig::SetConfig_Options() { * DEFAULT: -2 \ingroup Config*/ addDoubleOption("MG_STARTUP_CONVERGENCE", MGOptions.MG_Startup_Convergence, -2.0); /*!\brief MG_STARTUP_STAGNATION\n DESCRIPTION: Full-MG promotion on stagnation. If the active level's residual ratio - * between successive iterations exceeds this value for MG_STARTUP_STAGNATION_ITER consecutive iterations, promote to - * the next finer level without waiting out MG_STARTUP_ITER. 0 disables it. DEFAULT: 0.99 \ingroup Config*/ + * between successive iterations stays between this value and its inverse for MG_STARTUP_STAGNATION_ITER consecutive + * iterations, promote to the next finer level without waiting out MG_STARTUP_ITER. 0 disables it. + * DEFAULT: 0.99 \ingroup Config*/ addDoubleOption("MG_STARTUP_STAGNATION", MGOptions.MG_Startup_Stagnation, 0.99); /*!\brief MG_STARTUP_STAGNATION_ITER\n DESCRIPTION: Consecutive stalled iterations required before Full-MG promotes * on stagnation. 0 disables it, as MG_STARTUP_STAGNATION= 0 does. DEFAULT: 5 \ingroup Config*/ @@ -2146,6 +2152,9 @@ void CConfig::SetConfig_Options() { /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). * If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p); + /*!\brief MG_LINEAR_SOLVER_ITER\n DESCRIPTION: Maximum linear solver iterations on coarse MG levels. Entry i applies to level i+1; + * the last value is repeated, and the cap never exceeds LINEAR_SOLVER_ITER. DEFAULT: LINEAR_SOLVER_ITER \ingroup Config*/ + addUShortListOption("MG_LINEAR_SOLVER_ITER", nMG_LinSolIter_p, MG_LinSolIter_p); /*!\par CONFIG_CATEGORY: Spatial Discretization \ingroup Config*/ /*--- Options related to the spatial discretization ---*/ @@ -2176,6 +2185,10 @@ void CConfig::SetConfig_Options() { /*!\brief MUSCL_FLOW \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/ addBoolOption("MUSCL_FLOW", MUSCL_Flow, true); + /*!\brief MUSCL_FIRST_ORDER_WALL_ENDS \n DESCRIPTION: Reconstruct the flow to first order on the edges of points + * where a no-slip wall ends on a boundary continuing its surface, or has a sharp edge such as a trailing edge, and + * skip the multigrid correction of the flow there \ingroup Config*/ + addBoolOption("MUSCL_FIRST_ORDER_WALL_ENDS", MUSCL_FirstOrderWallEnds, false); /*!\brief MUSCL_KAPPA_FLOW \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/ addDoubleOption("MUSCL_KAPPA_FLOW", MUSCL_Kappa_Flow, 0.0); /*!\brief RAMP_MUSCL \n DESCRIPTION: Enable ramping of the MUSCL scheme from 1st to 2nd order using specified method*/ @@ -3115,6 +3128,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: ParMETIS load balancing weight for edges (equiv. to neighbors) */ addLongOption("PARMETIS_EDGE_WEIGHT", ParMETIS_edgeWgt, 1); + /* DESCRIPTION: Weight ParMETIS graph edges by cell anisotropy and keep anisotropic lines on one rank */ + addBoolOption("PARMETIS_ANISOTROPY_WEIGHT", ParMETIS_anisotropyWgt, false); + /*--- options that are used in the Hybrid RANS/LES Simulations ---*/ /*!\par CONFIG_CATEGORY:Hybrid_RANSLES Options\ingroup Config*/ @@ -4882,6 +4898,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } + if (RampMUSCL && !(MUSCL_Flow || MUSCL_Turb || MUSCL_Species || MUSCL_Heat)) { + if (rank == MASTER_NODE) cout << "WARNING: RAMP_MUSCL is ignored because no MUSCL reconstruction is active." << endl; + RampMUSCL = false; + } + if(RampMUSCL && !DiscreteAdjoint){ if (RampMUSCLParam.RampMUSCLPower <= 0.0) SU2_MPI::Error("RAMP_MUSCL_POWER cannot be less than or equal to zero!", CURRENT_FUNCTION); rampMUSCLValue = 0.0; @@ -5084,6 +5105,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Kappa_2nd_AdjFlow = jst_adj_coeff[0]; Kappa_4th_AdjFlow = jst_adj_coeff[1]; + /*--- The multigrid integration carries per-level arrays of this size. ---*/ + + if (nMGLevels > MAX_MGLEVELS) { + SU2_MPI::Error("MGLEVEL is larger than the supported maximum of " + std::to_string(MAX_MGLEVELS) + ".", + CURRENT_FUNCTION); + } + /*--- Fill MG smooth vectors to size nMGLevels+1. Use parsed values (truncating or extending by repeat) or defaults if not set. ---*/ @@ -5114,6 +5142,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i for (unsigned short i = 0; i < nMGLevels; ++i) MGOptions.MG_CflScaling[i] = 0.25; } + + /*--- Linear solver iteration cap per level; level 0 always uses LINEAR_SOLVER_ITER. ---*/ + MGOptions.MG_LinSolIter.assign(nMGLevels + 1, Linear_Solver_Iter); + for (unsigned short i = 1; i <= nMGLevels && nMG_LinSolIter_p != 0; ++i) { + const auto cap = MG_LinSolIter_p[min(i - 1, nMG_LinSolIter_p - 1)]; + MGOptions.MG_LinSolIter[i] = min(cap, Linear_Solver_Iter); + } } /*--- Override MG Smooth parameters ---*/ diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index acafc2b7c025..334572192364 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -627,6 +627,10 @@ void CGeometry::GetCommCountAndType(const CConfig* config, MPI_QUANTITIES commTy COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE::DOUBLE; break; + case MPI_QUANTITIES::WALL_DISTANCE: + COUNT_PER_POINT = 2; + MPI_TYPE = COMM_TYPE::DOUBLE; + break; case MPI_QUANTITIES::NEIGHBORS: COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT; @@ -718,6 +722,10 @@ void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, MPI_QU case MPI_QUANTITIES::MAX_LENGTH: bufDSend[buf_offset] = nodes->GetMaxLength(iPoint); break; + case MPI_QUANTITIES::WALL_DISTANCE: + bufDSend[buf_offset] = nodes->GetWall_Distance(iPoint); + bufDSend[buf_offset + 1] = nodes->GetRoughnessHeight(iPoint); + break; case MPI_QUANTITIES::NEIGHBORS: bufSSend[buf_offset] = geometry->nodes->GetnNeighbor(iPoint); break; @@ -811,6 +819,10 @@ void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, MPI_QU case MPI_QUANTITIES::MAX_LENGTH: nodes->SetMaxLength(iPoint, bufDRecv[buf_offset]); break; + case MPI_QUANTITIES::WALL_DISTANCE: + nodes->SetWall_Distance(iPoint, bufDRecv[buf_offset]); + nodes->SetRoughnessHeight(iPoint, bufDRecv[buf_offset + 1]); + break; case MPI_QUANTITIES::NEIGHBORS: nodes->SetnNeighbor(iPoint, bufSRecv[buf_offset]); break; @@ -2741,6 +2753,16 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE); geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE); geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]); + + /*--- SetCoord centred a halo agglomerate on the partial child list this rank holds. Take + the owner's coordinate, otherwise coarse stencils depend on the partitioning. Deformation + runs this inside a parallel region, so only one thread may reach the exchange. ---*/ + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + geometry_container[iMesh]->InitiateComms(geometry_container[iMesh], config, MPI_QUANTITIES::COORDINATES); + geometry_container[iMesh]->CompleteComms(geometry_container[iMesh], config, MPI_QUANTITIES::COORDINATES); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS } /*--- Compute the global surface areas for all markers. ---*/ @@ -4542,6 +4564,47 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); return 2 * Vol / GeometryToolbox::Norm(3, Normal); } + +/*--- Volume-averages the wall distance and roughness of the children onto a coarse grid, then sets the + * nearest-neighbor distance of its viscous wall vertices. ---*/ +void RestrictWallDistance(const CGeometry* geo_fine, CGeometry* geo_coarse, const CConfig* config) { + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geo_coarse->GetnPointDomain(); iPoint++) { + su2double dist = 0.0, roughness = 0.0, vol = 0.0; + for (auto iChild = 0u; iChild < geo_coarse->nodes->GetnChildren_CV(iPoint); iChild++) { + const auto jPoint = geo_coarse->nodes->GetChildren_CV(iPoint, iChild); + const su2double volChild = geo_fine->nodes->GetVolume(jPoint); + dist += geo_fine->nodes->GetWall_Distance(jPoint) * volChild; + roughness += geo_fine->nodes->GetRoughnessHeight(jPoint) * volChild; + vol += volChild; + } + geo_coarse->nodes->SetWall_Distance(iPoint, (vol > 0.0) ? su2double(dist / vol) : su2double(0.0)); + geo_coarse->nodes->SetRoughnessHeight(iPoint, (vol > 0.0) ? su2double(roughness / vol) : su2double(0.0)); + } + END_SU2_OMP_FOR + + /*--- A halo agglomerate only holds the children on this rank, so take the owner's values. ---*/ + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + geo_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::WALL_DISTANCE); + geo_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::WALL_DISTANCE); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) { + const auto viscous = config->GetViscous_Wall(iMarker); + + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->nVertex[iMarker], omp_get_num_threads())) + for (auto iVertex = 0u; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + const auto iPoint = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + const su2double dist = (viscous && geo_coarse->nodes->GetDomain(iPoint)) + ? NearestNeighborDistance(geo_coarse, config, iPoint) + : geo_coarse->nodes->GetWall_Distance(iPoint); + geo_coarse->vertex[iMarker][iVertex]->SetNearestNeighborDistance(dist); + } + END_SU2_OMP_FOR + } +} } // namespace void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container, @@ -4642,6 +4705,13 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo } END_SU2_OMP_FOR } + + /*--- The Full-MG startup solves the turbulence model on coarse grids. ---*/ + + for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { + RestrictWallDistance(geometry_container[iZone][iInst][iMesh - 1], geometry_container[iZone][iInst][iMesh], + config); + } } } } diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index b8caa00e9c94..e1ee13ecbca3 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -69,6 +69,39 @@ vector MarkerSetClasses(unsigned long nEntity, vectorGetnDim(); + if ((nDim != 3) || (config->GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE)) return false; + + auto unitNormal = [&](unsigned short jMarker, su2double* n) { + const long iVertex = grid->nodes->GetVertex(iPoint, jMarker); + if (iVertex < 0) return false; + grid->vertex[jMarker][iVertex]->GetNormal(n); + const su2double nrm = GeometryToolbox::Norm(nDim, n); + if (nrm <= 0.0) return false; + for (unsigned short d = 0; d < nDim; ++d) n[d] /= nrm; + return true; + }; + + su2double nSym[3] = {0.0}, nOther[3] = {0.0}; + if (!unitNormal(iMarker, nSym)) return false; + + bool hasOther = false; + for (auto jMarker = 0u; jMarker < grid->GetnMarker(); jMarker++) { + const auto bc = config->GetMarker_All_KindBC(jMarker); + if ((bc == SEND_RECEIVE) || (bc == SYMMETRY_PLANE) || (grid->nodes->GetVertex(iPoint, jMarker) < 0)) continue; + if ((bc == PERIODIC_BOUNDARY) || !unitNormal(jMarker, nOther)) return false; + if (fabs(GeometryToolbox::DotProduct(nDim, nSym, nOther)) >= sqrt(0.5)) return false; + hasOther = true; + } + return hasOther; +} + +bool HasElements(const CGeometry* grid); + } // namespace CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, unsigned short iMesh) : CGeometry() { @@ -387,6 +420,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un vector inCV(fine_grid->GetnPoint(), 0); vector isCandidate(fine_grid->GetnPoint(), 0); vector members, candidates; + vector> wholeCells; members.reserve(maxAgglomSize); /*--- A local frame at the seed: up to nDim of its incident edges, as mutually orthogonal as @@ -434,6 +468,56 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un addMember(iPoint); + /*--- A seed standing in a free quadrilateral or hexahedron takes that element whole, the one + * pressing hardest against what is already agglomerated, so leftover structured mesh is + * tiled by whole cells. A grid without elements, a coarse one, offers the 4-cycles of its + * connectivity as quadrilaterals. ---*/ + wholeCells.clear(); + if (HasElements(fine_grid)) { + for (auto iElem = 0u; iElem < fine_grid->nodes->GetnElem(iPoint); ++iElem) { + const auto* elem = fine_grid->elem[fine_grid->nodes->GetElem(iPoint, iElem)]; + if ((elem->GetVTK_Type() != QUADRILATERAL) && (elem->GetVTK_Type() != HEXAHEDRON)) continue; + wholeCells.emplace_back(); + for (auto iNode = 0u; iNode < elem->GetnNodes(); ++iNode) + if (elem->GetNode(iNode) != static_cast(iPoint)) + wholeCells.back().push_back(elem->GetNode(iNode)); + } + } else if (nDim == 2) { + const auto neighbours = fine_grid->nodes->GetPoints(iPoint); + auto touching = [&](unsigned long a, unsigned long b) { + const auto around = fine_grid->nodes->GetPoints(a); + return std::find(around.begin(), around.end(), b) != around.end(); + }; + for (auto a : neighbours) + for (auto b : neighbours) { + if ((a >= b) || touching(a, b)) continue; + for (auto m : fine_grid->nodes->GetPoints(a)) + if ((m != iPoint) && touching(m, b) && !touching(m, iPoint)) wholeCells.push_back({a, m, b}); + } + } + + auto bestCell = NO_CV; + unsigned long bestContact = 0, bestKey = NO_CV; + for (auto iCell = 0ul; iCell < wholeCells.size(); ++iCell) { + bool free = true; + unsigned long contact = 0, key = fine_grid->nodes->GetGlobalIndex(iPoint); + for (auto kPoint : wholeCells[iCell]) { + key = std::min(key, fine_grid->nodes->GetGlobalIndex(kPoint)); + free = free && !fine_grid->nodes->GetAgglomerate(kPoint) && fine_grid->nodes->GetDomain(kPoint) && + GeometricalCheck(kPoint, fine_grid, config); + if (!free) break; + for (auto mPoint : fine_grid->nodes->GetPoints(kPoint)) contact += fine_grid->nodes->GetAgglomerate(mPoint); + } + if (!free) continue; + if ((bestCell == NO_CV) || (contact > bestContact) || ((contact == bestContact) && (key < bestKey))) { + bestCell = iCell; + bestContact = contact; + bestKey = key; + } + } + if (bestCell != NO_CV) + for (auto kPoint : wholeCells[bestCell]) addMember(kPoint); + edgeDir.clear(); edgeLen.clear(); su2double seedShortest = std::numeric_limits::max(); @@ -494,7 +578,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } unsigned long best = NO_CV; - unsigned short best_shared = 0; + unsigned short best_shared = 0, best_flush = 0; su2double best_dist = std::numeric_limits::max(); for (auto CVPoint : candidates) { @@ -502,6 +586,27 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un unsigned short shared = 0; for (auto jPoint : fine_grid->nodes->GetPoints(CVPoint)) shared += inCV[jPoint]; + /*--- In 2D a CV of three or more nodes only takes nodes touching two members, so it cannot + * grow into a chain. ---*/ + if ((nDim == 2) && (nChildren >= 3) && (shared < 2)) continue; + + /*--- A pair does not grow into a line of three, the third node then seeds a pair of its own, + * so the cells of a row keep the phase of the cells beside them. ---*/ + if ((nDim == 2) && (nChildren == 2) && (shared == 1)) { + const auto around = fine_grid->nodes->GetPoints(CVPoint); + const bool firstTouches = std::find(around.begin(), around.end(), members[0]) != around.end(); + const auto tip = firstTouches ? members[0] : members[1], tail = firstTouches ? members[1] : members[0]; + su2double along[MAXNDIM] = {0.0}, across[MAXNDIM] = {0.0}; + for (unsigned short d = 0; d < nDim; ++d) { + along[d] = fine_grid->nodes->GetCoord(CVPoint, d) - fine_grid->nodes->GetCoord(tip, d); + across[d] = fine_grid->nodes->GetCoord(tip, d) - fine_grid->nodes->GetCoord(tail, d); + } + const su2double lenAlong = GeometryToolbox::Norm(nDim, along); + const su2double lenAcross = GeometryToolbox::Norm(nDim, across); + if ((lenAlong > 0.0) && (lenAcross > 0.0) && + (GeometryToolbox::DotProduct(nDim, along, across) > 0.5 * lenAlong * lenAcross)) + continue; + } /*--- Distance to the centroid in cells: the offset is scaled by the seed edge pointing * most nearly along it. ---*/ @@ -522,9 +627,24 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un dist = r * r; } - if ((shared > best_shared) || ((shared == best_shared) && (dist < best_dist))) { + /*--- Between equals, the node lying flush against the same neighbouring CVs as the members + * keeps the cells of a structured mesh in line with those already made. ---*/ + unsigned short flush = 0; + for (auto jPoint : fine_grid->nodes->GetPoints(CVPoint)) { + if (!fine_grid->nodes->GetAgglomerate(jPoint) || inCV[jPoint]) continue; + bool lies = false; + for (auto kPoint : members) + for (auto lPoint : fine_grid->nodes->GetPoints(kPoint)) + lies = lies || ((lPoint != jPoint) && fine_grid->nodes->GetAgglomerate(lPoint) && !inCV[lPoint] && + (fine_grid->nodes->GetParent_CV(lPoint) == fine_grid->nodes->GetParent_CV(jPoint))); + flush += lies; + } + + if ((shared > best_shared) || ((shared == best_shared) && (flush > best_flush)) || + ((shared == best_shared) && (flush == best_flush) && (dist < best_dist))) { best = CVPoint; best_shared = shared; + best_flush = flush; best_dist = dist; } } @@ -638,18 +758,42 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un nodes->SetnChildren_CV(fromCV, 0); }; - /*--- Neighbour with the fewest children that may take this CV, NO_CV if there is none. Neighbours - * an earlier merge in the same pass emptied are not revived. ---*/ - auto smallestPartner = [&](unsigned long iCV, unsigned short maxPartnerSize, bool keepOffPartition) { + /*--- Dual-face area and number of fine edges between the children of two coarse CVs. ---*/ + auto sharedFaces = [&](unsigned long iCV, unsigned long jCV) { + su2double area = 0.0; + unsigned short nEdges = 0; + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCV); iChildren++) { + const auto iFinePoint = nodes->GetChildren_CV(iCV, iChildren); + for (auto iNeigh = 0u; iNeigh < fine_grid->nodes->GetnPoint(iFinePoint); iNeigh++) { + if (fine_grid->nodes->GetParent_CV(fine_grid->nodes->GetPoint(iFinePoint, iNeigh)) != jCV) continue; + area += GeometryToolbox::Norm(nDim, fine_grid->edges->GetNormal(fine_grid->nodes->GetEdge(iFinePoint, iNeigh))); + nEdges++; + } + } + return std::make_pair(area, nEdges); + }; + + /*--- Neighbour sharing the most dual-face area that may take this CV, the one with fewer children on + * a tie, NO_CV if there is none. A compact merge must join the two through at least as many fine + * edges as the smaller one has children. Neighbours an earlier merge in the same pass emptied are + * not revived. ---*/ + auto bestPartner = [&](unsigned long iCV, unsigned short maxPartnerSize, bool keepOffPartition, bool compact) { auto best = NO_CV; auto best_nChildren = std::numeric_limits::max(); + su2double best_area = -1.0; for (auto jCV : nodes->GetPoints(iCV)) { const auto nChildren = nodes->GetnChildren_CV(jCV); - if ((nChildren == 0) || (nChildren > maxPartnerSize) || (nChildren >= best_nChildren)) continue; + if ((nChildren == 0) || (nChildren > maxPartnerSize)) continue; if (keepOffPartition && touchesPartition[jCV]) continue; if (!mayMerge(iCV, jCV)) continue; - best_nChildren = nChildren; - best = jCV; + const auto shared = sharedFaces(iCV, jCV); + if (compact && (shared.second < min(nChildren, nodes->GetnChildren_CV(iCV)))) continue; + if ((shared.first > best_area * (1.0 + 1e-6)) || + ((shared.first >= best_area * (1.0 - 1e-6)) && (nChildren < best_nChildren))) { + best_nChildren = nChildren; + best_area = shared.first; + best = jCV; + } } return best; }; @@ -686,8 +830,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un mergeInto(iCoarsePoint, iTarget); } - /*--- A CV left holding a single child joins the neighbour with the fewest, even where that takes - the neighbour one past maxAgglomSize. ---*/ + /*--- A CV left holding a single child joins the neighbour it shares the most face with, even where + that takes the neighbour one past maxAgglomSize. ---*/ for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { if (nodes->GetnChildren_CV(iCoarsePoint) != 1) continue; @@ -695,12 +839,12 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- Already handled above, or truly islanded. ---*/ if (nodes->GetnPoint(iCoarsePoint) <= 1) continue; - const auto iTarget = smallestPartner(iCoarsePoint, std::numeric_limits::max(), false); + const auto iTarget = bestPartner(iCoarsePoint, std::numeric_limits::max(), false, false); if (iTarget != NO_CV) mergeInto(iCoarsePoint, iTarget); } - /*--- Merge a paved CV of at most SMALL_STACK_CV children into an equally small neighbour, which - the two passes above do not reach. Two of them can never exceed maxAgglomSize together. ---*/ + /*--- Merge a paved CV of at most SMALL_STACK_CV children into an equally small neighbour it joins + compactly, which the two passes above do not reach. Two of them can never exceed maxAgglomSize together. ---*/ constexpr unsigned short SMALL_STACK_CV = 2; /*!< \brief The block size a narrow stack emits between flushes; see BlockFor. */ @@ -712,10 +856,74 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- Interior stack fragments only, a boundary row keeps its footprint. ---*/ if (cvOnBoundary[iCoarsePoint]) continue; - const auto iTarget = smallestPartner(iCoarsePoint, SMALL_STACK_CV, true); + const auto iTarget = bestPartner(iCoarsePoint, SMALL_STACK_CV, true, true); if (iTarget != NO_CV) mergeInto(iCoarsePoint, iTarget); } + /*--- In 2D, a CV of three nodes bent round a corner takes the fourth node of its quadrilateral from + * a neighbouring three-node CV that can spare it, which turns staggered seams into rectangles. ---*/ + if (nDim == 2) { + auto adjacent = [&](unsigned long iFinePoint, unsigned long jFinePoint) { + const auto neighbours = fine_grid->nodes->GetPoints(iFinePoint); + return std::find(neighbours.begin(), neighbours.end(), jFinePoint) != neighbours.end(); + }; + + for (auto iPass = 0u; iPass < 4; iPass++) { + bool changed = false; + + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { + if (nodes->GetnChildren_CV(iCoarsePoint) != 3) continue; + if (mustStayAlone[iCoarsePoint] || isStackBase(iCoarsePoint) || touchesPartition[iCoarsePoint]) continue; + if ((iCoarsePoint >= firstLineCV) && (iCoarsePoint < endLineCV)) continue; + + /*--- The middle node of the bend touches the other two, which do not touch each other. ---*/ + unsigned long child[3]; + for (auto k = 0u; k < 3; k++) child[k] = nodes->GetChildren_CV(iCoarsePoint, k); + long corner = -1; + for (auto k = 0u; k < 3; k++) { + const auto a = child[(k + 1) % 3], b = child[(k + 2) % 3]; + if (adjacent(child[k], a) && adjacent(child[k], b) && !adjacent(a, b)) corner = k; + } + if (corner < 0) continue; + const auto a = child[(corner + 1) % 3], b = child[(corner + 2) % 3]; + + /*--- The node closing the quadrilateral touches both ends of the bend, not its middle. ---*/ + auto donor = NO_CV; + unsigned long fourth = NO_CV; + for (auto m : fine_grid->nodes->GetPoints(a)) { + const auto jCoarsePoint = fine_grid->nodes->GetParent_CV(m); + if ((jCoarsePoint == NO_CV) || (jCoarsePoint >= nPointDomain) || (jCoarsePoint == iCoarsePoint)) continue; + if (!adjacent(m, b) || adjacent(m, child[corner])) continue; + if (nodes->GetnChildren_CV(jCoarsePoint) != 3) continue; + if (mustStayAlone[jCoarsePoint] || isStackBase(jCoarsePoint) || touchesPartition[jCoarsePoint]) continue; + if ((jCoarsePoint >= firstLineCV) && (jCoarsePoint < endLineCV)) continue; + if (cvMarkerClass[jCoarsePoint] != cvMarkerClass[iCoarsePoint]) continue; + + /*--- What the donor keeps must stay connected. ---*/ + vector kept; + for (auto k = 0u; k < 3; k++) + if (nodes->GetChildren_CV(jCoarsePoint, k) != m) kept.push_back(nodes->GetChildren_CV(jCoarsePoint, k)); + if (!adjacent(kept[0], kept[1])) continue; + + donor = jCoarsePoint; + fourth = m; + break; + } + if (donor == NO_CV) continue; + + const auto nDonor = nodes->GetnChildren_CV(donor); + for (auto k = 0u; k < nDonor; k++) + if (nodes->GetChildren_CV(donor, k) == fourth) + nodes->SetChildren_CV(donor, k, nodes->GetChildren_CV(donor, nDonor - 1)); + nodes->SetnChildren_CV(donor, nDonor - 1); + giveChild(fourth, iCoarsePoint); + changed = true; + } + + if (!changed) break; + } + } + /*--- Compact the coarse numbering, squeezing out the indices the repair passes emptied. The children lists, indirect-agglomeration flags and owned parent indices are remapped. ---*/ @@ -996,6 +1204,7 @@ vector CMultiGridGeometry::FindMixedBoundaryNodes(const CGeometry* fine_gr const auto bc = static_cast(config->GetMarker_All_KindBC(iMarker)); for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); + if (AbsorbedSymmetry(fine_grid, config, iPoint, iMarker)) continue; if (firstBC[iPoint] < 0) firstBC[iPoint] = bc; else if (firstBC[iPoint] != bc) @@ -1616,7 +1825,7 @@ CMultiGridGeometry::CFrontSeeds CMultiGridGeometry::SeedFrontNodes(const CGeomet /*--- Paving order: a no-slip wall first, a slip wall next, everything else last. ---*/ auto tierOfBC = [](unsigned short bc) -> char { if ((bc == HEAT_FLUX) || (bc == ISOTHERMAL) || (bc == CHT_WALL_INTERFACE) || (bc == SMOLUCHOWSKI_MAXWELL)) return 0; - return (bc == EULER_WALL) ? 1 : 2; + return ((bc == EULER_WALL) || (bc == SYMMETRY_PLANE)) ? 1 : 2; }; for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { @@ -1648,6 +1857,7 @@ CMultiGridGeometry::CFrontSeeds CMultiGridGeometry::SeedFrontNodes(const CGeomet std::array n0{}; for (unsigned short d = 0; d < nDim; ++d) n0[d] = Normal[d]; + seeds.node.push_back(iPoint); seeds.normal.push_back(n0); seeds.strength.push_back(strength); @@ -1894,9 +2104,15 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co const auto nPointFine = fine_grid->GetnPoint(); constexpr auto NO_COLUMN = std::numeric_limits::max(); const short int maxAgglomSize = MaxAgglomSize(nDim); + /*--- Largest CV a parity leftover may form by joining its neighbour. ---*/ + const size_t maxOversize = 3 * maxAgglomSize / 2; /*--- SeedFrontNodes returns the seeds ordered by anisotropy. ---*/ - const auto seeds = SeedFrontNodes(fine_grid, config); + auto seeds = SeedFrontNodes(fine_grid, config); + + /*--- A node where two conditions meet stays a CV of its own, and the column above it is paved + * before any other, so no column from either side reaches over to the other. ---*/ + for (auto si = 0ul; si < seeds.node.size(); ++si) seeds.tier[si] = mixedBC[seeds.node[si]] ? 0 : seeds.tier[si] + 1; vector walkOf; const auto patches = BuildFrontPatches(seeds, fine_grid, config, mixedBC, walkOf); @@ -2016,7 +2232,7 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co /*--- One tier is paved out everywhere before the next starts, so an inherited column belongs to * the tier that is running. ---*/ - constexpr char N_TIER = 3; + constexpr char N_TIER = 4; vector active; for (char tier = 0; tier < N_TIER; ++tier) { /*--- The columns of this tier still advancing, in seed order. A column that stops drops out at @@ -2129,7 +2345,7 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co if (admissible(best)) { candidates.push_back(best); - stepDir.push_back(bestStep); + stepDir.push_back(marching); } else if (!fine_grid->nodes->GetDomain(best) && (haloSlot[best] >= 0) && (haloClaim[best] == NO_COLUMN)) { haloWanted.push_back(best); @@ -2145,10 +2361,70 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co std::sort(distinct.begin(), distinct.end()); distinct.erase(std::unique(distinct.begin(), distinct.end()), distinct.end()); if (distinct.size() != width) return false; + /*--- The layer must stay connected, or its blocks split into strands. ---*/ + vector reached(1, distinct[0]); + for (size_t k = 0; k < reached.size(); ++k) + for (auto jPoint : fine_grid->nodes->GetPoints(reached[k])) + if (std::binary_search(distinct.begin(), distinct.end(), jPoint) && + (std::find(reached.begin(), reached.end(), jPoint) == reached.end())) + reached.push_back(jPoint); + return reached.size() == width; + }; + + /*--- A column follows the stiff direction of the mesh. A step whose edges couple far more + * weakly than the stiffest edge at their node runs along a thin layer, not across it, + * and is left to a column crossing that layer. A column rising from a no-slip wall + * stops only at a much thinner layer, so moderate stretching does not cut it short. ---*/ + auto stiffStep = [&](const vector& set) { + const passivedouble STIFF_FRACTION = (tierOf[iColumn] <= 1) ? 0.01 : 0.1; + su2double ratio = 0.0; + for (auto iPoint : layer[iColumn]) { + su2double wStep = 0.0, wMax = 0.0; + for (auto iNeigh = 0u; iNeigh < fine_grid->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = fine_grid->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = fine_grid->nodes->GetEdge(iPoint, iNeigh); + const su2double w = GeometryToolbox::Norm(nDim, fine_grid->edges->GetNormal(iEdge)) * + (1.0 / fine_grid->nodes->GetVolume(iPoint) + 1.0 / fine_grid->nodes->GetVolume(jPoint)); + wMax = std::max(wMax, w); + if (std::find(set.begin(), set.end(), jPoint) != set.end()) wStep = std::max(wStep, w); + } + ratio += (wMax > 0.0) ? wStep / wMax : 1.0; + } + return ratio >= STIFF_FRACTION * su2double(layer[iColumn].size()); + }; + + /*--- Whether every node of the step reaches a physical boundary within a few steps on the + * way it travels. A band two coarse cells thin is taken, the domain pass tiles it poorly. ---*/ + auto nearBoundary = [&](const vector& set, const vector>& dirs) { + constexpr unsigned short THIN_BAND = 4; + for (size_t k = 0; k < set.size(); ++k) { + auto iPoint = set[k]; + bool reached = false; + for (unsigned short iStep = 0; (iStep < THIN_BAND) && !reached; ++iStep) { + auto next = NO_COLUMN; + su2double bestAlign = 0.0; + for (auto jPoint : fine_grid->nodes->GetPoints(iPoint)) { + su2double step[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(jPoint), fine_grid->nodes->GetCoord(iPoint), + step); + const su2double len = GeometryToolbox::Norm(nDim, step); + if (len <= 0.0) continue; + const su2double align = GeometryToolbox::DotProduct(nDim, step, dirs[k].data()) / len; + if (align > bestAlign) { + bestAlign = align; + next = jPoint; + } + } + if (next == NO_COLUMN) return false; + reached = onPhysBoundary[next]; + iPoint = next; + } + if (!reached) return false; + } return true; }; - if (complete(candidates)) { + if (complete(candidates) && (stiffStep(candidates) || nearBoundary(candidates, stepDir))) { for (auto jPoint : candidates) { columnOf[jPoint] = iColumn; layerOf[jPoint] = iRound; @@ -2309,6 +2585,28 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co } } + /*--- Last-layer nodes whose next step, the way they travel, lands on free mesh. ---*/ + vector openAhead(nPointFine, 0); + for (auto iColumn = 0ul; iColumn < layer.size(); ++iColumn) { + for (size_t k = 0; k < layer[iColumn].size(); ++k) { + const auto iPoint = layer[iColumn][k]; + auto next = NO_COLUMN; + su2double bestAlign = 0.0; + for (auto jPoint : fine_grid->nodes->GetPoints(iPoint)) { + su2double step[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, fine_grid->nodes->GetCoord(jPoint), fine_grid->nodes->GetCoord(iPoint), step); + const su2double len = GeometryToolbox::Norm(nDim, step); + if (len <= 0.0) continue; + const su2double align = GeometryToolbox::DotProduct(nDim, step, dirOf[iColumn][k].data()) / len; + if (align > bestAlign) { + bestAlign = align; + next = jPoint; + } + } + openAhead[iPoint] = (next != NO_COLUMN) && admissible(next); + } + } + /*--- Gather each column by layer. ---*/ const auto nColumn = layer.size(); vector>> byLayer(nColumn); @@ -2319,6 +2617,75 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co byLayer[iColumn][layerOf[iPoint]].push_back(iPoint); } + /*--- A boundary run with an odd number of seeds, and every node where two boundaries meet, + * leaves one column a single node wide. Two such columns side by side are emitted together, + * so the domain gets square cells instead of a file of narrow ones. ---*/ + vector partnerOf(nColumn, -1); + { + auto isNarrow = [&](unsigned long iColumn) { + return isSeeded[iColumn] && (byLayer[iColumn].size() > 1) && (byLayer[iColumn][0].size() == 1); + }; + + /*--- Both columns must stand shoulder to shoulder over every layer they share, otherwise the + * cells the merge makes are not compact. ---*/ + auto aligned = [&](unsigned long iColumn, unsigned long jColumn) { + const auto nCommon = std::min(byLayer[iColumn].size(), byLayer[jColumn].size()); + for (auto iLayer = 1ul; iLayer < nCommon; ++iLayer) { + const auto iPoint = byLayer[iColumn][iLayer].front(); + const auto& jLayer = byLayer[jColumn][iLayer]; + bool touch = false; + for (auto kPoint : fine_grid->nodes->GetPoints(iPoint)) + touch = touch || (std::find(jLayer.begin(), jLayer.end(), kPoint) != jLayer.end()); + if (!touch) return false; + } + return nCommon > 1; + }; + + /*--- Taken in global index order, so the pairing does not depend on the partitioning. ---*/ + vector narrow; + for (auto iColumn = 0ul; iColumn < nColumn; ++iColumn) + if (isNarrow(iColumn)) narrow.push_back(iColumn); + std::sort(narrow.begin(), narrow.end(), [&](unsigned long a, unsigned long b) { + return fine_grid->nodes->GetGlobalIndex(byLayer[a][0].front()) < + fine_grid->nodes->GetGlobalIndex(byLayer[b][0].front()); + }); + + /*--- A column a narrow one may lean on when no narrow neighbour is left: seeded, grown, and + * narrow enough that the merged block stays within the oversize limit. ---*/ + auto isWide = [&](unsigned long iColumn) { + return isSeeded[iColumn] && (byLayer[iColumn].size() > 1) && (byLayer[iColumn][0].size() > 1) && + (2 * (byLayer[iColumn][0].size() + 1) <= maxOversize); + }; + + /*--- Narrow neighbours pair first, the narrow columns left over then join a wide one. ---*/ + for (const bool toWide : {false, true}) { + for (auto iColumn : narrow) { + if (partnerOf[iColumn] >= 0) continue; + auto best = NO_COLUMN; + auto bestKey = std::numeric_limits::max(); + + for (auto jPoint : fine_grid->nodes->GetPoints(byLayer[iColumn][1].front())) { + const auto jColumn = columnOf[jPoint]; + if ((jColumn == NO_COLUMN) || (jColumn == iColumn) || (partnerOf[jColumn] >= 0)) continue; + if (toWide ? !isWide(jColumn) : !isNarrow(jColumn)) continue; + /*--- The column above a corner, tier 0, may pair with either side. ---*/ + if ((tierOf[jColumn] != tierOf[iColumn]) && (tierOf[jColumn] != 0) && (tierOf[iColumn] != 0)) continue; + if (!aligned(iColumn, jColumn)) continue; + + const auto key = fine_grid->nodes->GetGlobalIndex(byLayer[jColumn][0].front()); + if (key < bestKey) { + bestKey = key; + best = jColumn; + } + } + + if (best == NO_COLUMN) continue; + partnerOf[iColumn] = static_cast(best); + partnerOf[best] = static_cast(iColumn); + } + } + } + auto emitGroup = [&](const vector& group) { nodes->SetChildren_CV(Index_CoarseCV, group); for (auto iPoint : group) { @@ -2331,11 +2698,12 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co auto minDepth = std::numeric_limits::max(), maxDepth = 0ul; vector group; + vector> paired; /*--- Hand a set of nodes out as coarse CVs, each connected and within the size limit. ---*/ vector inSet(nPointFine, 0); vector chunk, stack; - auto emitConnected = [&](const vector& set) { + auto emitConnected = [&](const vector& set, size_t maxSize) { for (auto iPoint : set) inSet[iPoint] = 1; for (auto iSeed : set) { @@ -2344,13 +2712,13 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co stack.assign(1, iSeed); inSet[iSeed] = 0; - while (!stack.empty() && (chunk.size() < static_cast(maxAgglomSize))) { + while (!stack.empty() && (chunk.size() < maxSize)) { const auto iPoint = stack.back(); stack.pop_back(); chunk.push_back(iPoint); for (auto jPoint : fine_grid->nodes->GetPoints(iPoint)) { if (!inSet[jPoint]) continue; - if (chunk.size() + stack.size() >= static_cast(maxAgglomSize)) break; + if (chunk.size() + stack.size() >= maxSize) break; inSet[jPoint] = 0; stack.push_back(jPoint); } @@ -2376,7 +2744,8 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co auto iLayer = 0ul; if (isSeeded[iColumn]) { - /*--- The boundary row is a coarse CV of its own, which fixes the footprint above it. ---*/ + /*--- The boundary row is a coarse CV of its own, which fixes the footprint above it. A + * paired column keeps its own row, which may hold a boundary condition of its own. ---*/ const auto baseCV = emitGroup(layers[0]); ct[P_CVS]++; ct[P_COVERED] += layers[0].size(); @@ -2384,22 +2753,58 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co iLayer = 1; } + /*--- Of a pair, the column of lower index emits what stands above both boundary rows. ---*/ + if ((partnerOf[iColumn] >= 0) && (partnerOf[iColumn] < static_cast(iColumn))) continue; + + const auto* emitted = &layers; + /*--- The block depth follows the widest of the columns, as it would for each on its own. ---*/ + auto blockWidth = layers.back().size(); + auto maxSize = static_cast(maxAgglomSize); + if (partnerOf[iColumn] > static_cast(iColumn)) { + const auto& other = byLayer[partnerOf[iColumn]]; + paired = layers; + for (auto k = 1ul; k < other.size(); ++k) { + if (k < paired.size()) + paired[k].insert(paired[k].end(), other[k].begin(), other[k].end()); + else + paired.push_back(other[k]); + } + emitted = &paired; + blockWidth = std::max(blockWidth, other.back().size()); + /*--- A narrow column leaning on a wide one is a parity strip, its cells may be oversized. ---*/ + if (blockWidth > 1) maxSize = maxOversize; + } + const auto block = BlockFor(maxAgglomSize, blockWidth); + /*--- Above it, consecutive layers are blocked so the coarse cell coarsens by the same ratio * along the column as the patch does across it. ---*/ - while (iLayer < layers.size()) { + while (iLayer < emitted->size()) { group.clear(); - const auto block = BlockFor(maxAgglomSize, layers[iLayer].size()); - for (auto k = 0ul; (k < block) && (iLayer < layers.size()); ++k) { - if (!group.empty() && (group.size() + layers[iLayer].size() > static_cast(maxAgglomSize))) break; - group.insert(group.end(), layers[iLayer].begin(), layers[iLayer].end()); + auto groupSize = maxSize; + for (auto k = 0ul; (k < block) && (iLayer < emitted->size()); ++k) { + if (!group.empty() && (group.size() + (*emitted)[iLayer].size() > maxSize)) break; + group.insert(group.end(), (*emitted)[iLayer].begin(), (*emitted)[iLayer].end()); iLayer++; } /*--- A single layer wider than the limit still has to go somewhere. ---*/ if (group.empty()) { - group = layers[iLayer]; + group = (*emitted)[iLayer]; iLayer++; } - emitConnected(group); + /*--- An odd last layer facing free mesh is left to the domain pass, which can pair it with + * that mesh. Otherwise it joins the block below instead of standing one layer thick. ---*/ + if ((block > 1) && (iLayer + 1 == emitted->size())) { + bool open = true; + for (auto iPoint : emitted->back()) open = open && openAhead[iPoint]; + if (open) { + iLayer++; + } else if (group.size() + emitted->back().size() <= static_cast(maxOversize)) { + group.insert(group.end(), emitted->back().begin(), emitted->back().end()); + groupSize = group.size(); + iLayer++; + } + } + emitConnected(group, groupSize); } } if (minDepth == std::numeric_limits::max()) minDepth = 0; diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 6928c2d922b5..ef7de89ec30f 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -2426,9 +2426,9 @@ void CPhysicalGeometry::LoadSurfaceElements(CConfig* config, CGeometry* geometry unsigned long Local_Nodes[N_POINTS_HEXAHEDRON]; - vector > Line_List; - vector > BoundTria_List; - vector > BoundQuad_List; + vector> Line_List; + vector> BoundTria_List; + vector> BoundQuad_List; vector Marker_Local; @@ -2967,14 +2967,14 @@ void CPhysicalGeometry::SetSendReceive(const CConfig* config) { unsigned short iNode, jNode; vector::iterator it; - vector > + vector> SendTransfLocal; /*!< \brief Vector to store the type of transformation for this send point. */ - vector > + vector> ReceivedTransfLocal; /*!< \brief Vector to store the type of transformation for this received point. */ - vector > SendDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the point - index of the node that must me sended. */ - vector > ReceivedDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the - point index of the node that must me sended. */ + vector> SendDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the point + index of the node that must me sended. */ + vector> ReceivedDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the + point index of the node that must me sended. */ unordered_map::const_iterator MI; @@ -4426,7 +4426,7 @@ void CPhysicalGeometry::SetPositive_ZArea(CConfig* config) { } void CPhysicalGeometry::SetPoint_Connectivity() { - vector > points(nPoint); + vector> points(nPoint); SU2_OMP_PARALLEL { unsigned short Node_Neighbor, iNode, iNeighbor; @@ -4434,7 +4434,7 @@ void CPhysicalGeometry::SetPoint_Connectivity() { /*--- Loop over all the elements ---*/ BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - vector > elems(nPoint); + vector> elems(nPoint); for (iElem = 0; iElem < nElem; iElem++) { /*--- Loop over all the nodes of an element ---*/ @@ -7361,6 +7361,277 @@ void CPhysicalGeometry::SetBoundTecPlot(char mesh_filename[MAX_STRING_SIZE], boo Tecplot_File.close(); } +#if defined(HAVE_MPI) && defined(HAVE_PARMETIS) +vector CPhysicalGeometry::FetchFromOwners(const vector& wanted, + const vector& localValue) const { + CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); + const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); + + vector nSend(size, 0), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (const auto iGlobal : wanted) nSend[pointPartitioner.GetRankContainingIndex(iGlobal)]++; + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + + vector sendIndex(sendDisp[size]), slotOf(wanted.size()), recvIndex(recvDisp[size]); + vector next(sendDisp.begin(), sendDisp.end() - 1); + for (auto k = 0ul; k < wanted.size(); k++) { + const auto slot = next[pointPartitioner.GetRankContainingIndex(wanted[k])]++; + sendIndex[slot] = wanted[k]; + slotOf[k] = slot; + } + SU2_MPI::Alltoallv(sendIndex.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvIndex.data(), nRecv.data(), + recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + vector reply(recvIndex.size()), answer(sendIndex.size()), result(wanted.size()); + for (auto k = 0ul; k < recvIndex.size(); k++) reply[k] = localValue[recvIndex[k] - firstIndex]; + SU2_MPI::Alltoallv(reply.data(), nRecv.data(), recvDisp.data(), MPI_UNSIGNED_LONG, answer.data(), nSend.data(), + sendDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + for (auto k = 0ul; k < wanted.size(); k++) result[k] = answer[slotOf[k]]; + return result; +} + +vector CPhysicalGeometry::LinkAnisotropicLines(const vector& adjwgt, + vector& lineEdge) const { + /*--- An edge at least this many times shorter than its element can join its points into a line. ---*/ + constexpr idx_t LINE_WEIGHT = 4; + + const unsigned long firstIndex = CLinearPartitioner(Global_nPointDomain, 0).GetFirstIndexOnRank(rank); + const vector neighbor(adjacency.begin(), adjacency.end()); + + vector strongest(nPoint, 0); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; k++) + strongest[iPoint] = max(strongest[iPoint], static_cast(adjwgt[k])); + const auto neighborStrongest = FetchFromOwners(neighbor, strongest); + + /*--- An edge joins a line if it is the strongest at one end, or nearly so at both. Equal or + * slowly varying spacing then keeps a line whole instead of splitting it where weights tie. ---*/ + lineEdge.assign(adjacency.size(), 0); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; k++) { + const auto w = static_cast(adjwgt[k]); + if (w < static_cast(LINE_WEIGHT)) continue; + const bool strongestAtEnd = (w == strongest[iPoint]) || (w == neighborStrongest[k]); + const bool strongAtBoth = (2 * w >= strongest[iPoint]) && (2 * w >= neighborStrongest[k]); + lineEdge[k] = strongestAtEnd || strongAtBoth; + } + } + + /*--- Each line is named by its lowest global index, found by propagating labels along the line + * edges and pointer jumping. ---*/ + vector root(nPoint); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) root[iPoint] = firstIndex + iPoint; + for (;;) { + unsigned long changed = 0, changedGlobal = 0; + const auto neighborRoot = FetchFromOwners(neighbor, root); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; k++) { + if (lineEdge[k] && (neighborRoot[k] < root[iPoint])) { + root[iPoint] = neighborRoot[k]; + changed++; + } + } + } + const auto rootOfRoot = FetchFromOwners(root, root); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + if (rootOfRoot[iPoint] < root[iPoint]) { + root[iPoint] = rootOfRoot[iPoint]; + changed++; + } + } + SU2_MPI::Allreduce(&changed, &changedGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + if (changedGlobal == 0) break; + } + return root; +} + +unsigned long CPhysicalGeometry::KeepLinesTogether(const vector& root, vector& part, + unsigned long& longestLine) const { + CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); + const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); + + /*--- The owner of each root counts in which partitions its line lies. ---*/ + vector nSend(size, 0), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) nSend[pointPartitioner.GetRankContainingIndex(root[iPoint])] += 2; + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + vector sendBuf(sendDisp[size]), recvBuf(recvDisp[size]); + vector next(sendDisp.begin(), sendDisp.end() - 1); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + auto& slot = next[pointPartitioner.GetRankContainingIndex(root[iPoint])]; + sendBuf[slot++] = root[iPoint]; + sendBuf[slot++] = static_cast(part[iPoint]); + } + SU2_MPI::Alltoallv(sendBuf.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvBuf.data(), nRecv.data(), + recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + map> votes; + for (auto k = 0ul; k < recvBuf.size(); k += 2) votes[recvBuf[k]][recvBuf[k + 1]]++; + + unsigned long longest = 0; + for (const auto& line : votes) { + unsigned long length = 0; + for (const auto& vote : line.second) length += vote.second; + longest = max(longest, length); + } + SU2_MPI::Allreduce(&longest, &longestLine, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm()); + + vector decision(nPoint); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) decision[iPoint] = static_cast(part[iPoint]); + for (const auto& line : votes) { + unsigned long best = 0, bestCount = 0; + for (const auto& vote : line.second) { + if (vote.second > bestCount) { + best = vote.first; + bestCount = vote.second; + } + } + decision[line.first - firstIndex] = best; + } + + /*--- Every point takes the partition of its line, unless that would leave a rank without points. ---*/ + const auto chosen = FetchFromOwners(root, decision); + vector count(size, 0), countGlobal(size); + for (const auto iPart : chosen) count[iPart]++; + SU2_MPI::Allreduce(count.data(), countGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + if (*min_element(countGlobal.begin(), countGlobal.end()) == 0) return 0; + + unsigned long moved = 0, movedGlobal = 0; + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + if (static_cast(part[iPoint]) != chosen[iPoint]) { + part[iPoint] = static_cast(chosen[iPoint]); + moved++; + } + } + SU2_MPI::Allreduce(&moved, &movedGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + return movedGlobal; +} + +/*! + * \brief Compute graph edge weights that encode mesh anisotropy for partitioning. + * + * The routine builds edge weights from geometric information of neighboring points. + * Since element connectivity can reference points owned by other MPI ranks, it first + * exchanges the remote coordinates that are required locally. The resulting weights are + * scaled/clamped to the user-provided \p maxWeight and returned in graph-edge order. + * + * \param[in] maxWeight Upper bound used when converting anisotropy metrics to integer + * partitioning weights. + * \return Edge weights suitable for graph partitioning (one value per graph edge entry). + */ +vector CPhysicalGeometry::ComputeAnisotropyEdgeWeights(long maxWeight) const { + /*--- Build local/global point-index helpers used throughout communication. ---*/ + CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); + const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); + const auto isLocal = [&](unsigned long iGlobal) { + return (iGlobal >= firstIndex) && (iGlobal - firstIndex < nPoint); + }; + + /*--- Identify remote element nodes whose coordinates must be requested. ---*/ + vector> request(size); + for (auto iElem = 0ul; iElem < nElem; iElem++) { + for (auto iNode = 0u; iNode < elem[iElem]->GetnNodes(); iNode++) { + const auto iGlobal = elem[iElem]->GetNode(iNode); + if (!isLocal(iGlobal)) request[pointPartitioner.GetRankContainingIndex(iGlobal)].push_back(iGlobal); + } + } + + /*--- Exchange request sizes (after deduplication) to prepare all-to-all transfers. ---*/ + vector nSend(size), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (int iRank = 0; iRank < size; iRank++) { + auto& list = request[iRank]; + sort(list.begin(), list.end()); + list.erase(unique(list.begin(), list.end()), list.end()); + nSend[iRank] = static_cast(list.size()); + } + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + + vector sendIndex, recvIndex(recvDisp[size]); + sendIndex.reserve(sendDisp[size]); + for (const auto& list : request) sendIndex.insert(sendIndex.end(), list.begin(), list.end()); + + SU2_MPI::Alltoallv(sendIndex.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvIndex.data(), nRecv.data(), + recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + vector replyCoord(recvIndex.size() * nDim), remoteCoord(sendIndex.size() * nDim); + for (auto i = 0ul; i < recvIndex.size(); i++) + for (auto iDim = 0u; iDim < nDim; iDim++) + replyCoord[i * nDim + iDim] = SU2_TYPE::GetValue(nodes->GetCoord(recvIndex[i] - firstIndex, iDim)); + + for (int iRank = 0; iRank < size; iRank++) { + nSend[iRank] *= nDim; + nRecv[iRank] *= nDim; + sendDisp[iRank] *= nDim; + recvDisp[iRank] *= nDim; + } + SU2_MPI::Alltoallv(replyCoord.data(), nRecv.data(), recvDisp.data(), MPI_DOUBLE, remoteCoord.data(), nSend.data(), + sendDisp.data(), MPI_DOUBLE, SU2_MPI::GetComm()); + + unordered_map remoteSlot; + for (auto i = 0ul; i < sendIndex.size(); i++) remoteSlot[sendIndex[i]] = i; + + auto getCoord = [&](unsigned long iGlobal, unsigned short iDim) -> passivedouble { + if (isLocal(iGlobal)) return SU2_TYPE::GetValue(nodes->GetCoord(iGlobal - firstIndex, iDim)); + return remoteCoord[remoteSlot.at(iGlobal) * nDim + iDim]; + }; + + /*--- Weight is the longest node distance of the element over the edge length, maximized over elements. ---*/ + + vector weight(adjacency.size(), 1); + + for (auto iElem = 0ul; iElem < nElem; iElem++) { + const auto nNodes = elem[iElem]->GetnNodes(); + passivedouble coord[N_POINTS_HEXAHEDRON][MAXNDIM] = {}; + for (auto iNode = 0u; iNode < nNodes; iNode++) + for (auto iDim = 0u; iDim < nDim; iDim++) coord[iNode][iDim] = getCoord(elem[iElem]->GetNode(iNode), iDim); + + auto distance2 = [&](unsigned short a, unsigned short b) { + passivedouble d2 = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) d2 += pow(coord[a][iDim] - coord[b][iDim], 2); + return d2; + }; + + passivedouble maxDist2 = 0.0; + for (auto a = 0u; a < nNodes; a++) + for (auto b = a + 1; b < nNodes; b++) maxDist2 = max(maxDist2, distance2(a, b)); + + for (auto a = 0u; a < nNodes; a++) { + const auto iGlobal = elem[iElem]->GetNode(a); + if (!isLocal(iGlobal)) continue; + const auto iLocal = iGlobal - firstIndex; + const auto begin = adjacency.begin() + xadj[iLocal]; + const auto end = adjacency.begin() + xadj[iLocal + 1]; + + for (auto b = 0u; b < nNodes; b++) { + const auto d2 = distance2(a, b); + if ((b == a) || (d2 <= 0.0)) continue; + + const auto jGlobal = static_cast(elem[iElem]->GetNode(b)); + const auto it = lower_bound(begin, end, jGlobal); + if ((it == end) || (*it != jGlobal)) continue; + + const auto ratio = static_cast(sqrt(maxDist2 / d2)); + auto& w = weight[it - adjacency.begin()]; + w = max(w, min(maxWeight, max(1, ratio))); + } + } + } + return weight; +} +#endif + void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { /*--- We need to have parallel support with MPI and have the ParMETIS library compiled and linked for parallel graph partitioning. ---*/ @@ -7412,6 +7683,23 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { vwgt[iPoint] = wp + we * (xadj[iPoint + 1] - xadj[iPoint]); } + /*--- Keep partition boundaries off the short edges of anisotropic cells. The edges that chain + * points into wall-normal lines get the largest weight, so the graph coarsening merges each + * line before it balances. The cap keeps the total edge weight within the range of idx_t. ---*/ + + const bool anisotropic = config->GetParMETIS_AnisotropyWeight(); + vector adjwgt; + vector lineRoot; + if (anisotropic) { + constexpr long maxWeight = (sizeof(idx_t) >= 8) ? 1000000 : 1000; + adjwgt = ComputeAnisotropyEdgeWeights(maxWeight); + vector lineEdge; + lineRoot = LinkAnisotropicLines(adjwgt, lineEdge); + for (auto k = 0ul; k < adjwgt.size(); k++) + if (lineEdge[k]) adjwgt[k] = maxWeight; + wgtflag = 3; + } + /*--- Create some structures that ParMETIS needs to output the partitioning. ---*/ idx_t edgecut; @@ -7420,12 +7708,30 @@ 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); + const int err = ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjacency.data(), vwgt.data(), + anisotropic ? adjwgt.data() : nullptr, &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; + cout << " graph partitioning complete (" << edgecut << (anisotropic ? " weighted" : "") << " edge cuts)." << endl; + } + + /*--- A line split between ranks loses its strongest couplings in the preconditioner, and its + * fragments give tangled coarse control volumes. ---*/ + + if (anisotropic) { + unsigned long longestLine = 0; + const auto moved = KeepLinesTogether(lineRoot, part, longestLine); + /*--- The imbalance is measured in the vertex weights ParMETIS balances. ---*/ + vector load(size, 0), loadGlobal(size); + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) load[part[iPoint]] += vwgt[iPoint]; + SU2_MPI::Allreduce(load.data(), loadGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, comm); + const auto maxLoad = *max_element(loadGlobal.begin(), loadGlobal.end()); + const auto sumLoad = accumulate(loadGlobal.begin(), loadGlobal.end(), 0ul); + if (rank == MASTER_NODE) { + cout << "Longest anisotropic line " << longestLine << " points, moved " << moved + << " points to keep lines on one rank, load imbalance " << double(maxLoad) * size / sumLoad << "." << endl; + } } /*--- Store the results of the partitioning (note that this is local diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 4e389fd0b3dc..19717ff37814 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -31,6 +31,9 @@ #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/toolboxes/allocation_toolbox.hpp" +#include + +#include #include #include #include @@ -350,9 +353,19 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi /*--- The GPU triangular solves are level-scheduled (exact, one pass per level), so they need * levels_ilu unconditionally; the host/OMP path only needs it when both multi-threaded and * requested via config. ---*/ - if (useCuda || (omp_get_max_threads() > 1 && config->GetIluOptions().LevelScheduling)) { + if (useCuda) { levels_ilu = computeLevels(pat_ilu.l); + } else if (omp_get_max_threads() > 1 && config->GetIluOptions().LevelScheduling) { + /*--- The levels only cover the domain points, like the thread partitions. ---*/ + const auto& pat_l = pat_ilu.l; + su2vector outerPtr(nPointDomain + 1); + for (auto i = 0ul; i <= nPointDomain; ++i) outerPtr(i) = pat_l.outerPtr()[i]; + su2vector innerIdx(outerPtr(nPointDomain)); + for (auto k = 0ul; k < innerIdx.size(); ++k) innerIdx(k) = pat_l.innerIdx()[k]; + levels_ilu = computeLevels(CCompressedSparsePatternUL(std::move(outerPtr), std::move(innerIdx))); } + level_waits_fwd.clear(); + level_waits_bwd.clear(); /*--- Coloring for the GPU iterative factorization, see IluFactorColorKernel. Colors are * true independent sets of the (domain-only, symmetric) dependency graph, computed the same @@ -521,6 +534,8 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi } } + if (ilu_needed && !levels_ilu.empty() && !useCuda && num_threads > 1) BuildLevelWaits(); + /*--- Generate MKL Kernels ---*/ #ifdef USE_MKL @@ -1060,6 +1075,102 @@ void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector +void CSysMatrix::BuildLevelWaits() { + const auto nThreads = static_cast(omp_get_max_threads()); + const auto nLevels = levels_ilu.getOuterSize(); + + /*--- Level and owning thread of each row, the rows of a level are split in equal contiguous slices. ---*/ + + std::vector level(nPointDomain, 0), owner(nPointDomain, 0); + for (unsigned long iLevel = 0; iLevel < nLevels; ++iLevel) { + const unsigned long n = levels_ilu.getNumNonZeros(iLevel); + for (unsigned long iThread = 0; iThread < nThreads; ++iThread) { + for (auto k = n * iThread / nThreads; k < n * (iThread + 1) / nThreads; ++k) { + const auto iPoint = levels_ilu.getInnerIdx(iLevel, k); + level[iPoint] = iLevel; + owner[iPoint] = iThread; + } + } + } + + std::vector fwd(nLevels * nThreads * nThreads, 0), bwd(fwd); + bool valid = true; + + for (unsigned long iLevel = 0; iLevel < nLevels; ++iLevel) { + const unsigned long n = levels_ilu.getNumNonZeros(iLevel); + const auto iStep = nLevels - 1 - iLevel; // Position in the backward sweep. + + for (unsigned long iThread = 0; iThread < nThreads; ++iThread) { + auto* wait_f = &fwd[(iLevel * nThreads + iThread) * nThreads]; + auto* wait_b = &bwd[(iStep * nThreads + iThread) * nThreads]; + + for (auto k = n * iThread / nThreads; k < n * (iThread + 1) / nThreads; ++k) { + const auto iPoint = levels_ilu.getInnerIdx(iLevel, k); + + for (auto kl = ilu.row_ptr_l[iPoint]; kl < ilu.row_ptr_l[iPoint + 1]; ++kl) { + const auto jPoint = ilu.col_ind_l[kl]; + if (level[jPoint] >= iLevel) valid = false; + wait_f[owner[jPoint]] = std::max(wait_f[owner[jPoint]], level[jPoint] + 1); + } + for (auto ku = ilu.row_ptr_u[iPoint]; ku < ilu.row_ptr_u[iPoint + 1]; ++ku) { + const auto jPoint = ilu.col_ind_u[ku]; + if (jPoint >= nPointDomain) break; + if (level[jPoint] <= iLevel) valid = false; + const auto jStep = nLevels - 1 - level[jPoint]; + wait_b[owner[jPoint]] = std::max(wait_b[owner[jPoint]], jStep + 1); + } + } + wait_f[iThread] = 0; + wait_b[iThread] = 0; + } + } + + /*--- The backward sweep needs a symmetric pattern, otherwise keep the barriers. ---*/ + if (!valid) return; + + level_waits_fwd = std::move(fwd); + level_waits_bwd = std::move(bwd); + level_num_threads = nThreads; + level_progress = std::vector(nThreads); +} + +template +template +void CSysMatrix::LevelSweep(bool forward, const F& row_func) const { + const auto nThreads = level_num_threads; + const auto iThread = static_cast(omp_get_thread_num()); + const auto nLevels = levels_ilu.getOuterSize(); + const auto& waits = forward ? level_waits_fwd : level_waits_bwd; + + /*--- All counters are equal at the start of a sweep. ---*/ + SU2_OMP_BARRIER + + auto& progress = level_progress[iThread].steps; + const uint64_t base = progress.load(std::memory_order_relaxed); + + for (unsigned long iStep = 0; iStep < nLevels; ++iStep) { + const auto* wait = &waits[(iStep * nThreads + iThread) * nThreads]; + + for (unsigned long jThread = 0; jThread < nThreads; ++jThread) { + if (wait[jThread] == 0) continue; + const uint64_t target = base + wait[jThread]; + for (unsigned spin = 0; level_progress[jThread].steps.load(std::memory_order_acquire) < target; ++spin) { + if (spin > 4096) sched_yield(); + } + } + + const auto iLevel = forward ? iStep : nLevels - 1 - iStep; + const unsigned long n = levels_ilu.getNumNonZeros(iLevel); + for (auto k = n * iThread / nThreads; k < n * (iThread + 1) / nThreads; ++k) { + row_func(levels_ilu.getInnerIdx(iLevel, k)); + } + progress.store(base + iStep + 1, std::memory_order_release); + } + + SU2_OMP_BARRIER +} + template void CSysMatrix::BuildILUPreconditioner() { SU2_ZONE_SCOPED @@ -1159,7 +1270,12 @@ void CSysMatrix::BuildILUPreconditioner() { InvertDiagonalBlockILUMatrix(iPoint); }; - if (levels_ilu.empty()) { + if (UseLevelWaits()) { + LevelSweep(true, [&](const auto iPoint) { + InitIluRow(iPoint); + BuildIluRow(iPoint, 0ul, nPointDomain); + }); + } else if (levels_ilu.empty()) { /*--- Each OMP thread will work on the submatrix defined from * row/col "begin" to row/col "end-1" (i.e. the range [begin,end[). * Which is exactly what the MPI-only implementation does. ---*/ @@ -1256,7 +1372,10 @@ void CSysMatrix::ComputeILUPreconditioner(const CSysVector::SetupInnerSolver(unsigned short kind_solver, const C template unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, const CSysVector& LinSysRes, CSysVector& LinSysSol, CGeometry* geometry, - const CConfig* config) { + const CConfig* config, unsigned long maxIter) { SU2_ZONE_SCOPED /*--- A word about the templated types. It is assumed that the residual and solution vectors are always of su2doubles, @@ -1458,7 +1458,7 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con default: { KindSolver = config->GetKind_Linear_Solver(); KindPrecond = config->GetKind_Linear_Solver_Prec(); - MaxIter = config->GetLinear_Solver_Iter(); + MaxIter = (maxIter != 0) ? maxIter : config->GetLinear_Solver_Iter(); SolverTol = SU2_TYPE::GetValue(config->GetLinear_Solver_Error()); ScreenOutput = false; break; diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index fbabb5c35867..85afd105ddeb 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -26,6 +26,9 @@ */ #include "CIntegration.hpp" +#include "../../../Common/include/containers/container_decorators.hpp" + +class CEulerVariable; /*! * \class CMultiGridIntegration @@ -103,10 +106,12 @@ class CMultiGridIntegration final : public CIntegration { /*! * \brief Set the value of the corrected fine grid solution. * \param[out] sol_fine - Pointer to the solution on the fine grid. + * \param[in] sol_flow - Flow solver on the fine grid, which flags the points that get no correction. * \param[in] geo_fine - Geometrical definition of the fine grid. * \param[in] config - Definition of the particular problem. */ - void SetProlongated_Correction(CSolver *sol_fine, CGeometry *geo_fine, CConfig *config, unsigned short iMesh); + void SetProlongated_Correction(CSolver *sol_fine, CSolver *sol_flow, CGeometry *geo_fine, CConfig *config, + unsigned short iMesh); /*! * \brief Compute the gradient in coarse grid using the fine grid information. @@ -141,9 +146,22 @@ class CMultiGridIntegration final : public CIntegration { * \param[in] geo_fine - Geometrical definition of the fine grid. * \param[in] geo_coarse - Geometrical definition of the coarse grid. * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the coarse mesh. */ void SetProlongated_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, unsigned short iMesh); + + /*! + * \brief Set the solution of every fine child from its coarse parent, constant or linearly reconstructed. + * \param[in,out] sol_fine - Solver on the fine grid. + * \param[in] sol_coarse - Solver on the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the coarse mesh. + */ + void InterpolateSolutionToChildren(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, + CGeometry *geo_coarse, CConfig *config, unsigned short iMesh); /*! @@ -190,7 +208,21 @@ class CMultiGridIntegration final : public CIntegration { * \param[in] config - Definition of the particular problem. */ void GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh); + + /*! + * \brief Least-squares gradient of the coarse-grid correction, limited so no child value + * leaves the range of the coarse stencil. + * \param[in] sol_coarse - Solver holding the correction in Solution_Old, or the solution itself. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the coarse mesh. + * \param[in] useSolution - Differentiate Solution instead of Solution_Old. + */ + void ComputeProlongationGradient(CSolver *sol_coarse, CGeometry *geo_coarse, CGeometry *geo_fine, + const CConfig *config, unsigned short iMesh, bool useSolution); /*! * \brief Do an implicit smoothing of the prolongated correction. @@ -296,9 +328,11 @@ class CMultiGridIntegration final : public CIntegration { * \param[in] geo_coarse - Geometry of the level handing over. * \param[in] config - Definition of the particular problem. * \param[in] eddy_viscosity - Carry the eddy viscosity across as well. + * \param[in] iMesh - Index of the coarse mesh. */ void SetProlongated_ScalarSolution(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, - CGeometry *geo_coarse, CConfig *config, bool eddy_viscosity); + CGeometry *geo_coarse, CConfig *config, bool eddy_viscosity, + unsigned short iMesh); /*! * \brief Helper function for early-exit logic during pre/post-smoothing. @@ -319,7 +353,8 @@ class CMultiGridIntegration final : public CIntegration { passivedouble lastRMS[2], char& exitReason, passivedouble& worstStepRatio, unsigned short& worstStep); - static constexpr int MAX_MG_LEVELS = 10; + /*--- CConfig rejects a larger MGLEVEL, so the per-level arrays below always fit. ---*/ + static constexpr int MAX_MG_LEVELS = MAX_MGLEVELS; /*--- Bounds the fixed-size stack buffers in the restriction and prolongation kernels, * independently of the CSysMatrix limit of the same name. ---*/ @@ -354,6 +389,14 @@ class CMultiGridIntegration final : public CIntegration { unsigned short lastPreSmoothWorstStep[MAX_MG_LEVELS+1] = {}; unsigned short lastPostSmoothWorstStep[MAX_MG_LEVELS+1] = {}; + /*--- Linear solver iterations and final relative residual of each smoothing step, per level. ---*/ + struct LinSolInfo { + unsigned long iters; + passivedouble res; + }; + std::vector lastPreLinSol[MAX_MG_LEVELS+1]; + std::vector lastPostLinSol[MAX_MG_LEVELS+1]; + /*! \brief FinestMesh observed on the previous call; 0 also serves as the sentinel forcing a * reset on the first call, since a FULL cycle starts at FinestMesh == nMGLevels > 0. */ unsigned short mg_ramp_last_FinestMesh = 0; @@ -368,6 +411,24 @@ class CMultiGridIntegration final : public CIntegration { enum class MGStartupPromote { NONE, BUDGET, CONVERGENCE, STAGNATION }; MGStartupPromote mg_startup_promote_reason = MGStartupPromote::NONE; + /*! \brief Limited least-squares gradient of the correction, indexed by coarse level. + * Allocated on first use, only when MG_LINEAR_PROLONGATION is on. */ + vector prolongGradient; + vector> correctionWeight; /*!< \brief Per level, weight of the correction near wall ends. */ + + /*! + * \brief Weight each level's prolongated correction from zero where a wall ends to one a few coarse CVs away. + * \param[in] geometry - Geometry of all levels of this zone and instance. + * \param[in] flowNodes - Flow variables of the finest grid, which flag the wall ends. + * \param[in] nLevels - Number of multigrid levels. + */ + void ComputeCorrectionWeights(CGeometry** geometry, const CEulerVariable* flowNodes, unsigned short nLevels); + + /*! \brief Coarse cells per level whose gradient fit is too ill-conditioned and keep the constant + * operator, and whether that level has been reported. */ + unsigned long prolongIllConditioned[MAX_MG_LEVELS + 1] = {}; + bool prolongReported[MAX_MG_LEVELS + 1] = {}; + vector mg_startup_conv_start; /*!< \brief Field values when the active level became active. */ vector mg_startup_conv_prev; /*!< \brief Field values on the previous iteration. */ unsigned long mg_startup_stall_count = 0; /*!< \brief Consecutive iterations without useful reduction. */ diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index 800167d4b0f5..c683d9445b11 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -119,9 +119,16 @@ FORCEINLINE CPair reconstructPrimitives(const Int& iEdge, fmax(solution.OutlierMitigation(iPoint[k]), solution.OutlierMitigation(jPoint[k])) / VariableType::MAX_OUTLIER_MITIGATION); } + /*--- Edges of flagged points drop to first order without the extra dissipation of a + * non-physical edge. ---*/ + Double firstOrder = nonPhysical; + for (size_t k = 0; k < Double::Size; ++k) { + firstOrder[k] = fmax(firstOrder[k], fmax(solution.FirstOrderReconstruction(iPoint[k]), + solution.FirstOrderReconstruction(jPoint[k]))); + } for (size_t iVar = 0; iVar < ReconVarType::nVar; ++iVar) { - V.i.all(iVar) = nonPhysical * V1st.i.all(iVar) + (1-nonPhysical) * V.i.all(iVar); - V.j.all(iVar) = nonPhysical * V1st.j.all(iVar) + (1-nonPhysical) * V.j.all(iVar); + V.i.all(iVar) = firstOrder * V1st.i.all(iVar) + (1-firstOrder) * V.i.all(iVar); + V.j.all(iVar) = firstOrder * V1st.j.all(iVar) + (1-firstOrder) * V.j.all(iVar); } } else { nonPhysical = 0; diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp index 1a8745b10bc2..3edc0232abd9 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -76,7 +76,6 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { const su2double gamma; const su2double gasConst; const su2double prandtlTurb; - const bool correct; const bool useSA_QCR; const bool wallFun; const bool uq; @@ -91,12 +90,11 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { * \brief Constructor, initialize constants and booleans. */ template - CCompressibleViscousFluxBase(const CConfig& config, int iMesh, + CCompressibleViscousFluxBase(const CConfig& config, int, const CVariable* turbVars_, Ts&...) : gamma(config.GetGamma()), gasConst(config.GetGas_ConstantND()), prandtlTurb(config.GetPrandtl_Turb()), - correct(iMesh == MESH_0), useSA_QCR(config.GetSAParsedOptions().qcr2000), wallFun(config.GetWall_Functions()), uq(config.GetSSTParsedOptions().uq), @@ -142,10 +140,10 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { Double mask = dist2_ij < EPS*EPS; dist2_ij += mask / (EPS*EPS); - /*--- Compute the corrected mean gradient. ---*/ + /*--- Compute the corrected mean gradient on all levels. ---*/ auto avgGrad = averageGradient(iPoint, jPoint, gradient); - if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad); + correctGradient(V, vector_ij, dist2_ij, avgGrad); /*--- Stress and heat flux tensors. ---*/ @@ -179,7 +177,10 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { /*--- Flux Jacobians. ---*/ - Double dist_ij = sqrt(dist2_ij); + /*--- Effective distance |d|^2/|d.n| matches the normal derivative of the corrected gradient. ---*/ + + const Double projDist = abs(dot(vector_ij, unitNormal)); + const Double dist_ij = dist2_ij / fmax(projDist, 0.1 * sqrt(dist2_ij)); auto dtau = stressTensorJacobian(avgV, unitNormal, dist_ij); /*--- Energy flux Jacobian. ---*/ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 3a46fc8c4db6..35f020b4ad40 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -629,7 +629,8 @@ void CFVMFlowSolverBase::ImplicitEuler_Iteration(CGeometry *geometry, CSol } END_SU2_OMP_FOR - auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); + auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config, + config->GetMGOptions().MG_LinSolIter[geometry->GetMGLevel()]); BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { SetIterLinSolver(iter); diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index af03659fe7f5..96c5f37e5cf6 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -586,7 +586,8 @@ void CScalarSolver::ImplicitEuler_Iteration(CGeometry* geometry, C } END_SU2_OMP_FOR - auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); + auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config, + config->GetMGOptions().MG_LinSolIter[geometry->GetMGLevel()]); BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { SetIterLinSolver(iter); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index e5cb9881077a..a7a62db8a964 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -82,6 +82,39 @@ class CSolver { vector NonLinRes_Series; /*!< \brief Vector holding the nonlinear residual indicator series. */ su2double Old_Func, /*!< \brief Old value of the nonlinear residual indicator. */ New_Func; /*!< \brief Current value of the nonlinear residual indicator. */ + /*--- State of the residual-driven CFL adaption, a negative ceiling or cap means none is active. ---*/ + vector ResTrend_History; /*!< \brief Recent mean log residuals. */ + vector> ResJump_History; /*!< \brief Recent log residuals of the flow variables. */ + vector> ResFlow_History; /*!< \brief Recent log residuals of the flow variables. */ + vector ResMean_History; /*!< \brief Recent mean log residuals of the flow variables. */ + passivedouble ResRise = 0.0; /*!< \brief Rise of the flow residual above its recent best. */ + unsigned long CFL_LastFlip = 0; /*!< \brief Iteration of the last flip-flop. */ + unsigned long CFL_LastRegress = 0; /*!< \brief Iteration of the last residual regression. */ + bool ResTrend_Increase = false; /*!< \brief The residual trend allows a larger CFL. */ + bool ResTrend_Reduce = false; /*!< \brief The residual trend requires a smaller CFL. */ + bool CFL_Cut = false; /*!< \brief Halve the CFL in this update. */ + passivedouble CFL_Ceiling = -1.0; /*!< \brief Temporary ceiling after a cut. */ + passivedouble CFL_Cap = -1.0; /*!< \brief Lasting limit learned from repeated events. */ + vector CFL_EventValues; /*!< \brief CFL at each cut. */ + unsigned long CFL_LastEvent = 0; /*!< \brief Iteration of the last cut. */ + unsigned long CFL_LastCapProbe = 0; /*!< \brief Iteration of the last cap increase. */ + vector> ResStall_History; /*!< \brief Log flow residuals over two stall windows. */ + vector CFLStall_History; /*!< \brief Mean CFL over two stall windows. */ + unsigned long Stall_Streak = 0; /*!< \brief Consecutive iterations without envelope progress. */ + unsigned long CFL_LastStall = 0; /*!< \brief Iteration of the last plateau. */ + passivedouble CFL_StallCap = -1.0; /*!< \brief Cap set at the last stall, which cap increases stay below. */ + passivedouble CFL_OscLevel = -1.0; /*!< \brief Lowest CFL at which the residuals flip-flopped. */ + unsigned long CFL_LastOscRelease = 0; /*!< \brief Iteration of the last flip-flop limit increase. */ + passivedouble Res_AtOsc = 0.0; /*!< \brief Mean log flow residual at that flip-flop. */ + bool Trial_Active = false; /*!< \brief CFL cuts at a residual plateau are being evaluated. */ + unsigned long Trial_Start = 0; /*!< \brief Iteration of the first cut of the trial. */ + unsigned long Trial_Hold = 0; /*!< \brief No trial starts before this iteration. */ + unsigned long Trial_Settle = 0; /*!< \brief Iteration after which the lowest CFL of the trial is judged. */ + passivedouble Trial_Res = 0.0; /*!< \brief Mean log flow residual when the trial started. */ + passivedouble Trial_CFL = 0.0; /*!< \brief Mean CFL when the trial started. */ + passivedouble Trial_Ceiling = -1.0; /*!< \brief CFL ceiling before the trial. */ + passivedouble Trial_Cap = -1.0; /*!< \brief CFL cap before the trial. */ + passivedouble Trial_StallCap = -1.0; /*!< \brief Plateau cap before the trial. */ unsigned short nVar, /*!< \brief Number of variables of the problem. */ nPrimVar, /*!< \brief Number of primitive variables of the problem. */ nPrimVarGrad, /*!< \brief Number of primitive variables of the problem in the gradient computation. */ diff --git a/SU2_CFD/include/variables/CAdjEulerVariable.hpp b/SU2_CFD/include/variables/CAdjEulerVariable.hpp index 1faf911816aa..cb039f291358 100644 --- a/SU2_CFD/include/variables/CAdjEulerVariable.hpp +++ b/SU2_CFD/include/variables/CAdjEulerVariable.hpp @@ -102,6 +102,11 @@ class CAdjEulerVariable : public CVariable { Solution_Old(iPoint,iDim+1) = val_velocity[iDim]*Solution(iPoint,0); } + /*! + * \brief Index of the adjoint momentum in the solution vector. + */ + inline short GetVelocityIndex() const final { return 1; } + /*! * \brief Set the momentum part of the truncation error to zero. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index abd63a3afd27..2466d6e61698 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -298,6 +298,11 @@ class CEulerVariable : public CFlowVariable { Solution_Old(iPoint,iDim+1) = val_velocity[iDim]*Solution(iPoint,0); } + /*! + * \brief Index of the momentum in the conservative solution vector. + */ + inline short GetVelocityIndex() const final { return 1; } + /*! * \brief Set the momentum part of the truncation error to zero. * \param[in] iPoint - Point index. @@ -369,4 +374,14 @@ class CEulerVariable : public CFlowVariable { */ su2vector OutlierMitigation; + /*! + * \brief Points whose edges are reconstructed to first order (1), others use MUSCL (0). + */ + su2vector FirstOrderReconstruction; + + /*! + * \brief Points, and their neighbours, where the multigrid correction is not prolongated (1), others (0). + */ + su2vector SkipMGCorrection; + }; diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index aaa6c437b665..7fc1860a826a 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -213,6 +213,11 @@ class CIncEulerVariable : public CFlowVariable { Solution_Old(iPoint,iDim+1) = val_velocity[iDim]; } + /*! + * \brief Index of the velocity in the solution vector. + */ + inline short GetVelocityIndex() const final { return 1; } + /*! * \brief Set the momentum part of the truncation error to zero. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index e0d4199144c4..37229a595945 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -298,6 +298,12 @@ class CNEMOEulerVariable : public CFlowVariable { } } + /*! + * \brief Index of the momentum in the conservative solution vector, which the primitive + * index Velocity() does not give. + */ + inline short GetVelocityIndex() const final { return static_cast(nSpecies); } + /*! * \brief A virtual member. * \return Value of the vibrational-electronic temperature. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 5fdad211e4dc..6459a409079b 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1621,6 +1621,12 @@ class CVariable { */ inline virtual void SetVelocity_Old(unsigned long iPoint, const su2double *val_velocity) {} + /*! + * \brief Index of the first velocity or momentum component in the solution vector. + * \return The index, or -1 when the solution carries no velocity. + */ + inline virtual short GetVelocityIndex() const { return -1; } + /*! * \brief A virtual member. * \param[in] laminarViscosity diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 965f8ba587a6..395ae8f4b18f 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -881,10 +881,6 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { geometry[iMGlevel]->SetCoord(geometry[iMGlevel-1]); - /*--- Find closest, most normal, neighbor to a surface point ---*/ - - geometry[iMGlevel]->FindNormal_Neighbor(config); - /*--- Store our multigrid index. ---*/ geometry[iMGlevel]->SetMGLevel(iMGlevel); @@ -934,6 +930,18 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + /*--- SetCoord centred a halo agglomerate on the partial child list this rank holds. Take + the owner's coordinate, otherwise coarse stencils depend on the partitioning. ---*/ + + if (iMGlevel > MESH_0) { + geometry[iMGlevel]->InitiateComms(geometry[iMGlevel], config, MPI_QUANTITIES::COORDINATES); + geometry[iMGlevel]->CompleteComms(geometry[iMGlevel], config, MPI_QUANTITIES::COORDINATES); + + /*--- Find closest, most normal, neighbor to a surface point ---*/ + + geometry[iMGlevel]->FindNormal_Neighbor(config); + } + /*--- Compute the max length. ---*/ if (!fea) { @@ -1809,7 +1817,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver /*--- Compressible flow Ideal gas ---*/ numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, true, config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, true, config); /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) @@ -1820,7 +1828,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver /*--- Compressible flow Real gas ---*/ numerics[MESH_0][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, true, config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, true, config); /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) @@ -1832,7 +1840,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver /*--- Incompressible flow, use preconditioning method ---*/ numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, true, config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, true, config); /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 607a0915d3c8..7d988cfe2a9d 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -28,7 +28,10 @@ #include "../../include/integration/CMultiGridIntegration.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include +#include +#include "../../include/variables/CEulerVariable.hpp" namespace { @@ -72,7 +75,9 @@ void CMultiGridIntegration::adaptDampingFactors(CConfig* config, passivedouble c applyGlobalTrend(config->GetDamp_Correc_Prolong(), crossCycleRatio)); } -CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { } +CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { + prolongGradient.resize(MAX_MG_LEVELS + 1); +} void CMultiGridIntegration::MonitorFullMG_Startup(const vector >& convFields, const CConfig *config) { @@ -101,14 +106,15 @@ void CMultiGridIntegration::MonitorFullMG_Startup(const vector 0.0 && !mg_startup_conv_prev.empty()) { + const passivedouble band = fabs(log10(stall_tol)); bool stalled = true; for (auto iField = 0ul; iField < nFields; iField++) - stalled = stalled && (convFields[iField].second - mg_startup_conv_prev[iField] >= log10(stall_tol)); + stalled = stalled && (fabs(convFields[iField].second - mg_startup_conv_prev[iField]) <= band); if (stalled) mg_startup_stall_count++; @@ -343,7 +349,7 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, solver_container[iZone][iInst][FinestMesh][Solver_Position], geometry[iZone][iInst][FinestMesh-1], geometry[iZone][iInst][FinestMesh], - config[iZone]); + config[iZone], FinestMesh); /*--- The scalar equations only restrict downward, so hand the new level their solution too, * otherwise they restart from their initial condition at every promotion. ---*/ @@ -357,7 +363,7 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, SetProlongated_ScalarSolution(scalar_fine, scalar_coarse, geometry[iZone][iInst][FinestMesh-1], geometry[iZone][iInst][FinestMesh], - config[iZone], Scalar_Position == TURB_SOL); + config[iZone], Scalar_Position == TURB_SOL, FinestMesh); } /*--- Report the promotion before the startup state is reset below. ---*/ @@ -451,6 +457,15 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, geometry[iZone][iInst][FinestMesh], config[iZone], true); } + /*--- The scalar solvers run on the active startup level and need its vorticity and strain rate. ---*/ + + if (fmg_warmup) { + solver_container[iZone][iInst][FinestMesh][Solver_Position]->Preprocessing(geometry[iZone][iInst][FinestMesh], + solver_container[iZone][iInst][FinestMesh], + config[iZone], FinestMesh, NO_RK_ITER, + RunTime_EqSystem, true); + } + /*--- Computes primitive variables and gradients in the finest mesh (useful for the next solver (turbulence) and output ---*/ solver_container[iZone][iInst][MESH_0][Solver_Position]->Preprocessing(geometry[iZone][iInst][MESH_0], @@ -537,6 +552,26 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, table << "-"; /*--- Corr.-smooth: defined on levels 0..nMGLevels-1; coarsest has none. ---*/ + /*--- Linear solver iterations and final relative residual of each smoothing step. ---*/ + auto linRows = [&](const char* name, const std::vector* info, unsigned short nLev) { + std::vector its(nMGLevels + 1, "-"), res(nMGLevels + 1, "-"); + for (unsigned short i = 0; i < nLev; ++i) { + std::ostringstream si, sr; + sr << std::scientific << std::setprecision(1); + for (std::size_t k = 0; k < info[i].size(); ++k) { + si << (k ? "," : "") << info[i][k].iters; + sr << (k ? "," : "") << info[i][k].res; + } + if (!info[i].empty()) { its[i] = si.str(); res[i] = sr.str(); } + } + table << std::string(name) + " its"; + for (const auto& str : its) table << str; + table << std::string(name) + " res"; + for (const auto& str : res) table << str; + }; + linRows("Pre-lin", lastPreLinSol, nMGLevels + 1); + linRows("Post-lin", lastPostLinSol, nMGLevels); + table << "Corr-smooth"; for (unsigned short i = 0; i < nMGLevels; ++i) table << cellStr(lastCorrecSmoothIters[i], mgOptsZone.MG_CorrecSmooth[i], @@ -599,6 +634,14 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, unsigned short iRKLimit = config->GetnRKStep(); + /*--- Wall end weights for every level, once the hierarchy exists. ---*/ + if ((iMesh == MESH_0) && correctionWeight.empty()) { + const auto* flowNodes = dynamic_cast(solver_container_fine[FLOW_SOL]->GetNodes()); + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + ComputeCorrectionWeights(geometry[iZone][iInst], flowNodes, config->GetnMGLevels()); + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + /*--- Do a presmoothing on the grid iMesh to be restricted to the grid iMesh+1 ---*/ PreSmoothing(RunTime_EqSystem, geometry, solver_container, config_container, solver_fine, numerics_fine, @@ -670,12 +713,13 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, /*--- Compute prolongated solution, and smooth the correction $u^(new)_k = u_k + Smooth(I^k_(k+1)(u_(k+1)-I^(k+1)_k u_k))$ ---*/ - GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, + iMesh+1); const auto& mgOpts = config->GetMGOptions(); SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, mgOpts.MG_CorrecSmooth[iMesh], mgOpts.MG_Smooth_Coeff, config, iMesh); - SetProlongated_Correction(solver_fine, geometry_fine, config, iMesh); + SetProlongated_Correction(solver_fine, solver_container_fine[FLOW_SOL], geometry_fine, config, iMesh); /*--- Solution post-smoothing in the prolongated grid. ---*/ @@ -743,7 +787,7 @@ void CMultiGridIntegration::PreSmoothing(unsigned short RunTime_EqSystem, ? SU2_TYPE::GetValue(mgOpts.MG_Smooth_StagnationTol) : passivedouble(1.0); /*--- Reset the shared early-exit flag (master only). ---*/ - SU2_OMP_SAFE_GLOBAL_ACCESS(mg_early_exit_flag = false;) + SU2_OMP_SAFE_GLOBAL_ACCESS(mg_early_exit_flag = false; lastPreLinSol[iMesh].clear();) for (unsigned short iPreSmooth = 0; iPreSmooth < nPreSmooth; iPreSmooth++) { /*--- Time and space integration ---*/ @@ -768,6 +812,14 @@ void CMultiGridIntegration::PreSmoothing(unsigned short RunTime_EqSystem, /*--- Time integration, update solution using the old solution plus the solution increment ---*/ Time_Integration(geometry_fine, solver_container_fine, config, iRKStep, RunTime_EqSystem); + if (mgOpts.MG_Smooth_Output) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + lastPreLinSol[iMesh].push_back({solver_fine->GetIterLinSolver(), + SU2_TYPE::GetValue(solver_fine->GetResLinSolver())}); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + /*--- At iRKStep==0 LinSysRes = R(u_k). ---*/ if (iRKStep == 0 && need_initial_rms) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { @@ -819,7 +871,7 @@ void CMultiGridIntegration::PostSmoothing(unsigned short RunTime_EqSystem, ? SU2_TYPE::GetValue(mgOpts.MG_Smooth_StagnationTol) : passivedouble(1.0); /*--- Reset the shared early-exit flag (master only). ---*/ - SU2_OMP_SAFE_GLOBAL_ACCESS(mg_early_exit_flag = false;) + SU2_OMP_SAFE_GLOBAL_ACCESS(mg_early_exit_flag = false; lastPostLinSol[iMesh].clear();) /*--- Do a postsmoothing on the grid iMesh after prolongation from the grid iMesh+1 ---*/ for (unsigned short iPostSmooth = 0; iPostSmooth < nPostSmooth; iPostSmooth++) { @@ -841,6 +893,14 @@ void CMultiGridIntegration::PostSmoothing(unsigned short RunTime_EqSystem, /*--- Time integration, update solution using the old solution plus the solution increment ---*/ Time_Integration(geometry_fine, solver_container_fine, config, iRKStep, RunTime_EqSystem); + if (mgOpts.MG_Smooth_Output) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + lastPostLinSol[iMesh].push_back({solver_fine->GetIterLinSolver(), + SU2_TYPE::GetValue(solver_fine->GetResLinSolver())}); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + /*--- At iRKStep==0 LinSysRes = R(u_k) ---*/ if (iRKStep == 0 && need_per_step_rms) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { @@ -873,7 +933,8 @@ void CMultiGridIntegration::PostSmoothing(unsigned short RunTime_EqSystem, void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh) { SU2_ZONE_SCOPED const unsigned short nVar = sol_coarse->GetnVar(); @@ -941,17 +1002,328 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS /*--- Interpolate the coarse-grid correction onto the fine * grid and store in LinSysRes. ---*/ - /*--- Halos too: the correction smoother reads them before its first exchange. ---*/ - SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) - for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { + if (!config->GetMGOptions().MG_Linear_Prolongation) { + + /*--- Halos too: the correction smoother reads them before its first exchange. ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { + const auto* Correction = sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + sol_fine->LinSysRes.SetBlock(Point_Fine, Correction); + } + } + END_SU2_OMP_FOR + + return; + } + + /*--- The coarse coordinate is the volume-weighted centroid of the children, so the + * reconstruction below averages back to the coarse correction exactly. ---*/ + + ComputeProlongationGradient(sol_coarse, geo_coarse, geo_fine, config, iMesh, false); + + const unsigned short nDim = geo_coarse->GetnDim(); + const auto& gradient = prolongGradient[iMesh]; + + /*--- A halo agglomerate holds a partial child list, so only owners write, and the + * exchange below fills the fine halos the correction smoother reads. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + const auto* Correction = sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse); + const auto* Coord_Coarse = geo_coarse->nodes->GetCoord(Point_Coarse); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - sol_fine->LinSysRes.SetBlock(Point_Fine, Correction); + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_fine->nodes->GetCoord(Point_Fine), Coord_Coarse, Distance); + + su2double Solution[MAXNVAR] = {0.0}; + for (auto iVar = 0u; iVar < nVar; iVar++) { + Solution[iVar] = Correction[iVar]; + for (auto iDim = 0u; iDim < nDim; iDim++) + Solution[iVar] += gradient(Point_Coarse, iVar, iDim) * Distance[iDim]; + } + sol_fine->LinSysRes.SetBlock(Point_Fine, Solution); } } END_SU2_OMP_FOR + SU2_OMP_BARRIER + CSysMatrixComms::Initiate(sol_fine->LinSysRes, geo_fine, config); + CSysMatrixComms::Complete(sol_fine->LinSysRes, geo_fine, config); + +} + +/*--- Largest condition number of the stencil directions for which the prolongation gradient is fitted. ---*/ +static constexpr passivedouble PROLONG_MAX_CONDITION = 1e3; + +/*--- Stencil aspect ratio up to which the fitted gradient is used in full, and from which the cell + * keeps the constant operator. In between the gradient fades out with the logarithm of the ratio. ---*/ +static constexpr passivedouble PROLONG_FULL_ASPECT = 10.0; +static constexpr passivedouble PROLONG_MAX_ASPECT = 100.0; + +namespace { +/*--- Condition number of a symmetric 2x2 or 3x3 matrix, infinite if it is not positive definite. ---*/ +passivedouble FitConditionNumber(unsigned short nDim, const su2double A[3][3]) { + constexpr passivedouble big = 1e300; + passivedouble lmin, lmax; + if (nDim == 2) { + const passivedouble a = SU2_TYPE::GetValue(A[0][0]), b = SU2_TYPE::GetValue(A[0][1]), d = SU2_TYPE::GetValue(A[1][1]); + const passivedouble mean = 0.5 * (a + d), radius = sqrt(0.25 * (a - d) * (a - d) + b * b); + lmin = mean - radius; + lmax = mean + radius; + } else { + passivedouble M[3][3]; + for (auto i = 0u; i < 3; ++i) + for (auto j = 0u; j < 3; ++j) M[i][j] = SU2_TYPE::GetValue(A[i][j]); + const passivedouble q = (M[0][0] + M[1][1] + M[2][2]) / 3.0; + const passivedouble p1 = M[0][1] * M[0][1] + M[0][2] * M[0][2] + M[1][2] * M[1][2]; + const passivedouble p2 = pow(M[0][0] - q, 2) + pow(M[1][1] - q, 2) + pow(M[2][2] - q, 2) + 2.0 * p1; + const passivedouble p = sqrt(p2 / 6.0); + if (p <= 0.0) return (q > 0.0) ? 1.0 : big; + passivedouble B[3][3]; + for (auto i = 0u; i < 3; ++i) + for (auto j = 0u; j < 3; ++j) B[i][j] = (M[i][j] - ((i == j) ? q : 0.0)) / p; + const passivedouble detB = B[0][0] * (B[1][1] * B[2][2] - B[1][2] * B[2][1]) - + B[0][1] * (B[1][0] * B[2][2] - B[1][2] * B[2][0]) + + B[0][2] * (B[1][0] * B[2][1] - B[1][1] * B[2][0]); + const passivedouble phi = acos(max(-1.0, min(1.0, 0.5 * detB))) / 3.0; + lmax = q + 2.0 * p * cos(phi); + lmin = q + 2.0 * p * cos(phi + 2.0 * M_PI / 3.0); + } + return (lmin > 0.0) ? lmax / lmin : big; +} + +/*--- Inverse of a symmetric positive definite matrix by Cholesky factorization, which stays accurate + * for the strong anisotropy of a boundary layer stencil. False if a pivot vanishes. ---*/ +bool CholeskyInverse(unsigned short n, const su2double A[3][3], su2double inv[3][3]) { + su2double L[3][3] = {{0.0}}; + for (auto i = 0u; i < n; ++i) { + for (auto j = 0u; j <= i; ++j) { + su2double sum = A[i][j]; + for (auto k = 0u; k < j; ++k) sum -= L[i][k] * L[j][k]; + if (i == j) { + if (sum <= 1e-14 * A[i][i]) return false; + L[i][i] = sqrt(sum); + } else { + L[i][j] = sum / L[j][j]; + } + } + } + for (auto c = 0u; c < n; ++c) { + su2double y[3] = {0.0}, x[3] = {0.0}; + for (auto i = 0u; i < n; ++i) { + su2double sum = (i == c) ? 1.0 : 0.0; + for (auto k = 0u; k < i; ++k) sum -= L[i][k] * y[k]; + y[i] = sum / L[i][i]; + } + for (int i = n - 1; i >= 0; --i) { + su2double sum = y[i]; + for (auto k = i + 1u; k < n; ++k) sum -= L[k][i] * x[k]; + x[i] = sum / L[i][i]; + } + for (auto i = 0u; i < n; ++i) inv[i][c] = x[i]; + } + return true; +} +} // namespace + +void CMultiGridIntegration::ComputeProlongationGradient(CSolver *sol_coarse, CGeometry *geo_coarse, + CGeometry *geo_fine, const CConfig *config, + unsigned short iMesh, bool useSolution) { + SU2_ZONE_SCOPED + + const unsigned short nVar = sol_coarse->GetnVar(); + const unsigned short nDim = geo_coarse->GetnDim(); + const auto nPointDomain = geo_coarse->GetnPointDomain(); + auto* nodes = sol_coarse->GetNodes(); + auto Field = [&](unsigned long iPoint) { + return useSolution ? nodes->GetSolution(iPoint) : nodes->GetSolution_Old(iPoint); + }; + + /*--- Every thread of this rank takes this branch together, and the loops below would all be + * empty anyway. It also keeps the sizing test off a container that was never resized. ---*/ + /*--- The coarse grid does not change, so each level is reported once. ---*/ + auto reportIllConditioned = [&]() { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + if (!prolongReported[iMesh]) { + prolongReported[iMesh] = true; + unsigned long local[2] = {prolongIllConditioned[iMesh], nPointDomain}, global[2] = {0, 0}; + SU2_MPI::Allreduce(local, global, 2, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + if ((global[0] > 0) && (SU2_MPI::GetRank() == MASTER_NODE)) + cout << "MG level " << iMesh << ": linear prolongation keeps the constant operator in " << global[0] + << " of " << global[1] << " coarse cells (stencil direction condition number above " + << PROLONG_MAX_CONDITION << " or aspect ratio above " << PROLONG_MAX_ASPECT << ")." << endl; + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + }; + + if (nPointDomain == 0) { + reportIllConditioned(); + return; + } + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { prolongIllConditioned[iMesh] = 0; } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + auto& grad = prolongGradient[iMesh]; + if ((grad.length() != nPointDomain) || (grad.rows() != nVar) || (grad.cols() != nDim)) + grad.resize(nPointDomain, nVar, nDim, 0.0); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + auto& gradient = prolongGradient[iMesh]; + + /*--- Inverse-distance-squared weights give every neighbour's directional difference the same + * weight, so the fit does not depend on how far apart the neighbours are. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain, omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) = 0.0; + + const auto* Coord_i = geo_coarse->nodes->GetCoord(iPoint); + const auto* Correction_i = Field(iPoint); + + /*--- The extent matrix measures the stencil's stretching, the fit matrix is the weighted normal matrix. ---*/ + su2double extentMat[MAXNDIM][MAXNDIM] = {{0.0}}, fitMat[MAXNDIM][MAXNDIM] = {{0.0}}; + su2double rhs[MAXNVAR][MAXNDIM] = {{0.0}}; + + for (auto jPoint : geo_coarse->nodes->GetPoints(iPoint)) { + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_coarse->nodes->GetCoord(jPoint), Coord_i, Distance); + const su2double dist2 = GeometryToolbox::SquaredNorm(nDim, Distance); + if (dist2 <= 0.0) continue; + + for (auto iDim = 0u; iDim < nDim; iDim++) + for (auto jDim = 0u; jDim < nDim; jDim++) { + extentMat[iDim][jDim] += Distance[iDim]*Distance[jDim]; + fitMat[iDim][jDim] += Distance[iDim]*Distance[jDim]/dist2; + } + + const auto* Correction_j = Field(jPoint); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double delta = Correction_j[iVar] - Correction_i[iVar]; + for (auto iDim = 0u; iDim < nDim; iDim++) rhs[iVar][iDim] += Distance[iDim]*delta/dist2; + } + } + + /*--- Neighbours that do not span every direction leave a slope undefined. A stencil much longer in + * one direction than another, from cell stretching, weakens the gradient or removes it. In both + * cases the cell falls back to the constant operator. ---*/ + su2double fitInv[MAXNDIM][MAXNDIM] = {{0.0}}; + const passivedouble aspect = sqrt(FitConditionNumber(nDim, extentMat)); + const passivedouble weight = + (aspect <= PROLONG_FULL_ASPECT) + ? 1.0 + : max(0.0, log(PROLONG_MAX_ASPECT / aspect) / log(PROLONG_MAX_ASPECT / PROLONG_FULL_ASPECT)); + if ((FitConditionNumber(nDim, fitMat) > PROLONG_MAX_CONDITION) || (weight <= 0.0) || + !CholeskyInverse(nDim, fitMat, fitInv)) { + SU2_OMP_ATOMIC + prolongIllConditioned[iMesh]++; + continue; + } + + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) { + su2double value = 0.0; + for (auto jDim = 0u; jDim < nDim; jDim++) value += fitInv[iDim][jDim]*rhs[iVar][jDim]; + gradient(iPoint, iVar, iDim) = weight * value; + } + } + END_SU2_OMP_FOR + + /*--- The children of a wall agglomerate are all wall nodes, so their velocity correction is + * the constrained one the caller already imposed and only its gradient has to go. The + * remaining variables carry no wall condition and reconstruct along the surface. ---*/ + + const short iVel = nodes->GetVelocityIndex(); + + if (iVel >= 0) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + + const auto kindBC = config->GetMarker_All_KindBC(iMarker); + if (!config->GetViscous_Wall(iMarker) && (kindBC != EULER_WALL) && (kindBC != SYMMETRY_PLANE)) continue; + + SU2_OMP_FOR_STAT(32) + for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + const auto iPoint = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + if (iPoint >= nPointDomain) continue; + for (auto iVar = iVel; iVar < min(iVel + nDim, nVar); iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) = 0.0; + } + END_SU2_OMP_FOR + } + } + + /*--- Barth-Jespersen limiter over the children. One factor per variable scales the whole + * gradient, which is what keeps the children averaging back to the coarse value. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain, omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + + const auto* Coord_i = geo_coarse->nodes->GetCoord(iPoint); + const auto* Correction_i = Field(iPoint); + + su2double Correction_Min[MAXNVAR], Correction_Max[MAXNVAR]; + for (auto iVar = 0u; iVar < nVar; iVar++) { + Correction_Min[iVar] = Correction_i[iVar]; + Correction_Max[iVar] = Correction_i[iVar]; + } + + for (auto jPoint : geo_coarse->nodes->GetPoints(iPoint)) { + const auto* Correction_j = Field(jPoint); + for (auto iVar = 0u; iVar < nVar; iVar++) { + Correction_Min[iVar] = min(Correction_Min[iVar], Correction_j[iVar]); + Correction_Max[iVar] = max(Correction_Max[iVar], Correction_j[iVar]); + } + } + + su2double limiter[MAXNVAR]; + for (auto iVar = 0u; iVar < nVar; iVar++) limiter[iVar] = 1.0; + + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(iPoint); iChildren++) { + + const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(iPoint, iChildren); + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_fine->nodes->GetCoord(Point_Fine), Coord_i, Distance); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + + su2double delta = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) delta += gradient(iPoint, iVar, iDim)*Distance[iDim]; + + /*--- Scaled with the correction so the cut-off tracks its magnitude. ---*/ + const su2double tol = EPS*max(fabs(Correction_Max[iVar]), fabs(Correction_Min[iVar])) + EPS*EPS; + + if (delta > tol) + limiter[iVar] = min(limiter[iVar], (Correction_Max[iVar] - Correction_i[iVar])/delta); + else if (delta < -tol) + limiter[iVar] = min(limiter[iVar], (Correction_Min[iVar] - Correction_i[iVar])/delta); + } + } + + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double factor = max(su2double(0.0), min(su2double(1.0), limiter[iVar])); + for (auto iDim = 0u; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) *= factor; + } + } + END_SU2_OMP_FOR + + reportIllConditioned(); } void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *solver, CGeometry *geometry, @@ -964,6 +1336,11 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ const unsigned short nVar = solver->GetnVar(); + if (config->GetMGOptions().MG_Smooth_Output) { + const su2double res = sqrt(solver->LinSysRes.squaredNorm() / (nVar * geometry->GetGlobal_nPointDomain())); + SU2_OMP_SAFE_GLOBAL_ACCESS(lastCorrecSmoothRMS[iMesh][0] = SU2_TYPE::GetValue(res);) + } + 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); @@ -1041,7 +1418,61 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } } -void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeometry *geo_fine, +void CMultiGridIntegration::ComputeCorrectionWeights(CGeometry** geometry, const CEulerVariable* flowNodes, + unsigned short nLevels) { + /*--- Coarse CVs over which the correction fades in from a wall end. ---*/ + constexpr unsigned long nBlend = 4; + + correctionWeight.assign(nLevels, {}); + if (!flowNodes) return; + + /*--- Points of each level holding a wall end, passed down from the finest grid. ---*/ + vector wallEnd(geometry[MESH_0]->GetnPoint(), 0); + for (auto iPoint = 0ul; iPoint < geometry[MESH_0]->GetnPointDomain(); iPoint++) + wallEnd[iPoint] = flowNodes->SkipMGCorrection(iPoint); + + constexpr auto FAR = std::numeric_limits::max(); + vector hops, front, next; + + for (unsigned short iMesh = 0; iMesh < nLevels; iMesh++) { + const auto* fine = geometry[iMesh]; + const auto* coarse = geometry[iMesh + 1]; + + vector coarseWallEnd(coarse->GetnPoint(), 0); + hops.assign(coarse->GetnPoint(), FAR); + front.clear(); + for (auto iCoarse = 0ul; iCoarse < coarse->GetnPointDomain(); iCoarse++) + for (auto iChild = 0u; iChild < coarse->nodes->GetnChildren_CV(iCoarse); iChild++) + if (wallEnd[coarse->nodes->GetChildren_CV(iCoarse, iChild)]) { + coarseWallEnd[iCoarse] = 1; + hops[iCoarse] = 0; + front.push_back(iCoarse); + break; + } + + for (unsigned long d = 1; (d < nBlend) && !front.empty(); d++) { + next.clear(); + for (auto iCoarse : front) + for (auto jCoarse : coarse->nodes->GetPoints(iCoarse)) + if (hops[jCoarse] == FAR) { + hops[jCoarse] = d; + next.push_back(jCoarse); + } + front.swap(next); + } + + auto& weight = correctionWeight[iMesh]; + weight.assign(fine->GetnPoint(), 1.0); + for (auto iCoarse = 0ul; iCoarse < coarse->GetnPointDomain(); iCoarse++) { + if (hops[iCoarse] == FAR) continue; + for (auto iChild = 0u; iChild < coarse->nodes->GetnChildren_CV(iCoarse); iChild++) + weight[coarse->nodes->GetChildren_CV(iCoarse, iChild)] = passivedouble(hops[iCoarse]) / passivedouble(nBlend); + } + wallEnd.swap(coarseWallEnd); + } +} + +void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CSolver *sol_flow, CGeometry *geo_fine, CConfig *config, unsigned short iMesh) { SU2_ZONE_SCOPED @@ -1049,20 +1480,23 @@ 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(); + const auto* flowNodes = sol_flow ? dynamic_cast(sol_flow->GetNodes()) : nullptr; + const bool skipHere = flowNodes && (sol_fine == sol_flow); SU2_OMP_FOR_STAT(roundUpDiv(geo_fine->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Fine = 0ul; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) { + if (skipHere && flowNodes->SkipMGCorrection(Point_Fine)) continue; auto* Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); auto* Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine); - 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]; + /*--- Prevent a fine grid divergence due to a coarse grid divergence ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) + if (Residual_Fine[iVar] != Residual_Fine[iVar]) Residual_Fine[iVar] = 0.0; + + su2double relax = factor; + if ((iMesh < correctionWeight.size()) && !correctionWeight[iMesh].empty()) relax *= correctionWeight[iMesh][Point_Fine]; - Solution_Fine[iVar] += correction; - } + for (auto iVar = 0u; iVar < nVar; iVar++) Solution_Fine[iVar] += relax * Residual_Fine[iVar]; } END_SU2_OMP_FOR @@ -1073,24 +1507,55 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet } -void CMultiGridIntegration::SetProlongated_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { +void CMultiGridIntegration::InterpolateSolutionToChildren(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, + CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh) { SU2_ZONE_SCOPED - const unsigned short Solver_Position = config->GetContainerPosition(RunTime_EqSystem); - const bool grid_movement = config->GetGrid_Movement(); + const unsigned short nVar = sol_coarse->GetnVar(); + const unsigned short nDim = geo_coarse->GetnDim(); + const bool linear = config->GetMGOptions().MG_Linear_Prolongation; - /*--- Constant injection: every fine child of a coarse CV takes its parent's value. ---*/ + if (linear) ComputeProlongationGradient(sol_coarse, geo_coarse, geo_fine, config, iMesh, true); SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { const auto* Solution_Coarse = sol_coarse->GetNodes()->GetSolution(Point_Coarse); + const auto* Coord_Coarse = geo_coarse->nodes->GetCoord(Point_Coarse); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - sol_fine->GetNodes()->SetSolution(Point_Fine, Solution_Coarse); + + if (!linear) { + sol_fine->GetNodes()->SetSolution(Point_Fine, Solution_Coarse); + continue; + } + + su2double Distance[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geo_fine->nodes->GetCoord(Point_Fine), Coord_Coarse, Distance); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + su2double value = Solution_Coarse[iVar]; + for (auto iDim = 0u; iDim < nDim; iDim++) + value += prolongGradient[iMesh](Point_Coarse, iVar, iDim) * Distance[iDim]; + sol_fine->GetNodes()->SetSolution(Point_Fine, iVar, value); + } } } END_SU2_OMP_FOR +} + +void CMultiGridIntegration::SetProlongated_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh) { + SU2_ZONE_SCOPED + + const unsigned short Solver_Position = config->GetContainerPosition(RunTime_EqSystem); + const bool grid_movement = config->GetGrid_Movement(); + + /*--- Constant injection, or the limited linear reconstruction when MG_Linear_Prolongation is set. ---*/ + + InterpolateSolutionToChildren(sol_fine, sol_coarse, geo_fine, geo_coarse, config, iMesh); /*--- The injected values do not satisfy the fine-grid wall conditions on their own. ---*/ @@ -1140,26 +1605,28 @@ void CMultiGridIntegration::SetProlongated_Solution(unsigned short RunTime_EqSys void CMultiGridIntegration::SetProlongated_ScalarSolution(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, - CConfig *config, bool eddy_viscosity) { + CConfig *config, bool eddy_viscosity, + unsigned short iMesh) { SU2_ZONE_SCOPED - /*--- Constant injection: every fine child of a coarse CV takes its parent's value. ---*/ + /*--- Constant injection, or the limited linear reconstruction when MG_Linear_Prolongation is set. ---*/ - SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + InterpolateSolutionToChildren(sol_fine, sol_coarse, geo_fine, geo_coarse, config, iMesh); - const auto* Solution_Coarse = sol_coarse->GetNodes()->GetSolution(Point_Coarse); - const su2double muT_Coarse = eddy_viscosity ? sol_coarse->GetNodes()->GetmuT(Point_Coarse) : su2double(0.0); + /*--- The eddy viscosity is injected as a constant. ---*/ - for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - const auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - sol_fine->GetNodes()->SetSolution(Point_Fine, Solution_Coarse); - if (eddy_viscosity) sol_fine->GetNodes()->SetmuT(Point_Fine, muT_Coarse); + if (eddy_viscosity) { + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + const su2double muT_Coarse = sol_coarse->GetNodes()->GetmuT(Point_Coarse); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) + sol_fine->GetNodes()->SetmuT(geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren), muT_Coarse); } + END_SU2_OMP_FOR } - END_SU2_OMP_FOR - /*--- The eddy viscosity is zero on a no-slip wall, the injected value is not. ---*/ + /*--- On a no-slip wall the eddy viscosity is zero, and the turbulence variables take the value + * of their wall parent, which holds the wall condition instead of an extrapolation. ---*/ if (eddy_viscosity) { for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { @@ -1168,6 +1635,8 @@ void CMultiGridIntegration::SetProlongated_ScalarSolution(CSolver *sol_fine, CSo SU2_OMP_FOR_STAT(32) for (auto iVertex = 0ul; iVertex < geo_fine->nVertex[iMarker]; iVertex++) { const auto Point_Fine = geo_fine->vertex[iMarker][iVertex]->GetNode(); + const auto Point_Coarse = geo_fine->nodes->GetParent_CV(Point_Fine); + sol_fine->GetNodes()->SetSolution(Point_Fine, sol_coarse->GetNodes()->GetSolution(Point_Coarse)); sol_fine->GetNodes()->SetmuT(Point_Fine, 0.0); } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index 7dca6641bce0..1e712b5079b7 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -87,6 +87,13 @@ void CSingleGridIntegration::SingleGrid_Iteration(CGeometry ****geometry, CSolve Time_Integration(geometry_fine, solvers_fine, config[iZone], NO_RK_ITER, RunTime_EqSystem); + /*--- Coarse-level residuals are per-rank sums of squares unless reduced, which the smoothing + * early exit already does. ---*/ + + if ((FinestMesh != MESH_0) && !config[iZone]->GetMGOptions().MG_Smooth_EarlyExit) { + solvers_fine[Solver_Position]->SetResidual_RMS(geometry_fine, config[iZone], true); + } + /*--- Postprocessing ---*/ solvers_fine[Solver_Position]->Postprocessing(geometry_fine, solvers_fine, config[iZone], FinestMesh); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index b1c61ad377de..c04ed191cd73 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -302,6 +302,110 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, nodes->NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0; } + /*--- The gradient is singular where a no-slip wall meets another kind of boundary. Unlimited + * reconstruction across it gives a defect the first-order multigrid levels cannot correct. ---*/ + + if ((iMesh == MESH_0) && config->GetMUSCL_Flow() && config->GetMUSCL_FirstOrderWallEnds()) { + /*--- Only where the other boundary continues the wall surface does the condition jump along a + * smooth surface. Where they meet at an angle, as at a wing root on a symmetry plane, the + * flow is smooth. Halfway between parallel and perpendicular separates the two. ---*/ + const su2double minAlignment = sqrt(0.5); + + auto unitNormal = [&](unsigned short iMarker, long iVertex, su2double* normal) { + geometry->vertex[iMarker][iVertex]->GetNormal(normal); + const su2double area = GeometryToolbox::Norm(nDim, normal); + for (unsigned short iDim = 0; iDim < nDim; iDim++) normal[iDim] /= max(area, EPS); + }; + + /*--- A sharp edge of the wall itself, such as a trailing edge, is singular as well. Its vertex + * normal is much shorter than the area of the faces around it, cos(half the angle between + * their normals); the same halfway threshold applies. ---*/ + vector wallArea(nPoint, 0.0); + for (unsigned short iWall = 0; iWall < config->GetnMarker_All(); iWall++) { + if (!config->GetViscous_Wall(iWall)) continue; + for (unsigned long iElem = 0; iElem < geometry->GetnElem_Bound(iWall); iElem++) { + const auto* elem = geometry->bound[iWall][iElem]; + const auto nNodes = elem->GetnNodes(); + su2double area = 0.0; + if (nDim == 2) { + su2double edge[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, geometry->nodes->GetCoord(elem->GetNode(1)), + geometry->nodes->GetCoord(elem->GetNode(0)), edge); + area = GeometryToolbox::Norm(nDim, edge); + } else { + for (unsigned short iNode = 1; iNode + 1 < nNodes; iNode++) { + su2double a[3] = {0.0}, b[3] = {0.0}, c[3] = {0.0}; + GeometryToolbox::Distance(3, geometry->nodes->GetCoord(elem->GetNode(iNode)), + geometry->nodes->GetCoord(elem->GetNode(0)), a); + GeometryToolbox::Distance(3, geometry->nodes->GetCoord(elem->GetNode(iNode + 1)), + geometry->nodes->GetCoord(elem->GetNode(0)), b); + GeometryToolbox::CrossProduct(a, b, c); + area += 0.5 * GeometryToolbox::Norm(3, c); + } + } + for (unsigned short iNode = 0; iNode < nNodes; iNode++) wallArea[elem->GetNode(iNode)] += area / nNodes; + } + } + const su2double minFlatness = cos(PI_NUMBER / 8); + + for (iPoint = 0; iPoint < nPoint; iPoint++) { + bool wallEnds = false; + for (unsigned short iWall = 0; iWall < config->GetnMarker_All(); iWall++) { + const auto wallVertex = geometry->nodes->GetVertex(iPoint, iWall); + if (!config->GetViscous_Wall(iWall) || (wallVertex < 0)) continue; + + const su2double vertexArea = GeometryToolbox::Norm(nDim, geometry->vertex[iWall][wallVertex]->GetNormal()); + if (vertexArea < minFlatness * wallArea[iPoint]) wallEnds = true; + + su2double wallNormal[MAXNDIM] = {0.0}; + unitNormal(iWall, wallVertex, wallNormal); + + for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + const auto kindBC = config->GetMarker_All_KindBC(iMarker); + if (config->GetViscous_Wall(iMarker) || (kindBC == SEND_RECEIVE) || (kindBC == PERIODIC_BOUNDARY) || + (kindBC == INTERNAL_BOUNDARY) || (kindBC == NEARFIELD_BOUNDARY) || (kindBC == FLUID_INTERFACE)) continue; + const auto iVertex = geometry->nodes->GetVertex(iPoint, iMarker); + if (iVertex < 0) continue; + su2double normal[MAXNDIM] = {0.0}; + unitNormal(iMarker, iVertex, normal); + if (fabs(GeometryToolbox::DotProduct(nDim, wallNormal, normal)) > minAlignment) wallEnds = true; + } + } + nodes->FirstOrderReconstruction(iPoint) = wallEnds; + } + + unsigned long nFlagged = 0, nFlaggedGlobal = 0; + for (iPoint = 0; iPoint < nPointDomain; iPoint++) nFlagged += nodes->FirstOrderReconstruction(iPoint); + SU2_MPI::Allreduce(&nFlagged, &nFlaggedGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + if (rank == MASTER_NODE) + cout << "First-order reconstruction at " << nFlaggedGlobal << " points where a no-slip wall ends or has a sharp edge." << endl; + + /*--- A halo may lack one of the markers its owner has, so it takes the owner's flag. ---*/ + for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) || (config->GetMarker_All_SendRecv(iMarker) <= 0)) + continue; + const auto markerS = iMarker, markerR = static_cast(iMarker + 1); + const int sendTo = config->GetMarker_All_SendRecv(markerS) - 1; + const int receiveFrom = abs(config->GetMarker_All_SendRecv(markerR)) - 1; + vector bufSend(geometry->nVertex[markerS]), bufRecv(geometry->nVertex[markerR]); + for (auto iVertex = 0ul; iVertex < geometry->nVertex[markerS]; iVertex++) + bufSend[iVertex] = nodes->FirstOrderReconstruction(geometry->vertex[markerS][iVertex]->GetNode()); + SU2_MPI::Sendrecv(bufSend.data(), bufSend.size(), MPI_UNSIGNED_SHORT, sendTo, 0, bufRecv.data(), bufRecv.size(), + MPI_UNSIGNED_SHORT, receiveFrom, 0, SU2_MPI::GetComm(), MPI_STATUS_IGNORE); + for (auto iVertex = 0ul; iVertex < geometry->nVertex[markerR]; iVertex++) + nodes->FirstOrderReconstruction(geometry->vertex[markerR][iVertex]->GetNode()) = bufRecv[iVertex]; + } + + /*--- The coarse grid cannot represent the singular flow there, its correction is not prolongated to + * those points and their neighbours. ---*/ + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + bool skip = nodes->FirstOrderReconstruction(iPoint); + for (auto iNeigh = 0u; !skip && (iNeigh < geometry->nodes->GetnPoint(iPoint)); iNeigh++) + skip = nodes->FirstOrderReconstruction(geometry->nodes->GetPoint(iPoint, iNeigh)); + nodes->SkipMGCorrection(iPoint) = skip; + } + } + /*--- Check that the initial solution is physical, report any non-physical nodes ---*/ counter_local = 0; @@ -1909,7 +2013,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set them with or without high order reconstruction using MUSCL strategy. ---*/ - if (!muscl) { + if (!muscl || nodes->FirstOrderReconstruction(iPoint) || nodes->FirstOrderReconstruction(jPoint)) { numerics->SetPrimitive(V_i, V_j); numerics->SetSecondary(S_i, S_j); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 643cb9a4f6c9..f76c03be4ccd 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -100,9 +100,17 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C } } - /*--- Compute gradient of the primitive variables ---*/ - - if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + /*--- Compute gradient of the primitive variables. Coarse levels take the viscous flux from the edge + * difference only, as their Jacobian does. ---*/ + + if (iMesh != MESH_0) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + for (unsigned short iVar = 0; iVar < nPrimVarGrad; iVar++) + for (unsigned short iDim = 0; iDim < nDim; iDim++) nodes->GetGradient_Primitive()(iPoint, iVar, iDim) = 0.0; + END_SU2_OMP_FOR + } + else if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { SetPrimitive_Gradient_GG(geometry, config); } else if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 407e1c9bc7a3..a47a4f899cb3 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -2035,6 +2035,282 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, } } } + + /*--- The finest-grid nonlinear residuals drive the CFL, so the decision does not depend on how many + * ranks split the preconditioner. The finest grid decides for all levels. ---*/ + + if (!config->OptionIsSet("OUTLIER_MITIGATION_PARAM")) { + constexpr size_t trendWindow = 5, jumpWindow = 3, spikeWindow = 5, frozenWindow = 30; + constexpr passivedouble growBelow = 0.01, shrinkAbove = 0.03, jumpAbove = 0.75, spikeAbove = 1.0; + constexpr passivedouble frozenSpread = 0.05, frozenProgress = 0.02; + constexpr passivedouble linearBrake = 0.5, linearFailure = 0.8; + constexpr unsigned long quietIters = 20, eventGap = 20, bestWindow = 100, regressGap = 5, regressEpisode = 30; + constexpr unsigned long capProbeQuiet = 50; + constexpr passivedouble capProbeFactor = 1.3, eventCapFactor = 0.8, eventMatch = 1.5, cutFactor = 0.7; + constexpr size_t flipWindow = 10, stuckWindow = 15; + constexpr passivedouble regressAbove = 1.0, regressCap = 0.3, growBelowBest = 0.5; + constexpr passivedouble stuckSpread = 0.15, flipAmplitude = 0.03, flipProgress = 0.15; + constexpr size_t stallWindow = 100; + constexpr unsigned long stallPersist = 20; + constexpr passivedouble stallProgress = 0.1, stallGrowth = 1.1; + constexpr size_t slopeWindow = 60; + constexpr unsigned long trialHold = 300, trialQuiet = 100, trialSettle = 150; + constexpr passivedouble trialCut = 0.75, trialBudget = 30.0, fallingSlope = 0.3, trialResponse = 0.5; + constexpr passivedouble oscCap = 0.5, oscReset = 3.0, oscRelease = 1.2; + constexpr unsigned long oscQuiet = 100; + + if (iMesh == MESH_0) { + const passivedouble avgCFL = SU2_TYPE::GetValue(Avg_CFL_Local); + + /*--- Short-window trend of the mean log residual. ---*/ + ResTrend_History.push_back(SU2_TYPE::GetValue(New_Func)); + if (ResTrend_History.size() > trendWindow + 1) ResTrend_History.erase(ResTrend_History.begin()); + const bool full = (ResTrend_History.size() > trendWindow); + const passivedouble trend = full ? (ResTrend_History.back() - ResTrend_History.front()) / trendWindow : 0.0; + ResTrend_Increase = full && (trend < growBelow); + ResTrend_Reduce = full && (trend > shrinkAbove); + + /*--- A divergence raises all flow residuals together, a transient usually only one, so the jump + * is the median rise over the flow variables. ---*/ + vector logRes(solverFlow->GetnVar()); + for (auto iVar = 0u; iVar < solverFlow->GetnVar(); iVar++) + logRes[iVar] = log10(SU2_TYPE::GetValue(solverFlow->GetRes_RMS(iVar))); + ResJump_History.push_back(logRes); + if (ResJump_History.size() > spikeWindow + 1) ResJump_History.erase(ResJump_History.begin()); + bool jump = false; + if (ResJump_History.size() > jumpWindow) { + const auto& past = ResJump_History[ResJump_History.size() - 1 - jumpWindow]; + vector rise(logRes.size()); + for (auto iVar = 0u; iVar < logRes.size(); iVar++) rise[iVar] = logRes[iVar] - past[iVar]; + sort(rise.begin(), rise.end()); + const auto n = rise.size(); + jump = (0.5 * (rise[(n - 1) / 2] + rise[n / 2]) > jumpAbove); + } + /*--- A spike in a single equation, such as the energy at a stagnation corner, is the first sign of + * a CFL that is too high even when the other residuals do not move. ---*/ + if (ResJump_History.size() > spikeWindow) { + for (auto iVar = 0u; iVar < logRes.size(); iVar++) + jump = jump || (logRes[iVar] - ResJump_History.front()[iVar] > spikeAbove); + } + + /*--- Flow residuals that all stop falling are frozen. The turbulence residual may still fall, so it + * is left out. The CFL cannot go below its lower limit, so nothing is done there. ---*/ + const passivedouble cflMin = SU2_TYPE::GetValue(CFLMin); + ResFlow_History.push_back(logRes); + if (ResFlow_History.size() > frozenWindow + 1) ResFlow_History.erase(ResFlow_History.begin()); + bool frozen = (ResFlow_History.size() > frozenWindow) && (iter >= startingIter + frozenWindow) && + (iter > CFL_LastEvent + frozenWindow) && (iter > CFL_LastStall + frozenWindow) && + (avgCFL > 1.05 * cflMin); + for (auto iVar = 0u; frozen && (iVar < logRes.size()); iVar++) { + passivedouble lowest = ResFlow_History.front()[iVar], highest = lowest; + for (const auto& past : ResFlow_History) { + lowest = min(lowest, past[iVar]); + highest = max(highest, past[iVar]); + } + frozen = (ResFlow_History.front()[iVar] - ResFlow_History.back()[iVar] < frozenProgress) && + (highest - lowest < frozenSpread); + } + + /*--- A divergence that grows slowly never looks like a jump from one iteration to the next, so the + * flow residual is also compared with its best value of the last bestWindow iterations. A residual + * far above it that no longer moves is stuck at a CFL that is too high. ---*/ + passivedouble flowMean = 0.0; + for (const auto value : logRes) flowMean += value / logRes.size(); + if (iter >= startingIter) { + ResMean_History.push_back(flowMean); + if (ResMean_History.size() > bestWindow + 1) ResMean_History.erase(ResMean_History.begin()); + } + passivedouble best = flowMean; + for (const auto value : ResMean_History) best = min(best, value); + ResRise = (iter >= startingIter) ? flowMean - best : 0.0; + bool regress = (ResRise > regressAbove) && (iter >= CFL_LastEvent + regressGap); + if (!regress && (ResMean_History.size() >= stuckWindow) && (ResRise > regressAbove) && + (iter >= CFL_LastEvent + stuckWindow)) { + passivedouble lowest = ResMean_History.back(), highest = lowest; + for (auto k = ResMean_History.size() - stuckWindow; k < ResMean_History.size(); ++k) { + lowest = min(lowest, ResMean_History[k]); + highest = max(highest, ResMean_History[k]); + } + regress = (highest - lowest < stuckSpread); + } + + /*--- Most flow residuals that alternate up and down at every iteration without net progress are a + * flip-flop, which a CFL that is too high produces. ---*/ + bool flip = (ResFlow_History.size() > flipWindow) && (iter >= startingIter + flipWindow) && + (iter >= CFL_LastEvent + flipWindow) && (iter >= CFL_LastFlip + flipWindow); + if (flip) { + const auto nHist = ResFlow_History.size(); + unsigned short nAlternating = 0; + passivedouble progress = 0.0; + for (auto iVar = 0u; iVar < logRes.size(); iVar++) { + unsigned short nChanges = 0; + passivedouble amplitude = 0.0; + for (auto k = nHist - flipWindow; k < nHist; ++k) { + const auto d1 = ResFlow_History[k][iVar] - ResFlow_History[k - 1][iVar]; + amplitude += fabs(d1) / flipWindow; + if (k > nHist - flipWindow) { + const auto d0 = ResFlow_History[k - 1][iVar] - ResFlow_History[k - 2][iVar]; + if (d0 * d1 < 0.0) nChanges++; + } + } + if ((nChanges + 1 >= flipWindow - 1) && (amplitude > flipAmplitude)) nAlternating++; + progress += (ResFlow_History[nHist - 1 - flipWindow][iVar] - ResFlow_History.back()[iVar]) / logRes.size(); + } + flip = (nAlternating * 4 >= logRes.size() * 3) && (progress < flipProgress); + } + if (flip) CFL_LastFlip = iter; + + /*--- Residuals that stop falling while the CFL no longer grows are stalled by a CFL that is too + * high, whatever the shape of their oscillation. The upper envelope of each residual over a + * window is compared with that of the window before. ---*/ + ResStall_History.push_back(logRes); + CFLStall_History.push_back(avgCFL); + if (ResStall_History.size() > 2 * stallWindow) { + ResStall_History.erase(ResStall_History.begin()); + CFLStall_History.erase(CFLStall_History.begin()); + } + bool stall = false; + if ((ResStall_History.size() == 2 * stallWindow) && (iter >= startingIter + 2 * stallWindow) && + (iter >= max(CFL_LastEvent, CFL_LastStall) + stallWindow) && (avgCFL > 1.05 * cflMin)) { + passivedouble progress = 0.0; + for (auto iVar = 0u; iVar < logRes.size(); iVar++) { + passivedouble envOld = ResStall_History.front()[iVar], envNew = ResStall_History.back()[iVar]; + for (size_t k = 0; k < stallWindow; k++) { + envOld = max(envOld, ResStall_History[k][iVar]); + envNew = max(envNew, ResStall_History[k + stallWindow][iVar]); + } + progress += (envOld - envNew) / logRes.size(); + } + const bool cflFlat = (avgCFL < stallGrowth * CFLStall_History[stallWindow - 1]); + Stall_Streak = ((progress < stallProgress) && cflFlat) ? Stall_Streak + 1 : 0; + stall = (Stall_Streak >= stallPersist); + } else { + Stall_Streak = 0; + } + /*--- Cut a plateau repeatedly unless the residual is falling, and undo the cuts if it does not respond. ---*/ + const auto nMean = ResMean_History.size(); + auto meanSlope = [&](size_t window) { /* Least-squares slope in decades per 100 iterations. */ + passivedouble sumX = 0.0, sumY = 0.0, sumXY = 0.0, sumXX = 0.0; + for (size_t k = 0; k <= window; k++) { + const passivedouble y = ResMean_History[nMean - 1 - window + k]; + sumX += k; sumY += y; sumXY += k * y; sumXX += passivedouble(k) * k; + } + const passivedouble n = window + 1; + return 100.0 * (n * sumXY - sumX * sumY) / (n * sumXX - sumX * sumX); + }; + if (Trial_Active && (CFL_LastEvent > Trial_Start)) Trial_Active = false; + if (Trial_Active && (iter >= CFL_LastStall + trialQuiet)) { + Trial_Active = false; + if (rank == MASTER_NODE) cout << "CFL_ADAPT: plateau ended, CFL cuts kept." << endl; + } + if ((frozen || stall) && (iter >= Trial_Hold) && (nMean > slopeWindow)) { + const passivedouble slope = meanSlope(slopeWindow); + CFL_LastStall = iter; + Stall_Streak = 0; + if (Trial_Active && (avgCFL < Trial_CFL / trialBudget) && (Trial_Settle == 0)) { + Trial_Settle = iter + trialSettle; + } else if (Trial_Active && (Trial_Settle > 0) && (iter < Trial_Settle)) { + /* The residual is given time to respond at the lowest CFL. */ + } else if (Trial_Active && (Trial_Settle > 0) && (Trial_Res - flowMean > trialResponse)) { + Trial_Active = false; + if (rank == MASTER_NODE) cout << "CFL_ADAPT: residuals responded to the CFL cuts, they are kept." << endl; + } else if (Trial_Active && (Trial_Settle > 0)) { + CFL_Cap = Trial_Cap; + CFL_Ceiling = Trial_Ceiling; + CFL_StallCap = Trial_StallCap; + Trial_Active = false; + Trial_Hold = iter + trialHold; + if (rank == MASTER_NODE) + cout << "CFL_ADAPT: residuals did not respond to the CFL cuts since iteration " << Trial_Start + << ", CFL cuts undone." << endl; + } else if (Trial_Active || (slope > -fallingSlope)) { + if (!Trial_Active) { + Trial_Active = true; + Trial_Start = iter; + Trial_CFL = avgCFL; + Trial_Res = flowMean; + Trial_Settle = 0; + Trial_Ceiling = CFL_Ceiling; + Trial_Cap = CFL_Cap; + Trial_StallCap = CFL_StallCap; + } + CFL_Cap = max(trialCut * ((CFL_Cap > 0.0) ? min(CFL_Cap, avgCFL) : avgCFL), cflMin); + CFL_Ceiling = CFL_Cap; + CFL_StallCap = CFL_Cap; + if (rank == MASTER_NODE) + cout << "CFL_ADAPT: residuals " << (frozen ? "frozen" : "stalled") << " at iteration " << iter + << ", CFL capped at " << CFL_Cap << "." << endl; + } + } + + /*--- Halve the CFL on a jump or a failed linear solve and keep it below a ceiling. A second event + * near the CFL of an earlier one sets a lasting cap below both. ---*/ + CFL_Cut = (jump || regress || flip || (linRes > linearFailure)) && (iter >= CFL_LastEvent + jumpWindow) && + (avgCFL > 1.05 * cflMin); + if (CFL_Cut) { + /*--- Cuts that follow each other closely belong to one event, only its first CFL is recorded. ---*/ + if (iter >= CFL_LastEvent + eventGap) { + for (const auto eventCFL : CFL_EventValues) + if ((avgCFL < eventMatch * eventCFL) && (eventCFL < eventMatch * avgCFL)) { + const passivedouble cap = max(eventCapFactor * min(eventCFL, avgCFL), cflMin); + CFL_Cap = (CFL_Cap > 0.0) ? min(CFL_Cap, cap) : cap; + } + CFL_EventValues.push_back(avgCFL); + } + if (regress && (iter >= CFL_LastRegress + regressEpisode)) { + const passivedouble cap = max(regressCap * avgCFL, cflMin); + CFL_Cap = (CFL_Cap > 0.0) ? min(CFL_Cap, cap) : cap; + } + if (regress) CFL_LastRegress = iter; + /*--- Flip-flops bound the CFL below the lowest level that produced one. ---*/ + if (flip) { + CFL_OscLevel = (CFL_OscLevel > 0.0) ? min(CFL_OscLevel, avgCFL) : avgCFL; + Res_AtOsc = flowMean; + const passivedouble cap = max(oscCap * CFL_OscLevel, cflMin); + CFL_Cap = (CFL_Cap > 0.0) ? min(CFL_Cap, cap) : cap; + } + CFL_Ceiling = (CFL_Cap > 0.0) ? min(cutFactor * avgCFL, CFL_Cap) : cutFactor * avgCFL; + CFL_LastEvent = iter; + if (rank == MASTER_NODE) + cout << "CFL_ADAPT: " + << (jump ? "residual jump" : regress ? "residual regression" : flip ? "residual flip-flop" : "linear solver failure") + << " at iteration " << iter << ", CFL halved from " << avgCFL << "." << endl; + } + + /*--- A cap that has held the CFL through a quiet, converging stretch is raised, so that it does not + * outlive the transient that caused it. A new event at the higher CFL brings it back down. It + * stays below the cap set at the last stall and below the flip-flop limit, which itself rises + * in small steps over quiet, converging stretches. ---*/ + if ((CFL_OscLevel > 0.0) && (flowMean < Res_AtOsc - oscReset)) CFL_OscLevel = -1.0; + if ((CFL_OscLevel > 0.0) && (iter > max(CFL_LastEvent, CFL_LastOscRelease) + oscQuiet) && ResTrend_Increase && + (ResRise < growBelowBest)) { + CFL_OscLevel *= oscRelease; + CFL_LastOscRelease = iter; + } + passivedouble probeLimit = (CFL_StallCap > 0.0) ? CFL_StallCap : 1e30; + if (CFL_OscLevel > 0.0) probeLimit = min(probeLimit, max(oscCap * CFL_OscLevel, cflMin)); + if ((CFL_Cap > 0.0) && (avgCFL > 0.95 * CFL_Cap) && (iter > max(CFL_LastEvent, CFL_LastStall) + capProbeQuiet) && + (iter > CFL_LastCapProbe + capProbeQuiet) && !Trial_Active && ResTrend_Increase && (ResRise < growBelowBest) && + (CFL_Cap < probeLimit)) { + CFL_Cap = min(CFL_Cap * capProbeFactor, probeLimit); + CFL_LastCapProbe = iter; + CFL_EventValues.clear(); + } + + /*--- Lift the ceiling again while the residuals fall, up to the cap. ---*/ + if ((CFL_Ceiling > 0.0) && (iter > CFL_LastEvent + quietIters) && !Trial_Active && ResTrend_Increase && (ResRise < growBelowBest)) { + CFL_Ceiling *= 1.05; + if (CFL_Cap > 0.0) CFL_Ceiling = min(CFL_Ceiling, CFL_Cap); + else if (CFL_Ceiling > SU2_TYPE::GetValue(CFLMax)) CFL_Ceiling = -1.0; + } + } + + /*--- These replace the linear-residual and oscillation rules above. A poor linear solve holds the + * CFL back gently, a failed one halves it above. ---*/ + canIncrease = ResTrend_Increase && !Trial_Active && (linRes < linearBrake) && (iter >= startingIter) && (ResRise < growBelowBest); + reduceCFL = (ResTrend_Reduce || (linRes > linearBrake)) && (iter >= startingIter); + resetCFL = false; + } } /* End safe global access, now all threads update the CFL number. */ END_SU2_OMP_SAFE_GLOBAL_ACCESS @@ -2086,10 +2362,10 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, /* Check if we are hitting the min or max and adjust. */ - if (CFL*CFLFactor <= CFLMin) { - CFL = CFLMin; + if (CFL * CFLFactor <= CFLMin) { + CFL = CFLMin; CFLFactor = MGFactor[iMesh]; - } else if (CFL*CFLFactor >= CFLMax) { + } else if (CFL * CFLFactor >= CFLMax) { CFL = CFLMax; CFLFactor = MGFactor[iMesh]; } @@ -2104,7 +2380,10 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, /* Apply the adjustment to the CFL and store local values. */ + if (CFL_Cut) CFLFactor = 0.5; CFL *= CFLFactor; + if (CFL_Cut) CFL = max(CFL, CFLMin * MGFactor[iMesh]); + if ((CFL_Ceiling > 0.0) && (CFL > CFL_Ceiling)) CFL = max(su2double(CFL_Ceiling), CFLMin); solverFlow->GetNodes()->SetLocalCFL(iPoint, CFL); if ((iMesh == MESH_0) && solverTurb) { solverTurb->GetNodes()->SetLocalCFL(iPoint, CFL * CFLTurbReduction); diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index f45c60060b81..19ecaba04314 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -107,6 +107,8 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 } OutlierMitigation.resize(nPoint) = 0; + FirstOrderReconstruction.resize(nPoint) = 0; + SkipMGCorrection.resize(nPoint) = 0; } bool CEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { diff --git a/config_template.cfg b/config_template.cfg index 0aa29398b90a..a4ff3cf86fa3 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1579,6 +1579,19 @@ NUM_METHOD_GRAD_RECON = LEAST_SQUARES CFL_NUMBER= 15.0 % % Adaptive CFL number (NO, YES) +% The CFL is adapted on the finest-grid nonlinear residuals, which makes it independent of the number +% of MPI ranks. It grows while the residuals do not rise and shrinks when they rise. The CFL is halved +% and kept below a ceiling that lifts while the residuals fall on: a sudden residual jump (also of a +% single flow equation), a flow residual more than one order of magnitude above its best value of the +% last 100 iterations (also when it stays there), most flow residuals alternating up and down at every +% iteration without progress (flip-flop), or a failed linear solve. The CFL does not grow while the +% residual is far above its best value. Repeated events, a regression, a flip-flop, or all flow +% residuals stalling for 30 iterations cap the CFL lasting below the level where that happened. Flow +% residuals whose upper envelope stops falling over 100 iterations while the CFL no longer grows cap +% the CFL at half its value, a quarter on a second such stall, an eighth on a third. A flip-flop caps +% the CFL at half the lowest CFL where one occurred, until the residuals fall three orders of magnitude. Uses +% factor up/down and the CFL limits of CFL_ADAPT_PARAM, and never goes below the lower limit. With +% OUTLIER_MITIGATION_PARAM the CFL is adapted on the linear solver convergence instead. CFL_ADAPT= NO % % Parameters of the adaptive CFL number (factor-down, factor-up, CFL min value, @@ -1631,6 +1644,13 @@ CUSTOM_OBJFUNC= 'DRAG + 10 * pow(fmax(0.4-LIFT, 0), 2)' % Required for 2nd order upwind schemes (NO, YES) MUSCL_FLOW= YES % +% Reconstruct to first order on the edges of points where a no-slip wall ends on a boundary that +% continues its surface (a leading edge after a symmetry plane) or has a sharp edge (a sharp +% trailing edge). The gradient is singular there, and MUSCL across it can stall or destabilize +% the solver, multigrid in particular. With multigrid, the coarse-grid correction of the flow +% solution is also not prolongated to these points and their neighbours (NO, YES) +MUSCL_FIRST_ORDER_WALL_ENDS= NO +% % Option to ramp the MUSCL reconstruction scheme gradually from first to second order RAMP_MUSCL= NO % @@ -1798,8 +1818,8 @@ MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) % Max. iterations for multi-grid post-smoothing level MG_POST_SMOOTH= ( 2, 2, 2, 2 ) % -% Jacobi implicit smoothing of the correction -MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +% Jacobi implicit smoothing of the correction (default 0) +MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 ) % % Damping factor for the residual restriction MG_DAMP_RESTRICTION= 0.5 @@ -1807,6 +1827,13 @@ MG_DAMP_RESTRICTION= 0.5 % Damping factor for the correction prolongation MG_DAMP_PROLONGATION= 0.5 % +% Prolong the correction with a limited least-squares gradient over the coarse control +% volume instead of injecting the parent value into every child. The reconstruction is +% exact for a linear correction and the children still average back to the coarse value. +% Control volumes on walls and symmetry planes, and those where the gradient fit has a condition number +% above 1e5 (strongly stretched agglomerates, reported on screen), keep the constant operator (NO, YES) +MG_LINEAR_PROLONGATION= NO +% % Enable early exit from multigrid smoothing when the residual RMS drops % below MG_SMOOTH_RES_THRESHOLD * initial_rms (NO, YES) MG_SMOOTH_EARLY_EXIT= YES @@ -1857,8 +1884,9 @@ MG_STARTUP_ITER= 100 MG_STARTUP_CONVERGENCE= -2 % % Full-MG promotion on stagnation: if the active level's residual ratio between -% successive iterations exceeds this for MG_STARTUP_STAGNATION_ITER consecutive -% iterations, promote to the next finer level early. 0 disables (default 0.99). +% successive iterations stays between this value and its inverse for +% MG_STARTUP_STAGNATION_ITER consecutive iterations, promote to the next finer level +% early. A growing residual does not count as stagnation. 0 disables (default 0.99). MG_STARTUP_STAGNATION= 0.99 % % Consecutive stalled iterations required before Full-MG promotes on stagnation. @@ -2559,6 +2587,13 @@ PARMETIS_TOLERANCE= 0.02 PARMETIS_EDGE_WEIGHT= 1 PARMETIS_POINT_WEIGHT= 0 % +% Make cutting a short edge of an anisotropic cell expensive, so that partition +% boundaries cross boundary layers in the wall-normal direction instead of running +% through them. Each graph edge is weighted by the longest node distance of its cells +% over its own length, and lines of strongly coupled points (wall-normal lines in a +% boundary layer) are kept whole on one rank (NO, YES) (default NO). +PARMETIS_ANISOTROPY_WEIGHT= NO +% % ----------------------- SOBOLEV GRADIENT SMOOTHING OPTIONS ----------------------% % % Activate the gradient smoothing solver for the discrete adjoint driver (NO, YES) diff --git a/preconfigure.py b/preconfigure.py index bfba6995c5bb..77772d1f1b3d 100755 --- a/preconfigure.py +++ b/preconfigure.py @@ -73,6 +73,11 @@ def build_ninja(): # Specify C++17 env = os.environ.copy() env["CXXFLAGS"] = env.get("CXXFLAGS", "") + " -std=c++17" + # ninja is a build tool and must not be instrumented by sanitizers + for var in ("CFLAGS", "CXXFLAGS", "LDFLAGS"): + env[var] = " ".join( + f for f in env.get(var, "").split() if not f.startswith("-fsanitize") + ) subprocess.run( ["python3", "configure.py", "--bootstrap"], cwd=ninjapath, env=env )