diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 8e625b719b6e..c477d884e6dc 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -762,6 +762,7 @@ class CConfig {
   Wrt_Performance,           /*!< \brief Write the performance summary at the end of a calculation.  */
   Wrt_AD_Statistics,         /*!< \brief Write the tape statistics (discrete adjoint).  */
   Wrt_MeshQuality,           /*!< \brief Write the mesh quality statistics to the visualization files.  */
+  Wrt_MultiGrid,             /*!< \brief Write the coarse grids to the visualization files.  */
   Wrt_Projected_Sensitivity, /*!< \brief Write projected sensitivities (dJ/dx) on surfaces to ASCII file. */
   Plot_Section_Forces;       /*!< \brief Write sectional forces for specified markers. */
   unsigned short
@@ -3043,6 +3044,11 @@ class CConfig {
    */
   bool GetWrt_MeshQuality(void) const { return Wrt_MeshQuality; }
 
+  /*!
+   * \brief Write coarse grids to the visualization files.
+   */
+  bool GetWrt_MultiGrid(void) const { return Wrt_MultiGrid; }
+
   /*!
    * \brief Get information about writing projected sensitivities on surfaces to an ASCII file with rows as x, y, z, dJ/dx, dJ/dy, dJ/dz for each vertex.
    * \return <code>TRUE</code> means that projected sensitivities on surfaces in an ASCII file with rows as x, y, z, dJ/dx, dJ/dy, dJ/dz for each vertex will be written.
diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp
index 93ed622fa815..8b78e9b8445c 100644
--- a/Common/include/geometry/CGeometry.hpp
+++ b/Common/include/geometry/CGeometry.hpp
@@ -185,6 +185,8 @@ class CGeometry {
   unsigned long edgeColorGroupSize{1};   /*!< \brief Size of the edge groups within each color. */
   unsigned long elemColorGroupSize{1};   /*!< \brief Size of the element groups within each color. */
 
+  ColMajorMatrix<uint8_t> CoarseGridColor_;  /*!< \brief Coarse grid levels, colorized. */
+
 public:
   /*--- Main geometric elements of the grid. ---*/
 
@@ -258,6 +260,8 @@ class CGeometry {
   vector<su2double> Aspect_Ratio;        /*!< \brief Measure of dual CV aspect ratio (max face area / min face area).  */
   vector<su2double> Volume_Ratio;        /*!< \brief Measure of dual CV volume ratio (max sub-element volume / min sub-element volume). */
 
+  const ColMajorMatrix<uint8_t>& CoarseGridColor = CoarseGridColor_;  /*!< \brief Coarse grid levels, colorized. */
+
   /*!
    * \brief Constructor of the class.
    */
@@ -1597,6 +1601,13 @@ class CGeometry {
    */
   inline virtual void ComputeMeshQualityStatistics(const CConfig *config) {}
 
+  /*!
+   * \brief Color multigrid levels for visualization.
+   * \param nMGLevels - Number of levels
+   * \param geometry - The levels
+   */
+  void ColorMGLevels(unsigned short nMGLevels, const CGeometry* const* geometry);
+
   /*!
    * \brief Get the sparse pattern of "type" with given level of fill.
    * \note This method builds the pattern if that has not been done yet.
diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp
index 7f79858541c4..75ab1a052c39 100644
--- a/Common/include/linear_algebra/CSysMatrix.hpp
+++ b/Common/include/linear_algebra/CSysMatrix.hpp
@@ -108,7 +108,7 @@ struct CSysMatrixComms {
 template<class ScalarType>
 class CSysMatrix {
 private:
-  friend class CSysMatrixComms;
+  friend struct CSysMatrixComms;
 
   const int rank;     /*!< \brief MPI Rank. */
   const int size;     /*!< \brief MPI Size. */
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 93b63f40f52b..3e1474ed409c 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -3207,6 +3207,14 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+  /*--- Check if MULTIGRID is requested in VOLUME_OUTPUT and set the config boolean accordingly. ---*/
+  Wrt_MultiGrid = false;
+  for (unsigned short iField = 0; iField < nVolumeOutput; iField++) {
+    if(VolumeOutput[iField].find("MULTIGRID") != string::npos) {
+      Wrt_MultiGrid = true;
+    }
+  }
+
   if (Kind_Solver == NAVIER_STOKES && Kind_Turb_Model != NONE){
     SU2_MPI::Error("KIND_TURB_MODEL must be NONE if SOLVER= NAVIER_STOKES", CURRENT_FUNCTION);
   }
@@ -4996,9 +5004,9 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   if (GetGasModel() == "ARGON") monoatomic = true;
 
   // This option is deprecated. After a grace period until 7.2.0 the usage warning should become an error.
-  if(OptionIsSet("CONV_CRITERIA")) {
-    cout << string("\n\nWARNING: CONV_CRITERIA is deprecated. SU2 will choose the criteria automatically based on the CONV_FIELD.\n") +
-            string("RESIDUAL for any RMS_* BGS_* value. CAUCHY for coefficients like DRAG etc.\n\n");
+  if(OptionIsSet("CONV_CRITERIA") && rank == MASTER_NODE) {
+    cout << "\n\nWARNING: CONV_CRITERIA is deprecated. SU2 will choose the criteria automatically based on the CONV_FIELD.\n"
+            "That is, RESIDUAL for any RMS_* BGS_* value, and CAUCHY for coefficients such as DRAG etc.\n" << endl;
   }
 }
 
diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp
index e238bfeac2ab..175d30bff821 100644
--- a/Common/src/geometry/CGeometry.cpp
+++ b/Common/src/geometry/CGeometry.cpp
@@ -3823,6 +3823,33 @@ void CGeometry::SetNaturalElementColoring()
   if (omp_get_max_threads() > 1) elemColorGroupSize = nElem;
 }
 
+void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const* geometry) {
+
+  using tColor = uint8_t;
+  constexpr auto nColor = numeric_limits<tColor>::max();
+
+  if (nMGLevels) CoarseGridColor_.resize(nPoint,nMGLevels) = 0;
+
+  for (auto iMesh = nMGLevels; iMesh >= 1; --iMesh) {
+    /*--- Color the coarse points. ---*/
+    vector<tColor> color;
+    const auto& adjacency = geometry[iMesh]->nodes->GetPoints();
+    if (colorSparsePattern<CCompressedSparsePatternUL, tColor, nColor>(adjacency, 1, false, &color).empty())
+      continue;
+
+    /*--- Propagate colors to fine mesh. ---*/
+    for (auto step = 0u; step < iMesh; ++step) {
+      auto coarseMesh = geometry[iMesh-1-step];
+      if (step)
+        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)
+          CoarseGridColor_(iPoint,step) = color[coarseMesh->nodes->GetParent_CV(iPoint)];
+    }
+  }
+}
+
 void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry ****geometry_container){
 
   int nZone = config_container[ZONE_0]->GetnZone();
diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp
index 694ff2843638..8e15d607ae1d 100644
--- a/SU2_CFD/include/output/CFlowOutput.hpp
+++ b/SU2_CFD/include/output/CFlowOutput.hpp
@@ -141,4 +141,15 @@ class CFlowOutput : public COutput{
    * \param node_flow
    */
   void LoadTimeAveragedData(unsigned long iPoint, CVariable *node_flow);
+
+  /*!
+   * \brief Add common FVM outputs.
+   */
+  void AddCommonFVMOutputs(const CConfig* config);
+
+  /*!
+   * \brief Load common FVM outputs.
+   */
+  void LoadCommonFVMOutputs(const CConfig* config, const CGeometry* geometry, unsigned long iPoint);
+
 };
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 87c3c2b58f23..c37c6283d924 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -898,6 +898,8 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet
 
   }
 
