[WIP] Feature full multigrid

Common/include/option_structure.hpp

@@ -1126,6 +1126,10 @@ struct CMGOptions {
   su2double MG_Smooth_StagnationTol{0.0}; /*!< \brief Stagnation early exit: stop if current_rms >= prev_rms * tol. 0 = disabled. */
   bool MG_Implicit_Lines{false};          /*!< \brief Enable implicit-lines agglomeration from walls. */
   unsigned long MG_Implicit_Lines_MaxLength{20}; /*!< \brief Maximum nodes on a wall-normal implicit line (including wall seed). */
+  bool MG_Implicit_Lines_Isotropic{false}; /*!< \brief Use isotropic (vs anisotropic) agglomeration along implicit lines. */
+  unsigned long MG_Startup_Iter{100};     /*!< \brief Number of iterations on coarsest mesh during FMG startup phase. */
+  bool TurbMG{false};                     /*!< \brief Run turbulence equations through a FAS MG V-cycle instead of single-grid. */
+  bool MG_Turb_Freeze_Source{false};      /*!< \brief Freeze turbulence source terms on coarse multigrid levels. */
 };
 
 /*!

Common/src/CConfig.cpp

@@ -2069,6 +2069,14 @@ void CConfig::SetConfig_Options() {
   addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false);
   /*!\brief MG_IMPLICIT_LINES_MAX_LENGTH\n DESCRIPTION: Maximum number of nodes on a wall-normal implicit agglomeration line (including the wall seed node). DEFAULT: 20 \ingroup Config*/
   addUnsignedLongOption("MG_IMPLICIT_LINES_MAX_LENGTH", MGOptions.MG_Implicit_Lines_MaxLength, 20);
