Use the eigenvalue estimation of PreconditionRelaxation (#1064)

Description of the problem The eigenvalue estimation for the smoother of the multigrid preconditioners was done through a PreconditionChebyshev. Description of the solution Due to recent changes in deal.II (dealii/dealii#16742), now we are able to do it using the PreconditionRelaxation directly. This allows us to remove the estimate omega function since it is done internally by the smoother if the relaxation parameter is set to 1. This allows also to remove the degree parameter that was only used in the case of Chebyshev.
chaos-polymtl · Mar 15, 2024 · 30f1a30 · 30f1a30
1 parent bc03310
commit 30f1a30
Show file tree

Hide file tree

Showing 16 changed files with 119 additions and 106 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -19,9 +19,15 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
 ## [Master] - 2024-03-11
 
+## Changed
+
+- MINOR The eigenvalue estimation for the multigrid preconditioners in the lethe-fluid-matrix-free solver is now done internally by the PreconditionRelaxation class instead of using a PreconditionChebyshev and an estimate omega function [#1064](https://github.com/lethe-cfd/lethe/pull/1064).
+
+## [Master] - 2024-03-11
+
 ### Added
 
-- MINOR Temperatue-dependent `stasis constraint` is now featured in the Melting Cavity heat transfer example. [#1061](https://github.com/lethe-cfd/lethe/pull/1061)
+- MINOR Temperature-dependent `stasis constraint` is now featured in the Melting Cavity heat transfer example. [#1061](https://github.com/lethe-cfd/lethe/pull/1061)
 
 ## [Master] - 2024-03-05
 

diff --git a/applications_tests/lethe-fluid-matrix-free/cylinder_kelly_estimator.prm b/applications_tests/lethe-fluid-matrix-free/cylinder_kelly_estimator.prm
@@ -166,7 +166,6 @@ subsection linear solver
     set mg smoother eig estimation = true #if set to true, previous parameter is not used
 
     # Eigenvalue estimation parameters
-    set eig estimation degree          = 3
     set eig estimation smoothing range = 5
     set eig estimation cg n iterations = 20
     set eig estimation verbosity       = quiet

diff --git a/.../cylinder_kelly_estimator.mpirun=1.output → ...mator.with_dealii.geq.9.6.mpirun=1.output b/.../cylinder_kelly_estimator.mpirun=1.output → ...mator.with_dealii.geq.9.6.mpirun=1.output
diff --git a/.../cylinder_kelly_estimator.mpirun=2.output → ...mator.with_dealii.geq.9.6.mpirun=2.output b/.../cylinder_kelly_estimator.mpirun=2.output → ...mator.with_dealii.geq.9.6.mpirun=2.output
diff --git a/...atrix-free/mms3d_fe1_gcmg.mpirun=2.output → ..._gcmg.with_dealii.geq.9.6.mpirun=2.output b/...atrix-free/mms3d_fe1_gcmg.mpirun=2.output → ..._gcmg.with_dealii.geq.9.6.mpirun=2.output
diff --git a/...e-fluid-matrix-free/mms3d_fe1_gcmg.output → ...mms3d_fe1_gcmg.with_dealii.geq.9.6.output b/...e-fluid-matrix-free/mms3d_fe1_gcmg.output → ...mms3d_fe1_gcmg.with_dealii.geq.9.6.output
diff --git a/applications_tests/lethe-fluid-matrix-free/tgv_restart_bdf1.prm b/applications_tests/lethe-fluid-matrix-free/tgv_restart_bdf1.prm
@@ -167,7 +167,6 @@ subsection linear solver
     set mg smoother eig estimation = true
 
     # Eigenvalue estimation parameters
-    set eig estimation degree          = 3
     set eig estimation smoothing range = 5
     set eig estimation cg n iterations = 20
     set eig estimation verbosity       = quiet

diff --git a/...rix-free/tgv_restart_bdf1.mpirun=1.output → ..._bdf1.with_dealii.geq.9.6.mpirun=1.output b/...rix-free/tgv_restart_bdf1.mpirun=1.output → ..._bdf1.with_dealii.geq.9.6.mpirun=1.output
diff --git a/...rix-free/tgv_restart_bdf1.mpirun=2.output → ..._bdf1.with_dealii.geq.9.6.mpirun=2.output b/...rix-free/tgv_restart_bdf1.mpirun=2.output → ..._bdf1.with_dealii.geq.9.6.mpirun=2.output
diff --git a/.../examples/incompressible-flow/3d-taylor-green-vortex/3d-taylor-green-vortex.rst b/.../examples/incompressible-flow/3d-taylor-green-vortex/3d-taylor-green-vortex.rst
@@ -242,7 +242,6 @@ The ``lethe-fluid-matrix-free`` has significantly more parameters for its linear
       set mg smoother eig estimation = true
       
       # Eigenvalue estimation parameters
-      set eig estimation degree          = 3
       set eig estimation smoothing range = 5
       set eig estimation cg n iterations = 20
       set eig estimation verbosity       = verbose

diff --git a/doc/source/parameters/cfd/linear_solver_control.rst b/doc/source/parameters/cfd/linear_solver_control.rst
@@ -198,7 +198,6 @@ Different parameters for the main components of the two geometric multigrid algo
     set mg smoother eig estimation = false #if set to true, previous parameter is not used
 
     # Eigenvalue estimation parameters
-    set eig estimation degree          = 3
     set eig estimation smoothing range = 10
     set eig estimation cg n iterations = 10
     set eig estimation verbosity       = quiet

diff --git a/examples/incompressible-flow/3d-taylor-green-vortex/tgv-matrix-free.prm b/examples/incompressible-flow/3d-taylor-green-vortex/tgv-matrix-free.prm
@@ -138,7 +138,6 @@ subsection linear solver
     set mg smoother eig estimation = true
 
     # Eigenvalue estimation parameters
-    set eig estimation degree          = 3
     set eig estimation smoothing range = 5
     set eig estimation cg n iterations = 10
     set eig estimation verbosity       = verbose

diff --git a/include/core/parameters.h b/include/core/parameters.h
@@ -1119,9 +1119,6 @@ namespace Parameters
     // MG eigenvalue estimation for smoother relaxation parameter
     bool mg_smoother_eig_estimation;
 
-    // MG degree of Chebyshev polynomial used for eigenvalue estimation
-    int eig_estimation_degree;
-
     // MG smoothing range to set range between eigenvalues
     int eig_estimation_smoothing_range;
 

diff --git a/include/solvers/mf_navier_stokes.h b/include/solvers/mf_navier_stokes.h
@@ -214,24 +214,6 @@ class MFNavierStokesSolver
   void
   solve_with_ILU(SolverGMRES<VectorType> &solver);
 
-  /**
-   * @brief Estimate the eigenvalues to obtain a relaxation parameter for the
-   * smoother used in all multigrid preconditioners.
-   *
-   * @param[in] operator Operator for which the estimation needs to be done.
-   *
-   * @param[in] level Corresponding multigrid level.
-   *
-   * @param[in] diagonal Pre-computed diagonal of level operator.
-   *
-   * @return Relaxation parameter omega.
-   */
-  double
-  estimate_omega(
-    std::shared_ptr<NavierStokesOperatorBase<dim, double>> &mg_operator,
-    const unsigned int                                     &level,
-    const VectorType                                       &diagonal);
-
 protected:
   /**
    * @brief Matrix-free operator in used for all the matrix-vector multiplications calls (vmult).

diff --git a/source/core/parameters.cc b/source/core/parameters.cc
@@ -2346,11 +2346,6 @@ namespace Parameters
                           Patterns::Bool(),
                           "estimate eigenvalues for relaxation parameter");
 
-        prm.declare_entry("eig estimation degree",
-                          "3",
-                          Patterns::Integer(),
-                          "degree used for the Chebyshev polynomial");
-
         prm.declare_entry("eig estimation smoothing range",
                           "10",
                           Patterns::Integer(),
@@ -2478,7 +2473,6 @@ namespace Parameters
         mg_smoother_iterations     = prm.get_integer("mg smoother iterations");
         mg_smoother_relaxation     = prm.get_double("mg smoother relaxation");
         mg_smoother_eig_estimation = prm.get_bool("mg smoother eig estimation");
-        eig_estimation_degree      = prm.get_integer("eig estimation degree");
         eig_estimation_smoothing_range =
           prm.get_integer("eig estimation smoothing range");
         eig_estimation_cg_n_iterations =

diff --git a/source/solvers/mf_navier_stokes.cc b/source/solvers/mf_navier_stokes.cc
@@ -447,74 +447,6 @@ MFNavierStokesSolver<dim>::calculate_time_derivative_previous_solutions()
     }
 }
 
-template <int dim>
-double
-MFNavierStokesSolver<dim>::estimate_omega(
-  std::shared_ptr<NavierStokesOperatorBase<dim, double>> &mg_operator,
-  const unsigned int                                     &level,
-  const VectorType                                       &diagonal)
-{
-  TimerOutput::Scope t(this->mg_computing_timer, "Estimate eigenvalues");
-
-  double omega = 0.0;
-
-  using OperatorType               = NavierStokesOperatorBase<dim, double>;
-  using SmootherPreconditionerType = DiagonalMatrix<VectorType>;
-  using ChebyshevPreconditionerType =
-    PreconditionChebyshev<OperatorType, VectorType, SmootherPreconditionerType>;
-  typename ChebyshevPreconditionerType::AdditionalData
-    chebyshev_additional_data;
-
-  chebyshev_additional_data.preconditioner =
-    std::make_shared<SmootherPreconditionerType>(diagonal);
-  chebyshev_additional_data.constraints.copy_from(this->zero_constraints);
-  chebyshev_additional_data.degree =
-    this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
-      .eig_estimation_degree;
-  chebyshev_additional_data.smoothing_range =
-    this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
-      .eig_estimation_smoothing_range;
-  chebyshev_additional_data.eig_cg_n_iterations =
-    this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
-      .eig_estimation_cg_n_iterations;
-  chebyshev_additional_data.eigenvalue_algorithm = ChebyshevPreconditionerType::
-    AdditionalData::EigenvalueAlgorithm::power_iteration;
-  chebyshev_additional_data.polynomial_type =
-    ChebyshevPreconditionerType::AdditionalData::PolynomialType::fourth_kind;
-
-  auto chebyshev = std::make_shared<ChebyshevPreconditionerType>();
-  chebyshev->initialize(*mg_operator, chebyshev_additional_data);
-
-  VectorType vec;
-  mg_operator->initialize_dof_vector(vec);
-
-  const auto evs = chebyshev->estimate_eigenvalues(vec);
-
-  const double alpha =
-    (chebyshev_additional_data.smoothing_range > 1. ?
-       evs.max_eigenvalue_estimate / chebyshev_additional_data.smoothing_range :
-       std::min(0.9 * evs.max_eigenvalue_estimate,
-                evs.min_eigenvalue_estimate));
-
-  omega = 2.0 / (alpha + evs.max_eigenvalue_estimate);
-
-  if (this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
-        .eig_estimation_verbose != Parameters::Verbosity::quiet)
-    {
-      this->pcout << std::endl;
-      this->pcout << "  -Eigenvalue estimation level " << level << ":"
-                  << std::endl;
-      this->pcout << "    Relaxation parameter: " << omega << std::endl;
-      this->pcout << "    Minimum eigenvalue: " << evs.min_eigenvalue_estimate
-                  << std::endl;
-      this->pcout << "    Maximum eigenvalue: " << evs.max_eigenvalue_estimate
-                  << std::endl;
-      this->pcout << std::endl;
-    }
-
-  return omega;
-}
-
 template <int dim>
 void
 MFNavierStokesSolver<dim>::solve_with_LSMG(SolverGMRES<VectorType> &solver)
@@ -740,11 +672,31 @@ MFNavierStokesSolver<dim>::solve_with_LSMG(SolverGMRES<VectorType> &solver)
       smoother_data[level].n_iterations =
         this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
           .mg_smoother_iterations;
+
       if (this->simulation_parameters.linear_solver
             .at(PhysicsID::fluid_dynamics)
             .mg_smoother_eig_estimation)
-        smoother_data[level].relaxation =
-          estimate_omega(mg_operators[level], level, diagonal_vector);
+        {
+#if DEAL_II_VERSION_GTE(9, 6, 0)
+          // Set relaxation to zero so that eigenvalues are estimated internally
+          smoother_data[level].relaxation = 0.0;
+          smoother_data[level].smoothing_range =
+            this->simulation_parameters.linear_solver
+              .at(PhysicsID::fluid_dynamics)
+              .eig_estimation_smoothing_range;
+          smoother_data[level].eig_cg_n_iterations =
+            this->simulation_parameters.linear_solver
+              .at(PhysicsID::fluid_dynamics)
+              .eig_estimation_cg_n_iterations;
+          smoother_data[level].eigenvalue_algorithm =
+            SmootherType::AdditionalData::EigenvalueAlgorithm::power_iteration;
+#else
+          AssertThrow(
+            false,
+            ExcMessage(
+              "The estimation of eigenvalues within LSMG requires a version of deal.II >= 9.6.0"));
+#endif
+        }
       else
         smoother_data[level].relaxation =
           this->simulation_parameters.linear_solver
@@ -754,8 +706,41 @@ MFNavierStokesSolver<dim>::solve_with_LSMG(SolverGMRES<VectorType> &solver)
 
   mg_smoother.initialize(mg_operators, smoother_data);
 
-  this->mg_computing_timer.leave_subsection("Set up and initialize smoother");
+#if DEAL_II_VERSION_GTE(9, 6, 0)
+  if (this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
+        .mg_smoother_eig_estimation &&
+      this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
+          .eig_estimation_verbose != Parameters::Verbosity::quiet)
+    {
+      // Print eigenvalue estimation for all levels
+      for (unsigned int level = minlevel; level <= maxlevel; ++level)
+        {
+          VectorType vec;
+          mg_operators[level]->initialize_dof_vector(vec);
+          const auto evs =
+            mg_smoother.smoothers[level].estimate_eigenvalues(vec);
+
+          this->pcout << std::endl;
+          this->pcout << "  -Eigenvalue estimation level " << level << ":"
+                      << std::endl;
+          this->pcout << "    Relaxation parameter: "
+                      << mg_smoother.smoothers[level].get_relaxation()
+                      << std::endl;
+          this->pcout << "    Minimum eigenvalue: "
+                      << evs.min_eigenvalue_estimate << std::endl;
+          this->pcout << "    Maximum eigenvalue: "
+                      << evs.max_eigenvalue_estimate << std::endl;
+          this->pcout << std::endl;
+        }
+    }
+#else
+  AssertThrow(
+    false,
+    ExcMessage(
+      "The estimation of eigenvalues within LSMG requires a version of deal.II >= 9.6.0"));
+#endif
 
+  this->mg_computing_timer.leave_subsection("Set up and initialize smoother");
 
   // If multigrid number of levels or minimum number of cells in level are
   // specified, change the min level for the coarse-grid solver and the
@@ -1262,11 +1247,31 @@ MFNavierStokesSolver<dim>::solve_with_GCMG(SolverGMRES<VectorType> &solver)
       smoother_data[level].n_iterations =
         this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
           .mg_smoother_iterations;
+
       if (this->simulation_parameters.linear_solver
             .at(PhysicsID::fluid_dynamics)
             .mg_smoother_eig_estimation)
-        smoother_data[level].relaxation =
-          estimate_omega(mg_operators[level], level, diagonal_vector);
+        {
+#if DEAL_II_VERSION_GTE(9, 6, 0)
+          // Set relaxation to zero so that eigenvalues are estimated internally
+          smoother_data[level].relaxation = 0.0;
+          smoother_data[level].smoothing_range =
+            this->simulation_parameters.linear_solver
+              .at(PhysicsID::fluid_dynamics)
+              .eig_estimation_smoothing_range;
+          smoother_data[level].eig_cg_n_iterations =
+            this->simulation_parameters.linear_solver
+              .at(PhysicsID::fluid_dynamics)
+              .eig_estimation_cg_n_iterations;
+          smoother_data[level].eigenvalue_algorithm =
+            SmootherType::AdditionalData::EigenvalueAlgorithm::power_iteration;
+#else
+          AssertThrow(
+            false,
+            ExcMessage(
+              "The estimation of eigenvalues within GCMG requires a version of deal.II >= 9.6.0"));
+#endif
+        }
       else
         smoother_data[level].relaxation =
           this->simulation_parameters.linear_solver
@@ -1276,6 +1281,40 @@ MFNavierStokesSolver<dim>::solve_with_GCMG(SolverGMRES<VectorType> &solver)
 
   mg_smoother.initialize(mg_operators, smoother_data);
 
+#if DEAL_II_VERSION_GTE(9, 6, 0)
+  if (this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
+        .mg_smoother_eig_estimation &&
+      this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics)
+          .eig_estimation_verbose != Parameters::Verbosity::quiet)
+    {
+      // Print eigenvalue estimation for all levels
+      for (unsigned int level = minlevel; level <= maxlevel; ++level)
+        {
+          VectorType vec;
+          mg_operators[level]->initialize_dof_vector(vec);
+          const auto evs =
+            mg_smoother.smoothers[level].estimate_eigenvalues(vec);
+
+          this->pcout << std::endl;
+          this->pcout << "  -Eigenvalue estimation level " << level << ":"
+                      << std::endl;
+          this->pcout << "    Relaxation parameter: "
+                      << mg_smoother.smoothers[level].get_relaxation()
+                      << std::endl;
+          this->pcout << "    Minimum eigenvalue: "
+                      << evs.min_eigenvalue_estimate << std::endl;
+          this->pcout << "    Maximum eigenvalue: "
+                      << evs.max_eigenvalue_estimate << std::endl;
+          this->pcout << std::endl;
+        }
+    }
+#else
+  AssertThrow(
+    false,
+    ExcMessage(
+      "The estimation of eigenvalues within GCMG requires a version of deal.II >= 9.6.0"));
+#endif
+
   this->mg_computing_timer.leave_subsection("Set up and initialize smoother");
 
   // Create coarse-grid GMRES solver and AMG preconditioner