+  if (config->GetWrt_MultiGrid()) geometry[MESH_0]->ColorMGLevels(config->GetnMGLevels(), geometry);
+
   /*--- For unsteady simulations, initialize the grid volumes
    and coordinates for previous solutions. Loop over all zones/grids ---*/
 
diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp
index 3e8f85e9dc86..077aa5b988f9 100644
--- a/SU2_CFD/src/output/CFlowCompOutput.cpp
+++ b/SU2_CFD/src/output/CFlowCompOutput.cpp
@@ -415,13 +415,7 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion");
   }
 
-  // Mesh quality metrics, computed in CPhysicalGeometry::ComputeMeshQualityStatistics.
-  AddVolumeOutput("ORTHOGONALITY", "Orthogonality", "MESH_QUALITY", "Orthogonality Angle (deg.)");
-  AddVolumeOutput("ASPECT_RATIO",  "Aspect_Ratio",  "MESH_QUALITY", "CV Face Area Aspect Ratio");
-  AddVolumeOutput("VOLUME_RATIO",  "Volume_Ratio",  "MESH_QUALITY", "CV Sub-Volume Ratio");
-
-  // MPI-Rank
-  AddVolumeOutput("RANK", "rank", "MPI", "Rank of the MPI-partition");
+  AddCommonFVMOutputs(config);
 
   if (config->GetTime_Domain()){
     SetTimeAveragedFields();
@@ -567,15 +561,7 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv
     SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(&(Node_Flow->GetGradient_Primitive(iPoint)[1])));
   }
 
