diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index dc96c21ad3b..b15c00e333b 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -488,7 +488,7 @@ 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,
@@ -2912,7 +2912,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.
    */
diff --git a/Common/include/geometry/CMultiGridGeometry.hpp b/Common/include/geometry/CMultiGridGeometry.hpp
index 00faa8126c6..692a85360d4 100644
--- a/Common/include/geometry/CMultiGridGeometry.hpp
+++ b/Common/include/geometry/CMultiGridGeometry.hpp
@@ -28,28 +28,29 @@
 #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.
@@ -57,7 +58,7 @@ class CMultiGridGeometry final : public CGeometry {
   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,21 @@ 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 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;
diff --git a/Common/include/geometry/CMultiGridQueue.hpp b/Common/include/geometry/CMultiGridQueue.hpp
index 296d4e8f73d..c16e23520e4 100644
--- a/Common/include/geometry/CMultiGridQueue.hpp
+++ b/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);
diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp
index db26d86930a..68bab580fe1 100644
--- a/Common/src/geometry/CGeometry.cpp
+++ b/Common/src/geometry/CGeometry.cpp
@@ -2591,16 +2591,14 @@ void CGeometry::ComputeSurfStraightness(const CConfig* config, bool print_on_scr
   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
+        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. ---*/
+        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);
   fill(boundIsStraight.begin(), boundIsStraight.end(), true);
 
@@ -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)];
+        }
     }
   }
 }
diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp
index a52d111237c..1fccf3fd240 100644
--- a/Common/src/geometry/CMultiGridGeometry.cpp
+++ b/Common/src/geometry/CMultiGridGeometry.cpp
@@ -29,20 +29,36 @@
 #include "../../include/geometry/CMultiGridQueue.hpp"
 #include "../../include/toolboxes/printing_toolbox.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
+#include <cmath>
+#include <unordered_map>
+#include <climits>
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+#include <map>
 
 CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, unsigned short iMesh) : CGeometry() {
   nDim = fine_grid->GetnDim();  // Write the number of dimensions of the coarse grid.
 
-  /*--- Create a queue system to do the agglomeration
+  /*--- Maximum agglomeration size in 2D is 4 nodes, in 3D is 8 nodes. ---*/
+  const short int maxAgglomSize = (nDim == 2) ? 4 : 8;
+
+  /*--- Inherit boundary properties from fine grid ---*/
+  boundIsStraight = fine_grid->boundIsStraight;
+
+  /*--- Agglomeration Scheme II (Nishikawa, Diskin, Thomas)
+        Create a queue system to do the agglomeration
    1st) More than two markers ---> Vertices (never agglomerate)
    2nd) Two markers ---> Edges (agglomerate if same BC, never agglomerate if different BC)
    3rd) One marker ---> Surface (always agglomerate)
    4th) No marker ---> Internal Volume (always agglomerate) ---*/
 
+  // note that for MPI, we introduce interfaces and we can choose to have agglomeration over
+  // the interface or not. Nishikawa chooses not to agglomerate over interfaces.
+
   /*--- Set a marker to indicate indirect agglomeration, for quads and hexs,
-   i.e. consider up to neighbors of neighbors of neighbors.
+   i.e. consider up to neighbors of neighbors.
    For other levels this information is propagated down during their construction. ---*/
-
   if (iMesh == MESH_1) {
     for (auto iPoint = 0ul; iPoint < fine_grid->GetnPoint(); iPoint++)
       fine_grid->nodes->SetAgglomerate_Indirect(iPoint, false);
@@ -59,7 +75,6 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
   }
 
   /*--- Create the coarse grid structure using as baseline the fine grid ---*/
-
   CMultiGridQueue MGQueue_InnerCV(fine_grid->GetnPoint());
   vector<unsigned long> Suitable_Indirect_Neighbors;
 
@@ -67,15 +82,26 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
 
   unsigned long Index_CoarseCV = 0;
 
-  /*--- The first step is the boundary agglomeration. ---*/
+  /*--- Statistics for Euler wall agglomeration ---*/
+  map<unsigned short, unsigned long> euler_wall_agglomerated, euler_wall_rejected_curvature,
+      euler_wall_rejected_straight;
+  for (unsigned short iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) {
+    if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) {
+      euler_wall_agglomerated[iMarker] = 0;
+      euler_wall_rejected_curvature[iMarker] = 0;
+      euler_wall_rejected_straight[iMarker] = 0;
+    }
+  }
 
+  /*--- STEP 1: The first step is the boundary agglomeration. ---*/
   for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) {
     for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) {
       const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode();
 
       /*--- If the element has not been previously agglomerated and it
-       belongs to this physical domain, and it meets the geometrical
-       criteria, the agglomeration is studied. ---*/
+      belongs to this physical domain, and it meets the geometrical
+      criteria, the agglomeration is studied. ---*/
+      vector<short> marker_seed;
 
       if ((!fine_grid->nodes->GetAgglomerate(iPoint)) && (fine_grid->nodes->GetDomain(iPoint)) &&
           (GeometricalCheck(iPoint, fine_grid, config))) {
@@ -89,94 +115,146 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
         /*--- We add the seed point (child) to the parent control volume ---*/
 
         nodes->SetChildren_CV(Index_CoarseCV, 0, iPoint);
-        bool agglomerate_seed = true;
+        bool agglomerate_seed = false;
         auto counter = 0;
         unsigned short copy_marker[3] = {};
-        const auto marker_seed = iMarker;
+        marker_seed.push_back(iMarker);
 
         /*--- For a particular point in the fine grid we save all the markers
          that are in that point ---*/
 
-        for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker() && counter < 3; jMarker++) {
+        for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker(); jMarker++) {
+          const string Marker_Tag = config->GetMarker_All_TagBound(iMarker);
           if (fine_grid->nodes->GetVertex(iPoint, jMarker) != -1) {
             copy_marker[counter] = jMarker;
             counter++;
+
+            if (jMarker != iMarker) {
+              marker_seed.push_back(jMarker);
+            }
           }
         }
 
-        /*--- To aglomerate a vertex it must have only one physical bc!!
-         This can be improved. If there is only a marker, it is a good
+        /*--- To agglomerate a vertex it must have only one physical bc.
+         This can be improved. If there is only one marker, it is a good
          candidate for agglomeration ---*/
 
+        /*--- 1 BC, so either an edge in 2D or the interior of a plane in 3D ---*/
+        /*--- Valley -> Valley : conditionally allowed when both points are on the same marker. ---*/
+        /*--- ! Note that in the case of MPI SEND_RECEIVE markers, we might need other conditions ---*/
         if (counter == 1) {
+          // The seed/parent is one valley, so we set this part to true
+          // if the child is on the same valley, we set it to true as well.
           agglomerate_seed = true;
-
-          /*--- Euler walls can be curved and agglomerating them leads to difficulties ---*/
-          if (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL) agglomerate_seed = false;
+          /*--- Euler walls: check curvature-based agglomeration criterion ---*/
+          if (config->GetMarker_All_KindBC(marker_seed[0]) == EULER_WALL) {
+            /*--- Allow agglomeration if marker is straight OR local curvature is small ---*/
+            if (!boundIsStraight[marker_seed[0]]) {
+              /*--- Compute local curvature at this point ---*/
+              su2double local_curvature = ComputeLocalCurvature(fine_grid, iPoint, marker_seed[0]);
+              // limit to 45 degrees
+              if (local_curvature >= 45.0) {
+                agglomerate_seed = false;  // High curvature: do not agglomerate
+                euler_wall_rejected_curvature[marker_seed[0]]++;
+              } else {
+                euler_wall_agglomerated[marker_seed[0]]++;
+              }
+            } else {
+              /*--- Straight wall: agglomerate ---*/
+              euler_wall_agglomerated[marker_seed[0]]++;
+            }
+          }
+          /*--- Note that if the (single) marker is a SEND_RECEIVE, then the node is actually an interior point.
+                In that case it can only be agglomerated with another interior point. ---*/
+          if (config->GetMarker_All_KindBC(marker_seed[0]) == SEND_RECEIVE) {
+            agglomerate_seed = true;
+          }
         }
-        /*--- If there are two markers, we will agglomerate if any of the
-         markers is SEND_RECEIVE ---*/
 
+        /*--- Note that in 2D, this is a corner and we do not agglomerate unless they both are SEND_RECEIVE. ---*/
+        /*--- In 3D, we agglomerate if the 2 markers are the same. ---*/
         if (counter == 2) {
-          agglomerate_seed = (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) ||
-                             (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE);
-
-          /* --- Euler walls can also not be agglomerated when the point has 2 markers ---*/
-          if ((config->GetMarker_All_KindBC(copy_marker[0]) == EULER_WALL) ||
-              (config->GetMarker_All_KindBC(copy_marker[1]) == EULER_WALL)) {
-            agglomerate_seed = false;
+          /*--- agglomerate if one of the 2 markers are MPI markers. ---*/
+          if (nDim == 2)
+            agglomerate_seed = (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) ||
+                               (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE);
+
+          /*--- agglomerate if both markers are the same. ---*/
+          if (nDim == 3) agglomerate_seed = (copy_marker[0] == copy_marker[1]);
+
+          /*--- Euler walls: check curvature-based agglomeration criterion for both markers ---*/
+          // only in 3d because in 2d it's a corner
+          bool euler_wall_rejected_here = false;
+          for (unsigned short i = 0; i < 2; i++) {
+            if ((nDim == 3) && (config->GetMarker_All_KindBC(copy_marker[i]) == EULER_WALL)) {
+              if (!boundIsStraight[copy_marker[i]]) {
+                /*--- Compute local curvature at this point ---*/
+                su2double local_curvature = ComputeLocalCurvature(fine_grid, iPoint, copy_marker[i]);
+                // limit to 45 degrees
+                if (local_curvature >= 45.0) {
+                  agglomerate_seed = false;  // High curvature: do not agglomerate
+                  euler_wall_rejected_curvature[copy_marker[i]]++;
+                  euler_wall_rejected_here = true;
+                }
+              }
+              /*--- Track agglomeration if not rejected ---*/
+              if (agglomerate_seed && !euler_wall_rejected_here) {
+                euler_wall_agglomerated[copy_marker[i]]++;
+              }
+            }
           }
         }
 
         /*--- If there are more than 2 markers, the aglomeration will be discarded ---*/
-
         if (counter > 2) agglomerate_seed = false;
 
-        /*--- If the seed can be agglomerated, we try to agglomerate more points ---*/
-
+        /*--- If the seed (parent) can be agglomerated, we try to agglomerate connected childs to the parent ---*/
+        /*--- Note that in 2D we allow a maximum of 4 nodes to be agglomerated ---*/
         if (agglomerate_seed) {
           /*--- Now we do a sweep over all the nodes that surround the seed point ---*/
-
           for (auto CVPoint : fine_grid->nodes->GetPoints(iPoint)) {
             /*--- The new point can be agglomerated ---*/
-
             if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) {
               /*--- We set the value of the parent ---*/
-
               fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV);
 
               /*--- We set the value of the child ---*/
-
               nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint);
               nChildren++;
+              /*--- In 2D, we agglomerate exactly 2 nodes if the nodes are on the line edge. ---*/
+              if ((nDim == 2) && (counter == 1)) break;
+              /*--- In 3D, we agglomerate exactly 2 nodes if the nodes are on the surface edge. ---*/
+              if ((nDim == 3) && (counter == 2)) break;
+              /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes (counter==1 in 3D). ---*/
+              if (nChildren == maxAgglomSize) break;
             }
           }
 