+  /*!\brief MG_IMPLICIT_LINES_ISOTROPIC\n DESCRIPTION: Use isotropic agglomeration along implicit lines (4 cells per coarse CV) instead of anisotropic (2 cells per coarse CV). DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_IMPLICIT_LINES_ISOTROPIC", MGOptions.MG_Implicit_Lines_Isotropic, false);
+  /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Number of iterations on the coarsest mesh during Full Multigrid (FMG) startup phase before advancing to finer meshes. DEFAULT: 100 \ingroup Config*/
+  addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100);
+  /*!\brief MG_TURB\n DESCRIPTION: Run turbulence equations through a FAS Multigrid V-cycle instead of single-grid. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_TURB", MGOptions.TurbMG, false);
+  /*!\brief MG_TURB_FREEZE_SOURCE\n DESCRIPTION: Freeze turbulence source terms on coarse multigrid levels using values from the fine grid. Reduces stiffness on coarse grids. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_TURB_FREEZE_SOURCE", MGOptions.MG_Turb_Freeze_Source, false);
   /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/
   addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p);
 

SU2_CFD/include/integration/CMultiGridIntegration.hpp

@@ -81,6 +81,18 @@ class CMultiGridIntegration final : public CIntegration {
   void SetForcing_Term(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine,
                        CGeometry *geo_coarse, CConfig *config, unsigned short iMesh);
 
+  /*!
+   * \brief Restrict the fine-grid residual defect to the coarse-grid FAS forcing term.
+   * \param[in] sol_fine - Pointer to the solution on the fine grid.
+   * \param[in] sol_coarse - Pointer to the solution on the coarse grid.
+   * \param[in] geo_fine - Geometrical definition of the fine grid.
+   * \param[in] geo_coarse - Geometrical definition of the coarse grid.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   */
+  void RestrictResidualToCoarseGrid(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine,
+                                    CGeometry *geo_coarse, CConfig *config, unsigned short iMesh);
+
   /*!
    * \brief Add the truncation error to the residual.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -180,6 +192,20 @@ class CMultiGridIntegration final : public CIntegration {
   void GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
                                  CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config);
 
+  /*!
+   * \brief Prolongate the coarse-grid state correction back to the fine-grid residual correction.
+   * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
+   * \param[in] sol_fine - Pointer to the solution on the fine grid.
+   * \param[in] sol_coarse - Pointer to the solution on the coarse grid.
+   * \param[in] geo_fine - Geometrical definition of the fine grid.
+   * \param[in] geo_coarse - Geometrical definition of the coarse grid.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   */
+  void ProlongateCorrectionToFineGrid(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
+                                      CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config,
+                                      unsigned short iMesh);
+
   /*!
    * \brief Do an implicit smoothing of the prolongated correction.
    * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
@@ -205,6 +231,18 @@ class CMultiGridIntegration final : public CIntegration {
   void SetRestricted_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
                               CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config);
 
+  /*!
+   * \brief Restrict frozen turbulence source terms from fine grid to coarse grid.
+   * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
+   * \param[in] sol_fine - Pointer to the solution on the fine grid.
+   * \param[out] sol_coarse - Pointer to the solution on the coarse grid.
+   * \param[in] geo_fine - Geometrical definition of the fine grid.
+   * \param[in] geo_coarse - Geometrical definition of the coarse grid.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetRestricted_FrozenSource(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
+                                  CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config);
+
   /*!
    * \brief Initialize the adjoint solution using the primal problem.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -248,6 +286,26 @@ class CMultiGridIntegration final : public CIntegration {
                         passivedouble lastRMS[2], char& exitReason,
                         passivedouble& worstStepRatio, unsigned short& worstStep);
 
+  /*!
+   * \brief Restrict turbulent eddy viscosity from fine to coarser mesh levels.
+   *
+   * After a turbulence FAS V-cycle completes, this function volume-weights restricts
+   * mu_t from the finest mesh down to all coarser levels. The flow solver on the next
+   * outer iteration uses these restricted mu_t values at every coarse level for the
+   * eddy-viscosity coupling. This ensures consistency between flow and turbulence
+   * solutions across the multigrid hierarchy.
+   *
+   * \param[in] geometry - Geometry hierarchy for one zone/instance (all levels).
+   * \param[in] solver - Solver hierarchy for one zone/instance (all levels).
+   * \param[in] config - Problem configuration.
+   * \param[in] FinestMesh - Current finest active mesh index.
+   * \param[in] nMGLevels - Total number of MG levels.
+   */
+  void RestrictTurbEddyViscToCoarseLevels(CGeometry** geometry, CSolver*** solver,
+                                          CConfig* config,
+                                          unsigned short FinestMesh,
+                                          unsigned short nMGLevels);
+
   static constexpr int MAX_MG_LEVELS = 10;
 
   /*--- Early-exit smoothing state (shared across OMP threads via master write + barrier). ---*/

SU2_CFD/include/integration/CSingleGridIntegration.hpp

