Skip to content

Commit

Permalink
Merge pull request #15812 from luca-heltai/fix_hanging_nodes_in_paral…
Browse files Browse the repository at this point in the history
…lel_nnmg

Add constrained entries to partitioner_coarse in NonNestedMG
  • Loading branch information
peterrum committed Aug 1, 2023
2 parents dd4023d + 0386188 commit 7c0fd73
Show file tree
Hide file tree
Showing 4 changed files with 268 additions and 49 deletions.
95 changes: 46 additions & 49 deletions include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
Original file line number Diff line number Diff line change
Expand Up @@ -1043,49 +1043,6 @@ namespace internal

class MGTwoLevelTransferImplementation
{
template <int dim, typename Number>
static std::shared_ptr<const Utilities::MPI::Partitioner>
create_coarse_partitioner(
const DoFHandler<dim> & dof_handler_coarse,
const dealii::AffineConstraints<Number> &constraints_coarse,
const unsigned int mg_level_coarse)
{
IndexSet locally_relevant_dofs =
(mg_level_coarse == numbers::invalid_unsigned_int) ?
DoFTools::extract_locally_active_dofs(dof_handler_coarse) :
DoFTools::extract_locally_active_level_dofs(dof_handler_coarse,
mg_level_coarse);

std::vector<types::global_dof_index> locally_relevant_dofs_temp;

for (const auto i : locally_relevant_dofs)
{
if (locally_relevant_dofs.is_element(i) == false)
locally_relevant_dofs_temp.emplace_back(i);

const auto constraints = constraints_coarse.get_constraint_entries(i);

if (constraints)
for (const auto &p : *constraints)
if (locally_relevant_dofs.is_element(p.first) == false)
locally_relevant_dofs_temp.emplace_back(p.first);
}

std::sort(locally_relevant_dofs_temp.begin(),
locally_relevant_dofs_temp.end());
locally_relevant_dofs.add_indices(locally_relevant_dofs_temp.begin(),
locally_relevant_dofs_temp.end());

return std::make_shared<Utilities::MPI::Partitioner>(
mg_level_coarse == numbers::invalid_unsigned_int ?
dof_handler_coarse.locally_owned_dofs() :
dof_handler_coarse.locally_owned_mg_dofs(mg_level_coarse),
locally_relevant_dofs,
dof_handler_coarse.get_communicator());
}



template <int dim, typename Number>
static void
compute_weights(
Expand Down Expand Up @@ -1212,6 +1169,49 @@ namespace internal


public:
template <int dim, typename Number>
static std::shared_ptr<const Utilities::MPI::Partitioner>
create_coarse_partitioner(
const DoFHandler<dim> & dof_handler_coarse,
const dealii::AffineConstraints<Number> &constraints_coarse,
const unsigned int mg_level_coarse)
{
IndexSet locally_relevant_dofs =
(mg_level_coarse == numbers::invalid_unsigned_int) ?
DoFTools::extract_locally_active_dofs(dof_handler_coarse) :
DoFTools::extract_locally_active_level_dofs(dof_handler_coarse,
mg_level_coarse);

std::vector<types::global_dof_index> locally_relevant_dofs_temp;

for (const auto i : locally_relevant_dofs)
{
if (locally_relevant_dofs.is_element(i) == false)
locally_relevant_dofs_temp.emplace_back(i);

const auto constraints = constraints_coarse.get_constraint_entries(i);

if (constraints)
for (const auto &p : *constraints)
if (locally_relevant_dofs.is_element(p.first) == false)
locally_relevant_dofs_temp.emplace_back(p.first);
}

std::sort(locally_relevant_dofs_temp.begin(),
locally_relevant_dofs_temp.end());
locally_relevant_dofs.add_indices(locally_relevant_dofs_temp.begin(),
locally_relevant_dofs_temp.end());

return std::make_shared<Utilities::MPI::Partitioner>(
mg_level_coarse == numbers::invalid_unsigned_int ?
dof_handler_coarse.locally_owned_dofs() :
dof_handler_coarse.locally_owned_mg_dofs(mg_level_coarse),
locally_relevant_dofs,
dof_handler_coarse.get_communicator());
}



template <int dim, typename Number>
static void
reinit_geometric_transfer(
Expand Down Expand Up @@ -3839,12 +3839,9 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::

// create partitioners and internal vectors
{
IndexSet locally_active_dofs =
DoFTools::extract_locally_active_dofs(dof_handler_coarse);
this->partitioner_coarse.reset(
new Utilities::MPI::Partitioner(dof_handler_coarse.locally_owned_dofs(),
locally_active_dofs,
dof_handler_coarse.get_communicator()));
this->partitioner_coarse =
internal::MGTwoLevelTransferImplementation::create_coarse_partitioner(
dof_handler_coarse, constraint_coarse, numbers::invalid_unsigned_int);

this->vec_coarse.reinit(this->partitioner_coarse);
}
Expand Down
209 changes: 209 additions & 0 deletions tests/multigrid-global-coarsening/non_nested_multigrid_04.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,209 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


/**
* Check that everything works also for locally-refined grids.
*
*/

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/quadrature_lib.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/mapping_q.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

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

#include "multigrid_util.h"


template <int dim>
void
create_quadrant(Triangulation<dim> &tria, const unsigned int n_refinements)
{
// Taken from @cite munch2022gc, according to the description in A FLEXIBLE,
// PARALLEL, ADAPTIVE GEOMETRIC MULTIGRID METHOD FOR FEM (Clevenger, Heister,
// Kanschat, Kronbichler): https://arxiv.org/pdf/1904.03317.pdf

GridGenerator::hyper_cube(tria, -1.0, +1.0);

if (n_refinements == 0)
return;

tria.refine_global(1);

for (unsigned int i = 1; i < n_refinements; i++)
{
for (auto cell : tria.active_cell_iterators())
if (cell->is_locally_owned())
{
bool flag = true;
for (int d = 0; d < dim; d++)
if (cell->center()[d] > 0.0)
flag = false;
if (flag)
cell->set_refine_flag();
}
tria.execute_coarsening_and_refinement();
}

AssertDimension(tria.n_global_levels() - 1, n_refinements);
}


template <int dim, typename Number = double>
void
test(const unsigned int n_refinements,
const unsigned int fe_degree_fine,
const bool do_simplex_mesh)
{
using VectorType = LinearAlgebra::distributed::Vector<Number>;

const unsigned int min_level = 0;
const unsigned int max_level = n_refinements;

MGLevelObject<std::shared_ptr<parallel::distributed::Triangulation<dim>>>
triangulations(min_level, max_level);
MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level);
MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<std::unique_ptr<Mapping<dim>>> mappings(min_level, max_level);
MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
transfers(min_level, max_level);
MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);


// set up levels
for (auto l = min_level; l <= max_level; ++l)
{
triangulations[l] =
std::make_shared<parallel::distributed::Triangulation<dim>>(
MPI_COMM_WORLD);
auto &dof_handler = dof_handlers[l];
auto &constraint = constraints[l];
auto &op = operators[l];

std::unique_ptr<FiniteElement<dim>> fe;
std::unique_ptr<Quadrature<dim>> quad;
std::unique_ptr<Mapping<dim>> mapping;

if (do_simplex_mesh)
{
fe = std::make_unique<FE_SimplexP<dim>>(fe_degree_fine);
quad = std::make_unique<QGaussSimplex<dim>>(fe_degree_fine + 1);
mapping = std::make_unique<MappingFE<dim>>(FE_SimplexP<dim>(1));
}
else
{
fe = std::make_unique<FE_Q<dim>>(fe_degree_fine);
quad = std::make_unique<QGauss<dim>>(fe_degree_fine + 1);
mapping = std::make_unique<MappingQ<dim>>(1);
}

mappings[l] = mapping->clone();

// set up triangulation

create_quadrant(*triangulations[l], 2 * l + 1);

// set up dofhandler
dof_handler.reinit(*triangulations[l]);
dof_handler.distribute_dofs(*fe);

// set up constraints
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler,
locally_relevant_dofs);
constraint.reinit(locally_relevant_dofs);
VectorTools::interpolate_boundary_values(
*mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraint);
DoFTools::make_hanging_node_constraints(dof_handler, constraint);
constraint.close();

// set up operator
op.reinit(*mapping, dof_handler, *quad, constraint);
}

// set up transfer operator
for (unsigned int l = min_level; l < max_level; ++l)
{
transfers[l + 1] =
std::make_shared<MGTwoLevelTransferNonNested<dim, VectorType>>();
transfers[l + 1]->reinit(dof_handlers[l + 1],
dof_handlers[l],
*mappings[l + 1],
*mappings[l],
constraints[l + 1],
constraints[l]);
}

MGTransferGlobalCoarsening<dim, VectorType> transfer(
transfers,
[&](const auto l, auto &vec) { operators[l].initialize_dof_vector(vec); });


GMGParameters mg_data; // TODO

VectorType dst, src;
operators[max_level].initialize_dof_vector(dst);
operators[max_level].initialize_dof_vector(src);

operators[max_level].rhs(src);

ReductionControl solver_control(
mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false);

mg_solve(solver_control,
dst,
src,
mg_data,
dof_handlers[max_level],
operators[max_level],
operators,
transfer);

deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
<< (do_simplex_mesh ? "tri " : "quad") << ' '
<< solver_control.last_step() << std::endl;
}

int
main(int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
MPILogInitAll all;

deallog.precision(8);

const unsigned int n_refinements = 3;
const unsigned int degree = 4;
test<2>(n_refinements, degree, false /*quadrilateral*/);
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@

DEAL:0::2 4 3 quad 3
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@

DEAL:0::2 4 3 quad 3

DEAL:1::2 4 3 quad 3


DEAL:2::2 4 3 quad 3


DEAL:3::2 4 3 quad 3

0 comments on commit 7c0fd73

Please sign in to comment.