-          Suitable_Indirect_Neighbors.clear();
-
-          if (fine_grid->nodes->GetAgglomerate_Indirect(iPoint))
-            SetSuitableNeighbors(Suitable_Indirect_Neighbors, iPoint, Index_CoarseCV, fine_grid);
-
-          /*--- Now we do a sweep over all the indirect nodes that can be added ---*/
-
-          for (auto CVPoint : Suitable_Indirect_Neighbors) {
-            /*--- The new point can be agglomerated ---*/
-
-            if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) {
-              /*--- We set the value of the parent ---*/
-
-              fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV);
-
-              /*--- We set the indirect agglomeration information of the corse point
-               based on its children in the fine grid. ---*/
-
-              if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint))
-                nodes->SetAgglomerate_Indirect(Index_CoarseCV, true);
-
-              /*--- We set the value of the child ---*/
-
-              nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint);
-              nChildren++;
+          /*--- Only take into account indirect neighbors for 3D faces, not 2D. ---*/
+          if (nDim == 3) {
+            Suitable_Indirect_Neighbors.clear();
+
+            if (fine_grid->nodes->GetAgglomerate_Indirect(iPoint))
+              SetSuitableNeighbors(Suitable_Indirect_Neighbors, iPoint, Index_CoarseCV, fine_grid);
+
+            /*--- Now we do a sweep over all the indirect nodes that can be added ---*/
+            for (auto CVPoint : Suitable_Indirect_Neighbors) {
+              /*--- The new point can be agglomerated ---*/
+              if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) {
+                /*--- We set the value of the parent ---*/
+                fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV);
+
+                /*--- We set the indirect agglomeration information of the corse point
+                based on its children in the fine grid. ---*/
+                if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint))
+                  nodes->SetAgglomerate_Indirect(Index_CoarseCV, true);
+
+                /*--- We set the value of the child ---*/
+                nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint);
+                nChildren++;
+                /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes. ---*/
+                if (nChildren == maxAgglomSize) break;
+              }
             }
           }
         }
@@ -223,8 +301,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
     }
   }
 
-  /*--- Agglomerate the domain points. ---*/
-
+  /*--- STEP 2: Agglomerate the domain points. ---*/
   auto iteration = 0ul;
   while (!MGQueue_InnerCV.EmptyQueue() && (iteration < fine_grid->GetnPoint())) {
     const auto iPoint = MGQueue_InnerCV.NextCV();
@@ -271,6 +348,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
 
           MGQueue_InnerCV.Update(CVPoint, fine_grid);
         }
+        if (nChildren == maxAgglomSize) break;
       }
 
       /*--- Identify the indirect neighbors ---*/
@@ -282,11 +360,11 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
       /*--- Now we do a sweep over all the indirect nodes that can be added ---*/
 
       for (auto CVPoint : Suitable_Indirect_Neighbors) {
+        // if we have reached the maximum, get out.
+        if (nChildren == maxAgglomSize) break;
         /*--- The new point can be agglomerated ---*/
-
         if ((!fine_grid->nodes->GetAgglomerate(CVPoint)) && (fine_grid->nodes->GetDomain(CVPoint))) {
           /*--- We set the value of the parent ---*/
-
           fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV);
 
           /*--- We set the indirect agglomeration information ---*/
@@ -340,16 +418,56 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
   for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) {
     if (nodes->GetnPoint(iCoarsePoint) == 1) {
       /*--- Find the neighbor of the isolated point. This neighbor is the right control volume ---*/
-
       const auto iCoarsePoint_Complete = nodes->GetPoint(iCoarsePoint, 0);
 
-      /*--- Add the children to the connected control volume (and modify its parent indexing).
-       Identify the child CV from the finest grid and add it to the correct control volume.
-       Set the parent CV of iFinePoint. Instead of using the original one
-       (iCoarsePoint), use the new one (iCoarsePoint_Complete) ---*/
+      /*--- Check if merging would exceed the maximum agglomeration size ---*/
+      auto nChildren_Target = nodes->GetnChildren_CV(iCoarsePoint_Complete);
+      auto nChildren_Isolated = nodes->GetnChildren_CV(iCoarsePoint);
+      auto nChildren_Total = nChildren_Target + nChildren_Isolated;
+
+      /*--- If the total would exceed maxAgglomSize, try to redistribute children to neighbors ---*/
+      if (nChildren_Total > maxAgglomSize) {
+        /*--- Find neighbors of the target coarse point that have room ---*/
+        unsigned short nChildrenToRedistribute = nChildren_Total - maxAgglomSize;
+
+        for (auto jCoarsePoint : nodes->GetPoints(iCoarsePoint_Complete)) {
+          if (nChildrenToRedistribute == 0) break;
+
+          auto nChildren_Neighbor = nodes->GetnChildren_CV(jCoarsePoint);
+          if (nChildren_Neighbor < maxAgglomSize) {
+            unsigned short nCanTransfer =
+                min(nChildrenToRedistribute, static_cast<unsigned short>(maxAgglomSize - nChildren_Neighbor));
+
+            /*--- Transfer children from target to neighbor ---*/
+            for (unsigned short iTransfer = 0; iTransfer < nCanTransfer; iTransfer++) {
+              /*--- Take from the end of the target's children list ---*/
+              auto nChildren_Current = nodes->GetnChildren_CV(iCoarsePoint_Complete);
+              if (nChildren_Current > 0) {
+                auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint_Complete, nChildren_Current - 1);
 
-      auto nChildren = nodes->GetnChildren_CV(iCoarsePoint_Complete);
+                /*--- Add to neighbor ---*/
+                auto nChildren_Neighbor_Current = nodes->GetnChildren_CV(jCoarsePoint);
+                nodes->SetChildren_CV(jCoarsePoint, nChildren_Neighbor_Current, iFinePoint);
+                nodes->SetnChildren_CV(jCoarsePoint, nChildren_Neighbor_Current + 1);
+
+                /*--- Update parent ---*/
+                fine_grid->nodes->SetParent_CV(iFinePoint, jCoarsePoint);
+
+                /*--- Remove from target (by reducing count) ---*/
+                nodes->SetnChildren_CV(iCoarsePoint_Complete, nChildren_Current - 1);
+
+                nChildrenToRedistribute--;
+              }
+            }
+          }
+        }
 
+        /*--- Update the target's child count after redistribution ---*/
+        nChildren_Target = nodes->GetnChildren_CV(iCoarsePoint_Complete);
+      }
+
+      /*--- Add the isolated point's children to the target control volume ---*/
+      auto nChildren = nChildren_Target;
       for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
         const auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
         nodes->SetChildren_CV(iCoarsePoint_Complete, nChildren, iFinePoint);
@@ -358,7 +476,6 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
       }
 
       /*--- Update the number of children control volumes ---*/
-
       nodes->SetnChildren_CV(iCoarsePoint_Complete, nChildren);
       nodes->SetnChildren_CV(iCoarsePoint, 0);
     }
@@ -369,15 +486,28 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
   nodes->ResetPoints();
 
 #ifdef HAVE_MPI
+  /*--- Reset halo point parents before MPI agglomeration.
+   This is critical for multi-level multigrid: when creating level N from level N-1,
+   the fine grid (level N-1) already has Parent_CV set from when it was created from level N-2.
+   Those parent indices point to level N, but when creating level N+1, they would be
+   incorrectly interpreted as level N+1 indices. Resetting ensures clean agglomeration. ---*/
+
+  for (auto iPoint = fine_grid->GetnPointDomain(); iPoint < fine_grid->GetnPoint(); iPoint++) {
+    fine_grid->nodes->SetParent_CV(iPoint, ULONG_MAX);
+  }
+
   /*--- Dealing with MPI parallelization, the objective is that the received nodes must be agglomerated
    in the same way as the donor (send) nodes. Send the node agglomeration information of the donor
    (parent and children). The agglomerated halos of this rank are set according to the rank where
    they are domain points. ---*/
 
   for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+    cout << " marker name = " << config->GetMarker_All_TagBound(iMarker) << endl;
