Skip to content

Commit

Permalink
Fix MGTransferGlobalCoarsening::interpolate_to_mg()
Browse files Browse the repository at this point in the history
  • Loading branch information
@peterrum
peterrum committed Apr 8, 2022
1 parent 8197bc4 commit 2dd1181
Show file tree
Hide file tree
Showing 4 changed files with 466 additions and 7 deletions.
30 changes: 23 additions & 7 deletions include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1481,7 +1481,8 @@ namespace internal

// indices
{
transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back());
transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back(),
numbers::invalid_dof_index);

std::vector<types::global_dof_index> local_dof_indices(
transfer.schemes[0].n_dofs_per_cell_coarse);
Expand Down Expand Up @@ -1525,7 +1526,9 @@ namespace internal
dof_handler_coarse,
transfer.schemes[0].n_coarse_cells +
transfer.schemes[1].n_coarse_cells,
use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse));
constraints_coarse.n_constraints() > 0 &&
use_fast_hanging_node_algorithm(dof_handler_coarse,
mg_level_coarse));

process_cells(
[&](const auto &cell_coarse, const auto &cell_fine) {
Expand Down Expand Up @@ -1574,9 +1577,20 @@ namespace internal
for (unsigned int i = 0;
i < transfer.schemes[1].n_dofs_per_cell_coarse;
++i)
level_dof_indices_fine_1[cell_local_children_indices[c][i]] =
transfer.partitioner_fine->global_to_local(
{
const auto index = transfer.partitioner_fine->global_to_local(
local_dof_indices[lexicographic_numbering_fine[i]]);

Assert(level_dof_indices_fine_1
[cell_local_children_indices[c][i]] ==
numbers::invalid_dof_index ||
level_dof_indices_fine_1
[cell_local_children_indices[c][i]] == index,
ExcInternalError());

level_dof_indices_fine_1[cell_local_children_indices[c][i]] =
index;
}
}

// move pointers (only once at the end)
Expand Down Expand Up @@ -1653,7 +1667,7 @@ namespace internal
for (unsigned int j = 0; j < fe->n_dofs_per_cell(); ++j)
transfer.schemes[1]
.restriction_matrix_1d[i * n_child_dofs_1d + j +
c * shift] =
c * shift] +=
matrix(renumbering[i], renumbering[j]);
}
}
Expand Down Expand Up @@ -1694,7 +1708,7 @@ namespace internal
transfer.schemes[1].restriction_matrix
[i * n_dofs_per_cell *
GeometryInfo<dim>::max_children_per_cell +
j + c * n_dofs_per_cell] = matrix(i, j);
j + c * n_dofs_per_cell] += matrix(i, j);
}
}
}
Expand Down Expand Up @@ -2014,7 +2028,9 @@ namespace internal
transfer.constraint_info.reinit(
dof_handler_coarse,
cell_no.back(),
use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse));
constraints_coarse.n_constraints() > 0 &&
use_fast_hanging_node_algorithm(dof_handler_coarse,
mg_level_coarse));

