add convergence check and early exit for pre,post and correction smoothing

Author: bigfooted

Common/include/CConfig.hpp

@@ -496,6 +496,10 @@ class CConfig {
   unsigned short *MG_PreSmooth,       /*!< \brief Multigrid Pre smoothing. */
   *MG_PostSmooth,                     /*!< \brief Multigrid Post smoothing. */
   *MG_CorrecSmooth;                   /*!< \brief Multigrid Jacobi implicit smoothing of the correction. */
+  bool MG_Smooth_EarlyExit; /*!< \brief Enable early exit for MG smoothing. */
+  su2double MG_Smooth_Res_Threshold; /*!< \brief Residual reduction threshold for early exit. */
+  bool MG_Smooth_Output; /*!< \brief Output per-iteration multigrid smoothing info. */
+  su2double MG_Smooth_Coeff; /*!< \brief Smoothing coefficient for multigrid correction smoothing. */
   su2double *LocationStations;        /*!< \brief Airfoil sections in wing slicing subroutine. */
 
   ENUM_MULTIZONE Kind_MZSolver;    /*!< \brief Kind of multizone solver.  */
@@ -3837,6 +3841,30 @@ class CConfig {
     return MG_CorrecSmooth[val_mesh];
   }
 
+  /*!
+   * \brief Get whether early exit is enabled for MG smoothing.
+   * \return True if early exit is enabled.
+   */
+  bool GetMG_Smooth_EarlyExit() const { return MG_Smooth_EarlyExit; }
+
+  /*!
+   * \brief Get the residual threshold for early exit in MG smoothing.
+   * \return Residual threshold.
+   */
+  su2double GetMG_Smooth_Res_Threshold() const { return MG_Smooth_Res_Threshold; }
+
+  /*!
+   * \brief Get whether per-iteration output is enabled for MG smoothing.
+   * \return True if output is enabled.
+   */
+  bool GetMG_Smooth_Output() const { return MG_Smooth_Output; }
+
+  /*!
+   * \brief Get the smoothing coefficient for MG correction smoothing.
+   * \return Smoothing coefficient.
+   */
+  su2double GetMG_Smooth_Coeff() const { return MG_Smooth_Coeff; }
+
   /*!
    * \brief plane of the FFD (I axis) that should be fixed.
    * \param[in] val_index - Index of the arrray with all the planes in the I direction that should be fixed.

Common/src/CConfig.cpp

@@ -1955,6 +1955,17 @@ void CConfig::SetConfig_Options() {
   /*!\brief MG_DAMP_PROLONGATION\n DESCRIPTION: Damping factor for the correction prolongation. DEFAULT 0.75 \ingroup Config*/
   addDoubleOption("MG_DAMP_PROLONGATION", Damp_Correc_Prolong, 0.75);
 
+  /*!\brief MG_SMOOTH_EARLY_EXIT\n DESCRIPTION: Enable early exit for MG smoothing iterations based on RMS residual. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_SMOOTH_EARLY_EXIT", MG_Smooth_EarlyExit, false);
+  /*!\brief MG_SMOOTH_RES_THRESHOLD\n DESCRIPTION: RMS residual threshold for early exit in MG smoothing. DEFAULT: 1e-2 \ingroup Config*/
+  addDoubleOption("MG_SMOOTH_RES_THRESHOLD", MG_Smooth_Res_Threshold, 1e-2);
+
+  /*!\brief MG_SMOOTH_OUTPUT\n DESCRIPTION: Output per-iteration RMS for MG smoothing. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_SMOOTH_OUTPUT", MG_Smooth_Output, false);
+
+  /*!\brief MG_SMOOTH_COEFF\n DESCRIPTION: Smoothing coefficient for MG correction smoothing. DEFAULT: 1.25 \ingroup Config*/
+  addDoubleOption("MG_SMOOTH_COEFF", MG_Smooth_Coeff, 1.25);
+
   /*!\par CONFIG_CATEGORY: Spatial Discretization \ingroup Config*/
   /*--- Options related to the spatial discretization ---*/
 

Common/src/geometry/CMultiGridGeometry.cpp

@@ -610,9 +610,9 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
 
   if (iMesh != MESH_0) {
     //const su2double factor = 1.5; //nijso: too high
-    const su2double factor = 2.0;
+    const su2double factor = 1.1;
     const su2double Coeff = pow(su2double(Global_nPointFine) / Global_nPointCoarse, 1.0 / nDim);
-    const su2double CFL = factor * config->GetCFL(iMesh - 1) / Coeff;
+    const su2double CFL = factor * config->GetCFL(iMesh - 1);// / Coeff;
     config->SetCFL(iMesh, CFL);
   }
 

QuickStart/inv_NACA0012.cfg

@@ -105,7 +105,7 @@ CFL_ADAPT= NO
 CFL_ADAPT_PARAM= ( 0.1, 2.0, 10.0, 1e10 )
 %
 % Number of total iterations
-ITER= 250
+ITER= 5
 
 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
 %
@@ -136,13 +136,22 @@ MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
 MG_POST_SMOOTH= ( 0, 0, 0, 0 )
 %
 % Jacobi implicit smoothing of the correction
-MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
+MG_CORRECTION_SMOOTH= ( 0, 3, 3, 3 )
 %
 % Damping factor for the residual restriction
 MG_DAMP_RESTRICTION= 1.0
 %
 % Damping factor for the correction prolongation
 MG_DAMP_PROLONGATION= 1.0
+%
+% Enable early exit for multigrid smoothing iterations
+MG_SMOOTH_EARLY_EXIT= YES
+%
+% RMS residual threshold for early exit (fraction of initial RMS)
+MG_SMOOTH_RES_THRESHOLD= 0.99
+%
+% Enable output of RMS residual during smoothing iterations
+MG_SMOOTH_OUTPUT= YES
 
 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
 %

SU2_CFD/include/integration/CMultiGridIntegration.hpp

@@ -177,7 +177,7 @@ class CMultiGridIntegration final : public CIntegration {
    * \param[in] config - Definition of the particular problem.
    */
   void SmoothProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *solver, CGeometry *geometry,
-                                    unsigned short val_nSmooth, su2double val_smooth_coeff, CConfig *config);
+                                    unsigned short val_nSmooth, su2double val_smooth_coeff, CConfig *config, unsigned short iMesh);
 
   /*!
    * \brief Restrict solution from fine grid to a coarse grid.

SU2_CFD/src/integration/CMultiGridIntegration.cpp

@@ -28,6 +28,8 @@
 #include "../../include/integration/CMultiGridIntegration.hpp"
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 
+#include "ComputeLinSysResRMS.hpp"
+
 
 CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { }
 
@@ -217,7 +219,7 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
     /*--- Compute prolongated solution, and smooth the correction $u^(new)_k = u_k +  Smooth(I^k_(k+1)(u_(k+1)-I^(k+1)_k u_k))$ ---*/
 
     GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config);
