multigrid linear prolongation with tigher full multigrid turbulence integration

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;

Common/include/option_structure.hpp

@@ -80,6 +80,7 @@ const unsigned int MAX_NUMBER_FFD = 15;       /*!< \brief Maximum number of FFDB
 enum: unsigned int{MAX_SOLS = 14};            /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */
 const unsigned int MAX_TERMS = 7;             /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */
 const unsigned int MAX_ZONES = 3;             /*!< \brief Maximum number of zones. */
+const unsigned short MAX_MGLEVELS = 10;       /*!< \brief Maximum number of coarse multigrid levels, which bounds the per-level arrays of the multigrid integration. */
 const unsigned int MAX_FE_KINDS = 7;          /*!< \brief Maximum number of Finite Elements. */
 const unsigned int NO_RK_ITER = 0;            /*!< \brief No Runge-Kutta iteration. */
 
@@ -1234,6 +1235,8 @@ struct CMGOptions {
   bool MG_Smooth_Output{false};           /*!< \brief Output compact per-cycle smoothing summary. */
   su2double MG_Smooth_StagnationTol{0.0}; /*!< \brief Stagnation early exit: stop if current_rms >= prev_rms * tol. 0 = disabled. */
   bool MG_Implicit_Lines{false};          /*!< \brief Enable implicit-lines agglomeration from walls. */
+  bool MG_Linear_Prolongation{false};     /*!< \brief Prolong the correction with a limited least-squares
+                                                 gradient instead of piecewise-constant injection. */
   unsigned long MG_Startup_Iter{100};     /*!< \brief Iterations per mesh during FMG startup, and the length of each level's CFL ramp. 0 = no iteration budget. */
   su2double MG_Startup_Convergence{-2.0}; /*!< \brief FMG: orders of magnitude (log10) that CONV_FIELD must drop on the
                                                  active level before promoting to the next finer one. Negative is a
@@ -2875,6 +2878,7 @@ enum class MPI_QUANTITIES {
   COORDINATES          ,  /*!< \brief Vertex coordinates communication. */
   COORDINATES_OLD      ,  /*!< \brief Old vertex coordinates communication. */
   MAX_LENGTH           ,  /*!< \brief Maximum length communication. */
+  WALL_DISTANCE        ,  /*!< \brief Wall distance and roughness of the nearest wall communication. */
   GRID_VELOCITY        ,  /*!< \brief Grid velocity communication. */
   SOLUTION_EDDY        ,  /*!< \brief Turbulent solution plus eddy viscosity communication. */
   STOCH_SOURCE_LANG    ,  /*!< \brief Stochastic source term for Langevin equations communication. */

Common/src/CConfig.cpp

@@ -2104,6 +2104,10 @@ void CConfig::SetConfig_Options() {
   /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Pave the coarse grid with advancing fronts raised from boundaries
    * that carry a stretched layer normal to themselves. DEFAULT: NO \ingroup Config*/
   addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false);
+  /*!\brief MG_LINEAR_PROLONGATION\n DESCRIPTION: Prolong the multigrid correction with a limited least-squares
+   * gradient over the coarse control volume instead of injecting the parent value into every child. Coarse CVs on
+   * walls and symmetry planes keep the constant operator. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_LINEAR_PROLONGATION", MGOptions.MG_Linear_Prolongation, false);
   /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Max number of iterations spent on each mesh during the Full
    * Multigrid (FMG) startup phase. DEFAULT: 100 \ingroup Config*/
   addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100);