@@ -63,6 +63,18 @@ class CSingleGridIntegration final : public CIntegration {
   void SetRestricted_EddyVisc(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
                               CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config);
 
+  /*!
+   * \brief Restrict frozen turbulence source terms from fine grid to coarse grid.
+   * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
+   * \param[in] sol_fine - Pointer to the solution on the fine grid.
+   * \param[out] sol_coarse - Pointer to the solution on the coarse grid.
+   * \param[in] geo_fine - Geometrical definition of the fine grid.
+   * \param[in] geo_coarse - Geometrical definition of the coarse grid.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetRestricted_FrozenSource(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse,
+                                  CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config);
+
 public:
   /*!
    * \brief Constructor of the class.

SU2_CFD/include/solvers/CSolver.hpp

@@ -389,6 +389,16 @@ class CSolver {
    */
   inline su2double GetAvg_CFL_Local(void) const { return Avg_CFL_Local; }
 
+  /*!\
+   * \brief Set min/max/avg local CFL summary statistics.
+   * \param[in] val_cfl - Uniform CFL value to report.
+   */
+  inline void SetCFL_Local_Stats(su2double val_cfl) {
+    Min_CFL_Local = val_cfl;
+    Max_CFL_Local = val_cfl;
+    Avg_CFL_Local = val_cfl;
+  }
+
   /*!
    * \brief Get the number of variables of the problem.
    */

SU2_CFD/include/solvers/CTurbSASolver.hpp

@@ -43,6 +43,25 @@ class CTurbSASolver final : public CTurbSolver {
   su2double nu_tilde_Engine[4] = {0.0};
   su2double nu_tilde_ActDisk[4] = {0.0};
 
+  /*!
+   * \brief Override SetTime_Step to include source term stiffness in the turbulence time step.
+   * \details The base CScalarSolver implementation scales the flow time step by the CFL ratio,
+   *          which only captures convective physics. For the SA model the source term destruction
+   *          dS/d(nu_tilde) ~ c_w*nu_tilde/d^2 dominates near walls (d->0), adding large negative
+   *          diagonal contributions. The frozen source Jacobian (cached each fine-grid iteration)
+   *          captures this stiffness. Limiting dt so that the diagonal is source-dominated rather
+   *          than time-dominated prevents the implicit update from being O(R*dt/V->0), which
+   *          produces r=1.000 (no residual reduction) in the multigrid pre-smoother.
+   *
+   * \param[in] geometry - Geometrical definition.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   * \param[in] Iteration - External iteration number.
+   */
+  void SetTime_Step(CGeometry* geometry, CSolver** solver_container, CConfig* config,
+                    unsigned short iMesh, unsigned long Iteration) override;
+
   /*!
    * \brief A virtual member.
    * \param[in] solver - Solver container

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -38,6 +38,8 @@
 class CTurbVariable : public CScalarVariable {
 protected:
   VectorType muT; /*!< \brief Eddy viscosity. */
+  VectorType frozen_source; /*!< \brief Frozen source term from fine grid for multigrid coarse levels. */
+  VectorType frozen_source_jacobian; /*!< \brief Frozen source Jacobian diagonal for implicit coupling on coarse grids. */
 
 public:
   static constexpr size_t MAXNVAR = 4;
@@ -106,6 +108,34 @@ class CTurbVariable : public CScalarVariable {
    * \param[in] val_DC_kw - diffusion coefficient value
    */
 
+  /*!
+   * \brief Get the frozen source term for multigrid coarse levels.
+   * \param[in] iPoint - Point index.
+   * \return Frozen source term.
+   */
+  inline su2double GetFrozenSource(unsigned long iPoint) const { return frozen_source(iPoint); }
+
+  /*!
+   * \brief Set the frozen source term for multigrid coarse levels.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_source - Frozen source term value.
+   */
+  inline void SetFrozenSource(unsigned long iPoint, su2double val_source) { frozen_source(iPoint) = val_source; }
+
+  /*!
+   * \brief Get the frozen source Jacobian diagonal for multigrid coarse levels.
+   * \param[in] iPoint - Point index.
+   * \return Frozen source Jacobian diagonal.
+   */
+  inline su2double GetFrozenSourceJacobian(unsigned long iPoint) const { return frozen_source_jacobian(iPoint); }
+
+  /*!
+   * \brief Set the frozen source Jacobian diagonal for multigrid coarse levels.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_jac - Frozen source Jacobian diagonal value.
+   */
+  inline void SetFrozenSourceJacobian(unsigned long iPoint, su2double val_jac) { frozen_source_jacobian(iPoint) = val_jac; }
+
   /*!
    * \brief Register eddy viscosity (muT) as Input or Output of an AD recording.
    * \param[in] input - Boolean whether In- or Output should be registered.

SU2_CFD/include/variables/CVariable.hpp

@@ -1778,6 +1778,18 @@ class CVariable {
    */
   inline virtual su2double GetTurbIndex(unsigned long iPoint) const {return 0.0;}
 
+  /*!
+   * \brief Get the frozen turbulence source term density for multigrid.
+   * \return Frozen source density (source per unit volume).
+   */
+  inline virtual su2double GetFrozenSource(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief Set the frozen turbulence source term density for multigrid.
+   * \param[in] val_source - Frozen source density value.
+   */
+  inline virtual void SetFrozenSource(unsigned long iPoint, su2double val_source) {}
+
   /*!
    * \brief A virtual member.
    * \param[in] iVar - Index of the variable.