From 19cf04c0d3a756f6d050029e00a7b3eb5c2bd9da Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 8 Apr 2022 14:05:29 +0200 Subject: [PATCH] Fix MGTransferGlobalCoarsening::interpolate_to_mg() --- include/deal.II/matrix_free/constraint_info.h | 7 +- .../mg_transfer_global_coarsening.templates.h | 30 ++- .../interpolate_01.cc | 201 ++++++++++++++++++ ...rpolate_01.with_p4est=true.mpirun=1.output | 121 +++++++++++ ...rpolate_01.with_p4est=true.mpirun=4.output | 121 +++++++++++ 5 files changed, 471 insertions(+), 9 deletions(-) create mode 100644 tests/multigrid-global-coarsening/interpolate_01.cc create mode 100644 tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output create mode 100644 tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h index dd5c97a483..fdf56f4429 100644 --- a/include/deal.II/matrix_free/constraint_info.h +++ b/include/deal.II/matrix_free/constraint_info.h @@ -125,8 +125,11 @@ namespace internal this->plain_dof_indices_per_cell.resize(n_cells); this->constraint_indicator_per_cell.resize(n_cells); - if (use_fast_hanging_node_algorithm && - dof_handler.get_triangulation().has_hanging_nodes()) + // note: has_hanging_nodes() is a global operatrion + const bool has_hanging_nodes = + dof_handler.get_triangulation().has_hanging_nodes(); + + if (use_fast_hanging_node_algorithm && has_hanging_nodes) { hanging_nodes = std::make_unique>( dof_handler.get_triangulation()); diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index 62a10e09f2..3a14658609 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -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_unsigned_int); std::vector local_dof_indices( transfer.schemes[0].n_dofs_per_cell_coarse); @@ -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) { @@ -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_unsigned_int || + 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) @@ -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]); } } @@ -1694,7 +1708,7 @@ namespace internal transfer.schemes[1].restriction_matrix [i * n_dofs_per_cell * GeometryInfo::max_children_per_cell + - j + c * n_dofs_per_cell] = matrix(i, j); + j + c * n_dofs_per_cell] += matrix(i, j); } } } @@ -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 = diff --git a/tests/multigrid-global-coarsening/interpolate_01.cc b/tests/multigrid-global-coarsening/interpolate_01.cc new file mode 100644 index 0000000000..9b9de8ec0f --- /dev/null +++ b/tests/multigrid-global-coarsening/interpolate_01.cc @@ -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 +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include + +#include "../tests.h" + + +template +void +create_quadrant(Triangulation &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 +std::shared_ptr +create_partitioner(const DoFHandler &dof_handler) +{ + return std::make_shared( + dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler), + dof_handler.get_communicator()); +} + + + +template +class RightHandSideFunction : public Function +{ +public: + RightHandSideFunction() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const + { + (void)component; + + return p[0]; + } +}; + + + +template +void +test(const unsigned int n_refinements, + const unsigned int fe_degree_fine, + const bool do_dg) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + const unsigned int min_level = 0; + const unsigned int max_level = n_refinements; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + create_quadrant(tria, n_refinements); + + const auto trias = + MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(tria); + + MGLevelObject> dof_handlers(min_level, max_level); + MGLevelObject> 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{fe_degree_fine}); + else + dof_handler.distribute_dofs(FE_Q{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 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 fu; + + VectorTools::interpolate(dof_handlers[max_level], fu, vec); + + MGLevelObject 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 norm_per_cell(trias[l]->n_active_cells()); + VectorTools::integrate_difference(dof_handlers[l], + results[l], + fu, + norm_per_cell, + QGauss(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); + } +} diff --git a/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..ec37f09d9c --- /dev/null +++ b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output @@ -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:: diff --git a/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..ec37f09d9c --- /dev/null +++ b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output @@ -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:: -- 2.39.5