From 7bc367614755f439e1bed5a3c04a68a650344c66 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 21 Aug 2021 14:16:18 +0200 Subject: [PATCH] Extend MGTwoLevelTransfer::reinit_polynomial_transfer() so that it works on any refinement level without HN --- .../multigrid/mg_transfer_global_coarsening.h | 4 +- .../mg_transfer_global_coarsening.templates.h | 36 +++- .../multigrid_p_03.cc | 179 ++++++++++++++++++ ...id_p_03.mpirun=1.with_trilinos=true.output | 37 ++++ 4 files changed, 244 insertions(+), 12 deletions(-) create mode 100644 tests/multigrid-global-coarsening/multigrid_p_03.cc create mode 100644 tests/multigrid-global-coarsening/multigrid_p_03.mpirun=1.with_trilinos=true.output diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index f2ca2ebcf4..c226cef185 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -205,7 +205,7 @@ public: * Set up polynomial coarsening between the given DoFHandler objects ( * @p dof_handler_fine and @p dof_handler_coarse). Polynomial transfers * can be only performed on active levels (`numbers::invalid_unsigned_int`) - * or on coarse-grid levels. + * or on coarse-grid levels, i.e., levels without hanging nodes. * * @note The function polynomial_transfer_supported() can be used to * check if the given polynomial coarsening strategy is supported. @@ -227,7 +227,7 @@ public: * * @note While geometric transfer can be only performed on active levels * (`numbers::invalid_unsigned_int`), polynomial transfers can also be - * performed on coarse-grid levels. + * performed on coarse-grid levels, i.e., levels without hanging nodes. * * @note The function polynomial_transfer_supported() can be used to * check if the given polynomial coarsening strategy is supported. 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 3ec986903f..1aff175162 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -1763,16 +1763,21 @@ namespace internal &transfer) { Assert( - mg_level_fine == 0 || mg_level_fine == numbers::invalid_unsigned_int, + mg_level_fine == numbers::invalid_unsigned_int || + mg_level_fine <= MGTools::max_level_for_coarse_mesh( + dof_handler_fine.get_triangulation()), ExcMessage( "Polynomial transfer is only allowed on the active level " - "(numbers::invalid_unsigned_int) or the coarse-grid level (0).")); + "(numbers::invalid_unsigned_int) or on refinement levels without " + "hanging nodes.")); Assert( - mg_level_coarse == 0 || - mg_level_coarse == numbers::invalid_unsigned_int, + mg_level_coarse == numbers::invalid_unsigned_int || + mg_level_coarse <= MGTools::max_level_for_coarse_mesh( + dof_handler_coarse.get_triangulation()), ExcMessage( "Polynomial transfer is only allowed on the active level " - "(numbers::invalid_unsigned_int) or the coarse-grid level (0).")); + "(numbers::invalid_unsigned_int) or on refinement levels without " + "hanging nodes.")); const PermutationFineDoFHandlerView view(dof_handler_fine, dof_handler_coarse, @@ -1888,12 +1893,20 @@ namespace internal } { IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs(dof_handler_coarse, - locally_relevant_dofs); + + if (mg_level_coarse == numbers::invalid_unsigned_int) + DoFTools::extract_locally_relevant_dofs(dof_handler_coarse, + locally_relevant_dofs); + else + DoFTools::extract_locally_relevant_level_dofs(dof_handler_coarse, + mg_level_coarse, + locally_relevant_dofs); transfer.partitioner_coarse = std::make_shared( - dof_handler_coarse.locally_owned_dofs(), + 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, comm); transfer.vec_coarse.reinit(transfer.partitioner_coarse); @@ -1993,8 +2006,11 @@ namespace internal transfer.schemes[i].level_dof_indices_fine.data(); } - bool fine_indices_touch_remote_dofs = false; - IndexSet locally_owned_dofs = dof_handler_fine.locally_owned_dofs(); + bool fine_indices_touch_remote_dofs = false; + const IndexSet locally_owned_dofs = + mg_level_fine == numbers::invalid_unsigned_int ? + dof_handler_fine.locally_owned_dofs() : + dof_handler_fine.locally_owned_mg_dofs(mg_level_fine); process_cells([&](const auto &cell_coarse, const auto &cell_fine) { const auto fe_pair_no = diff --git a/tests/multigrid-global-coarsening/multigrid_p_03.cc b/tests/multigrid-global-coarsening/multigrid_p_03.cc new file mode 100644 index 0000000000..3af75fdf7e --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_03.cc @@ -0,0 +1,179 @@ +// --------------------------------------------------------------------- +// +// 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 p-multigrid so that it also works on any refinement level without + * hanging nodes. + */ + +#include "multigrid_util.h" + +template +void +test(const unsigned int n_refinements, + const unsigned int level, + const unsigned int fe_degree_fine) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::MeshSmoothing::none, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + GridGenerator::subdivided_hyper_cube(tria, 2); + tria.refine_global(n_refinements); + + const auto level_degrees = + MGTransferGlobalCoarseningTools::create_polynomial_coarsening_sequence( + fe_degree_fine, + MGTransferGlobalCoarseningTools::PolynomialCoarseningSequenceType:: + bisect); + + const unsigned int min_level = 0; + const unsigned int max_level = level_degrees.size() - 1; + + MGLevelObject> dof_handlers(min_level, max_level, tria); + MGLevelObject> constraints(min_level, max_level); + MGLevelObject> transfers(min_level, + max_level); + MGLevelObject> operators(min_level, max_level); + + std::unique_ptr> mapping_; + + // set up levels + for (auto l = min_level; l <= max_level; ++l) + { + auto &dof_handler = dof_handlers[l]; + auto &constraint = constraints[l]; + auto &op = operators[l]; + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + fe = std::make_unique>(level_degrees[l]); + quad = std::make_unique>(level_degrees[l] + 1); + mapping = std::make_unique>(FE_Q(1)); + + if (l == max_level) + mapping_ = mapping->clone(); + + // set up dofhandler + dof_handler.distribute_dofs(*fe); + dof_handler.distribute_mg_dofs(); + + // set up constraints + std::set dirichlet_boundary; + dirichlet_boundary.insert(0); + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, + dirichlet_boundary); + + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + level, + relevant_dofs); + constraint.reinit(relevant_dofs); + constraint.add_lines(mg_constrained_dofs.get_boundary_indices(level)); + constraint.close(); + + constraint.print(std::cout); + + // set up operator + op.reinit(*mapping, dof_handler, *quad, constraint, level); + } + + // set up transfer operator + for (unsigned int l = min_level; l < max_level; ++l) + transfers[l + 1].reinit(dof_handlers[l + 1], + dof_handlers[l], + constraints[l + 1], + constraints[l], + level, + level); + + MGTransferGlobalCoarsening 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); + + constraints[max_level].distribute(dst); + + deallog << dim << " " << fe_degree_fine << " " << n_refinements << " " + << level << " " << dst.size() << " " << solver_control.last_step() + << std::endl; + + return; + + static unsigned int counter = 0; + + MGLevelObject results(min_level, max_level); + + transfer.interpolate_to_mg(dof_handlers[max_level], results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handlers[l]); + data_out.add_data_vector( + results[l], + "solution", + DataOut_DoFData::DataVectorType::type_dof_data); + data_out.build_patches(*mapping_, 2); + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter) + "." + std::to_string(l) + + ".vtk"); + data_out.write_vtk(output); + } + + counter++; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + deallog.precision(8); + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 4; ++degree) + for (unsigned int level = 0; level <= n_refinements; ++level) + test<2>(n_refinements, level, degree); +} diff --git a/tests/multigrid-global-coarsening/multigrid_p_03.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_p_03.mpirun=1.with_trilinos=true.output new file mode 100644 index 0000000000..170074986e --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_03.mpirun=1.with_trilinos=true.output @@ -0,0 +1,37 @@ + +DEAL:0::2 2 2 0 25 2 +DEAL:0::2 2 2 1 81 2 +DEAL:0::2 2 2 2 289 3 +DEAL:0::2 3 2 0 49 2 +DEAL:0::2 3 2 1 169 3 +DEAL:0::2 3 2 2 625 3 +DEAL:0::2 4 2 0 81 3 +DEAL:0::2 4 2 1 289 3 +DEAL:0::2 4 2 2 1089 3 +DEAL:0::2 2 3 0 25 2 +DEAL:0::2 2 3 1 81 2 +DEAL:0::2 2 3 2 289 3 +DEAL:0::2 2 3 3 1089 3 +DEAL:0::2 3 3 0 49 2 +DEAL:0::2 3 3 1 169 3 +DEAL:0::2 3 3 2 625 3 +DEAL:0::2 3 3 3 2401 3 +DEAL:0::2 4 3 0 81 3 +DEAL:0::2 4 3 1 289 3 +DEAL:0::2 4 3 2 1089 3 +DEAL:0::2 4 3 3 4225 3 +DEAL:0::2 2 4 0 25 2 +DEAL:0::2 2 4 1 81 2 +DEAL:0::2 2 4 2 289 3 +DEAL:0::2 2 4 3 1089 3 +DEAL:0::2 2 4 4 4225 3 +DEAL:0::2 3 4 0 49 2 +DEAL:0::2 3 4 1 169 3 +DEAL:0::2 3 4 2 625 3 +DEAL:0::2 3 4 3 2401 3 +DEAL:0::2 3 4 4 9409 3 +DEAL:0::2 4 4 0 81 3 +DEAL:0::2 4 4 1 289 3 +DEAL:0::2 4 4 2 1089 3 +DEAL:0::2 4 4 3 4225 3 +DEAL:0::2 4 4 4 16641 3 -- 2.39.5