@@ -2112,8 +2116,9 @@ void CConfig::SetConfig_Options() {
    * DEFAULT: -2 \ingroup Config*/
   addDoubleOption("MG_STARTUP_CONVERGENCE", MGOptions.MG_Startup_Convergence, -2.0);
   /*!\brief MG_STARTUP_STAGNATION\n DESCRIPTION: Full-MG promotion on stagnation. If the active level's residual ratio
-   * between successive iterations exceeds this value for MG_STARTUP_STAGNATION_ITER consecutive iterations, promote to
-   * the next finer level without waiting out MG_STARTUP_ITER. 0 disables it. DEFAULT: 0.99 \ingroup Config*/
+   * between successive iterations stays between this value and its inverse for MG_STARTUP_STAGNATION_ITER consecutive
+   * iterations, promote to the next finer level without waiting out MG_STARTUP_ITER. 0 disables it.
+   * DEFAULT: 0.99 \ingroup Config*/
   addDoubleOption("MG_STARTUP_STAGNATION", MGOptions.MG_Startup_Stagnation, 0.99);
   /*!\brief MG_STARTUP_STAGNATION_ITER\n DESCRIPTION: Consecutive stalled iterations required before Full-MG promotes
    * on stagnation. 0 disables it, as MG_STARTUP_STAGNATION= 0 does. DEFAULT: 5 \ingroup Config*/
@@ -4979,6 +4984,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   Kappa_2nd_AdjFlow = jst_adj_coeff[0];
   Kappa_4th_AdjFlow = jst_adj_coeff[1];
 
+  /*--- The multigrid integration carries per-level arrays of this size. ---*/
+
+  if (nMGLevels > MAX_MGLEVELS) {
+    SU2_MPI::Error("MGLEVEL is larger than the supported maximum of " + std::to_string(MAX_MGLEVELS) + ".",
+                   CURRENT_FUNCTION);
+  }
+
   /*--- Fill MG smooth vectors to size nMGLevels+1.
    Use parsed values (truncating or extending by repeat) or defaults if not set. ---*/
 

Common/src/geometry/CGeometry.cpp

@@ -627,6 +627,10 @@ void CGeometry::GetCommCountAndType(const CConfig* config, MPI_QUANTITIES commTy
       COUNT_PER_POINT = 1;
       MPI_TYPE = COMM_TYPE::DOUBLE;
       break;
+    case MPI_QUANTITIES::WALL_DISTANCE:
+      COUNT_PER_POINT = 2;
+      MPI_TYPE = COMM_TYPE::DOUBLE;
+      break;
     case MPI_QUANTITIES::NEIGHBORS:
       COUNT_PER_POINT = 1;
       MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT;
@@ -718,6 +722,10 @@ void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, MPI_QU
         case MPI_QUANTITIES::MAX_LENGTH:
           bufDSend[buf_offset] = nodes->GetMaxLength(iPoint);
           break;
+        case MPI_QUANTITIES::WALL_DISTANCE:
+          bufDSend[buf_offset] = nodes->GetWall_Distance(iPoint);
+          bufDSend[buf_offset + 1] = nodes->GetRoughnessHeight(iPoint);
+          break;
         case MPI_QUANTITIES::NEIGHBORS:
           bufSSend[buf_offset] = geometry->nodes->GetnNeighbor(iPoint);
           break;
@@ -811,6 +819,10 @@ void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, MPI_QU
         case MPI_QUANTITIES::MAX_LENGTH:
           nodes->SetMaxLength(iPoint, bufDRecv[buf_offset]);
           break;
+        case MPI_QUANTITIES::WALL_DISTANCE:
+          nodes->SetWall_Distance(iPoint, bufDRecv[buf_offset]);
+          nodes->SetRoughnessHeight(iPoint, bufDRecv[buf_offset + 1]);
+          break;
         case MPI_QUANTITIES::NEIGHBORS:
           nodes->SetnNeighbor(iPoint, bufSRecv[buf_offset]);
           break;
@@ -2741,6 +2753,16 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config)
     geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE);
     geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE);
     geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]);
+
+    /*--- SetCoord centred a halo agglomerate on the partial child list this rank holds. Take
+     the owner's coordinate, otherwise coarse stencils depend on the partitioning. Deformation
+     runs this inside a parallel region, so only one thread may reach the exchange. ---*/
+
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+      geometry_container[iMesh]->InitiateComms(geometry_container[iMesh], config, MPI_QUANTITIES::COORDINATES);
+      geometry_container[iMesh]->CompleteComms(geometry_container[iMesh], config, MPI_QUANTITIES::COORDINATES);
+    }
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
   /*--- Compute the global surface areas for all markers. ---*/