-    SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, config->GetMG_CorrecSmooth(iMesh), 1.25, config);
+    SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, config->GetMG_CorrecSmooth(iMesh), config->GetMG_Smooth_Coeff(), config, iMesh);
     SetProlongated_Correction(solver_fine, geometry_fine, config, iMesh);
 
 
@@ -245,43 +247,52 @@ unsigned short iRKLimit) {
   const unsigned short nPreSmooth = config->GetMG_PreSmooth(iMesh);
   const unsigned long timeIter = config->GetTimeIter();
 
+  // Early exit settings from config
+  const bool early_exit_enabled = config->GetMG_Smooth_EarlyExit();
+  const su2double early_exit_threshold = config->GetMG_Smooth_Res_Threshold();
+  const bool output_enabled = config->GetMG_Smooth_Output();
+
+  su2double initial_rms = 0.0;
+  if (early_exit_enabled || output_enabled) {
+    initial_rms = ComputeLinSysResRMS(solver_fine, geometry_fine);
+    if (output_enabled) {
+      cout << "MG Pre-Smoothing Level " << iMesh << " Initial RMS: " << initial_rms << endl;
+    }
+  }
+
   /*--- Do a presmoothing on the grid iMesh to be restricted to the grid iMesh+1 ---*/
   for (unsigned short iPreSmooth = 0; iPreSmooth < nPreSmooth; iPreSmooth++) {
-
     /*--- Time and space integration ---*/
     for (unsigned short iRKStep = 0; iRKStep < iRKLimit; iRKStep++) {
-
       /*--- Send-Receive boundary conditions, and preprocessing ---*/
       solver_fine->Preprocessing(geometry_fine, solver_container_fine, config, iMesh, iRKStep, RunTime_EqSystem, false);
-
       if (iRKStep == 0) {
-
         /*--- Set the old solution ---*/
         solver_fine->Set_OldSolution();
-
         if (classical_rk4) solver_fine->Set_NewSolution();
-        // nijso asks: only call when classical_rk4?
-        // this copies solution to old solution, which we already did above.
-        //solver_fine->Set_OldSolution();
-
-        /*--- Compute time step, max eigenvalue, and integration scheme (steady and unsteady problems) ---*/
         solver_fine->SetTime_Step(geometry_fine, solver_container_fine, config, iMesh, timeIter);
-
-        /*--- Restrict the solution and gradient for the adjoint problem ---*/
-        // nijso asks: why do we call this here but not in the post-smoothing?
         Adjoint_Setup(geometry, solver_container, config_container, RunTime_EqSystem, timeIter, iZone);
-
       }
-      //cout << "pre-smoothing mesh " << iMesh << " rkstep " << iRKStep << endl;
       /*--- Space integration ---*/
       Space_Integration(geometry_fine, solver_container_fine, numerics_fine, config, iMesh, iRKStep, RunTime_EqSystem);
-
       /*--- Time integration, update solution using the old solution plus the solution increment ---*/
       Time_Integration(geometry_fine, solver_container_fine, config, iRKStep, RunTime_EqSystem);
-
       /*--- Send-Receive boundary conditions, and postprocessing ---*/
       solver_fine->Postprocessing(geometry_fine, solver_container_fine, config, iMesh);
+    }
 
+    // Early exit check and output
+    if (early_exit_enabled || output_enabled) {
+      su2double current_rms = ComputeLinSysResRMS(solver_fine, geometry_fine);
+      if (output_enabled) {
+        cout << "MG Pre-Smoothing Level " << iMesh << " Iteration " << iPreSmooth + 1 << "/" << nPreSmooth << " RMS: " << current_rms << endl;
+      }
+      if (early_exit_enabled && current_rms < early_exit_threshold * initial_rms) {
+        if (output_enabled) {
+          cout << "MG Pre-Smoothing Level " << iMesh << " Early exit at iteration " << iPreSmooth + 1 << endl;
+        }
+        break;
+      }
     }
   }
 }