+    cout << " marker type = " << config->GetMarker_All_KindBC(iMarker) << endl;
+    cout << " send/recv = " << config->GetMarker_All_SendRecv(iMarker) << endl;
     if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0)) {
-      const auto MarkerS = iMarker;
-      const auto MarkerR = iMarker + 1;
+      const auto MarkerS = iMarker;      // sending marker
+      const auto MarkerR = iMarker + 1;  // receiving marker
 
       const auto send_to = config->GetMarker_All_SendRecv(MarkerS) - 1;
       const auto receive_from = abs(config->GetMarker_All_SendRecv(MarkerR)) - 1;
@@ -424,38 +554,151 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
       vector<unsigned long> Parent_Local(nVertexR);
       vector<unsigned long> Children_Local(nVertexR);
 
+      /*--- First pass: Determine which parents will actually be used (have non-skipped children).
+       This prevents creating orphaned halo CVs that have coordinates (0,0,0). ---*/
+      vector<bool> parent_used(Aux_Parent.size(), false);
+      vector<unsigned long> parent_local_index(Aux_Parent.size(), ULONG_MAX);
+
       for (auto iVertex = 0ul; iVertex < nVertexR; iVertex++) {
-        /*--- We use the same sorting as in the donor domain, i.e. the local parents
-         are numbered according to their order in the remote rank. ---*/
+        const auto iPoint_Fine = fine_grid->vertex[MarkerR][iVertex]->GetNode();
+        auto existing_parent = fine_grid->nodes->GetParent_CV(iPoint_Fine);
+
+        /*--- Skip if already agglomerated (first-wins policy) ---*/
+        if (existing_parent != ULONG_MAX) continue;
+
+        /*--- Find which parent this vertex maps to ---*/
+        for (auto jVertex = 0ul; jVertex < Aux_Parent.size(); jVertex++) {
+          if (Parent_Remote[iVertex] == Aux_Parent[jVertex]) {
+            parent_used[jVertex] = true;
+            break;
+          }
+        }
+      }
+
+      /*--- Assign local indices only to used parents ---*/
+      unsigned long nUsedParents = 0;
+      for (auto jVertex = 0ul; jVertex < Aux_Parent.size(); jVertex++) {
+        if (parent_used[jVertex]) {
+          parent_local_index[jVertex] = Index_CoarseCV + nUsedParents;
+          nUsedParents++;
+        }
+      }
 
+      /*--- Now map each received vertex to its local parent ---*/
+      for (auto iVertex = 0ul; iVertex < nVertexR; iVertex++) {
+        Parent_Local[iVertex] = ULONG_MAX;
         for (auto jVertex = 0ul; jVertex < Aux_Parent.size(); jVertex++) {
           if (Parent_Remote[iVertex] == Aux_Parent[jVertex]) {
-            Parent_Local[iVertex] = jVertex + Index_CoarseCV;
+            Parent_Local[iVertex] = parent_local_index[jVertex];
             break;
           }
         }
+
+        /*--- Validate that parent mapping was found (only matters if not skipped later) ---*/
+        if (Parent_Local[iVertex] == ULONG_MAX) {
+          SU2_MPI::Error(string("MPI agglomeration failed to map parent index ") + to_string(Parent_Remote[iVertex]) +
+                             string(" for vertex ") + to_string(iVertex),
+                         CURRENT_FUNCTION);
+        }
+
         Children_Local[iVertex] = fine_grid->vertex[MarkerR][iVertex]->GetNode();
       }
 
-      Index_CoarseCV += Aux_Parent.size();
+      /*--- Only increment by the number of parents that will actually be used ---*/
+      Index_CoarseCV += nUsedParents;
 
-      vector<unsigned short> nChildren_MPI(Index_CoarseCV, 0);
+      /*--- Debug counters ---*/
+      unsigned long nConflicts = 0, nSkipped = 0, nOutOfBounds = 0, nSuccess = 0;
 
       /*--- Create the final structure ---*/
       for (auto iVertex = 0ul; iVertex < nVertexR; iVertex++) {
         const auto iPoint_Coarse = Parent_Local[iVertex];
         const auto iPoint_Fine = Children_Local[iVertex];
 
-        /*--- Be careful, it is possible that a node changes the agglomeration configuration,
-         the priority is always when receiving the information. ---*/
+        /*--- Debug: Check for out-of-bounds access ---*/
+        if (iPoint_Coarse >= Index_CoarseCV) {
+          cout << "ERROR [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR
+               << "]: Out-of-bounds coarse CV index " << iPoint_Coarse << " >= " << Index_CoarseCV << " (vertex "
+               << iVertex << ", fine point " << iPoint_Fine << ")" << endl;
+          nOutOfBounds++;
+          continue;
+        }
+
+        /*--- Solution 1: Skip if this halo point was already agglomerated ---*/
+        auto existing_parent = fine_grid->nodes->GetParent_CV(iPoint_Fine);
+        if (existing_parent != ULONG_MAX) {
+          if (existing_parent != iPoint_Coarse) {
+            /*--- Conflict detected: different parent from different interface ---*/
+            nConflicts++;
+
+            /*--- Only print detailed info for first few conflicts or if suspicious ---*/
+            if (nConflicts <= 5 || existing_parent < nPointDomain) {
+              cout << "INFO [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR << "]: Halo point "
+                   << iPoint_Fine << " already agglomerated to parent " << existing_parent
+                   << (existing_parent < nPointDomain ? " (DOMAIN CV!)" : " (halo CV)") << ", skipping reassignment to "
+                   << iPoint_Coarse << " (from rank " << receive_from << ")" << endl;
+            }
+          } else {
+            /*--- Same parent from different interface (duplicate) - just skip silently ---*/
+            nSkipped++;
+          }
+          continue;  // First-wins: keep existing assignment
+        }
+
+        /*--- First assignment for this halo point - proceed with agglomeration ---*/
+
+        /*--- Critical fix: Append to existing children, don't overwrite ---*/
+        auto existing_children_count = nodes->GetnChildren_CV(iPoint_Coarse);
+
         fine_grid->nodes->SetParent_CV(iPoint_Fine, iPoint_Coarse);
-        nodes->SetChildren_CV(iPoint_Coarse, nChildren_MPI[iPoint_Coarse], iPoint_Fine);
-        nChildren_MPI[iPoint_Coarse]++;
-        nodes->SetnChildren_CV(iPoint_Coarse, nChildren_MPI[iPoint_Coarse]);
+        nodes->SetChildren_CV(iPoint_Coarse, existing_children_count, iPoint_Fine);
+        nodes->SetnChildren_CV(iPoint_Coarse, existing_children_count + 1);
         nodes->SetDomain(iPoint_Coarse, false);
+        nSuccess++;
+      }
+
+      /*--- Debug: Report statistics for this marker pair ---*/
+      cout << "MPI Agglomeration [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR << " (rank " << send_to
+           << " <-> " << receive_from << ")]: " << nSuccess << " assigned, " << nSkipped << " duplicates, "
+           << nConflicts << " conflicts";
+      if (nOutOfBounds > 0) {
+        cout << ", " << nOutOfBounds << " OUT-OF-BOUNDS (CRITICAL!)";
+      }
+      cout << endl;
+
+      if (nConflicts > 5) {
+        cout << "  Note: Only first 5 conflicts shown in detail, total conflicts = " << nConflicts << endl;
+      }
+
+      /*--- Debug: Validate buffer size assumption ---*/
+      if (nVertexS != nVertexR) {
+        cout << "WARNING [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR
+             << "]: Asymmetric interface - nVertexS=" << nVertexS << " != nVertexR=" << nVertexR << endl;
       }
     }
   }
+
+  /*--- Post-process: Check for orphaned coarse CVs (should be none with new logic) ---*/
+
+  if (size > SINGLE_NODE) {
+    unsigned long nOrphaned = 0;
+
+    /*--- Count orphaned CVs for reporting ---*/
+    for (auto iCoarse = nPointDomain; iCoarse < Index_CoarseCV; iCoarse++) {
+      if (nodes->GetnChildren_CV(iCoarse) == 0) {
+        nOrphaned++;
+        /*--- This shouldn't happen with the new parent prefiltering logic ---*/
+        cout << "WARNING [Rank " << rank << "]: Orphaned halo CV " << iCoarse
+             << " detected (should not occur with current logic)" << endl;
+      }
+    }
+
+    if (nOrphaned > 0) {
+      cout << "WARNING [Rank " << rank << "]: " << nOrphaned
+           << " orphaned halo coarse CVs found - this indicates a logic error!" << endl;
+    }
+  }
+
 #endif  // HAVE_MPI
 
   /*--- Update the number of points after the MPI agglomeration ---*/
@@ -463,25 +706,28 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
   nPoint = Index_CoarseCV;
 
   /*--- Console output with the summary of the agglomeration ---*/
-
-  unsigned long nPointFine = fine_grid->GetnPoint();
+  // nijso: do not include halo points in the count
+  unsigned long nPointFine = fine_grid->GetnPointDomain();
   unsigned long Global_nPointCoarse, Global_nPointFine;
 
-  SU2_MPI::Allreduce(&nPoint, &Global_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&nPointDomain, &Global_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
   SU2_MPI::Allreduce(&nPointFine, &Global_nPointFine, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
 
   SetGlobal_nPointDomain(Global_nPointCoarse);
 
   if (iMesh != MESH_0) {
-    const su2double factor = 1.5;
-    const su2double Coeff = pow(su2double(Global_nPointFine) / Global_nPointCoarse, 1.0 / nDim);
+    /*--- Note: CFL at the coarse levels have a large impact on convergence,
+          this should be rewritten to use adaptive CFL. ---*/
+    const su2double factor = 1.0;
+    const su2double Coeff = 1.1;  // pow(su2double(Global_nPointFine) / Global_nPointCoarse, 1.0 / nDim);
     const su2double CFL = factor * config->GetCFL(iMesh - 1) / Coeff;
     config->SetCFL(iMesh, CFL);
   }
 
   const su2double ratio = su2double(Global_nPointFine) / su2double(Global_nPointCoarse);
-
-  if (((nDim == 2) && (ratio < 2.5)) || ((nDim == 3) && (ratio < 2.5))) {
+  cout << "********** ratio = " << ratio << endl;
+  // lower value leads to more levels being accepted.
+  if (((nDim == 2) && (ratio < 1.5)) || ((nDim == 3) && (ratio < 1.5))) {
     config->SetMGLevels(iMesh - 1);
   } else if (rank == MASTER_NODE) {
     PrintingToolbox::CTablePrinter MGTable(&std::cout);
@@ -503,11 +749,70 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
     }
   }
 
+  /*--- Output Euler wall agglomeration statistics ---*/
+  if (rank == MASTER_NODE) {
+    /*--- Gather global statistics for Euler walls ---*/
+    bool has_euler_walls = false;
+    for (unsigned short iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) {
+      if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) {
+        has_euler_walls = true;
+        break;
+      }
+    }
+
+    if (has_euler_walls) {
+      cout << endl;
+      cout << "Euler Wall Agglomeration Statistics (45° curvature threshold):" << endl;
+      cout << "----------------------------------------------------------------" << endl;
+
+      for (unsigned short iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) {
+        if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) {
+          string marker_name = config->GetMarker_All_TagBound(iMarker);
+          unsigned long agglomerated = euler_wall_agglomerated[iMarker];
+          unsigned long rejected = euler_wall_rejected_curvature[iMarker];
+          unsigned long total = agglomerated + rejected;
+
+          if (total > 0) {
+            su2double accept_rate = 100.0 * su2double(agglomerated) / su2double(total);
+            cout << "  Marker: " << marker_name << endl;
+            cout << "    Seeds agglomerated:       " << agglomerated << " (" << std::setprecision(1) << std::fixed
+                 << accept_rate << "%)" << endl;
+            cout << "    Seeds rejected (>45° curv): " << rejected << " (" << std::setprecision(1) << std::fixed
+                 << (100.0 - accept_rate) << "%)" << endl;
+          }
+        }
+      }
+      cout << "----------------------------------------------------------------" << endl;
+    }
+  }
+
   edgeColorGroupSize = config->GetEdgeColoringGroupSize();
 }
 
-bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, short marker_seed, const CGeometry* fine_grid,
-                                               const CConfig* config) const {
+bool CMultiGridGeometry::GeometricalCheck(unsigned long iPoint, const CGeometry* fine_grid,
+                                          const CConfig* config) const {
+  su2double max_dimension = 1.2;
+
+  /*--- Evaluate the total size of the element ---*/
+
+  bool Volume = true;
+  su2double ratio = pow(fine_grid->nodes->GetVolume(iPoint), 1.0 / su2double(nDim)) * max_dimension;
+  su2double limit = pow(config->GetDomainVolume(), 1.0 / su2double(nDim));
+  if (ratio > limit) {
+    Volume = false;
+    cout << "Volume limit reached!" << endl;
+  }
+  /*--- Evaluate the stretching of the element ---*/
+
+  bool Stretching = true;
+
+  return (Stretching && Volume);
+}
+
+void CMultiGridGeometry::ComputeSurfStraightness(CConfig* config) { CGeometry::ComputeSurfStraightness(config, true); }
+
+bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, vector<short> marker_seed,
+                                               const CGeometry* fine_grid, const CConfig* config) const {
   bool agglomerate_CV = false;
 
   /*--- Basic condition, the point has not been previously agglomerated, it belongs to the domain,
@@ -517,12 +822,13 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, short mark
       (GeometricalCheck(CVPoint, fine_grid, config))) {
     /*--- If the point belongs to a boundary, its type must be compatible with the seed marker. ---*/
 
+    int counter = 0;
+    unsigned short copy_marker[3] = {};
+
     if (fine_grid->nodes->GetBoundary(CVPoint)) {
       /*--- Identify the markers of the vertex that we want to agglomerate ---*/
 
       // count number of markers on the agglomeration candidate
-      int counter = 0;
-      unsigned short copy_marker[3] = {};
       for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker() && counter < 3; jMarker++) {
         if (fine_grid->nodes->GetVertex(CVPoint, jMarker) != -1) {
           copy_marker[counter] = jMarker;
@@ -530,94 +836,56 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, short mark
         }
       }
 
-      /*--- The basic condition is that the aglomerated vertex must have the same physical marker,
+      /*--- The basic condition is that the agglomerated vertex must have the same physical marker,
        but eventually a send-receive condition ---*/
 
-      /*--- Only one marker in the vertex that is going to be aglomerated ---*/
+      /*--- Only one marker in the vertex that is going to be agglomerated ---*/
 
+      /*--- Valley -> Valley: only if of the same type---*/
       if (counter == 1) {
         /*--- We agglomerate if there is only one marker and it is the same marker as the seed marker ---*/
-        // note that this should be the same marker id, not just the same marker type
-        if (copy_marker[0] == marker_seed) agglomerate_CV = true;
-
-        /*--- If there is only one marker, but the marker is the SEND_RECEIVE ---*/
-
-        if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) {
-          agglomerate_CV = true;
+        // So this is the case when in 2D we are on an edge, and in 3D we are in the interior of a surface.
+        // note that this should be the same marker id, not just the same marker type.
+        // also note that the seed point can have 2 markers, one of them may be a send-receive.
+        if ((marker_seed.size() == 1) && (copy_marker[0] == marker_seed[0])) agglomerate_CV = true;
+        if ((marker_seed.size() == 2) && (config->GetMarker_All_KindBC(marker_seed[0]) == SEND_RECEIVE)) {
+          if (copy_marker[0] == marker_seed[1]) {
+            agglomerate_CV = true;
+          }
         }
-
-        if ((config->GetMarker_All_KindBC(marker_seed) == SYMMETRY_PLANE) ||
-            (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL)) {
-          if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) {
-            agglomerate_CV = false;
+        if ((marker_seed.size() == 2) && (config->GetMarker_All_KindBC(marker_seed[1]) == SEND_RECEIVE)) {
+          if (copy_marker[0] == marker_seed[0]) {
+            agglomerate_CV = true;
           }
         }
+
+        /*--- Note: If there is only one marker, but the marker is the SEND_RECEIVE, then the point is actually an
+              interior point and we do not agglomerate.  ---*/
       }
 
       /*--- If there are two markers in the vertex that is going to be aglomerated ---*/
 
       if (counter == 2) {
-        /*--- First we verify that the seed is a physical boundary ---*/
-
-        if (config->GetMarker_All_KindBC(marker_seed) != SEND_RECEIVE) {
-          /*--- Then we check that one of the markers is equal to the seed marker, and the other is send/receive ---*/
+        /*--- Both markers have to be the same. ---*/
 
-          if (((copy_marker[0] == marker_seed) && (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE)) ||
-              ((config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) && (copy_marker[1] == marker_seed))) {
+        if (marker_seed.size() == 2) {
+          if (((copy_marker[0] == marker_seed[0]) && (copy_marker[1] == marker_seed[1])) ||
+              ((copy_marker[0] == marker_seed[1]) && (copy_marker[1] == marker_seed[0]))) {
             agglomerate_CV = true;
           }
         }
       }
     }
-    /*--- If the element belongs to the domain, it is always agglomerated. ---*/
+    /*--- If the element belongs to the domain, it is never agglomerated with a boundary node. ---*/
     else {
-      agglomerate_CV = true;
-
-      // actually, for symmetry (and possibly other cells) we only agglomerate cells that are on the marker
-      // at this point, the seed was on the boundary and the CV was not. so we check if the seed is a symmetry
-      if ((config->GetMarker_All_KindBC(marker_seed) == SYMMETRY_PLANE) ||
-          (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL)) {
-        agglomerate_CV = false;
-      }
+      agglomerate_CV = false;
     }
   }
 
   return agglomerate_CV;
 }
 
-bool CMultiGridGeometry::GeometricalCheck(unsigned long iPoint, const CGeometry* fine_grid,
-                                          const CConfig* config) const {
-  su2double max_dimension = 1.2;
-
-  /*--- Evaluate the total size of the element ---*/
-
-  bool Volume = true;
-  su2double ratio = pow(fine_grid->nodes->GetVolume(iPoint), 1.0 / su2double(nDim)) * max_dimension;
-  su2double limit = pow(config->GetDomainVolume(), 1.0 / su2double(nDim));
-  if (ratio > limit) Volume = false;
-
-  /*--- Evaluate the stretching of the element ---*/
-
-  bool Stretching = true;
-
-  /* unsigned short iNode, iDim;
-   unsigned long jPoint;
-   su2double *Coord_i = fine_grid->nodes->GetCoord(iPoint);
-   su2double max_dist = 0.0 ; su2double min_dist = 1E20;
-   for (iNode = 0; iNode < fine_grid->nodes->GetnPoint(iPoint); iNode ++) {
-   jPoint = fine_grid->nodes->GetPoint(iPoint, iNode);
-   su2double *Coord_j = fine_grid->nodes->GetCoord(jPoint);
-   su2double distance = 0.0;
-   for (iDim = 0; iDim < nDim; iDim++)
-   distance += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]);
-   distance = sqrt(distance);
-   max_dist = max(distance, max_dist);
-   min_dist = min(distance, min_dist);
-   }
-   if ( max_dist/min_dist > 100.0 ) Stretching = false;*/
-
-  return (Stretching && Volume);
-}
+/*--- ---*/
 
 void CMultiGridGeometry::SetSuitableNeighbors(vector<unsigned long>& Suitable_Indirect_Neighbors, unsigned long iPoint,
                                               unsigned long Index_CoarseCV, const CGeometry* fine_grid) const {
@@ -637,8 +905,8 @@ void CMultiGridGeometry::SetSuitableNeighbors(vector<unsigned long>& Suitable_In
 
       auto end = First_Neighbor_Points.end();
       if (find(First_Neighbor_Points.begin(), end, kPoint) == end) {
-        Second_Neighbor_Points.push_back(kPoint);
-        Second_Origin_Points.push_back(jPoint);
+        Second_Neighbor_Points.push_back(kPoint);  // neighbor of a neighbor, not connected to original ipoint
+        Second_Origin_Points.push_back(jPoint);    // the neighbor that is connected to ipoint
       }
     }
   }
@@ -668,53 +936,11 @@ void CMultiGridGeometry::SetSuitableNeighbors(vector<unsigned long>& Suitable_In
   sort(Suitable_Second_Neighbors.begin(), Suitable_Second_Neighbors.end());
   auto it1 = unique(Suitable_Second_Neighbors.begin(), Suitable_Second_Neighbors.end());
   Suitable_Second_Neighbors.resize(it1 - Suitable_Second_Neighbors.begin());
-
-  /*--- Create a list with the third neighbors, without first, second, and seed neighbors. ---*/
-
-  /// TODO: This repeats the process above but I doubt it catches any more points.
-
-  vector<unsigned long> Third_Neighbor_Points, Third_Origin_Points;
-
-  for (auto kPoint : Suitable_Second_Neighbors) {
-    for (auto lPoint : fine_grid->nodes->GetPoints(kPoint)) {
-      /*--- Check that the third neighbor does not belong to the first neighbors or the seed ---*/
-
-      auto end1 = First_Neighbor_Points.end();
-      if (find(First_Neighbor_Points.begin(), end1, lPoint) != end1) continue;
-
-      /*--- Check that the third neighbor does not belong to the second neighbors ---*/
-
-      auto end2 = Suitable_Second_Neighbors.end();
-      if (find(Suitable_Second_Neighbors.begin(), end2, lPoint) != end2) continue;
-
-      Third_Neighbor_Points.push_back(lPoint);
-      Third_Origin_Points.push_back(kPoint);
-    }
-  }
-
-  /*--- Identify those third neighbors that are repeated (candidate to be added). ---*/
-
-  for (auto iNeighbor = 0ul; iNeighbor < Third_Neighbor_Points.size(); iNeighbor++) {
-    for (auto jNeighbor = iNeighbor + 1; jNeighbor < Third_Neighbor_Points.size(); jNeighbor++) {
-      /*--- Repeated third neighbor with different origin ---*/
-
-      if ((Third_Neighbor_Points[iNeighbor] == Third_Neighbor_Points[jNeighbor]) &&
-          (Third_Origin_Points[iNeighbor] != Third_Origin_Points[jNeighbor])) {
-        Suitable_Indirect_Neighbors.push_back(Third_Neighbor_Points[iNeighbor]);
-      }
-    }
-  }
-
-  /*--- Remove duplicates from the final list of Suitable Indirect Neighbors. ---*/
-
-  sort(Suitable_Indirect_Neighbors.begin(), Suitable_Indirect_Neighbors.end());
-  auto it2 = unique(Suitable_Indirect_Neighbors.begin(), Suitable_Indirect_Neighbors.end());
-  Suitable_Indirect_Neighbors.resize(it2 - Suitable_Indirect_Neighbors.begin());
 }
 
 void CMultiGridGeometry::SetPoint_Connectivity(const CGeometry* fine_grid) {
   /*--- Temporary, CPoint (nodes) then compresses this structure. ---*/
-  vector<vector<unsigned long> > points(nPoint);
+  vector<vector<unsigned long>> points(nPoint);
 
   for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) {
     /*--- For each child CV (of the fine grid), ---*/
@@ -741,17 +967,14 @@ void CMultiGridGeometry::SetPoint_Connectivity(const CGeometry* fine_grid) {
 }
 
 void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* config) {
-  unsigned long iVertex, iFinePoint, iCoarsePoint;
-  unsigned short iMarker, iMarker_Tag, iChildren;
-
   nMarker = fine_grid->GetnMarker();
   unsigned short nMarker_Max = config->GetnMarker_Max();
 
   /*--- If any children node belong to the boundary then the entire control
    volume will belong to the boundary ---*/
-  for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++)
-    for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
-      iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
+  for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++)
+    for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
+      auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
       if (fine_grid->nodes->GetBoundary(iFinePoint)) {
         nodes->SetBoundary(iCoarsePoint, nMarker);
         break;
@@ -762,20 +985,20 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co
   nVertex = new unsigned long[nMarker];
 
   Tag_to_Marker = new string[nMarker_Max];
-  for (iMarker_Tag = 0; iMarker_Tag < nMarker_Max; iMarker_Tag++)
+  for (auto iMarker_Tag = 0u; iMarker_Tag < nMarker_Max; iMarker_Tag++)
     Tag_to_Marker[iMarker_Tag] = fine_grid->GetMarker_Tag(iMarker_Tag);
 
   /*--- Compute the number of vertices to do the dimensionalization ---*/
-  for (iMarker = 0; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0;
+  for (auto iMarker = 0u; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0;
 
-  for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) {
+  for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) {
     if (nodes->GetBoundary(iCoarsePoint)) {
-      for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
-        iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
-        for (iMarker = 0; iMarker < nMarker; iMarker++) {
+      for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
+        auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
+        for (auto iMarker = 0u; iMarker < nMarker; iMarker++) {
           if ((fine_grid->nodes->GetVertex(iFinePoint, iMarker) != -1) &&
               (nodes->GetVertex(iCoarsePoint, iMarker) == -1)) {
-            iVertex = nVertex[iMarker];
+            auto iVertex = nVertex[iMarker];
             nodes->SetVertex(iCoarsePoint, iVertex, iMarker);
             nVertex[iMarker]++;
           }
@@ -784,25 +1007,25 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co
     }
   }
 
-  for (iMarker = 0; iMarker < nMarker; iMarker++) {
+  for (auto iMarker = 0u; iMarker < nMarker; iMarker++) {
     vertex[iMarker] = new CVertex*[fine_grid->GetnVertex(iMarker) + 1];
     nVertex[iMarker] = 0;
   }
 
-  for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++)
+  for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++)
     if (nodes->GetBoundary(iCoarsePoint))
-      for (iMarker = 0; iMarker < nMarker; iMarker++) nodes->SetVertex(iCoarsePoint, -1, iMarker);
+      for (auto iMarker = 0u; iMarker < nMarker; iMarker++) nodes->SetVertex(iCoarsePoint, -1, iMarker);
 
-  for (iMarker = 0; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0;
+  for (auto iMarker = 0u; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0;
 
-  for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) {
+  for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) {
     if (nodes->GetBoundary(iCoarsePoint)) {
-      for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
-        iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
-        for (iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) {
+      for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
+        auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
+        for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) {
           if ((fine_grid->nodes->GetVertex(iFinePoint, iMarker) != -1) &&
               (nodes->GetVertex(iCoarsePoint, iMarker) == -1)) {
-            iVertex = nVertex[iMarker];
+            auto iVertex = nVertex[iMarker];
             vertex[iMarker][iVertex] = new CVertex(iCoarsePoint, nDim);
             nodes->SetVertex(iCoarsePoint, iVertex, iMarker);
 
@@ -819,15 +1042,13 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co
 }
 
 void CMultiGridGeometry::MatchActuator_Disk(const CConfig* config) {
-  unsigned short iMarker;
-  unsigned long iVertex, iPoint;
   int iProcessor = size;
 
-  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if ((config->GetMarker_All_KindBC(iMarker) == ACTDISK_INLET) ||
         (config->GetMarker_All_KindBC(iMarker) == ACTDISK_OUTLET)) {
-      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
-        iPoint = vertex[iMarker][iVertex]->GetNode();
+      for (auto iVertex = 0u; iVertex < nVertex[iMarker]; iVertex++) {
+        auto iPoint = vertex[iMarker][iVertex]->GetNode();
         if (nodes->GetDomain(iPoint)) {
           vertex[iMarker][iVertex]->SetDonorPoint(iPoint, nodes->GetGlobalIndex(iPoint), iVertex, iMarker, iProcessor);
         }
@@ -837,20 +1058,18 @@ void CMultiGridGeometry::MatchActuator_Disk(const CConfig* config) {
 }
 
 void CMultiGridGeometry::MatchPeriodic(const CConfig* config, unsigned short val_periodic) {
-  unsigned short iMarker, iPeriodic, nPeriodic;
-  unsigned long iVertex, iPoint;
   int iProcessor = rank;
 
   /*--- Evaluate the number of periodic boundary conditions ---*/
 
-  nPeriodic = config->GetnMarker_Periodic();
+  auto nPeriodic = config->GetnMarker_Periodic();
 
-  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
-      iPeriodic = config->GetMarker_All_PerBound(iMarker);
+      auto iPeriodic = config->GetMarker_All_PerBound(iMarker);
       if ((iPeriodic == val_periodic) || (iPeriodic == val_periodic + nPeriodic / 2)) {
-        for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
-          iPoint = vertex[iMarker][iVertex]->GetNode();
+        for (auto iVertex = 0u; iVertex < nVertex[iMarker]; iVertex++) {
+          auto iPoint = vertex[iMarker][iVertex]->GetNode();
           if (nodes->GetDomain(iPoint)) {
             vertex[iMarker][iVertex]->SetDonorPoint(iPoint, nodes->GetGlobalIndex(iPoint), iVertex, iMarker,
                                                     iProcessor);
@@ -863,18 +1082,12 @@ void CMultiGridGeometry::MatchPeriodic(const CConfig* config, unsigned short val
 
 void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned short action) {
   BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-    unsigned long iFinePoint, iCoarsePoint, iEdge, iParent;
-    long FineEdge, CoarseEdge;
-    unsigned short iChildren;
-    bool change_face_orientation;
-    su2double Coarse_Volume, Area;
-
     /*--- Compute the area of the coarse volume ---*/
-    for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) {
+    for (auto iCoarsePoint = 0u; iCoarsePoint < nPoint; iCoarsePoint++) {
       nodes->SetVolume(iCoarsePoint, 0.0);
-      Coarse_Volume = 0.0;
-      for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
-        iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
+      su2double Coarse_Volume = 0.0;
+      for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
+        auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
         Coarse_Volume += fine_grid->nodes->GetVolume(iFinePoint);
       }
       nodes->SetVolume(iCoarsePoint, Coarse_Volume);
@@ -885,19 +1098,19 @@ void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned s
       edges->SetZeroValues();
     }
 
-    for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++)
-      for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
-        iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
+    for (auto iCoarsePoint = 0u; iCoarsePoint < nPoint; iCoarsePoint++)
+      for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
+        auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
 
         for (auto iFinePoint_Neighbor : fine_grid->nodes->GetPoints(iFinePoint)) {
-          iParent = fine_grid->nodes->GetParent_CV(iFinePoint_Neighbor);
+          auto iParent = fine_grid->nodes->GetParent_CV(iFinePoint_Neighbor);
           if ((iParent != iCoarsePoint) && (iParent < iCoarsePoint)) {
-            FineEdge = fine_grid->FindEdge(iFinePoint, iFinePoint_Neighbor);
+            auto FineEdge = fine_grid->FindEdge(iFinePoint, iFinePoint_Neighbor);
 
-            change_face_orientation = false;
+            bool change_face_orientation = false;
             if (iFinePoint < iFinePoint_Neighbor) change_face_orientation = true;
 
-            CoarseEdge = FindEdge(iParent, iCoarsePoint);
+            auto CoarseEdge = FindEdge(iParent, iCoarsePoint);
 
             const auto Normal = fine_grid->edges->GetNormal(FineEdge);
 
@@ -912,9 +1125,9 @@ void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned s
 
     /*--- Check if there is a normal with null area ---*/
 
-    for (iEdge = 0; iEdge < nEdge; iEdge++) {
+    for (auto iEdge = 0u; iEdge < nEdge; iEdge++) {
       const auto NormalFace = edges->GetNormal(iEdge);
-      Area = GeometryToolbox::Norm(nDim, NormalFace);
+      su2double Area = GeometryToolbox::Norm(nDim, NormalFace);
       if (Area == 0.0) {
         su2double DefaultNormal[3] = {EPS * EPS};
         edges->SetNormal(iEdge, DefaultNormal);
@@ -926,24 +1139,23 @@ void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned s
 
 void CMultiGridGeometry::SetBoundControlVolume(const CGeometry* fine_grid, const CConfig* config,
                                                unsigned short action) {
-  unsigned long iCoarsePoint, iFinePoint, FineVertex, iVertex;
-  su2double Normal[MAXNDIM] = {0.0}, Area, *NormalFace = nullptr;
+  su2double Normal[MAXNDIM] = {0.0}, *NormalFace = nullptr;
 
   if (action != ALLOCATE) {
     SU2_OMP_FOR_DYN(1)
-    for (auto iMarker = 0; iMarker < nMarker; iMarker++)
-      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) vertex[iMarker][iVertex]->SetZeroValues();
+    for (auto iMarker = 0u; iMarker < nMarker; iMarker++)
+      for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) vertex[iMarker][iVertex]->SetZeroValues();
     END_SU2_OMP_FOR
   }
 
   SU2_OMP_FOR_DYN(1)
-  for (auto iMarker = 0; iMarker < nMarker; iMarker++) {
-    for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
-      iCoarsePoint = vertex[iMarker][iVertex]->GetNode();
+  for (auto iMarker = 0u; iMarker < nMarker; iMarker++) {
+    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
+      auto iCoarsePoint = vertex[iMarker][iVertex]->GetNode();
       for (auto iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) {
-        iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
+        auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
         if (fine_grid->nodes->GetVertex(iFinePoint, iMarker) != -1) {
-          FineVertex = fine_grid->nodes->GetVertex(iFinePoint, iMarker);
+          auto FineVertex = fine_grid->nodes->GetVertex(iFinePoint, iMarker);
           fine_grid->vertex[iMarker][FineVertex]->GetNormal(Normal);
           vertex[iMarker][iVertex]->AddNormal(Normal);
         }
@@ -954,10 +1166,10 @@ void CMultiGridGeometry::SetBoundControlVolume(const CGeometry* fine_grid, const
 
   /*--- Check if there is a normal with null area ---*/
   SU2_OMP_FOR_DYN(1)
-  for (auto iMarker = 0; iMarker < nMarker; iMarker++) {
-    for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
+  for (auto iMarker = 0u; iMarker < nMarker; iMarker++) {
+    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
       NormalFace = vertex[iMarker][iVertex]->GetNormal();
-      Area = GeometryToolbox::Norm(nDim, NormalFace);
+      su2double Area = GeometryToolbox::Norm(nDim, NormalFace);
       if (Area == 0.0)
         for (auto iDim = 0; iDim < nDim; iDim++) NormalFace[iDim] = EPS * EPS;
     }
@@ -1046,36 +1258,32 @@ void CMultiGridGeometry::SetRestricted_GridVelocity(const CGeometry* fine_grid)
 }
 
 void CMultiGridGeometry::FindNormal_Neighbor(const CConfig* config) {
-  unsigned short iMarker, iDim;
-  unsigned long iPoint, iVertex;
-
-  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE &&
         config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY &&
         config->GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) {
-      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
-        iPoint = vertex[iMarker][iVertex]->GetNode();
+      for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
+        auto iPoint = vertex[iMarker][iVertex]->GetNode();
 
         /*--- If the node belong to the domain ---*/
         if (nodes->GetDomain(iPoint)) {
           /*--- Compute closest normal neighbor ---*/
-          su2double cos_max, scalar_prod, norm_vect, norm_Normal, cos_alpha, diff_coord;
           unsigned long Point_Normal = 0;
           su2double* Normal = vertex[iMarker][iVertex]->GetNormal();
-          cos_max = -1.0;
+          su2double cos_max = -1.0;
           for (auto jPoint : nodes->GetPoints(iPoint)) {
-            scalar_prod = 0.0;
-            norm_vect = 0.0;
-            norm_Normal = 0.0;
-            for (iDim = 0; iDim < nDim; iDim++) {
-              diff_coord = nodes->GetCoord(jPoint, iDim) - nodes->GetCoord(iPoint, iDim);
+            su2double scalar_prod = 0.0;
+            su2double norm_vect = 0.0;
+            su2double norm_Normal = 0.0;
+            for (auto iDim = 0u; iDim < nDim; iDim++) {
+              su2double diff_coord = nodes->GetCoord(jPoint, iDim) - nodes->GetCoord(iPoint, iDim);
               scalar_prod += diff_coord * Normal[iDim];
               norm_vect += diff_coord * diff_coord;
               norm_Normal += Normal[iDim] * Normal[iDim];
             }
             norm_vect = sqrt(norm_vect);
             norm_Normal = sqrt(norm_Normal);
-            cos_alpha = scalar_prod / (norm_vect * norm_Normal);
+            su2double cos_alpha = scalar_prod / (norm_vect * norm_Normal);
 
             /*--- Get maximum cosine (not minimum because normals are oriented inwards) ---*/
             if (cos_alpha >= cos_max) {
@@ -1089,3 +1297,67 @@ void CMultiGridGeometry::FindNormal_Neighbor(const CConfig* config) {
     }
   }
 }
+
+su2double CMultiGridGeometry::ComputeLocalCurvature(const CGeometry* fine_grid, unsigned long iPoint,
+                                                    unsigned short iMarker) const {
+  /*--- Compute local curvature (maximum angle between adjacent face normals) at a boundary vertex.
+        This is used to determine if agglomeration is safe based on a curvature threshold. ---*/
+
+  /*--- Get the vertex index for this point on this marker ---*/
+  long iVertex = fine_grid->nodes->GetVertex(iPoint, iMarker);
+  if (iVertex < 0) return 0.0;  // Point not on this marker
+
+  /*--- Get the normal at this vertex ---*/
+  su2double Normal_i[MAXNDIM] = {0.0};
+  fine_grid->vertex[iMarker][iVertex]->GetNormal(Normal_i);
+  su2double Area_i = GeometryToolbox::Norm(int(nDim), Normal_i);
+
+  if (Area_i < EPS) return 0.0;  // Skip degenerate vertices
+
+  /*--- Normalize the normal ---*/
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+    Normal_i[iDim] /= Area_i;
+  }
+
+  /*--- Find maximum angle with neighboring vertices on the same marker ---*/
+  su2double max_angle = 0.0;
+
+  /*--- Loop over edges connected to this point ---*/
+  for (unsigned short iEdge = 0; iEdge < fine_grid->nodes->GetnPoint(iPoint); iEdge++) {
+    unsigned long jPoint = fine_grid->nodes->GetPoint(iPoint, iEdge);
+
+    /*--- Check if neighbor is also on this marker ---*/
+    long jVertex = fine_grid->nodes->GetVertex(jPoint, iMarker);
+    if (jVertex < 0) continue;  // Not on this marker
+
+    /*--- Get normal at neighbor vertex ---*/
+    su2double Normal_j[MAXNDIM] = {0.0};
+    fine_grid->vertex[iMarker][jVertex]->GetNormal(Normal_j);
+    su2double Area_j = GeometryToolbox::Norm(int(nDim), Normal_j);
+
+    if (Area_j < EPS) continue;  // Skip degenerate neighbor
+
+    /*--- Normalize the neighbor normal ---*/
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      Normal_j[iDim] /= Area_j;
+    }
+
+    /*--- Compute dot product: cos(angle) = n_i · n_j ---*/
+    su2double dot_product = 0.0;
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      dot_product += Normal_i[iDim] * Normal_j[iDim];
+    }
+
+    /*--- Clamp to [-1, 1] to avoid numerical issues with acos ---*/
+    dot_product = max(-1.0, min(1.0, dot_product));
+
+    /*--- Compute angle in degrees ---*/
+    su2double angle_rad = acos(dot_product);
+    su2double angle_deg = angle_rad * 180.0 / PI_NUMBER;
+
+    /*--- Track maximum angle ---*/
+    max_angle = max(max_angle, angle_deg);
+  }
+
+  return max_angle;
+}
diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp
index 361792eac5e..b5bc0dfdf35 100644
--- a/Common/src/geometry/dual_grid/CPoint.cpp
+++ b/Common/src/geometry/dual_grid/CPoint.cpp
@@ -25,6 +25,7 @@
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <climits>
 #include "../../../include/geometry/dual_grid/CPoint.hpp"
 #include "../../../include/CConfig.hpp"
 #include "../../../include/parallelization/omp_structure.hpp"
@@ -73,7 +74,7 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) {
 
   /*--- Multigrid structures. ---*/
   if (config->GetnMGLevels() > 0) {
-    Parent_CV.resize(npoint) = 0;
+    Parent_CV.resize(npoint) = ULONG_MAX;  // Initialize to ULONG_MAX to indicate unagglomerated
     Agglomerate.resize(npoint) = false;
     Agglomerate_Indirect.resize(npoint) = false;
     /*--- The finest grid does not have children CV's. ---*/
diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp
index eb2b7ded797..39e99f27bd5 100644
--- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp
+++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp
@@ -168,7 +168,7 @@ class CMultiGridIntegration final : public CIntegration {
    * \param[in] InclSharedDomain - Include the shared domain in the interpolation.
    */
   void SetRestricted_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 = 999);
 
   /*!
    * \brief Initialize the adjoint solution using the primal problem.
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
index b61d3d61f33..26c4f8be33c 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
@@ -578,7 +578,7 @@ class CFVMFlowSolverBase : public CSolver {
 
     }
 
-    /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is diferent from 0.
+    /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is different from 0.
      * This is only done once because in dual time the time step cannot be variable. ---*/
 
     if (dual_time && (Iteration == config->GetRestart_Iter()) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) {
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 961ce04832f..77c093f6c10 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -810,7 +810,7 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) {
 
   geometry[MESH_0]->SetMGLevel(MESH_0);
   if ((config->GetnMGLevels() != 0) && (rank == MASTER_NODE))
-    cout << "Setting the multigrid structure." << endl;
+    cout << "Setting the multigrid structure. NMG="<< config->GetnMGLevels() << endl;
 
   /*--- Loop over all the new grid ---*/
 
diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp
index b616f0e6649..8378bde18d8 100644
--- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp
+++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp
@@ -28,6 +28,60 @@
 #include "../../include/integration/CMultiGridIntegration.hpp"
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 
+namespace {
+/*!
+ * \brief Helper function to enforce Euler wall BC by projecting momentum to tangent plane.
+ * \param[in] geo_coarse - Coarse grid geometry.
+ * \param[in] config - Problem configuration.
+ * \param[in,out] sol_coarse - Coarse grid solver (to access and modify solution/correction).
+ * \param[in] use_solution_old - If true, project Solution_Old (corrections); if false, project Solution.
+ */
+void ProjectEulerWallToTangentPlane(CGeometry* geo_coarse, const CConfig* config, CSolver* sol_coarse,
+                                    bool use_solution_old) {
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+    if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) {
+
+      SU2_OMP_FOR_STAT(32)
+      for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) {
+        const auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode();
+
+        if (!geo_coarse->nodes->GetDomain(Point_Coarse)) continue;
+
+        /*--- Get coarse grid normal ---*/
+        su2double Normal[3] = {0.0};
+        geo_coarse->vertex[iMarker][iVertex]->GetNormal(Normal);
+        const auto nDim = geo_coarse->GetnDim();
+        su2double Area = GeometryToolbox::Norm(nDim, Normal);
+
+        if (Area < EPS) continue;
+
+        /*--- Normalize normal vector ---*/
+        su2double UnitNormal[3] = {0.0};
+        for (auto iDim = 0u; iDim < nDim; iDim++) {
+          UnitNormal[iDim] = Normal[iDim] / Area;
+        }
+
+        /*--- Get current solution or correction ---*/
+        su2double* solution_coarse = use_solution_old ?
+          sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse) :
+          sol_coarse->GetNodes()->GetSolution(Point_Coarse);
+
+        /*--- Compute normal component of momentum (v·n or correction·n) ---*/
+        su2double momentum_n = 0.0;
+        for (auto iDim = 0u; iDim < nDim; iDim++) {
+          momentum_n += solution_coarse[iDim + 1] * UnitNormal[iDim];
+        }
+
+        /*--- Project to tangent plane: solution_coarse -= (solution_coarse·n)n ---*/
+        for (auto iDim = 0u; iDim < nDim; iDim++) {
+          solution_coarse[iDim + 1] -= momentum_n * UnitNormal[iDim];
+        }
+      }
+      END_SU2_OMP_FOR
+    }
+  }
+}
+}  // anonymous namespace
 
 CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { }
 
@@ -240,12 +294,27 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
 
     /*--- Compute $r_(k+1) = F_(k+1)(I^(k+1)_k u_k)$ ---*/
 
-    SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config);
+    SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh);
 
     solver_coarse->Preprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem, false);
 
     Space_Integration(geometry_coarse, solver_container_coarse, numerics_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem);
 
+    /*--- Store initial coarse grid residual before MG cycle ---*/
+    su2double max_res_coarse_before = 0.0;
+    for (unsigned long iPoint = 0; iPoint < geometry_coarse->GetnPointDomain(); iPoint++) {
+      const su2double* res = solver_coarse->LinSysRes.GetBlock(iPoint);
+      for (unsigned short iVar = 0; iVar < solver_coarse->GetnVar(); iVar++) {
+        max_res_coarse_before = max(max_res_coarse_before, fabs(res[iVar]));
+      }
+    }
+
+#ifdef HAVE_MPI
+    su2double global_max_res_coarse_before = 0.0;
+    SU2_MPI::Allreduce(&max_res_coarse_before, &global_max_res_coarse_before, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+    max_res_coarse_before = global_max_res_coarse_before;
+#endif
+
     /*--- Compute $P_(k+1) = I^(k+1)_k(r_k) - r_(k+1) ---*/
 
     SetForcing_Term(solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh+1);
@@ -272,6 +341,126 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
 
     GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config);
 
+    /*--- Measure coarse grid residual after MG cycle to assess convergence rate ---*/
+    su2double max_res_coarse_after = 0.0;
+    for (unsigned long iPoint = 0; iPoint < geometry_coarse->GetnPointDomain(); iPoint++) {
+      const su2double* res = solver_coarse->LinSysRes.GetBlock(iPoint);
+      for (unsigned short iVar = 0; iVar < solver_coarse->GetnVar(); iVar++) {
+        max_res_coarse_after = max(max_res_coarse_after, fabs(res[iVar]));
+      }
+    }
+
+#ifdef HAVE_MPI
+    su2double global_max_res_coarse_after = 0.0;
+    SU2_MPI::Allreduce(&max_res_coarse_after, &global_max_res_coarse_after, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+    max_res_coarse_after = global_max_res_coarse_after;
+#endif
+
+    /*--- Adaptive CFL based on coarse grid convergence rate ---*/
+    if (max_res_coarse_before > 1e-16) {
+      su2double convergence_rate = max_res_coarse_after / max_res_coarse_before;
+
+      /*--- Get current CFL values ---*/
+      su2double CFL_fine = config->GetCFL(iMesh);
+      su2double CFL_coarse_current = config->GetCFL(iMesh+1);
+      su2double current_coeff = CFL_fine / CFL_coarse_current;
+
+      /*--- Adaptive CFL for coarse multigrid levels ---*/
+      constexpr int MAX_MG_LEVELS = 10;
+      static int consecutive_decreasing[MAX_MG_LEVELS] = {};
+      static su2double prev_fine_res[MAX_MG_LEVELS] = {};
+      static su2double min_fine_res[MAX_MG_LEVELS] = {};
+      static int iterations_since_improvement[MAX_MG_LEVELS] = {};
+      static bool initialized = false;
+
+      const unsigned short nMGLevels = config->GetnMGLevels();
+      const unsigned short max_lvl = min(nMGLevels, (unsigned short)MAX_MG_LEVELS);
+
+      if (!initialized) {
+        for (int lvl = 0; lvl < MAX_MG_LEVELS; lvl++) {
+          consecutive_decreasing[lvl] = 0;
+          prev_fine_res[lvl] = 0.0;
+          min_fine_res[lvl] = 1e10;
+          iterations_since_improvement[lvl] = 0;
+        }
+        initialized = true;
+      }
+
+      unsigned short lvl = min(iMesh, (unsigned short)(max_lvl - 1));
+
+      /*--- Track fine grid RMS residual ---*/
+      su2double rms_res_fine = solver_fine->GetRes_RMS(0);
+      static su2double smoothed_fine_res[MAX_MG_LEVELS] = {};
+
+      if (smoothed_fine_res[lvl] < 1e-30) {
+        smoothed_fine_res[lvl] = rms_res_fine;
+        prev_fine_res[lvl] = rms_res_fine;
+        min_fine_res[lvl] = rms_res_fine;
+      } else {
+        smoothed_fine_res[lvl] = 0.9 * rms_res_fine + 0.1 * smoothed_fine_res[lvl];
+      }
+
+      /*--- Track consecutive decreasing iterations ---*/
+      if (smoothed_fine_res[lvl] < prev_fine_res[lvl]) {
+        consecutive_decreasing[lvl] = min(consecutive_decreasing[lvl] + 1, 10);
+      } else {
+        consecutive_decreasing[lvl] = 0;
+      }
+      prev_fine_res[lvl] = smoothed_fine_res[lvl];
+
+      /*--- Track best residual for divergence detection only ---*/
+      if (smoothed_fine_res[lvl] < min_fine_res[lvl]) {
+        min_fine_res[lvl] = smoothed_fine_res[lvl];
+        iterations_since_improvement[lvl] = 0;
+      } else {
+        iterations_since_improvement[lvl]++;
+      }
+
+      /*--- Detect actual divergence (not stall during slow convergence) ---*/
+      bool fine_grid_diverging = (smoothed_fine_res[lvl] > 1.5 * min_fine_res[lvl] && min_fine_res[lvl] > 1e-30);
+
+      /*--- Sustained convergence: residual decreasing for 5 consecutive iterations ---*/
+      bool sustained_convergence = (consecutive_decreasing[lvl] >= 5);
+
+      /*--- Adapt coefficient based on convergence state ---*/
+      su2double new_coeff = current_coeff;
+
+      if (fine_grid_diverging) {
+        new_coeff = current_coeff * 1.1;
+      } else if (convergence_rate > 1.0) {
+        new_coeff = current_coeff * (1.0 + 0.10 * min(convergence_rate - 1.0, 1.0));
+      } else if (sustained_convergence) {
+        new_coeff = current_coeff * 0.99;
+      }
+
+      /*--- Clamp and enforce hierarchy ---*/
+      new_coeff = min(3.0, max(1.0, new_coeff));
+      su2double min_safe_coeff = CFL_fine / (0.95 * CFL_fine);
+      new_coeff = max(new_coeff, min_safe_coeff);
+
+      /*--- Update coarse grid CFL if changed significantly ---*/
+      su2double CFL_coarse_new = CFL_fine / new_coeff;
+
+      if (fabs(CFL_coarse_new - CFL_coarse_current) > 1e-8) {
+        config->SetCFL(iMesh+1, CFL_coarse_new);
+
+        /*--- Update LocalCFL at each coarse grid point ---*/
+        SU2_OMP_FOR_STAT(roundUpDiv(geometry_coarse->GetnPoint(), omp_get_num_threads()))
+        for (auto iPoint = 0ul; iPoint < geometry_coarse->GetnPoint(); iPoint++) {
+          solver_coarse->GetNodes()->SetLocalCFL(iPoint, CFL_coarse_new);
+        }
+        END_SU2_OMP_FOR
+      }
+
+      /*--- Output monitoring information periodically ---*/
+      if (SU2_MPI::GetRank() == 0 && config->GetTimeIter() % 50 == 0) {
+        cout << "  MG Level " << iMesh+1 << ": conv_rate = " << convergence_rate
+             << ", consecutive = " << consecutive_decreasing[lvl]
+             << ", sustained = " << (sustained_convergence ? "YES" : "NO")
+             << ", coeff = " << new_coeff << ", CFL = " << CFL_coarse_new << endl;
+      }
+    }
+
     SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, config->GetMG_CorrecSmooth(iMesh), 1.25, config);
 
     SetProlongated_Correction(solver_fine, geometry_fine, config, iMesh);
@@ -307,9 +496,6 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
 
 void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
                                                       CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) {
-  unsigned long Point_Fine, Point_Coarse, iVertex;
-  unsigned short iMarker, iChildren, iVar;
-  su2double Area_Parent, Area_Children;
   const su2double *Solution_Fine = nullptr, *Solution_Coarse = nullptr;
 
   const unsigned short nVar = sol_coarse->GetnVar();
@@ -317,41 +503,43 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS
   auto *Solution = new su2double[nVar];
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads()))
-  for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
+  for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
 
-    Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse);
+    su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse);
 
-    for (iVar = 0; iVar < nVar; iVar++) Solution[iVar] = 0.0;
+    for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] = 0.0;
 
-    for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
-      Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
-      Area_Children = geo_fine->nodes->GetVolume(Point_Fine);
+    for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+      auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
+      su2double Area_Children = geo_fine->nodes->GetVolume(Point_Fine);
       Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine);
-      for (iVar = 0; iVar < nVar; iVar++)
+      for (auto iVar = 0u; iVar < nVar; iVar++)
         Solution[iVar] -= Solution_Fine[iVar]*Area_Children/Area_Parent;
     }
 
     Solution_Coarse = sol_coarse->GetNodes()->GetSolution(Point_Coarse);
 
-    for (iVar = 0; iVar < nVar; iVar++)
+    for (auto iVar = 0u; iVar < nVar; iVar++)
       Solution[iVar] += Solution_Coarse[iVar];
 
-    for (iVar = 0; iVar < nVar; iVar++)
+    for (auto iVar = 0u; iVar < nVar; iVar++)
       sol_coarse->GetNodes()->SetSolution_Old(Point_Coarse,Solution);
   }
   END_SU2_OMP_FOR
 
   delete [] Solution;
 
+  /*--- Enforce Euler wall BC on corrections by projecting to tangent plane ---*/
+  ProjectEulerWallToTangentPlane(geo_coarse, config, sol_coarse, true);
+
   /*--- Remove any contributions from no-slip walls. ---*/
 
-  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if (config->GetViscous_Wall(iMarker)) {
 
       SU2_OMP_FOR_STAT(32)
-      for (iVertex = 0; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) {
-
-        Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode();
+      for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) {
+        auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode();
 
         /*--- For dirichlet boundary condtions, set the correction to zero.
          Note that Solution_Old stores the correction not the actual value ---*/
@@ -370,9 +558,9 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS
   sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD);
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads()))
-  for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
-    for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
-      Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
+  for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
+    for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+      auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
       sol_fine->LinSysRes.SetBlock(Point_Fine, sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse));
     }
   }
@@ -387,13 +575,11 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_
   if (val_nSmooth == 0) return;
 
   const su2double *Residual_Old, *Residual_Sum, *Residual_j;
-  unsigned short iVar, iSmooth, iMarker, iNeigh;
-  unsigned long iPoint, jPoint, iVertex;
 
   const unsigned short nVar = solver->GetnVar();
 
   SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads()))
-  for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) {
+  for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); iPoint++) {
     Residual_Old = solver->LinSysRes.GetBlock(iPoint);
     solver->GetNodes()->SetResidual_Old(iPoint,Residual_Old);
   }
