Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add different variants of p-multgrid #1209

Merged
merged 4 commits into from
Jul 31, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions include/core/parameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/parsed_function.h>

#include <deal.II/multigrid/mg_transfer_global_coarsening.h>


using namespace dealii;

Expand Down Expand Up @@ -1187,6 +1189,20 @@ namespace Parameters
/// MG enable hessians in jacobian
bool mg_enable_hessians_jacobian;

/// Type of multigrid
enum class MultigridCoarseningSequenceType
peterrum marked this conversation as resolved.
Show resolved Hide resolved
{
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;

Expand Down
47 changes: 47 additions & 0 deletions source/core/parameters.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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"),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is bisect?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From the deal.II documentation:

Screenshot from 2024-07-30 19-47-19

"mg p coarsening type for gcmg");

prm.declare_entry("mg gmres max iterations",
"2000",
Patterns::Integer(),
Expand Down Expand Up @@ -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");
Expand Down
88 changes: 81 additions & 7 deletions source/solvers/mf_navier_stokes.cc
Original file line number Diff line number Diff line change
Expand Up @@ -535,6 +535,10 @@ MFNavierStokesPreconditionGMG<dim>::MFNavierStokesPreconditionGMG(
.preconditioner ==
Parameters::LinearSolver::PreconditionerType::gcmg)
{
AssertThrow(*cell_quadrature ==
QGauss<dim>(this->dof_handler.get_fe().degree + 1),
ExcNotImplemented());

// Create triangulations
this->mg_setup_timer.enter_subsection("Create level triangulations");
this->coarse_grid_triangulations =
Expand Down Expand Up @@ -599,10 +603,68 @@ MFNavierStokesPreconditionGMG<dim>::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<std::pair<unsigned int, unsigned int>> 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<AffineConstraints<typename VectorType::value_type>>
Expand All @@ -620,7 +682,8 @@ MFNavierStokesPreconditionGMG<dim>::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
Expand All @@ -630,6 +693,11 @@ MFNavierStokesPreconditionGMG<dim>::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();
Expand All @@ -638,7 +706,8 @@ MFNavierStokesPreconditionGMG<dim>::MFNavierStokesPreconditionGMG(
FESystem<dim>(FE_Q_iso_Q1<dim>(points), dim + 1));
}
else
this->dof_handlers[l].distribute_dofs(this->dof_handler.get_fe());
this->dof_handlers[l].distribute_dofs(
FESystem<dim>(FE_Q<dim>(levels[l].second), dim + 1));
}

this->mg_setup_timer.leave_subsection(
Expand All @@ -654,7 +723,7 @@ MFNavierStokesPreconditionGMG<dim>::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;
Expand Down Expand Up @@ -802,12 +871,17 @@ MFNavierStokesPreconditionGMG<dim>::MFNavierStokesPreconditionGMG(

// Provide appropriate quadrature depending on the type of elements of
// the level
auto quadrature_mg = *cell_quadrature;
Quadrature<dim> quadrature_mg = QGauss<dim>(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();
Expand Down
Loading