@@ -293,32 +304,49 @@ void CMultiGridIntegration::PostSmoothing(unsigned short RunTime_EqSystem, CSolv
   const unsigned short nPostSmooth = config->GetMG_PostSmooth(iMesh);
   const unsigned long timeIter = config->GetTimeIter();
 
+  // Early exit settings from config
+  const bool early_exit_enabled = config->GetMG_Smooth_EarlyExit();
+  const su2double early_exit_threshold = config->GetMG_Smooth_Res_Threshold();
+  const bool output_enabled = config->GetMG_Smooth_Output();
+
+  su2double initial_rms = 0.0;
+  if (early_exit_enabled || output_enabled) {
+    initial_rms = ComputeLinSysResRMS(solver_fine, geometry_fine);
+    if (output_enabled) {
+      cout << "MG Post-Smoothing Level " << iMesh << " Initial RMS: " << initial_rms << endl;
+    }
+  }
+
   /*--- Do a postsmoothing on the grid iMesh after prolongation from the grid iMesh+1 ---*/
   for (unsigned short iPostSmooth = 0; iPostSmooth < nPostSmooth; iPostSmooth++) {
-
     for (unsigned short iRKStep = 0; iRKStep < iRKLimit; iRKStep++) {
-
       solver_fine->Preprocessing(geometry_fine, solver_container_fine, config, iMesh, iRKStep, RunTime_EqSystem, false);
-
       if (iRKStep == 0) {
-
         /*--- Set the old solution ---*/
         solver_fine->Set_OldSolution();
-
         if (classical_rk4) solver_fine->Set_NewSolution();
-          solver_fine->SetTime_Step(geometry_fine, solver_container_fine, config, iMesh,  timeIter);
-
+        solver_fine->SetTime_Step(geometry_fine, solver_container_fine, config, iMesh,  timeIter);
       }
-
       /*--- Space integration ---*/
       Space_Integration(geometry_fine, solver_container_fine, numerics_fine, config, iMesh, iRKStep, RunTime_EqSystem);
-      //cout << "post-smoothing mesh " << iMesh << " rkstep " << iRKStep << endl;
       /*--- Time integration, update solution using the old solution plus the solution increment ---*/
       Time_Integration(geometry_fine, solver_container_fine, config, iRKStep, RunTime_EqSystem);
-
       /*--- Send-Receive boundary conditions, and postprocessing ---*/
       solver_fine->Postprocessing(geometry_fine, solver_container_fine, config, iMesh);
+    }
 
+    // Early exit check and output
+    if (early_exit_enabled || output_enabled) {
+      su2double current_rms = ComputeLinSysResRMS(solver_fine, geometry_fine);
+      if (output_enabled) {
+        cout << "MG Post-Smoothing Level " << iMesh << " Iteration " << iPostSmooth + 1 << "/" << nPostSmooth << " RMS: " << current_rms << endl;
+      }
+      if (early_exit_enabled && current_rms < early_exit_threshold * initial_rms) {
+        if (output_enabled) {
+          cout << "MG Post-Smoothing Level " << iMesh << " Early exit at iteration " << iPostSmooth + 1 << endl;
+        }
+        break;
+      }
     }
   }
 }
