[WIP] Multigrid with 3D implicit line agglomeration

Common/include/CConfig.hpp

@@ -639,7 +639,8 @@ class CConfig {
   su2double Linear_Solver_Error;   /*!< \brief Min error of the linear solver for the implicit formulation. */
   su2double Deform_Linear_Solver_Error;          /*!< \brief Min error of the linear solver for the implicit formulation. */
   su2double Linear_Solver_Smoother_Relaxation;   /*!< \brief Relaxation factor for iterative linear smoothers. */
-  unsigned long Linear_Solver_Iter;              /*!< \brief Max iterations of the linear solver for the implicit formulation. */
+  unsigned long Linear_Solver_Iter;
+  unsigned long Linear_Solver_Prec_Freeze;  /*!< \brief Reuse the finest-grid preconditioner for this many solves. */              /*!< \brief Max iterations of the linear solver for the implicit formulation. */
   unsigned long Deform_Linear_Solver_Iter;       /*!< \brief Max iterations of the linear solver for the implicit formulation. */
   unsigned long Linear_Solver_Restart_Frequency; /*!< \brief Restart frequency of the linear solver for the implicit formulation. */
   unsigned long Linear_Solver_Restart_Deflation; /*!< \brief Number of vectors used for deflated restarts. */
@@ -1098,6 +1099,7 @@ class CConfig {
   su2double ParMETIS_tolerance;     /*!< \brief Load balancing tolerance for ParMETIS. */
   long ParMETIS_pointWgt;           /*!< \brief Load balancing weight given to points. */
   long ParMETIS_edgeWgt;            /*!< \brief Load balancing weight given to edges. */
+  su2double ParMETIS_anisoWgt;      /*!< \brief Strength of the anisotropy-aware ParMETIS edge weights. 0 disables them. */
   unsigned short DirectDiff;        /*!< \brief Direct Differentation mode. */
   bool DiscreteAdjoint,                /*!< \brief AD-based discrete adjoint mode. */
   DiscreteAdjointDebug;                /*!< \brief Discrete adjoint debug mode using tags. */
@@ -4380,6 +4382,12 @@ class CConfig {
    */
   unsigned long GetLinear_Solver_Iter(void) const { return Linear_Solver_Iter; }
 
+  /*!
+   * \brief Number of consecutive linear solves that reuse one finest-grid preconditioner.
+   * \return Freeze period, 1 meaning rebuild on every solve.
+   */
+  unsigned long GetLinear_Solver_Prec_Freeze(void) const { return Linear_Solver_Prec_Freeze; }
+
   /*!
    * \brief Get max number of iterations of the linear solver for the implicit formulation.
    * \return Max number of iterations of the linear solver for the implicit formulation.
@@ -10166,6 +10174,11 @@ class CConfig {
    */
   long GetParMETIS_EdgeWeight() const { return ParMETIS_edgeWgt; }
 
+  /*!
+   * \brief Get the strength of the anisotropy-aware ParMETIS edge weights (0 disables them).
+   */
+  passivedouble GetParMETIS_AnisoWeight() const { return SU2_TYPE::GetValue(ParMETIS_anisoWgt); }
+
   /*!
    * \brief Find the marker index (if any) that is part of a given interface pair.
    * \param[in] iInterface - Number of the interface pair being tested, starting at 0.

Common/include/linear_algebra/CSysSolve.hpp

@@ -136,6 +136,13 @@ class CSysSolve {
   /*!< \brief Inner solver for nested preconditioning. */
   std::unique_ptr<CSysSolve<ScalarType>> inner_solver;
 
+  /*--- Preconditioner freezing on coarse multigrid levels. The factorization lives in the
+   *    CSysMatrix (not in the short-lived CPreconditioner object built in Solve), so simply
+   *    skipping Build() reuses the previous one. This instance belongs to one solver on one
+   *    grid level, so the counter is naturally per-level. See MG_COARSE_PREC_FREEZE. ---*/
+  mutable unsigned long precSolveCount = 0; /*!< \brief Linear solves done by this instance. */
+  mutable bool buildPrecThisSolve = true;   /*!< \brief Decision for the current solve, shared by all threads. */
+
   /*!
    * \brief sign transfer function
    * \param[in] x - value having sign prescribed

Common/include/option_structure.hpp

@@ -1116,7 +1116,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
 struct CMGOptions {
   su2double MG_Smooth_Res_Threshold{0.0}; /*!< \brief RMS reduction threshold for MG smoothing early exit. */
   su2double MG_Smooth_Coeff{0.0};         /*!< \brief Jacobi smoother coefficient for coarse-grid correction. */
-  unsigned long MG_Min_MeshSize{0};       /*!< \brief Minimum CVs on coarsest MG level. */
+  unsigned long MG_Min_MeshSize{0};       /*!< \brief Minimum CVs on coarsest MG level, per MPI rank. */
   std::vector<unsigned short> MG_PreSmooth;    /*!< \brief Multigrid pre-smoothing iterations per level. */
   std::vector<unsigned short> MG_PostSmooth;   /*!< \brief Multigrid post-smoothing iterations per level. */
   std::vector<unsigned short> MG_CorrecSmooth; /*!< \brief Multigrid Jacobi correction-smoothing per level. */
@@ -1126,6 +1126,16 @@ 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_Implicit_Lines_Max_Group{0}; /*!< \brief Max number of parallel implicit lines merged into one coarse
+                                                      CV tangential to the wall. 0 = dimension-appropriate default
+                                                      (2 in 2D, 4 in 3D). See CMultiGridGeometry::AgglomerateImplicitLines. */
+  su2double MG_Implicit_Lines_Min_AR{2.0}; /*!< \brief Smallest local cell aspect ratio for which a node still counts as
+                                                  part of a stretched layer. Ends a line where the mesh stops being
+                                                  stretched along it, and decides which boundaries carry a layer normal
+                                                  to them. See CMultiGridGeometry::AgglomerateImplicitLines. */
+  unsigned long MG_Coarse_Prec_Freeze{1}; /*!< \brief On MG levels > 0, reuse the linear-solver preconditioner for this many consecutive solves. 1 = rebuild every solve. */
+  su2double MG_Correction_Limit{0.0};     /*!< \brief Max relative change of any solution component from one prolongated FAS correction. 0 = no limit. */
   unsigned long MG_Startup_Iter{100};     /*!< \brief Iterations per mesh during FMG startup, and the length of each level's CFL ramp. 0 = no iteration budget. */
   su2double MG_Startup_Convergence{-2.0}; /*!< \brief FMG: orders of magnitude (log10) that CONV_FIELD must drop on the
                                                  active level before promoting to the next finer one. Negative is a

Common/src/CConfig.cpp

@@ -2065,12 +2065,36 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("MG_SMOOTH_STAGNATION_TOL", MGOptions.MG_Smooth_StagnationTol, 0.99);
   /*!\brief MG_SMOOTH_COEFF\n DESCRIPTION: Smoothing coefficient for the correction prolongation Jacobi smoother. DEFAULT: 1.25 \ingroup Config*/
   addDoubleOption("MG_SMOOTH_COEFF", MGOptions.MG_Smooth_Coeff, 1.25);
-  /*!\brief MG_MIN_MESHSIZE\n DESCRIPTION: Minimum number of CVs on the coarsest multigrid level. Levels that would produce fewer CVs are not created. DEFAULT: 50 \ingroup Config*/
+  /*!\brief MG_COARSE_PREC_FREEZE\n DESCRIPTION: On multigrid levels above MESH_0, reuse the linear-solver preconditioner
+   * (e.g. the ILU factorization) for this many consecutive linear solves instead of rebuilding it every time.
+   * 1 reproduces the previous behaviour exactly. DEFAULT: 1 \ingroup Config*/
+  addUnsignedLongOption("MG_COARSE_PREC_FREEZE", MGOptions.MG_Coarse_Prec_Freeze, 1);
+  /*!\brief LINEAR_SOLVER_PREC_FREEZE\n DESCRIPTION: Reuse the linear-solver preconditioner (e.g. the ILU factorization)
+   * for this many consecutive solves on the finest grid, instead of rebuilding it every time. Applies to MESH_0 and to
+   * single-grid runs; MG_COARSE_PREC_FREEZE covers the coarse levels. The fine-grid preconditioner drives the outer
+   * nonlinear convergence, so raise this one with more care. 1 rebuilds every solve. DEFAULT: 1 \ingroup Config*/
+  addUnsignedLongOption("LINEAR_SOLVER_PREC_FREEZE", Linear_Solver_Prec_Freeze, 1);
+  /*!\brief MG_CORRECTION_LIMIT\n DESCRIPTION: Largest relative change any solution component may undergo from a single
+   * prolongated multigrid correction, e.g. 0.1 caps it at 10%. The whole correction vector at a point is scaled by one
+   * factor so its direction is preserved. 0 disables the limiter (previous behaviour). DEFAULT: 0 \ingroup Config*/
+  addDoubleOption("MG_CORRECTION_LIMIT", MGOptions.MG_Correction_Limit, 0.0);
+  /*!\brief MG_MIN_MESHSIZE\n DESCRIPTION: Minimum number of CVs on the coarsest multigrid level, checked per MPI rank (i.e. on the smallest partition). Levels that would produce fewer CVs on any rank are not created. DEFAULT: 50 \ingroup Config*/
   addUnsignedLongOption("MG_MIN_MESHSIZE", MGOptions.MG_Min_MeshSize, 500);
   /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Enable agglomeration along implicit lines from wall seeds. DEFAULT: NO \ingroup Config*/
   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_IMPLICIT_LINES_MAX_GROUP\n DESCRIPTION: Maximum number of parallel implicit lines merged tangential to
+   * the wall into one coarse CV (2D: always 2; 3D: e.g. 4 for a wall quad/hex corner, 3 for a triangular prism apex).
+   * 0 uses the dimension-appropriate default (2 in 2D, 4 in 3D). DEFAULT: 0 \ingroup Config*/
+  addUnsignedLongOption("MG_IMPLICIT_LINES_MAX_GROUP", MGOptions.MG_Implicit_Lines_Max_Group, 0);
+  /*!\brief MG_IMPLICIT_LINES_MIN_AR\n DESCRIPTION: Smallest local cell aspect ratio for which a node still counts as part
+   * of a stretched layer, measured from the ratio of dual-grid edge weights. Ends an implicit line where the mesh stops
+   * being stretched along it, instead of letting it run to the far field, and decides which boundaries carry a layer
+   * normal to them and may therefore seed lines. 1.0 disables both tests. DEFAULT: 2.0 \ingroup Config*/
+  addDoubleOption("MG_IMPLICIT_LINES_MIN_AR", MGOptions.MG_Implicit_Lines_Min_AR, 2.0);
   /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Max number of iterations spent on each mesh during the Full
    * Multigrid (FMG) startup phase. DEFAULT: 100 \ingroup Config*/
   addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100);
@@ -3053,6 +3077,14 @@ void CConfig::SetConfig_Options() {
 
   /* DESCRIPTION: ParMETIS load balancing weight for edges (equiv. to neighbors) */
   addLongOption("PARMETIS_EDGE_WEIGHT", ParMETIS_edgeWgt, 1);
+  /* DESCRIPTION: Strength of the anisotropy-aware ParMETIS edge weights. ParMETIS is otherwise given no edge weights at
+   * all, so every edge is equally cheap to cut and partition boundaries slice straight through the stretched cells of a
+   * boundary layer, splitting the wall-normal columns that implicit-line agglomeration and line-implicit smoothing rely
+   * on. Weighting an edge by the inverse of its length makes the short wall-normal edges expensive to cut and pushes the
+   * cuts into the tangential direction instead. On a mesh without stretching all edges are of similar length, the
+   * weights come out uniform, and the partitioning is the same as with no weights at all. 0 disables the weights.
+   * DEFAULT: 0 */
+  addDoubleOption("PARMETIS_ANISO_WEIGHT", ParMETIS_anisoWgt, 0.0);
 
   /*--- options that are used in the Hybrid RANS/LES Simulations  ---*/
   /*!\par CONFIG_CATEGORY:Hybrid_RANSLES Options\ingroup Config*/