[WIP] fix multigrid agglomeration

Common/include/CConfig.hpp

@@ -483,19 +483,28 @@ class CConfig {
   unsigned short **DegreeFFDBox;      /*!< \brief Degree of the FFD boxes. */
   string *FFDTag;                     /*!< \brief Parameters of the design variable. */
   string *TagFFDBox;                  /*!< \brief Tag of the FFD box. */
-  unsigned short GeometryMode;        /*!< \brief Gemoetry mode (analysis or gradient computation). */
+  unsigned short GeometryMode;        /*!< \brief Geometry mode (analysis or gradient computation). */
   unsigned short MGCycle;             /*!< \brief Kind of multigrid cycle. */
   unsigned short FinestMesh;          /*!< \brief Finest mesh for the full multigrid approach. */
   unsigned short nFFD_Fix_IDir,
   nFFD_Fix_JDir, nFFD_Fix_KDir;       /*!< \brief Number of planes fixed in the FFD. */
   unsigned short nMG_PreSmooth,       /*!< \brief Number of MG pre-smooth parameters found in config file. */
   nMG_PostSmooth,                     /*!< \brief Number of MG post-smooth parameters found in config file. */
   nMG_CorrecSmooth;                   /*!< \brief Number of MG correct-smooth parameters found in config file. */
+  unsigned long MG_Min_MeshSize;      /*!< \brief Minimum mesh size per coarsest level to allow another MG level. */
+
   short *FFD_Fix_IDir,
   *FFD_Fix_JDir, *FFD_Fix_KDir;       /*!< \brief Exact sections. */
   unsigned short *MG_PreSmooth,       /*!< \brief Multigrid Pre smoothing. */
   *MG_PostSmooth,                     /*!< \brief Multigrid Post smoothing. */
   *MG_CorrecSmooth;                   /*!< \brief Multigrid Jacobi implicit smoothing of the correction. */
+  bool MG_Smooth_EarlyExit; /*!< \brief Enable early exit for MG smoothing. */
+  su2double MG_Smooth_Res_Threshold; /*!< \brief Residual reduction threshold for early exit. */
+  bool MG_Smooth_Output; /*!< \brief Output per-iteration multigrid smoothing info. */
+  bool MG_Implicit_Lines; /*!< \\brief Enable implicit-lines agglomeration from walls. */
+  bool MG_Implicit_Debug; /*!< \brief Enable debug output for implicit-lines agglomeration. */
+  bool MG_DebugHaloCoordinates; /*!< \brief Enable halo CV coordinate validation for multigrid. */
+  su2double MG_Smooth_Coeff; /*!< \brief Smoothing coefficient for multigrid correction smoothing. */
   su2double *LocationStations;        /*!< \brief Airfoil sections in wing slicing subroutine. */
 
   ENUM_MULTIZONE Kind_MZSolver;    /*!< \brief Kind of multizone solver.  */
@@ -2888,7 +2897,7 @@ class CConfig {
   unsigned short GetFinestMesh(void) const { return FinestMesh; }
 
   /*!
-   * \brief Get the kind of multigrid (V or W).
+   * \brief Get the kind of multigrid (V, W or FULLMG).
    * \note This variable is used in a recursive way to perform the different kind of cycles
    * \return 0 or 1 depending of we are dealing with a V or W cycle.
    */
@@ -3031,7 +3040,28 @@ class CConfig {
    * \brief Get the number of Runge-Kutta steps.
    * \return Number of Runge-Kutta steps.
    */
-  unsigned short GetnRKStep(void) const { return nRKStep; }
+ unsigned short GetnRKStep(void) const {
+
+    unsigned short iRKLimit = 1;
+
+    switch (GetKind_TimeIntScheme()) {
+      case RUNGE_KUTTA_EXPLICIT:
+        iRKLimit = GetnRKStep();
+        break;
+      case CLASSICAL_RK4_EXPLICIT:
+        iRKLimit = 4;
+        break;
+      case EULER_EXPLICIT:
+      case EULER_IMPLICIT:
+        iRKLimit = 1;
+        break;
+      default:
+        iRKLimit = 1;
+        break;
+    }
+  //return nRKStep;
+  return iRKLimit;
+}
 
   /*!
    * \brief Get the number of time levels for time accurate local time stepping.
@@ -3816,6 +3846,55 @@ class CConfig {
     return MG_CorrecSmooth[val_mesh];
   }
 
+  /*!
+   * \brief Get whether early exit is enabled for MG smoothing.
+   * \return True if early exit is enabled.
+   */
+  bool GetMG_Smooth_EarlyExit() const { return MG_Smooth_EarlyExit; }
+
+  /*!
+   * \brief Get the residual threshold for early exit in MG smoothing.
+   * \return Residual threshold.
+   */
+  su2double GetMG_Smooth_Res_Threshold() const { return MG_Smooth_Res_Threshold; }
+
+  /*!
+   * \brief Get whether per-iteration output is enabled for MG smoothing.
+   * \return True if output is enabled.
+   */
+  bool GetMG_Smooth_Output() const { return MG_Smooth_Output; }
+
+  /*!\
+   * \brief Get whether implicit-lines agglomeration is enabled.
+   * \return True if implicit-lines agglomeration is enabled.
+   */
+  bool GetMG_Implicit_Lines() const { return MG_Implicit_Lines; }
+
+  /*!\
+   * \brief Get whether implicit-lines debug output is enabled.
+   * \return True if implicit-lines debug output is enabled.
+   */
+  bool GetMG_Implicit_Debug() const { return MG_Implicit_Debug; }
+
+  /*!\
+   * \brief Get whether halo CV coordinate validation is enabled for multigrid.
+   * \return True if halo coordinate validation is enabled.
+   */
+  bool GetMG_DebugHaloCoordinates() const { return MG_DebugHaloCoordinates; }
+
+  /*!\
+   * \brief Get the minimum mesh size threshold used to compute effective MG levels.
+   * \return Minimum mesh size per coarsest level.
+   */
+  unsigned long GetMG_Min_MeshSize() const { return MG_Min_MeshSize; }
+
+
+  /*!
+   * \brief Get the smoothing coefficient for MG correction smoothing.
+   * \return Smoothing coefficient.
+   */
+  su2double GetMG_Smooth_Coeff() const { return MG_Smooth_Coeff; }
+
   /*!
    * \brief plane of the FFD (I axis) that should be fixed.
    * \param[in] val_index - Index of the arrray with all the planes in the I direction that should be fixed.
@@ -6736,6 +6815,9 @@ class CConfig {
    */
   su2double GetDamp_Res_Restric(void) const { return Damp_Res_Restric; }
 
+  /*!\n+   * \brief Set the damping factor for the residual restriction.\n+   * \param[in] val New value for the damping factor.\n+   */
+  void SetDamp_Res_Restric(su2double val) { Damp_Res_Restric = val; }
+
   /*!
    * \brief Value of the damping factor for the correction prolongation.
    * \return Value of the damping factor.

Common/include/geometry/CGeometry.hpp

@@ -245,6 +245,7 @@ class CGeometry {
 
   /*!< \brief Bool if boundary-marker is straight(2D)/plane(3D) for each local marker. */
   vector<bool> boundIsStraight;
+
   vector<su2double> SurfaceAreaCfgFile; /*!< \brief Total Surface area for all markers. */
 
   /*--- Partitioning-specific variables ---*/

Common/include/geometry/CMultiGridGeometry.hpp

@@ -28,36 +28,37 @@
 #pragma once
 
 #include "CGeometry.hpp"
+class CMultiGridQueue;
 
 /*!
  * \class CMultiGridGeometry
- * \brief Class for defining the multigrid geometry, the main delicated part is the
+ * \brief Class for defining the multigrid geometry, the main delegated part is the
  *        agglomeration stage, which is done in the declaration.
  * \author F. Palacios
  */
 class CMultiGridGeometry final : public CGeometry {
  private:
   /*!
-   * \brief Determine if a CVPoint van be agglomerated, if it have the same marker point as the seed.
+   * \brief Determine if a CVPoint can be agglomerated, if it has the same marker point as the seed.
    * \param[in] CVPoint - Control volume to be agglomerated.
    * \param[in] marker_seed - Marker of the seed.
    * \param[in] fine_grid - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    * \return <code>TRUE</code> or <code>FALSE</code> depending if the control volume can be agglomerated.
    */
-  bool SetBoundAgglomeration(unsigned long CVPoint, short marker_seed, const CGeometry* fine_grid,
+  bool SetBoundAgglomeration(unsigned long CVPoint, vector<short> marker_seed, const CGeometry* fine_grid,
                              const CConfig* config) const;
 
   /*!
-   * \brief Determine if a can be agglomerated using geometrical criteria.
+   * \brief Determine if a Point can be agglomerated using geometrical criteria.
    * \param[in] iPoint - Seed point.
    * \param[in] fine_grid - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
   bool GeometricalCheck(unsigned long iPoint, const CGeometry* fine_grid, const CConfig* config) const;
 
   /*!
-   * \brief Determine if a CVPoint van be agglomerated, if it have the same marker point as the seed.
+   * \brief Determine if a CVPoint can be agglomerated, if it has the same marker point as the seed.
    * \param[out] Suitable_Indirect_Neighbors - List of Indirect Neighbours that can be agglomerated.
    * \param[in] iPoint - Seed point.
    * \param[in] Index_CoarseCV - Index of agglomerated point.
@@ -66,6 +67,32 @@ class CMultiGridGeometry final : public CGeometry {
   void SetSuitableNeighbors(vector<unsigned long>& Suitable_Indirect_Neighbors, unsigned long iPoint,
                             unsigned long Index_CoarseCV, const CGeometry* fine_grid) const;
 
+  /*!
+   * \brief Agglomerate high-aspect ratio interior cells along implicit lines starting from wall vertices.
+   * \param[in,out] Index_CoarseCV - Reference to the current coarse control volume index that will be
+   *                                  incremented as new coarse CVs are created.
+   * \param[in] fine_grid - Geometrical definition of the problem (fine grid).
+   * \param[in] config - Definition of the particular problem.
+   * \param[in,out] MGQueue_InnerCV - Queue used for STEP 2; agglomerated points will be removed from it.
+   */
+  void AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const CGeometry* fine_grid, const CConfig* config,
+                                CMultiGridQueue& MGQueue_InnerCV);
+
+  /*!
+   * \brief Compute surface straightness for multigrid geometry.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeSurfStraightness(CConfig* config);
+
+  /*!
+   * \brief Compute local curvature at a boundary vertex on Euler wall.
+   * \param[in] fine_grid - Fine grid geometry.
+   * \param[in] iPoint - Point index.
+   * \param[in] iMarker - Marker index.
+   * \return Maximum angle (in degrees) between this vertex normal and adjacent vertex normals.
+   */
+  su2double ComputeLocalCurvature(const CGeometry* fine_grid, unsigned long iPoint, unsigned short iMarker) const;
+
  public:
   /*--- This is to suppress Woverloaded-virtual, omitting it has no negative impact. ---*/
   using CGeometry::SetBoundControlVolume;
@@ -154,4 +181,11 @@ class CMultiGridGeometry final : public CGeometry {
    * \param[in] val_marker - Index of the boundary marker.
    */
   void SetMultiGridWallTemperature(const CGeometry* fine_grid, unsigned short val_marker) override;
+
+  /*!
+   * \brief Validate that halo CV coordinates match corresponding domain CVs on remote ranks (debug feature).
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Multigrid level for reporting.
+   */
+  void ValidateHaloCoordinates(const CConfig* config, unsigned short iMesh) const;
 };

Common/include/geometry/CMultiGridQueue.hpp

@@ -93,7 +93,7 @@ class CMultiGridQueue {
   void IncrPriorityCV(unsigned long incrPoint);
 
   /*!
-   * \brief Increase the priority of the CV.
+   * \brief Reduce the priority of the CV.
    * \param[in] redPoint - Index of the control volume.
    */
   void RedPriorityCV(unsigned long redPoint);

Common/src/CConfig.cpp

@@ -1955,6 +1955,31 @@ void CConfig::SetConfig_Options() {
   /*!\brief MG_DAMP_PROLONGATION\n DESCRIPTION: Damping factor for the correction prolongation. DEFAULT 0.75 \ingroup Config*/
   addDoubleOption("MG_DAMP_PROLONGATION", Damp_Correc_Prolong, 0.75);
 
+  /*!\brief MG_SMOOTH_EARLY_EXIT\n DESCRIPTION: Enable early exit for MG smoothing iterations based on RMS residual. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_SMOOTH_EARLY_EXIT", MG_Smooth_EarlyExit, false);
+  /*!\brief MG_SMOOTH_RES_THRESHOLD\n DESCRIPTION: RMS residual threshold for early exit in MG smoothing. DEFAULT: 1e-2 \ingroup Config*/
+  addDoubleOption("MG_SMOOTH_RES_THRESHOLD", MG_Smooth_Res_Threshold, 1e-2);
+
+  /*!\brief MG_SMOOTH_OUTPUT\n DESCRIPTION: Output per-iteration RMS for MG smoothing. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_SMOOTH_OUTPUT", MG_Smooth_Output, false);
+
+  /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Enable agglomeration along implicit lines from wall seeds. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_IMPLICIT_LINES", MG_Implicit_Lines, false);
+
+  /*!\brief MG_IMPLICIT_DEBUG\n DESCRIPTION: Enable debug output for implicit-lines agglomeration. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_IMPLICIT_DEBUG", MG_Implicit_Debug, false);
+
+  /*!\brief MG_DEBUG_HALO_COORDINATES\n DESCRIPTION: Enable halo CV coordinate validation for multigrid (expensive MPI check). DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_DEBUG_HALO_COORDINATES", MG_DebugHaloCoordinates, false);
+
+  /*!\brief MG_MIN_MESHSIZE
+   \ DESCRIPTION: Minimum global mesh size (points) to allow another multigrid level. DEFAULT: 1000 \ingroup Config*/
+  addUnsignedLongOption("MG_MIN_MESHSIZE", MG_Min_MeshSize, 1000);
+
+
+  /*!\brief MG_SMOOTH_COEFF\n DESCRIPTION: Smoothing coefficient for MG correction smoothing. DEFAULT: 1.25 \ingroup Config*/
+  addDoubleOption("MG_SMOOTH_COEFF", MG_Smooth_Coeff, 1.25);
+
   /*!\par CONFIG_CATEGORY: Spatial Discretization \ingroup Config*/
   /*--- Options related to the spatial discretization ---*/
 
@@ -3538,7 +3563,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   }
 
   /*--- Check if MULTIGRID is requested in VOLUME_OUTPUT and set the config boolean accordingly. ---*/
-  Wrt_MultiGrid = false;
+  Wrt_MultiGrid = true;
   for (unsigned short iField = 0; iField < nVolumeOutput; iField++) {
     if(VolumeOutput[iField].find("MULTIGRID") != string::npos) {
       Wrt_MultiGrid = true;

Common/src/geometry/CGeometry.cpp

@@ -2472,8 +2472,7 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
    * All nodes that are shared by multiple symmetries have to get a corrected normal. */
 
   /*--- Compute if markers are straight lines or planes. ---*/
-  ComputeSurfStraightness(config, false);
-
+  ComputeSurfStraightness(config, true);
   symmetryNormals.clear();
   symmetryNormals.resize(nMarker);
   std::vector<unsigned short> symMarkers, curvedSymMarkers;
@@ -2589,19 +2588,18 @@ void CGeometry::ComputeSurfStraightness(const CConfig* config, bool print_on_scr
   constexpr passivedouble epsilon = 1.0e-6;
   su2double Area;
   string Local_TagBound, Global_TagBound;
-
-  vector<su2double> Normal(nDim), UnitNormal(nDim), RefUnitNormal(nDim);
-
-  /*--- Assume now that this boundary marker is straight. As soon as one
-        AreaElement is found that is not aligend with a Reference then it is
-        certain that the boundary marker is not straight and one can stop
-        searching. Another possibility is that this process doesn't own
-        any nodes of that boundary, in that case we also have to assume the
-        boundary is straight.
-        Any boundary type other than SYMMETRY_PLANE or EULER_WALL gets
-        the value false (or see cases specified in the conditional below)
-        which could be wrong. ---*/
+  cout << "Computing surface straightness for symmetry and Euler wall boundary markers..." << endl;
+
+  vector<su2double> Normal(nDim), UnitNormal(nDim), RefUnitNormal(nDim); /*--- Assume now that this boundary marker is
+       straight. As soon as one AreaElement is found that is not aligned with a Reference then it is certain that the
+       boundary marker is not straight and one can stop searching. Another possibility is that this process doesn't own
+       any nodes of that boundary, in that case we also have to assume the
+       boundary is straight.
+       Any boundary type other than SYMMETRY_PLANE or EULER_WALL gets
+       the value false (or see cases specified in the conditional below)
+       which could be wrong. ---*/
   boundIsStraight.resize(nMarker);
+  cout << "boundisstraight size = " << boundIsStraight.size() << endl;
   fill(boundIsStraight.begin(), boundIsStraight.end(), true);
 
   /*--- Loop over all local markers ---*/
@@ -2685,7 +2683,7 @@ void CGeometry::ComputeSurfStraightness(const CConfig* config, bool print_on_scr
     /*--- Product of type <int>(bool) is equivalnt to a 'logical and' ---*/
     SU2_MPI::Allreduce(Buff_Send_isStraight.data(), Buff_Recv_isStraight.data(), nMarker_Global, MPI_INT, MPI_PROD,
                        SU2_MPI::GetComm());
-
+    cout << "global markers = " << nMarker_Global << endl;
     /*--- Print results on screen. ---*/
     if (rank == MASTER_NODE) {
       for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
@@ -3902,11 +3900,13 @@ void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const*
     for (auto step = 0u; step < iMesh; ++step) {
       auto coarseMesh = geometry[iMesh - 1 - step];
       if (step)
-        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint)
+        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint) {
           CoarseGridColor_(iPoint, step) = CoarseGridColor_(coarseMesh->nodes->GetParent_CV(iPoint), step - 1);
+        }
       else
-        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint)
+        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint) {
           CoarseGridColor_(iPoint, step) = color[coarseMesh->nodes->GetParent_CV(iPoint)];
+        }
     }
   }
 }