@@ -401,17 +587,17 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_
 
   /*--- Jacobi iterations. ---*/
 
-  for (iSmooth = 0; iSmooth < val_nSmooth; iSmooth++) {
+  for (auto iSmooth = 0u; iSmooth < val_nSmooth; iSmooth++) {
 
     /*--- Loop over all mesh points (sum the residuals of direct neighbors). ---*/
 
     SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads()))
-    for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) {
+    for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) {
 
       solver->GetNodes()->SetResidualSumZero(iPoint);
 
-      for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
-        jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+      for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+        auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
         Residual_j = solver->LinSysRes.GetBlock(jPoint);
         solver->GetNodes()->AddResidual_Sum(iPoint, Residual_j);
       }
@@ -422,28 +608,28 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_
     /*--- Loop over all mesh points (update residuals with the neighbor averages). ---*/
 
     SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads()))
-    for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) {
+    for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) {
 
       su2double factor = 1.0/(1.0+val_smooth_coeff*su2double(geometry->nodes->GetnPoint(iPoint)));
 
       Residual_Sum = solver->GetNodes()->GetResidual_Sum(iPoint);
       Residual_Old = solver->GetNodes()->GetResidual_Old(iPoint);
 
-      for (iVar = 0; iVar < nVar; iVar++)
+      for (auto iVar = 0u; iVar < nVar; iVar++)
         solver->LinSysRes(iPoint,iVar) = (Residual_Old[iVar] + val_smooth_coeff*Residual_Sum[iVar])*factor;
     }
     END_SU2_OMP_FOR
 
     /*--- Restore original residuals (without average) at boundary points. ---*/
 
