From 64de8b9dc4eb3dd05d48de204b61c6a1da3e6bd2 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 29 Oct 2022 13:27:27 +0200 Subject: [PATCH] Introduce 4th-kind Chebyshev preconditioner --- doc/doxygen/references.bib | 14 ++ .../changes/minor/20221112KronbichlerMunch | 4 + include/deal.II/lac/precondition.h | 139 ++++++++++++++---- tests/lac/precondition_chebyshev_06.cc | 12 +- tests/lac/precondition_chebyshev_06.output | 18 +-- 5 files changed, 152 insertions(+), 35 deletions(-) create mode 100644 doc/news/changes/minor/20221112KronbichlerMunch diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 4a7e7a3c63..15424db560 100644 --- 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 index 0000000000..ac4e283a61 --- /dev/null +++ b/doc/news/changes/minor/20221112KronbichlerMunch @@ -0,0 +1,4 @@ +New: PreconditionChebyshev now also supports +James Lottes’s fourth-kind Chebyshev. +
+(Martin Kronbichler, Peter Munch, 2022/11/12) diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index bcf4f473a7..6ff12b63e1 100644 --- 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:: 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:: 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:: const_cast< PreconditionChebyshev *>(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 *>(this) ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5; @@ -3712,7 +3741,9 @@ PreconditionChebyshev::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::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::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::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::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::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::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::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 53cfb9bb55..445fddc67a 100644 --- 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 v(dim); Vector 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 0143cd8cf0..cedd2585cb 100644 --- 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 -- 2.39.5