@@ -4538,6 +4560,47 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co
   const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
   return 2 * Vol / GeometryToolbox::Norm(3, Normal);
 }
+
+/*--- Volume-averages the wall distance and roughness of the children onto a coarse grid, then sets the
+ *    nearest-neighbor distance of its viscous wall vertices. ---*/
+void RestrictWallDistance(const CGeometry* geo_fine, CGeometry* geo_coarse, const CConfig* config) {
+  SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads()))
+  for (auto iPoint = 0ul; iPoint < geo_coarse->GetnPointDomain(); iPoint++) {
+    su2double dist = 0.0, roughness = 0.0, vol = 0.0;
+    for (auto iChild = 0u; iChild < geo_coarse->nodes->GetnChildren_CV(iPoint); iChild++) {
+      const auto jPoint = geo_coarse->nodes->GetChildren_CV(iPoint, iChild);
+      const su2double volChild = geo_fine->nodes->GetVolume(jPoint);
+      dist += geo_fine->nodes->GetWall_Distance(jPoint) * volChild;
+      roughness += geo_fine->nodes->GetRoughnessHeight(jPoint) * volChild;
+      vol += volChild;
+    }
+    geo_coarse->nodes->SetWall_Distance(iPoint, (vol > 0.0) ? dist / vol : 0.0);
+    geo_coarse->nodes->SetRoughnessHeight(iPoint, (vol > 0.0) ? roughness / vol : 0.0);
+  }
+  END_SU2_OMP_FOR
+
+  /*--- A halo agglomerate only holds the children on this rank, so take the owner's values. ---*/
+
+  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+    geo_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::WALL_DISTANCE);
+    geo_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::WALL_DISTANCE);
+  }
+  END_SU2_OMP_SAFE_GLOBAL_ACCESS
+
+  for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) {
+    const auto viscous = config->GetViscous_Wall(iMarker);
+
+    SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
+    for (auto iVertex = 0u; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) {
+      const auto iPoint = geo_coarse->vertex[iMarker][iVertex]->GetNode();
+      const su2double dist = (viscous && geo_coarse->nodes->GetDomain(iPoint))
+                                 ? NearestNeighborDistance(geo_coarse, config, iPoint)
+                                 : geo_coarse->nodes->GetWall_Distance(iPoint);
+      geo_coarse->vertex[iMarker][iVertex]->SetNearestNeighborDistance(dist);
+    }
+    END_SU2_OMP_FOR
+  }
+}
 }  // namespace
 
 void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container,
@@ -4638,6 +4701,13 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo
         }
         END_SU2_OMP_FOR
       }
+
+      /*--- The Full-MG startup solves the turbulence model on coarse grids. ---*/
+
+      for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) {
+        RestrictWallDistance(geometry_container[iZone][iInst][iMesh - 1], geometry_container[iZone][iInst][iMesh],
+                             config);
+      }
     }
   }
 }

Common/src/geometry/CMultiGridGeometry.cpp

@@ -1619,7 +1619,7 @@ CMultiGridGeometry::CFrontSeeds CMultiGridGeometry::SeedFrontNodes(const CGeomet
   /*--- Paving order: a no-slip wall first, a slip wall next, everything else last. ---*/
   auto tierOfBC = [](unsigned short bc) -> char {
     if ((bc == HEAT_FLUX) || (bc == ISOTHERMAL) || (bc == CHT_WALL_INTERFACE) || (bc == SMOLUCHOWSKI_MAXWELL)) return 0;
-    return (bc == EULER_WALL) ? 1 : 2;
+    return ((bc == EULER_WALL) || (bc == SYMMETRY_PLANE)) ? 1 : 2;
   };
 
   for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) {
@@ -2322,6 +2322,63 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co
     byLayer[iColumn][layerOf[iPoint]].push_back(iPoint);
   }
 