-    for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) {
+    for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++) {
       if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) &&
           (config->GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) &&
           (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) {
 
         SU2_OMP_FOR_STAT(32)
-        for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
-          iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+        for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
+          auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
           Residual_Old = solver->GetNodes()->GetResidual_Old(iPoint);
           solver->LinSysRes.SetBlock(iPoint, Residual_Old);
         }
@@ -457,26 +643,49 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_
 
 void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeometry *geo_fine,
                                                       CConfig *config, unsigned short iMesh) {
-  unsigned long Point_Fine;
-  unsigned short iVar;
   su2double *Solution_Fine, *Residual_Fine;
 
   const unsigned short nVar = sol_fine->GetnVar();
-  const su2double factor = config->GetDamp_Correc_Prolong();
+  //const su2double factor = config->GetDamp_Correc_Prolong();
+
+  /*--- Apply level-dependent damping: more aggressive damping on deeper coarse levels ---*/
+  const su2double base_damp = config->GetDamp_Correc_Prolong();
+  const su2double level_factor = pow(0.75, iMesh);  // 0.75^iMesh reduces damping progressively
+  const su2double factor = base_damp; // * level_factor;
+
+  /*--- Track maximum correction magnitude for monitoring ---*/
+  su2double max_correction_local = 0.0;
+  su2double max_correction_global = 0.0;
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_fine->GetnPointDomain(), omp_get_num_threads()))
-  for (Point_Fine = 0; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) {
+  for (auto Point_Fine = 0ul; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) {
     Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine);
     Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine);