@@ -400,7 +428,7 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS
 }
 
 void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *solver, CGeometry *geometry,
-                                                         unsigned short val_nSmooth, su2double val_smooth_coeff, CConfig *config) {
+                                                         unsigned short val_nSmooth, su2double val_smooth_coeff, CConfig *config, unsigned short iMesh) {
 
   /*--- Check if there is work to do. ---*/
   if (val_nSmooth == 0) return;
@@ -418,7 +446,19 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_
   }
   END_SU2_OMP_FOR
 
-  /*--- Jacobi iterations. ---*/
+  /*--- Compute initial RMS for adaptive smoothing ---*/
+  su2double initial_rms = 0.0;
+  if (config->GetMG_Smooth_EarlyExit()) {
+    initial_rms = ComputeLinSysResRMS(solver, geometry);
+  }
+
+  /*--- Output initial RMS if enabled ---*/
+  if (config->GetMG_Smooth_Output()) {
+    cout << "MG Correction-Smoothing Level " << iMesh << " Initial RMS: " << initial_rms << endl;
+  }
+
+  /*--- Jacobi iterations with adaptive early exit ---*/
+  unsigned short actual_iterations = val_nSmooth;
 
   for (iSmooth = 0; iSmooth < val_nSmooth; iSmooth++) {
 
@@ -470,6 +510,35 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_
       }
     }
 
