Introduce 4th-kind Chebyshev preconditioner

diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib
index 4a7e7a3c63b1c82e739b6d8517fe4ad50f2a8985..15424db560ab8fe74633b8cf12e4ed739b3549df 100644 (file)
--- a/doc/doxygen/references.bib
+++ b/doc/doxygen/references.bib
@@ -1950,3 +1950,17 @@ url = {https://doi.org/10.1016/0045-7930(73)90027-3}
   publisher = {arXiv},
   year = {2022},
 }
+
+@article{lottes2022optimal,
+  title={Optimal polynomial smoothers for multigrid V-cycles},
+  author={Lottes, James},
+  journal={arXiv preprint arXiv:2202.08830},
+  year={2022}
+}
+
+@article{phillips2022optimal,
+  title={Optimal Chebyshev Smoothers and One-sided V-cycles},
+  author={Phillips, Malachi and Fischer, Paul},
+  journal={arXiv preprint arXiv:2210.03179},
+  year={2022}
+}

diff --git a/doc/news/changes/minor/20221112KronbichlerMunch b/doc/news/changes/minor/20221112KronbichlerMunch
new file mode 100644 (file)
index 0000000..ac4e283
--- /dev/null
+++ b/doc/news/changes/minor/20221112KronbichlerMunch
@@ -0,0 +1,4 @@
+New: PreconditionChebyshev now also supports
+James Lottes’s fourth-kind Chebyshev.
+<br>
+(Martin Kronbichler, Peter Munch, 2022/11/12)

diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h
index bcf4f473a7e171c06c715c8edf2e02c072f965de..6ff12b63e1426c4880e367338c6da5d1c947ed70 100644 (file)
--- a/include/deal.II/lac/precondition.h
+++ b/include/deal.II/lac/precondition.h
@@ -1966,16 +1966,34 @@ public:
       power_iteration
     };
 
+    /**
+     * An enum to define the available types of polynomial types.
+     */
+    enum class PolynomialType
+    {
+      /**
+       * First-kind Chebyshev polynomials.
+       */
+      first_kind,
+      /**
+       * Fourth-kind Chebyshev polynomials according to @cite lottes2022optimal
+       * and @cite phillips2022optimal.
+       */
+      fourth_kind
+    };
+
     /**
      * Constructor.
      */
-    AdditionalData(const unsigned int        degree              = 1,
-                   const double              smoothing_range     = 0.,
-                   const unsigned int        eig_cg_n_iterations = 8,
-                   const double              eig_cg_residual     = 1e-2,
-                   const double              max_eigenvalue      = 1,
-                   const EigenvalueAlgorithm eigenvalue_algorithm =
-                     EigenvalueAlgorithm::lanczos);
+    AdditionalData(
+      const unsigned int        degree              = 1,
+      const double              smoothing_range     = 0.,
+      const unsigned int        eig_cg_n_iterations = 8,
+      const double              eig_cg_residual     = 1e-2,
+      const double              max_eigenvalue      = 1,
+      const EigenvalueAlgorithm eigenvalue_algorithm =
+        EigenvalueAlgorithm::lanczos,
+      const PolynomialType polynomial_type = PolynomialType::first_kind);
 
     /**
      * Copy assignment operator.
@@ -2046,6 +2064,11 @@ public:
      * Specifies the underlying eigenvalue estimation algorithm.
      */
     EigenvalueAlgorithm eigenvalue_algorithm;
+
+    /**
+     * Specifies the polynomial type to be used.
+     */
+    PolynomialType polynomial_type;
   };
 
 
@@ -3462,13 +3485,15 @@ inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
                                  const unsigned int        eig_cg_n_iterations,
                                  const double              eig_cg_residual,
                                  const double              max_eigenvalue,
-                                 const EigenvalueAlgorithm eigenvalue_algorithm)
+                                 const EigenvalueAlgorithm eigenvalue_algorithm,
+                                 const PolynomialType      polynomial_type)
   : degree(degree)
   , smoothing_range(smoothing_range)
   , eig_cg_n_iterations(eig_cg_n_iterations)
   , eig_cg_residual(eig_cg_residual)
   , max_eigenvalue(max_eigenvalue)
   , eigenvalue_algorithm(eigenvalue_algorithm)
+  , polynomial_type(polynomial_type)
 {}
 
 
@@ -3487,6 +3512,7 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
   max_eigenvalue       = other_data.max_eigenvalue;
   preconditioner       = other_data.preconditioner;
   eigenvalue_algorithm = other_data.eigenvalue_algorithm;