process_cells([&](const auto &cell_coarse, const auto &cell_fine) {
const auto fe_pair_no =
Expand Down
201 changes: 201 additions & 0 deletions tests/multigrid-global-coarsening/interpolate_01.cc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,201 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 2021 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.
//
// ---------------------------------------------------------------------


/**
* Test MGTransferGlobalCoarsening::interpolate_to_mg() for FE_Q and FE_DGQ.
*/

#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_system.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/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 <deal.II/numerics/vector_tools.h>

#include "../tests.h"


template <int dim>
void
create_quadrant(Triangulation<dim> &tria, const unsigned int n_refinements)
{
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>
std::shared_ptr<Utilities::MPI::Partitioner>
create_partitioner(const DoFHandler<dim> &dof_handler)
{
return std::make_shared<Utilities::MPI::Partitioner>(
dof_handler.locally_owned_dofs(),
DoFTools::extract_locally_active_dofs(dof_handler),
dof_handler.get_communicator());
}



template <int dim>
class RightHandSideFunction : public Function<dim>
{
public:
RightHandSideFunction()
: Function<dim>(1)
{}

virtual double
value(const Point<dim> &p, const unsigned int component = 0) const
{
(void)component;

return p[0];
}
};



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

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

parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
create_quadrant(tria, n_refinements);

const auto trias =
MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(tria);

MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level);
MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> transfers(min_level,
max_level);

for (auto l = min_level; l <= max_level; ++l)
{
auto &dof_handler = dof_handlers[l];

dof_handler.reinit(*trias[l]);

if (do_dg)
dof_handler.distribute_dofs(FE_DGQ<dim>{fe_degree_fine});
else
dof_handler.distribute_dofs(FE_Q<dim>{fe_degree_fine});
}

for (unsigned int l = min_level; l < max_level; ++l)
transfers[l + 1].reinit(dof_handlers[l + 1], dof_handlers[l]);

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

VectorType vec;
vec.reinit(create_partitioner(dof_handlers[max_level]));

RightHandSideFunction<dim> fu;

VectorTools::interpolate(dof_handlers[max_level], fu, vec);

MGLevelObject<VectorType> results(min_level, max_level);

transfer.interpolate_to_mg(dof_handlers[max_level], results, vec);

for (unsigned int l = min_level; l <= max_level; ++l)
{
results[l].update_ghost_values();
Vector<float> norm_per_cell(trias[l]->n_active_cells());
VectorTools::integrate_difference(dof_handlers[l],
results[l],
fu,
norm_per_cell,
QGauss<dim>(fe_degree_fine + 2),
VectorTools::L2_norm);
const double error_L2_norm =
VectorTools::compute_global_error(*trias[l],
norm_per_cell,
VectorTools::L2_norm);


deallog << l << " " << results[l].l2_norm() << " " << error_L2_norm
<< std::endl;
}

deallog << std::endl;
}

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

deallog.precision(8);

for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
for (unsigned int degree = 1; degree <= 4; ++degree)
{
test<2>(n_refinements, degree, true);
test<2>(n_refinements, degree, false);
}
}
121 changes: 121 additions & 0 deletions tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@