+    /*--- Output RMS residual if enabled ---*/
+    if (config->GetMG_Smooth_Output()) {
+      const su2double RMS_Res = ComputeLinSysResRMS(solver, geometry);
+      cout << "MG Correction-Smoothing Level " << iMesh << " Iteration " << iSmooth+1 << "/" << val_nSmooth << " RMS: " << RMS_Res << endl;
+    }
+
+    /*--- Adaptive early exit check after first iteration ---*/
+    if (config->GetMG_Smooth_EarlyExit() && iSmooth == 0 && val_nSmooth > 1) {
+      su2double current_rms = ComputeLinSysResRMS(solver, geometry);
+      su2double reduction_ratio = current_rms / initial_rms;
+
+      // If RMS reduction is sufficient (ratio <= threshold), additional iterations may not be necessary
+      if (reduction_ratio <= config->GetMG_Smooth_Res_Threshold()) {
+        if (config->GetMG_Smooth_Output()) {
+          cout << "MG Correction-Smoothing Level " << iMesh << " Early exit: sufficient RMS reduction ("
+               << reduction_ratio << " <= " << config->GetMG_Smooth_Res_Threshold() << ")" << endl;
+        }
+        actual_iterations = 1;  // Only do this one iteration
+        break;
+      }
+      // If reduction is insufficient (ratio > threshold), continue with remaining iterations
+    }
+
+  }
+
+  /*--- Log if iterations were skipped ---*/
+  if (config->GetMG_Smooth_Output() && actual_iterations < val_nSmooth) {
+    cout << "MG Correction-Smoothing Level " << iMesh << " completed " << actual_iterations
+         << "/" << val_nSmooth << " iterations (adaptive early exit)" << endl;
   }
 
 }

SU2_CFD/src/integration/ComputeLinSysResRMS.hpp

@@ -0,0 +1,26 @@
+#pragma once
+#include "../../../Common/include/geometry/CGeometry.hpp"
+#include "../../include/solvers/CSolver.hpp"
+#include <vector>
+#include <cmath>
+#include "../../../Common/include/CConfig.hpp"
+
+inline su2double ComputeLinSysResRMS(const CSolver* solver, const CGeometry* geometry) {
+    unsigned short nVar = solver->GetnVar();
+    unsigned long nPointDomain = geometry->GetnPointDomain();
+    std::vector<su2double> sumRes(nVar, 0.0);
+    su2double localSum = 0.0;
+    for (unsigned long i = 0; i < nPointDomain; ++i) {
+        const su2double* res = solver->LinSysRes.GetBlock(i);
+        for (unsigned short v = 0; v < nVar; ++v) {
+            sumRes[v] += res[v] * res[v];
+        }
+    }
+    for (unsigned short v = 0; v < nVar; ++v) localSum += sumRes[v];
+    su2double globalSum = 0.0;
+    SU2_MPI::Allreduce(&localSum, &globalSum, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    unsigned long globalNPointDomain = 0;
+    SU2_MPI::Allreduce(&nPointDomain, &globalNPointDomain, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+    if (globalNPointDomain == 0) return 0.0;
+    return std::sqrt(globalSum / (globalNPointDomain * nVar));
+}

SU2_CFD/src/solvers/CSolver.cpp

@@ -1955,11 +1955,12 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
           const su2double nPointFine = su2double(geometry[iMGLevel-1]->GetGlobal_nPointDomain());
 
           /*--- Apply the same scaling factor as used during geometry construction. ---*/
-          const su2double factor = 2.0;
+          const su2double factor = 1.1;
           const su2double Coeff = pow(nPointFine / nPointCoarse, 1.0 / nDim);
-          const su2double CFL_Coarse = factor * config->GetCFL(iMGLevel - 1) / Coeff;
-
-          config->SetCFL(iMGLevel, CFL_Coarse);
+          const su2double CFL_Coarse = factor * config->GetCFL(iMGLevel - 1);// / Coeff;
+          cout << "Adapted CFL for MG Level " << iMGLevel << ": " << CFL_Coarse << endl;
+          // nijso: comment this if we do not adapt the coarse level CFLs
+          //config->SetCFL(iMGLevel, CFL_Coarse);
         }
       }
       END_SU2_OMP_SAFE_GLOBAL_ACCESS

config_template.cfg

@@ -1597,6 +1597,9 @@ MG_POST_SMOOTH= ( 0, 0, 0, 0 )
 % Jacobi implicit smoothing of the correction
 MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
 %
+% Smoothing coefficient for multigrid correction smoothing (1.25 by default)
+MG_SMOOTH_COEFF= 1.25
+%
 % Damping factor for the residual restriction
 MG_DAMP_RESTRICTION= 0.75
 %