+  polynomial_type      = other_data.polynomial_type;
   constraints.copy_from(other_data.constraints);
 
   return *this;
@@ -3661,7 +3687,10 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 
   const_cast<
     PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
-    ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
+    ->delta =
+    (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
+      (info.max_eigenvalue_estimate) :
+      ((info.max_eigenvalue_estimate - alpha) * 0.5);
   const_cast<
     PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
     ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
@@ -3712,7 +3741,9 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
     rhs,
     0,
     0.,
-    1. / theta,
+    (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
+      (4. / (3. * delta)) :
+      (1. / theta),
     solution,
     solution_old,
     temp_vector1,
@@ -3726,9 +3757,22 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
   double rhok = delta / theta, sigma = theta / delta;
   for (unsigned int k = 0; k < data.degree - 1; ++k)
     {
-      const double rhokp   = 1. / (2. * sigma - rhok);
-      const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
-      rhok = rhokp;
+      double factor1 = 0.0;
+      double factor2 = 0.0;
+
+      if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
+        {
+          factor1 = (2 * k + 1.) / (2 * k + 5.);
+          factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
+        }
+      else
+        {
+          const double rhokp = 1. / (2. * sigma - rhok);
+          factor1            = rhokp * rhok;
+          factor2            = 2. * rhokp / delta;
+          rhok               = rhokp;
+        }
+
       internal::PreconditionChebyshevImplementation::vmult_and_update(
         *matrix_ptr,
         *data.preconditioner,
@@ -3760,7 +3804,9 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
     *data.preconditioner,
     0,
     0.,
-    1. / theta,
+    (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
+      (4. / (3. * delta)) :
+      (1. / theta),
     solution_old,
     temp_vector1,
     temp_vector2,
@@ -3772,9 +3818,22 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
   double rhok = delta / theta, sigma = theta / delta;
   for (unsigned int k = 0; k < data.degree - 1; ++k)
     {
-      const double rhokp   = 1. / (2. * sigma - rhok);
-      const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
-      rhok = rhokp;
+      double factor1 = 0.0;
+      double factor2 = 0.0;
+
+      if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
+        {
+          factor1 = (2 * k + 1.) / (2 * k + 5.);
+          factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
+        }
+      else
+        {
+          const double rhokp = 1. / (2. * sigma - rhok);
+          factor1            = rhokp * rhok;
+          factor2            = 2. * rhokp / delta;
+          rhok               = rhokp;
+        }
+
       matrix_ptr->Tvmult(temp_vector1, solution);
       internal::PreconditionChebyshevImplementation::vector_updates(
         rhs,
@@ -3807,7 +3866,9 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
     rhs,
     1,
     0.,
-    1. / theta,
+    (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
+      (4. / (3. * delta)) :
+      (1. / theta),
     solution,
     solution_old,
     temp_vector1,
@@ -3819,9 +3880,22 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
   double rhok = delta / theta, sigma = theta / delta;
   for (unsigned int k = 0; k < data.degree - 1; ++k)
     {
-      const double rhokp   = 1. / (2. * sigma - rhok);
-      const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
-      rhok = rhokp;
+      double factor1 = 0.0;
+      double factor2 = 0.0;
+
+      if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
+        {
+          factor1 = (2 * k + 1.) / (2 * k + 5.);
+          factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
+        }
+      else
+        {
+          const double rhokp = 1. / (2. * sigma - rhok);
+          factor1            = rhokp * rhok;
+          factor2            = 2. * rhokp / delta;
+          rhok               = rhokp;
+        }
+
       internal::PreconditionChebyshevImplementation::vmult_and_update(
         *matrix_ptr,
         *data.preconditioner,
@@ -3854,7 +3928,9 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
     *data.preconditioner,
     1,
     0.,
-    1. / theta,
+    (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
+      (4. / (3. * delta)) :
+      (1. / theta),
     solution_old,
     temp_vector1,
     temp_vector2,
@@ -3866,9 +3942,22 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
   double rhok = delta / theta, sigma = theta / delta;
   for (unsigned int k = 0; k < data.degree - 1; ++k)
     {
-      const double rhokp   = 1. / (2. * sigma - rhok);
-      const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
-      rhok = rhokp;
+      double factor1 = 0.0;
+      double factor2 = 0.0;
+
+      if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
+        {
+          factor1 = (2 * k + 1.) / (2 * k + 5.);
+          factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
+        }
+      else
+        {
+          const double rhokp = 1. / (2. * sigma - rhok);
+          factor1            = rhokp * rhok;
+          factor2            = 2. * rhokp / delta;
+          rhok               = rhokp;
+        }
+
       matrix_ptr->Tvmult(temp_vector1, solution);
       internal::PreconditionChebyshevImplementation::vector_updates(
         rhs,

diff --git a/tests/lac/precondition_chebyshev_06.cc b/tests/lac/precondition_chebyshev_06.cc
index 53cfb9bb55a360a872426ab18a9801efa018fe23..445fddc67a3b957ab822456071a0d68d866b6cb1 100644 (file)
--- a/tests/lac/precondition_chebyshev_06.cc
+++ b/tests/lac/precondition_chebyshev_06.cc
@@ -62,6 +62,11 @@ main()
         Chebyshev::AdditionalData::EigenvalueAlgorithm::power_iteration;
       cheby.initialize(A, cheby_data);
 
+      Chebyshev cheby_4;
+      cheby_data.polynomial_type =
+        Chebyshev::AdditionalData::PolynomialType::fourth_kind;
+      cheby_4.initialize(A, cheby_data);
+
       Vector<double> v(dim);
       Vector<double> tmp1(dim), tmp2(dim);
       for (unsigned int i = 0; i < 3; ++i)
@@ -79,9 +84,14 @@ main()
           tmp2 -= v;
           const double cheby_residual = tmp2.l2_norm();
 
+          A.vmult(tmp1, v);
+          cheby_4.vmult(tmp2, tmp1);
+          tmp2 -= v;
+          const double cheby_4_residual = tmp2.l2_norm();
+
           deallog << "Residual step i=" << i << ":  "
                   << " ilu=" << ilu_residual << ", cheby=" << cheby_residual
-                  << std::endl;
+                  << ", 4th-kind cheby=" << cheby_4_residual << std::endl;
         }
     }
 

diff --git a/tests/lac/precondition_chebyshev_06.output b/tests/lac/precondition_chebyshev_06.output
index 0143cd8cf0a8e82dccea8c350bffc74bd1353469..cedd2585cb7f02612b79d88d27030f24c9af4171 100644 (file)
--- a/tests/lac/precondition_chebyshev_06.output
+++ b/tests/lac/precondition_chebyshev_06.output
@@ -1,13 +1,13 @@
 
 DEAL::Size 4 Unknowns 9
-DEAL::Residual step i=0:   ilu=0.3321067622, cheby=0.05099432420
-DEAL::Residual step i=1:   ilu=0.3592209788, cheby=0.04600439733
-DEAL::Residual step i=2:   ilu=0.1716913688, cheby=0.02710123523
+DEAL::Residual step i=0:   ilu=0.3321067622, cheby=0.05099432420, 4th-kind cheby=0.03534224422
+DEAL::Residual step i=1:   ilu=0.3592209788, cheby=0.04600439733, 4th-kind cheby=0.03663547728
+DEAL::Residual step i=2:   ilu=0.1716913688, cheby=0.02710123523, 4th-kind cheby=0.01763749678
 DEAL::Size 8 Unknowns 49
-DEAL::Residual step i=0:   ilu=1.943632434, cheby=0.1724610774
-DEAL::Residual step i=1:   ilu=2.148727169, cheby=0.1731956717
-DEAL::Residual step i=2:   ilu=1.756318039, cheby=0.1687993258
+DEAL::Residual step i=0:   ilu=1.943632434, cheby=0.1724610774, 4th-kind cheby=0.09757580669
+DEAL::Residual step i=1:   ilu=2.148727169, cheby=0.1731956717, 4th-kind cheby=0.1005156479
+DEAL::Residual step i=2:   ilu=1.756318039, cheby=0.1687993258, 4th-kind cheby=0.1079517431
 DEAL::Size 16 Unknowns 225
-DEAL::Residual step i=0:   ilu=6.146258348, cheby=0.4251547823
-DEAL::Residual step i=1:   ilu=6.252364591, cheby=0.4316860132
-DEAL::Residual step i=2:   ilu=5.653977586, cheby=0.4063348761
+DEAL::Residual step i=0:   ilu=6.146258348, cheby=0.4251547823, 4th-kind cheby=0.6729443511
+DEAL::Residual step i=1:   ilu=6.252364591, cheby=0.4316860132, 4th-kind cheby=0.6896778102
+DEAL::Residual step i=2:   ilu=5.653977586, cheby=0.4063348761, 4th-kind cheby=0.6284682272