DEAL::0 2.0000000 8.1041682e-16
DEAL::1 2.8284271 3.6480780e-16
DEAL::2 3.4641016 0.0000000
DEAL::
DEAL::0 2.0000000 0.0000000
DEAL::1 2.4494897 0.0000000
DEAL::2 2.7838822 0.0000000
DEAL::
DEAL::0 2.4494897 4.3469376e-16
DEAL::1 3.8729833 1.9846630e-16
DEAL::2 4.8989795 1.2272899e-16
DEAL::
DEAL::0 2.4494897 1.1453553e-16
DEAL::1 3.5355339 9.9806058e-17
DEAL::2 4.2866070 1.1488487e-16
DEAL::
DEAL::0 3.0983867 6.6306930e-16
DEAL::1 5.0596443 2.7087997e-16
DEAL::2 6.4498062 2.3247867e-16
DEAL::
DEAL::0 3.0983867 3.0349490e-16
DEAL::1 4.7328638 2.0099536e-16
DEAL::2 5.8694122 2.1659599e-16
DEAL::
DEAL::0 3.7796447 1.1027361e-15
DEAL::1 6.2678317 3.9148216e-16
DEAL::2 8.0178373 3.1129099e-16
DEAL::
DEAL::0 3.7796447 2.3002617e-16
DEAL::1 5.9461873 2.9830841e-16
DEAL::2 7.4558223 3.2071973e-16
DEAL::
DEAL::0 2.0000000 1.1994317e-15
DEAL::1 2.8284271 8.0596747e-16
DEAL::2 3.4641016 6.0282265e-16
DEAL::3 6.3245553 0.0000000
DEAL::
DEAL::0 2.0000000 0.0000000
DEAL::1 2.4494897 0.0000000
DEAL::2 2.7838822 0.0000000
DEAL::3 4.2573466 0.0000000
DEAL::
DEAL::0 2.4494897 6.6147021e-16
DEAL::1 3.8729833 4.1436381e-16
DEAL::2 4.8989795 3.4330648e-16
DEAL::3 9.3273791 1.4568899e-16
DEAL::
DEAL::0 2.4494897 1.1453553e-16
DEAL::1 3.5355339 9.9806058e-17
DEAL::2 4.2866070 1.1488487e-16
DEAL::3 7.3271754 1.3584636e-16
DEAL::
DEAL::0 3.0983867 9.8309852e-16
DEAL::1 5.0596443 5.6517728e-16
DEAL::2 6.4498062 3.8488399e-16
DEAL::3 12.393547 2.3245655e-16
DEAL::
DEAL::0 3.0983867 4.3584763e-16
DEAL::1 4.7328638 2.7306050e-16
DEAL::2 5.8694122 2.3526684e-16
DEAL::3 10.419933 2.2576644e-16
DEAL::
DEAL::0 3.7796447 1.4273981e-15
DEAL::1 6.2678317 7.7585762e-16
DEAL::2 8.0178373 7.4056572e-16
DEAL::3 15.468863 3.7969912e-16
DEAL::
DEAL::0 3.7796447 2.8778492e-16
DEAL::1 5.9461873 3.3252973e-16
DEAL::2 7.4558223 3.4411320e-16
DEAL::3 13.509587 3.6867808e-16
DEAL::
DEAL::0 2.0000000 1.4131724e-15
DEAL::1 2.8284271 1.0606982e-15
DEAL::2 3.4641016 8.5294693e-16
DEAL::3 6.3245553 4.3812845e-16
DEAL::4 10.723805 0.0000000
DEAL::
DEAL::0 2.0000000 0.0000000
DEAL::1 2.4494897 0.0000000
DEAL::2 2.7838822 0.0000000
DEAL::3 4.2573466 0.0000000
DEAL::4 6.4128777 0.0000000
DEAL::
DEAL::0 2.4494897 8.1236828e-16
DEAL::1 3.8729833 5.3562441e-16
DEAL::2 4.8989795 4.2374356e-16
DEAL::3 9.3273791 2.4185602e-16
DEAL::4 15.992186 1.4728562e-16
DEAL::
DEAL::0 2.4494897 1.1453553e-16
DEAL::1 3.5355339 9.9806058e-17
DEAL::2 4.2866070 1.1488487e-16
DEAL::3 7.3271754 1.3584636e-16
DEAL::4 11.805719 1.3167380e-16
DEAL::
DEAL::0 3.0983867 1.4470137e-15
DEAL::1 5.0596443 8.9797785e-16
DEAL::2 6.4498062 6.6972640e-16
DEAL::3 12.393547 3.5145888e-16
DEAL::4 21.297887 2.3268951e-16
DEAL::
DEAL::0 3.0983867 6.3814779e-16
DEAL::1 4.7328638 4.5675034e-16
DEAL::2 5.8694122 3.7189323e-16
DEAL::3 10.419933 2.4820397e-16
DEAL::4 17.160274 2.2156718e-16
DEAL::
DEAL::0 3.7796447 1.4543742e-15
DEAL::1 6.2678317 8.5375061e-16
DEAL::2 8.0178373 8.4327874e-16
DEAL::3 15.468863 5.4672827e-16
DEAL::4 26.608940 3.8089812e-16
DEAL::
DEAL::0 3.7796447 3.6858565e-16
DEAL::1 5.9461873 3.5838440e-16
DEAL::2 7.4558223 3.6691834e-16
DEAL::3 13.509587 3.6351703e-16
DEAL::4 22.497222 3.6587400e-16
DEAL::

0 comments on commit 2dd1181

Please sign in to comment.