Use the eigenvalue estimation of PreconditionRelaxation

CHANGELOG.md

@@ -5,9 +5,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
 

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

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

doc/source/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

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

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

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;
 

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).

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 =

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