-    for (iVar = 0; iVar < nVar; iVar++) {
+    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;
-      Solution_Fine[iVar] += factor*Residual_Fine[iVar];
+
+      su2double correction = factor * Residual_Fine[iVar];
+      Solution_Fine[iVar] += correction;
+
+     /*--- Track maximum correction ---*/
+      SU2_OMP_CRITICAL
+      max_correction_local = max(max_correction_local, fabs(correction));
     }
   }
   END_SU2_OMP_FOR
 
+  /*--- Reduce maximum correction across all ranks ---*/
+  SU2_MPI::Allreduce(&max_correction_local, &max_correction_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+
+  /*--- Output correction magnitude for monitoring ---*/
+  if (SU2_MPI::GetRank() == 0) {
+    cout << "  MG Level " << iMesh << ": Max correction = " << max_correction_global
+         << ", Damping factor = " << factor << " (base=" << base_damp
+         << " × level_factor=" << level_factor << ")" << endl;
+  }
+
   /*--- MPI the new interpolated solution ---*/
 
   sol_fine->InitiateComms(geo_fine, config, MPI_QUANTITIES::SOLUTION);
@@ -486,13 +695,11 @@ 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) {
-  unsigned long Point_Fine, Point_Coarse;
-  unsigned short iChildren;
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads()))
-  for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
-    for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
-      Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
+  for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
+    for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+      auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
       sol_fine->GetNodes()->SetSolution(Point_Fine, sol_coarse->GetNodes()->GetSolution(Point_Coarse));
     }
   }
