From 91608d7dd8265c31a13d8772a46a7fa35d2879f9 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 31 Jul 2024 20:47:25 +0200 Subject: [PATCH] Add different variants of p-multgrid (#1209) Co-authored-by: Laura Prieto Saavedra Co-authored-by: Laura Prieto Saavedra <40216050+lpsaavedra@users.noreply.github.com> Co-authored-by: Bruno Blais --- CHANGELOG.md | 7 ++ .../parameters/cfd/linear_solver_control.rst | 35 +++++++- include/core/parameters.h | 16 ++++ source/core/parameters.cc | 47 ++++++++++ source/solvers/mf_navier_stokes.cc | 88 +++++++++++++++++-- 5 files changed, 182 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0fa7391bbc..68f25860a6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,13 @@ All notable changes to the Lethe project will be documented in this file. The format is based on [Keep a Changelog](http://keepachangelog.com/). +## [Master] - 2024-07-31 + +### Added + +- MINOR P-multigrid was added to the gcmg preconditioner of the lethe-fluid-matrix-free application. It supports three different coarsening strategies to define the degree p of the different levels. It also allows to use hybrid hp- and ph-multigrid strategies. [#1209](https://github.com/chaos-polymtl/lethe/pull/1209) + + ## [Master] - 2024-07-24 ### Changed diff --git a/doc/source/parameters/cfd/linear_solver_control.rst b/doc/source/parameters/cfd/linear_solver_control.rst index 9054af02e5..347636886a 100644 --- a/doc/source/parameters/cfd/linear_solver_control.rst +++ b/doc/source/parameters/cfd/linear_solver_control.rst @@ -190,9 +190,15 @@ AMG preconditioner .. seealso:: For more information about the ``amg`` preconditioner parameters, the reader is referred to the deal.II documentation for the `AMG preconditioner `_ and its `Additional Data `_. ------------------------------- -LSMG and GCMG preconditioners ------------------------------- +-------------------------- +Multigrid preconditioners +-------------------------- + +Lethe supports two types of geometric multigrid preconditioners that only differ when dealing with locally-refined meshes: + +* Global coarsening ``gcmg``: coarsens all cells simultaneously, i.e., each level contains all the cells at their most refined state. + +* Local smoothing ``lsmg``: uses the refinement hierarchy to create the multigrid levels and to perform smoothing refinement level by refinement level, i.e., cells of less refined parts of the mesh are skipped. Different parameters for the main components of the two geometric multigrid algorithms can be specified. The parameters can be general or can belong to either the smoother or the coarse-grid solver. Lethe supports different coarse-grid solvers: ``gmres``, ``amg``, ``ilu`` and ``direct``. The ``gmres`` coarse-grid solver supports two preconditioners ``amg`` and ``ilu``. @@ -251,4 +257,25 @@ Different parameters for the main components of the two geometric multigrid algo If your coarse-grid level is small enough, it might be worth it for some problems to set ``mg amg use default parameters = true`` to use a direct solver. On the other hand, if high order elements are used, it might be useful to set ``set mg coarse grid use fe q iso q1 = true`` to solve the coarse grid problem using `FE_Q_iso_Q1 elements `_. .. tip:: - Evaluating terms involving the hessian is expensive. Therefore, one can turn on or off those terms in the mg level operators to improve performance by setting ``mg enable hessians in jacobian`` to ``false``. This is useful for certain problems and must be used carefully. \ No newline at end of file + Evaluating terms involving the hessian is expensive. Therefore, one can turn on or off those terms in the mg level operators to improve performance by setting ``mg enable hessians in jacobian`` to ``false``. This is useful for certain problems and must be used carefully. + +In addition, Lethe supports `p-multigrid` through the ``gcmg`` preconditioner. It can be used by specifying two additional parameters: + +.. code-block:: text + + set mg coarsening type = p + set mg p coarsening type = decrease by one + +This multigrid preconditioner creates the different multigrid levels by keeping the mesh constant but reducing the polynomial degree `p` of the shape functions. Three strategies to create the `p-multigrid` levels can be used by specifying the ``mg p coarsening type`` parameter: + +* ``bisect``: half polynomial degree. + +* ``decrease by one``: decrease the polynomial degree by one for every level. + +* ``go to one``: decrease the polynomial degree to one directly. + +In addition, Lethe supports hybrid strategies that combine h- and p-multigrid, and can be specified through the ``mg coarsening type`` parameter: + +* ``hp``: first levels with different mesh and then levels with different degree `p`. + +* ``ph``: first levels with different degree `p` and then levels with different mesh. \ No newline at end of file diff --git a/include/core/parameters.h b/include/core/parameters.h index d761eeb4fa..a4a190a78f 100644 --- a/include/core/parameters.h +++ b/include/core/parameters.h @@ -35,6 +35,8 @@ #include #include +#include + using namespace dealii; @@ -1184,6 +1186,20 @@ namespace Parameters /// MG enable hessians in jacobian bool mg_enable_hessians_jacobian; + /// Type of multigrid + enum class MultigridCoarseningSequenceType + { + h, + p, + hp, + ph + }; + MultigridCoarseningSequenceType mg_coarsening_type; + + /// Type of p coarsening sequence + MGTransferGlobalCoarseningTools::PolynomialCoarseningSequenceType + mg_p_coarsening_type; + /// MG smoother number of iterations int mg_smoother_iterations; diff --git a/source/core/parameters.cc b/source/core/parameters.cc index 7d3825b0e1..32900db1ca 100644 --- a/source/core/parameters.cc +++ b/source/core/parameters.cc @@ -2550,6 +2550,17 @@ namespace Parameters Patterns::Bool(), "use elements with linear interpolation for coarse grid"); + prm.declare_entry("mg coarsening type", + "h", + Patterns::Selection("h|p|hp|ph"), + "mg coarsening type for gcmg"); + + prm.declare_entry("mg p coarsening type", + "decrease by one", + Patterns::Selection( + "decrease by one|bisect|go to one"), + "mg p coarsening type for gcmg"); + prm.declare_entry("mg gmres max iterations", "2000", Patterns::Integer(), @@ -2706,6 +2717,42 @@ namespace Parameters mg_use_fe_q_iso_q1 = prm.get_bool("mg coarse grid use fe q iso q1"); + const std::string mg_coarsening_type = prm.get("mg coarsening type"); + if (mg_coarsening_type == "h") + this->mg_coarsening_type = MultigridCoarseningSequenceType::h; + else if (mg_coarsening_type == "p") + this->mg_coarsening_type = MultigridCoarseningSequenceType::p; + else if (mg_coarsening_type == "hp") + this->mg_coarsening_type = MultigridCoarseningSequenceType::hp; + else if (mg_coarsening_type == "ph") + this->mg_coarsening_type = MultigridCoarseningSequenceType::ph; + else + AssertThrow(false, ExcNotImplemented()); + + const std::string mg_p_coarsening_type = + prm.get("mg p coarsening type"); + if (mg_p_coarsening_type == "decrease by one") + this->mg_p_coarsening_type = MGTransferGlobalCoarseningTools:: + PolynomialCoarseningSequenceType::decrease_by_one; + else if (mg_p_coarsening_type == "bisect") + this->mg_p_coarsening_type = MGTransferGlobalCoarseningTools:: + PolynomialCoarseningSequenceType::bisect; + else if (mg_p_coarsening_type == "go to one") + this->mg_p_coarsening_type = MGTransferGlobalCoarseningTools:: + PolynomialCoarseningSequenceType::go_to_one; + else + AssertThrow(false, ExcNotImplemented()); + + AssertThrow((!mg_use_fe_q_iso_q1) || + (this->mg_coarsening_type == + MultigridCoarseningSequenceType::h), + ExcNotImplemented()); + + AssertThrow((preconditioner != PreconditionerType::lsmg) || + (this->mg_coarsening_type == + MultigridCoarseningSequenceType::h), + ExcNotImplemented()); + mg_gmres_max_iterations = prm.get_integer("mg gmres max iterations"); mg_gmres_tolerance = prm.get_double("mg gmres tolerance"); mg_gmres_reduce = prm.get_double("mg gmres reduce"); diff --git a/source/solvers/mf_navier_stokes.cc b/source/solvers/mf_navier_stokes.cc index dc0b25005e..b86c10468f 100644 --- a/source/solvers/mf_navier_stokes.cc +++ b/source/solvers/mf_navier_stokes.cc @@ -535,6 +535,10 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( .preconditioner == Parameters::LinearSolver::PreconditionerType::gcmg) { + AssertThrow(*cell_quadrature == + QGauss(this->dof_handler.get_fe().degree + 1), + ExcNotImplemented()); + // Create triangulations this->mg_setup_timer.enter_subsection("Create level triangulations"); this->coarse_grid_triangulations = @@ -599,10 +603,68 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( this->coarse_grid_triangulations = temp; + // p-multigrid + const auto mg_coarsening_type = + this->simulation_parameters.linear_solver.at(PhysicsID::fluid_dynamics) + .mg_coarsening_type; + + const auto polynomial_coarsening_sequence = + MGTransferGlobalCoarseningTools::create_polynomial_coarsening_sequence( + this->dof_handler.get_fe().degree, + this->simulation_parameters.linear_solver + .at(PhysicsID::fluid_dynamics) + .mg_p_coarsening_type); + + std::vector> levels; + + if (mg_coarsening_type == + Parameters::LinearSolver::MultigridCoarseningSequenceType::hp) + { + // p + for (const auto i : polynomial_coarsening_sequence) + levels.emplace_back(0, i); + + // h + for (unsigned int i = 1; i < this->coarse_grid_triangulations.size(); + ++i) + levels.emplace_back(i, polynomial_coarsening_sequence.back()); + } + else if (mg_coarsening_type == + Parameters::LinearSolver::MultigridCoarseningSequenceType::ph) + { + // h + for (unsigned int i = 0; + i < this->coarse_grid_triangulations.size() - 1; + ++i) + levels.emplace_back(i, polynomial_coarsening_sequence.front()); + + // p + for (const auto i : polynomial_coarsening_sequence) + levels.emplace_back(this->coarse_grid_triangulations.size() - 1, i); + } + else if (mg_coarsening_type == + Parameters::LinearSolver::MultigridCoarseningSequenceType::p) + { + // p + for (const auto i : polynomial_coarsening_sequence) + levels.emplace_back(this->coarse_grid_triangulations.size() - 1, i); + } + else if (mg_coarsening_type == + Parameters::LinearSolver::MultigridCoarseningSequenceType::h) + { + // h + for (unsigned int i = 0; i < this->coarse_grid_triangulations.size(); + ++i) + levels.emplace_back(i, polynomial_coarsening_sequence.back()); + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + // Define maximum and minimum level according to triangulations - const unsigned int n_h_levels = this->coarse_grid_triangulations.size(); - this->minlevel = 0; - this->maxlevel = n_h_levels - 1; + this->minlevel = 0; + this->maxlevel = levels.size() - 1; // Local object for constraints of the different levels MGLevelObject> @@ -620,7 +682,8 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( for (unsigned int l = this->minlevel; l <= this->maxlevel; ++l) { - this->dof_handlers[l].reinit(*this->coarse_grid_triangulations[l]); + this->dof_handlers[l].reinit( + *this->coarse_grid_triangulations[levels[l].first]); // To use elements with linear interpolation for coarse-grid we need // to create the min level dof handler with the appropriate element @@ -630,6 +693,11 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( .mg_use_fe_q_iso_q1 && l == this->minlevel) { + AssertThrow( + mg_coarsening_type == + Parameters::LinearSolver::MultigridCoarseningSequenceType::h, + ExcNotImplemented()); + const auto points = QGaussLobatto<1>(this->dof_handler.get_fe().degree + 1) .get_points(); @@ -638,7 +706,8 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( FESystem(FE_Q_iso_Q1(points), dim + 1)); } else - this->dof_handlers[l].distribute_dofs(this->dof_handler.get_fe()); + this->dof_handlers[l].distribute_dofs( + FESystem(FE_Q(levels[l].second), dim + 1)); } this->mg_setup_timer.leave_subsection( @@ -654,7 +723,7 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( ++level) this->pcout << " Level " << level << ": " << this->dof_handlers[level].n_dofs() << " DoFs, " - << this->coarse_grid_triangulations[level] + << this->coarse_grid_triangulations[levels[level].first] ->n_global_active_cells() << " cells" << std::endl; this->pcout << std::endl; @@ -802,12 +871,17 @@ MFNavierStokesPreconditionGMG::MFNavierStokesPreconditionGMG( // Provide appropriate quadrature depending on the type of elements of // the level - auto quadrature_mg = *cell_quadrature; + Quadrature quadrature_mg = QGauss(levels[level].second + 1); if (this->simulation_parameters.linear_solver .at(PhysicsID::fluid_dynamics) .mg_use_fe_q_iso_q1 && level == this->minlevel) { + AssertThrow( + mg_coarsening_type == + Parameters::LinearSolver::MultigridCoarseningSequenceType::h, + ExcNotImplemented()); + const auto points = QGaussLobatto<1>(this->dof_handler.get_fe().degree + 1) .get_points();