+  /*--- A boundary run with an odd number of seeds, and every node where two boundaries meet,
+   *    leaves one column a single node wide. Two such columns side by side are emitted together,
+   *    so the domain gets square cells instead of a file of narrow ones. ---*/
+  vector<long> partnerOf(nColumn, -1);
+  {
+    auto isNarrow = [&](unsigned long iColumn) {
+      return isSeeded[iColumn] && (byLayer[iColumn].size() > 1) && (byLayer[iColumn][0].size() == 1);
+    };
+
+    /*--- Both columns must stand shoulder to shoulder over every layer they share, otherwise the
+     *    cells the merge makes are not compact. ---*/
+    auto aligned = [&](unsigned long iColumn, unsigned long jColumn) {
+      const auto nCommon = std::min(byLayer[iColumn].size(), byLayer[jColumn].size());
+      for (auto iLayer = 1ul; iLayer < nCommon; ++iLayer) {
+        const auto iPoint = byLayer[iColumn][iLayer].front();
+        const auto jPoint = byLayer[jColumn][iLayer].front();
+        bool touch = false;
+        for (auto kPoint : fine_grid->nodes->GetPoints(iPoint)) touch = touch || (kPoint == jPoint);
+        if (!touch) return false;
+      }
+      return nCommon > 1;
+    };
+
+    /*--- Taken in global index order, so the pairing does not depend on the partitioning. ---*/
+    vector<unsigned long> narrow;
+    for (auto iColumn = 0ul; iColumn < nColumn; ++iColumn)
+      if (isNarrow(iColumn)) narrow.push_back(iColumn);
+    std::sort(narrow.begin(), narrow.end(), [&](unsigned long a, unsigned long b) {
+      return fine_grid->nodes->GetGlobalIndex(byLayer[a][0].front()) <
+             fine_grid->nodes->GetGlobalIndex(byLayer[b][0].front());
+    });
+
+    for (auto iColumn : narrow) {
+      if (partnerOf[iColumn] >= 0) continue;
+      auto best = NO_COLUMN;
+      auto bestKey = std::numeric_limits<unsigned long>::max();
+
+      for (auto jPoint : fine_grid->nodes->GetPoints(byLayer[iColumn][1].front())) {
+        const auto jColumn = columnOf[jPoint];
+        if ((jColumn == NO_COLUMN) || (jColumn == iColumn)) continue;
+        if ((partnerOf[jColumn] >= 0) || !isNarrow(jColumn)) continue;
+        if (tierOf[jColumn] != tierOf[iColumn]) continue;
+        if (!aligned(iColumn, jColumn)) continue;
+
+        const auto key = fine_grid->nodes->GetGlobalIndex(byLayer[jColumn][0].front());
+        if (key < bestKey) {
+          bestKey = key;
+          best = jColumn;
+        }
+      }
+
+      if (best == NO_COLUMN) continue;
+      partnerOf[iColumn] = static_cast<long>(best);
+      partnerOf[best] = static_cast<long>(iColumn);
+    }
+  }
+
   auto emitGroup = [&](const vector<unsigned long>& group) {
     nodes->SetChildren_CV(Index_CoarseCV, group);
     for (auto iPoint : group) {
@@ -2334,6 +2391,7 @@ string CMultiGridGeometry::PaveAdvancingFronts(unsigned long& Index_CoarseCV, co
 
   auto minDepth = std::numeric_limits<unsigned long>::max(), maxDepth = 0ul;
   vector<unsigned long> group;
+  vector<vector<unsigned long>> paired;
 
   /*--- Hand a set of nodes out as coarse CVs, each connected and within the size limit. ---*/
   vector<char> inSet(nPointFine, 0);
@@ -2379,27 +2437,44 @@ 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();
       if (depthOf[iColumn] == 0) neverGrewCV.push_back(baseCV);
       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<long>(iColumn))) continue;
+
+    const auto* emitted = &layers;
+    if (partnerOf[iColumn] > static_cast<long>(iColumn)) {
+      const auto& other = byLayer[partnerOf[iColumn]];
+      paired = layers;
+      for (auto k = 1ul; k < other.size(); ++k) {
+        if (k < paired.size())
+          paired[k].insert(paired[k].end(), other[k].begin(), other[k].end());
+        else
+          paired.push_back(other[k]);
+      }
+      emitted = &paired;
+    }
+
     /*--- Above it, consecutive layers are blocked so the coarse cell coarsens by the same ratio
      *    along the column as the patch does across it. ---*/
-    while (iLayer < layers.size()) {
+    while (iLayer < emitted->size()) {
       group.clear();
-      const auto block = BlockFor(maxAgglomSize, layers[iLayer].size());
-      for (auto k = 0ul; (k < block) && (iLayer < layers.size()); ++k) {
-        if (!group.empty() && (group.size() + layers[iLayer].size() > static_cast<size_t>(maxAgglomSize))) break;
-        group.insert(group.end(), layers[iLayer].begin(), layers[iLayer].end());
+      const auto block = BlockFor(maxAgglomSize, (*emitted)[iLayer].size());
+      for (auto k = 0ul; (k < block) && (iLayer < emitted->size()); ++k) {
+        if (!group.empty() && (group.size() + (*emitted)[iLayer].size() > static_cast<size_t>(maxAgglomSize))) break;
+        group.insert(group.end(), (*emitted)[iLayer].begin(), (*emitted)[iLayer].end());
         iLayer++;
       }
       /*--- A single layer wider than the limit still has to go somewhere. ---*/
       if (group.empty()) {
-        group = layers[iLayer];
+        group = (*emitted)[iLayer];
         iLayer++;
       }
       emitConnected(group);

SU2_CFD/include/integration/CMultiGridIntegration.hpp

@@ -26,6 +26,7 @@
  */
 
 #include "CIntegration.hpp"
+#include "../../../Common/include/containers/container_decorators.hpp"
 
 /*!
  * \class CMultiGridIntegration
@@ -190,7 +191,20 @@ class CMultiGridIntegration final : public CIntegration {
    * \param[in] config - Definition of the particular problem.
    */
   void GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
-                                 CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config);
+                                 CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config,
+                                 unsigned short iMesh);
+
+  /*!
+   * \brief Least-squares gradient of the coarse-grid correction, limited so no child value
+   * leaves the range of the coarse stencil.
+   * \param[in] sol_coarse - Solver holding the correction in Solution_Old.
+   * \param[in] geo_coarse - Geometrical definition of the coarse grid.
+   * \param[in] geo_fine - Geometrical definition of the fine grid.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the coarse mesh.
+   */
+  void ComputeProlongationGradient(CSolver *sol_coarse, CGeometry *geo_coarse, CGeometry *geo_fine,
+                                   const CConfig *config, unsigned short iMesh);
 
   /*!
    * \brief Do an implicit smoothing of the prolongated correction.
@@ -319,7 +333,8 @@ class CMultiGridIntegration final : public CIntegration {
                         passivedouble lastRMS[2], char& exitReason,
                         passivedouble& worstStepRatio, unsigned short& worstStep);
 
-  static constexpr int MAX_MG_LEVELS = 10;
+  /*--- CConfig rejects a larger MGLEVEL, so the per-level arrays below always fit. ---*/
+  static constexpr int MAX_MG_LEVELS = MAX_MGLEVELS;
 
   /*--- Bounds the fixed-size stack buffers in the restriction and prolongation kernels,
    *    independently of the CSysMatrix limit of the same name. ---*/
@@ -367,6 +382,10 @@ class CMultiGridIntegration final : public CIntegration {
   enum class MGStartupPromote { NONE, BUDGET, CONVERGENCE, STAGNATION };
   MGStartupPromote mg_startup_promote_reason = MGStartupPromote::NONE;
 
+  /*! \brief Limited least-squares gradient of the correction, indexed by coarse level.
+   *  Allocated on first use, only when MG_LINEAR_PROLONGATION is on. */
+  vector<CVectorOfMatrix> prolongGradient;
+
   vector<passivedouble> mg_startup_conv_start; /*!< \brief Field values when the active level became active. */
   vector<passivedouble> mg_startup_conv_prev;  /*!< \brief Field values on the previous iteration. */
   unsigned long mg_startup_stall_count = 0;    /*!< \brief Consecutive iterations without useful reduction. */