@@ -502,8 +709,6 @@ void CMultiGridIntegration::SetProlongated_Solution(unsigned short RunTime_EqSys
 void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine,
                                             CGeometry *geo_coarse, CConfig *config, unsigned short iMesh) {
 
-  unsigned long Point_Fine, Point_Coarse, iVertex;
-  unsigned short iMarker, iVar, iChildren;
   const su2double *Residual_Fine;
 
   const unsigned short nVar = sol_coarse->GetnVar();
@@ -512,16 +717,15 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar
   auto *Residual = new su2double[nVar];
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads()))
-  for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
+  for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
 
     sol_coarse->GetNodes()->SetRes_TruncErrorZero(Point_Coarse);
 
-    for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0;
-
-    for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
-      Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
+    for (auto iVar = 0u; iVar < nVar; iVar++) Residual[iVar] = 0.0;
+    for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+      auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
       Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine);
-      for (iVar = 0; iVar < nVar; iVar++)
+      for (auto iVar = 0u; iVar < nVar; iVar++)
         Residual[iVar] += factor*Residual_Fine[iVar];
     }
     sol_coarse->GetNodes()->AddRes_TruncError(Point_Coarse, Residual);
@@ -530,11 +734,11 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar
 
   delete [] Residual;
 
-  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if (config->GetViscous_Wall(iMarker)) {
       SU2_OMP_FOR_STAT(32)
-      for (iVertex = 0; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) {
-        Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode();
+      for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) {
+        auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode();
         sol_coarse->GetNodes()->SetVel_ResTruncError_Zero(Point_Coarse);
       }
       END_SU2_OMP_FOR
@@ -542,7 +746,7 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar
   }
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads()))
-  for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
+  for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) {
     sol_coarse->GetNodes()->SubtractRes_TruncError(Point_Coarse, sol_coarse->LinSysRes.GetBlock(Point_Coarse));
   }
   END_SU2_OMP_FOR
@@ -553,7 +757,7 @@ void CMultiGridIntegration::SetResidual_Term(CGeometry *geometry, CSolver *solve
 
   AD::StartNoSharedReading();
   SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPointDomain(), omp_get_num_threads()))
-  for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++)
+  for (auto iPoint = 0ul; iPoint < geometry->GetnPointDomain(); iPoint++)
     solver->LinSysRes.AddBlock(iPoint, solver->GetNodes()->GetResTruncError(iPoint));
   END_SU2_OMP_FOR
   AD::EndNoSharedReading();
@@ -561,7 +765,7 @@ void CMultiGridIntegration::SetResidual_Term(CGeometry *geometry, CSolver *solve
 }
 
 void CMultiGridIntegration::SetRestricted_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) {
 
   const unsigned short Solver_Position = config->GetContainerPosition(RunTime_EqSystem);
   const bool grid_movement = config->GetGrid_Movement();
@@ -606,6 +810,9 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst
     }
   }
 
+  /*--- Enforce Euler wall BC by projecting velocity to tangent plane ---*/
+  ProjectEulerWallToTangentPlane(geo_coarse, config, sol_coarse, false);
+
   /*--- MPI the new interpolated solution ---*/
 
   sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION);
@@ -615,39 +822,36 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst
 
 void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
                                                    CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) {
-  unsigned long Point_Fine, Point_Coarse;
-  unsigned short iVar, iDim, iChildren;
-  su2double Area_Parent, Area_Children;
 
   const unsigned short nDim = geo_coarse->GetnDim();
   const unsigned short nVar = sol_coarse->GetnVar();
 
   auto **Gradient = new su2double* [nVar];
-  for (iVar = 0; iVar < nVar; iVar++)
+  for (auto iVar = 0u; iVar < nVar; iVar++)
     Gradient[iVar] = new su2double [nDim];
 
   SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads()))
-  for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) {
-    Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse);
+  for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) {
+    su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse);
 
-    for (iVar = 0; iVar < nVar; iVar++)
-      for (iDim = 0; iDim < nDim; iDim++)
+    for (auto iVar = 0u; iVar < nVar; iVar++)
+      for (auto iDim = 0u; iDim < nDim; iDim++)
         Gradient[iVar][iDim] = 0.0;
 
-    for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
-      Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
-      Area_Children = geo_fine->nodes->GetVolume(Point_Fine);
+    for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+      unsigned long Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren);
+      su2double Area_Children = geo_fine->nodes->GetVolume(Point_Fine);
       auto Gradient_fine = sol_fine->GetNodes()->GetGradient(Point_Fine);
 
-      for (iVar = 0; iVar < nVar; iVar++)
-        for (iDim = 0; iDim < nDim; iDim++)
+      for (auto iVar = 0u; iVar < nVar; iVar++)
+        for (auto iDim = 0u; iDim < nDim; iDim++)
           Gradient[iVar][iDim] += Gradient_fine[iVar][iDim]*Area_Children/Area_Parent;
     }
     sol_coarse->GetNodes()->SetGradient(Point_Coarse,Gradient);
   }
   END_SU2_OMP_FOR
 
-  for (iVar = 0; iVar < nVar; iVar++)
+  for (auto iVar = 0u; iVar < nVar; iVar++)
     delete [] Gradient[iVar];
   delete [] Gradient;
 
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index 6477fa6f361..7113ba92a51 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -739,6 +739,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
 
       const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
       const auto Point_Normal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
+
       /*--- On the finest mesh compute also on halo nodes to avoid communication of tau wall. ---*/
       if ((!geometry->nodes->GetDomain(iPoint)) && !(MGLevel==MESH_0)) continue;
 
diff --git a/TestCases/euler/channel/inv_channel_RK.cfg b/TestCases/euler/channel/inv_channel_RK.cfg
index a719b413d73..bf374054e77 100644
--- a/TestCases/euler/channel/inv_channel_RK.cfg
+++ b/TestCases/euler/channel/inv_channel_RK.cfg
@@ -44,24 +44,24 @@ MARKER_MONITORING= ( upper_wall, lower_wall )
 % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
 %
 NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
-CFL_NUMBER= 2.5
+CFL_NUMBER= 1.0
 CFL_ADAPT= NO
 CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
 RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
 ITER= 110
 LINEAR_SOLVER= BCGSTAB
-LINEAR_SOLVER_ERROR= 1E-6
-LINEAR_SOLVER_ITER= 5
+LINEAR_SOLVER_ERROR= 1E-1
+LINEAR_SOLVER_ITER= 10
 
 % -------------------------- MULTIGRID PARAMETERS -----------------------------%
 %
 MGLEVEL= 3
 MGCYCLE= W_CYCLE
-MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
-MG_POST_SMOOTH= ( 0, 0, 0, 0 )
-MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
-MG_DAMP_RESTRICTION= 0.9
-MG_DAMP_PROLONGATION= 0.9
+MG_PRE_SMOOTH= ( 4, 4, 4, 4 )
+MG_POST_SMOOTH= ( 4, 4, 4, 4 )
+MG_CORRECTION_SMOOTH= ( 4, 4, 4, 4 )
+MG_DAMP_RESTRICTION= 0.5
+MG_DAMP_PROLONGATION= 0.5
 
 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
 %