-  // Mesh quality metrics
-  if (config->GetWrt_MeshQuality()) {
-    SetVolumeOutputValue("ORTHOGONALITY", iPoint, geometry->Orthogonality[iPoint]);
-    SetVolumeOutputValue("ASPECT_RATIO",  iPoint, geometry->Aspect_Ratio[iPoint]);
-    SetVolumeOutputValue("VOLUME_RATIO",  iPoint, geometry->Volume_Ratio[iPoint]);
-  }
-
-  // MPI-Rank
-  SetVolumeOutputValue("RANK", iPoint, rank);
+  LoadCommonFVMOutputs(config, geometry, iPoint);
 
   if (config->GetTime_Domain()){
     LoadTimeAveragedData(iPoint, Node_Flow);
diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp
index 1a46865d4dbd..8b279613a501 100644
--- a/SU2_CFD/src/output/CFlowIncOutput.cpp
+++ b/SU2_CFD/src/output/CFlowIncOutput.cpp
@@ -498,11 +498,6 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion");
   }
 
-  // Mesh quality metrics, computed in CPhysicalGeometry::ComputeMeshQualityStatistics.
-  AddVolumeOutput("ORTHOGONALITY", "Orthogonality", "MESH_QUALITY", "Orthogonality Angle (deg.)");
-  AddVolumeOutput("ASPECT_RATIO",  "Aspect_Ratio",  "MESH_QUALITY", "CV Face Area Aspect Ratio");
-  AddVolumeOutput("VOLUME_RATIO",  "Volume_Ratio",  "MESH_QUALITY", "CV Sub-Volume Ratio");
-
   // Streamwise Periodicity
   if(streamwisePeriodic) {
     AddVolumeOutput("RECOVERED_PRESSURE", "Recovered_Pressure", "SOLUTION", "Recovered physical pressure");
@@ -510,8 +505,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
       AddVolumeOutput("RECOVERED_TEMPERATURE", "Recovered_Temperature", "SOLUTION", "Recovered physical temperature");
   }
 
-  // MPI-Rank
-  AddVolumeOutput("RANK", "Rank", "MPI", "Rank of the MPI-partition");
+  AddCommonFVMOutputs(config);
 }
 
 void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
@@ -669,15 +663,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
       SetVolumeOutputValue("RECOVERED_TEMPERATURE", iPoint, Node_Flow->GetStreamwise_Periodic_RecoveredTemperature(iPoint));
   }
 
-  // Mesh quality metrics
-  if (config->GetWrt_MeshQuality()) {
-    SetVolumeOutputValue("ORTHOGONALITY", iPoint, geometry->Orthogonality[iPoint]);
-    SetVolumeOutputValue("ASPECT_RATIO",  iPoint, geometry->Aspect_Ratio[iPoint]);
-    SetVolumeOutputValue("VOLUME_RATIO",  iPoint, geometry->Volume_Ratio[iPoint]);
-  }
-
-  // MPI-Rank
-  SetVolumeOutputValue("RANK", iPoint, rank);
+  LoadCommonFVMOutputs(config, geometry, iPoint);
 }
 
 void CFlowIncOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, unsigned short iMarker, unsigned long iVertex){
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 7a6f2992a4b4..815b400b13a4 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -2735,6 +2735,41 @@ bool CFlowOutput::WriteVolume_Output(CConfig *config, unsigned long Iter, bool f
   return false || force_writing;
 }
 
+void CFlowOutput::AddCommonFVMOutputs(const CConfig *config) {
+
+  AddVolumeOutput("ORTHOGONALITY", "Orthogonality", "MESH_QUALITY", "Orthogonality Angle (deg.)");
+  AddVolumeOutput("ASPECT_RATIO",  "Aspect_Ratio",  "MESH_QUALITY", "CV Face Area Aspect Ratio");
+  AddVolumeOutput("VOLUME_RATIO",  "Volume_Ratio",  "MESH_QUALITY", "CV Sub-Volume Ratio");
+
+  AddVolumeOutput("RANK", "rank", "MPI", "Rank of the MPI-partition");
+
+  for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); ++iMesh) {
+    stringstream key, name;
+    key << "MG_" << iMesh;
+    name << "Coarse_Grid_" << iMesh;
+    AddVolumeOutput(key.str(), name.str(), "MULTIGRID", "Coarse mesh");
+  }
+}
+
+void CFlowOutput::LoadCommonFVMOutputs(const CConfig* config, const CGeometry* geometry, unsigned long iPoint) {
+
+  if (config->GetWrt_MeshQuality()) {
+    SetVolumeOutputValue("ORTHOGONALITY", iPoint, geometry->Orthogonality[iPoint]);
+    SetVolumeOutputValue("ASPECT_RATIO",  iPoint, geometry->Aspect_Ratio[iPoint]);
+    SetVolumeOutputValue("VOLUME_RATIO",  iPoint, geometry->Volume_Ratio[iPoint]);
+  }
+
+  SetVolumeOutputValue("RANK", iPoint, rank);
+
+  if (config->GetWrt_MultiGrid()) {
+    for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); ++iMesh) {
+      stringstream key;
+      key << "MG_" << iMesh;
+      SetVolumeOutputValue(key.str(), iPoint, geometry->CoarseGridColor(iPoint,iMesh-1));
+    }
+  }
+}
+
 void CFlowOutput::SetTimeAveragedFields(){
   AddVolumeOutput("MEAN_DENSITY", "MeanDensity", "TIME_AVERAGE", "Mean density");
   AddVolumeOutput("MEAN_VELOCITY-X", "MeanVelocity_x", "TIME_AVERAGE", "Mean velocity x-component");
diff --git a/SU2_CFD/src/output/CNEMOCompOutput.cpp b/SU2_CFD/src/output/CNEMOCompOutput.cpp
index b4709a3de6a3..1aae0b381ac0 100644
--- a/SU2_CFD/src/output/CNEMOCompOutput.cpp
+++ b/SU2_CFD/src/output/CNEMOCompOutput.cpp
@@ -419,6 +419,8 @@ void CNEMOCompOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("VORTICITY_Z", "Vorticity_z", "VORTEX_IDENTIFICATION", "z-component of the vorticity vector");
   }
 
+  AddCommonFVMOutputs(config);
+
   if (config->GetTime_Domain()){
     SetTimeAveragedFields();
   }
@@ -555,6 +557,8 @@ void CNEMOCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv
     SetVolumeOutputValue("ROE_DISSIPATION", iPoint, Node_Flow->GetRoe_Dissipation(iPoint));
   }
 
+  LoadCommonFVMOutputs(config, geometry, iPoint);
+
   if (config->GetTime_Domain()){
     LoadTimeAveragedData(iPoint, Node_Flow);
   }