From 0a55d007400eebf8c17f96f0ca20a36f5ea8cc2d Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 15 Jul 2023 14:34:58 +0200 Subject: [PATCH] MGTransferGlobalCoarsening: initialize with MGConstrainedDoFs --- .../deal.II/multigrid/mg_constrained_dofs.h | 72 ++++++ .../multigrid/mg_transfer_global_coarsening.h | 232 +++++++++++++++++- .../multigrid_a_03.cc | 153 ++++++++++++ ...id_a_03.mpirun=1.with_trilinos=true.output | 10 + 4 files changed, 464 insertions(+), 3 deletions(-) create mode 100644 tests/multigrid-global-coarsening/multigrid_a_03.cc create mode 100644 tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index 9d0f2511d7..82f2d2711b 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -225,6 +225,27 @@ public: const AffineConstraints & get_user_constraint_matrix(const unsigned int level) const; + /** + * Merge selected constraints of a specifiedlevel into a given single + * AffineConstraints object. + * + * @param constraints AffineConstraints object to be filled. + * @param level Refinement to be considered. + * @param add_boundary_indices Add boundary indices. + * @param add_refinement_edge_indices Add refinement-edge indices. + * @param add_level_constraints Add level constraints including the one passed + * during initialize() and periodicy constraints. + * @param add_user_constraints Add user constraints. + */ + template + void + merge_constraints(AffineConstraints &constraints, + const unsigned int level, + const bool add_boundary_indices, + const bool add_refinement_edge_indices, + const bool add_level_constraints, + const bool add_user_constraints) const; + private: /** * The indices of boundary dofs for each level. @@ -565,6 +586,57 @@ MGConstrainedDoFs::get_user_constraint_matrix(const unsigned int level) const } +template +inline void +MGConstrainedDoFs::merge_constraints(AffineConstraints &constraints, + const unsigned int level, + const bool add_boundary_indices, + const bool add_refinement_edge_indices, + const bool add_level_constraints, + const bool add_user_constraints) const +{ + constraints.clear(); + + // determine local lines + IndexSet index_set(this->get_refinement_edge_indices(level).size()); + + if (add_boundary_indices && this->have_boundary_indices()) + index_set.add_indices(this->get_boundary_indices(level)); + + if (add_refinement_edge_indices) + index_set.add_indices(this->get_refinement_edge_indices(level)); + + if (add_level_constraints) + index_set.add_indices(this->get_level_constraints(level).get_local_lines()); + + if (add_user_constraints) + index_set.add_indices( + this->get_user_constraint_matrix(level).get_local_lines()); + + constraints.reinit(index_set); + + // merge constraints + if (add_boundary_indices && this->have_boundary_indices()) + constraints.add_lines(this->get_boundary_indices(level)); + + if (add_refinement_edge_indices) + constraints.add_lines(this->get_refinement_edge_indices(level)); + + if (add_level_constraints) + constraints.merge(this->get_level_constraints(level), + AffineConstraints::left_object_wins, + true); + + if (add_user_constraints) + constraints.merge(this->get_user_constraint_matrix(level), + AffineConstraints::left_object_wins, + true); + + // finalize setup + constraints.close(); +} + + DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index 0ef32204b4..3a72e69c9c 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -806,7 +806,23 @@ private: /** * Implementation of the MGTransferBase. In contrast to * other multigrid transfer operators, the user can provide separate - * transfer operators of type MGTwoLevelTransfer between each level. + * transfer operators of type MGTwoLevelTransfer between each level. The + * sequence of functions calls for setup is: + * @code + * MGTransferGlobalCoarsening mg_transfer; + * mg_transfer.intitialize_two_level_transfers(two_level_transfers); + * mg_transfer.build(partitioners); + * @endcode + * + * Alternatively, this class can also be set up as in the case of + * MGTransferMatrixFree: + * @code + * MGTransferGlobalCoarsening mg_transfer; + * mg_transfer.initialize_constraints(mg_constrained_dofs); + * mg_transfer.build(dof_handler, partitioners); + * @endcode + * However, this is way to set up is currently only working for globally + * refined meshes. * * This class currently only works for the tensor-product finite elements * FE_Q and FE_DGQ and simplex elements FE_SimplexP and FE_SimplexDGP as well as @@ -824,6 +840,16 @@ public: */ using Number = typename VectorType::value_type; + /** + * Default constructor. + */ + MGTransferGlobalCoarsening() = default; + + /** + * @name Global coarsening. + */ + /** @{ */ + /** * Constructor taking a collection of transfer operators (with the coarsest * level kept empty in @p transfer) and an optional function that initializes the @@ -839,6 +865,14 @@ public: const std::function &initialize_dof_vector = {}); + /** + * Set two-level transfers. + */ + template + void + intitialize_two_level_transfers( + const MGLevelObject &transfer); + /** * Similar function to MGTransferMatrixFree::build() with the difference that * the information for the prolongation for each level has already been built. @@ -859,6 +893,57 @@ public: build(const std::function &initialize_dof_vector); + /** @} */ + + /** + * @name Local smoothing. + */ + /** @{ */ + + /** + * Constructor with constraints. Equivalent to the default constructor + * followed by initialize_constraints(). + * + * @note See also MGTransferMatrixFree. + */ + MGTransferGlobalCoarsening(const MGConstrainedDoFs &mg_constrained_dofs); + + /** + * Initialize the constraints to be used in build(). + * + * @note See also MGTransferMatrixFree. + */ + void + initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs); + + /** + * Actually build the information for the prolongation for each level. + * + * @note See also MGTransferMatrixFree. + */ + void + build(const DoFHandler &dof_handler, + const std::vector> + &external_partitioners = {}); + + /** + * Same as above but taking a lambda for initializing vector instead of + * partitioners. + * + * @note See also MGTransferMatrixFree. + */ + void + build(const DoFHandler &dof_handler, + const std::function + &initialize_dof_vector); + + /** @} */ + + /** + * @name Tranfer functions. + */ + /** @{ */ + /** * Perform prolongation. */ @@ -935,6 +1020,13 @@ public: MGLevelObject &dst, const InVector & src) const; + /** @} */ + + /** + * @name Utility functions. + */ + /** @{ */ + /** * Return the memory consumption of the allocated memory in this class. * @@ -956,13 +1048,47 @@ public: unsigned int max_level() const; + /** @} */ + private: + /** + * Initial internal transfer operator. + * + * @note See also MGTransferMatrixFree. + */ + void + intitialize_internal_transfer( + const DoFHandler & dof_handler, + const SmartPointer &mg_constrained_dofs); + + /** + * Set references to two-level transfer operators to be used. + */ + template + void + intitialize_transfer_references( + const MGLevelObject &transfer); + /** * Function to initialize internal level vectors. */ void initialize_dof_vector(const unsigned int level, VectorType &vector) const; + /** + * MGConstrainedDoFs passed during build(). + * + * @note See also MGTransferMatrixFree. + */ + SmartPointer mg_constrained_dofs; + + /** + * Internal transfer operator. + * + * @note See also MGTransferMatrixFree. + */ + MGLevelObject> internal_transfer; + /** * Collection of the two-level transfer operators. */ @@ -1021,6 +1147,80 @@ MGTransferGlobalCoarsening::MGTransferGlobalCoarsening( const MGLevelObject &transfer, const std::function &initialize_dof_vector) +{ + this->intitialize_transfer_references(transfer); + this->build(initialize_dof_vector); +} + + + +template +MGTransferGlobalCoarsening::MGTransferGlobalCoarsening( + const MGConstrainedDoFs &mg_constrained_dofs) +{ + this->initialize_constraints(mg_constrained_dofs); +} + + + +template +template +void +MGTransferGlobalCoarsening::intitialize_two_level_transfers( + const MGLevelObject &transfer) +{ + this->intitialize_transfer_references(transfer); +} + + + +template +void +MGTransferGlobalCoarsening::initialize_constraints( + const MGConstrainedDoFs &mg_constrained_dofs) +{ + this->mg_constrained_dofs = &mg_constrained_dofs; +} + + + +template +void +MGTransferGlobalCoarsening::intitialize_internal_transfer( + const DoFHandler & dof_handler, + const SmartPointer &mg_constrained_dofs) +{ + const unsigned int min_level = 0; + const unsigned int max_level = + dof_handler.get_triangulation().n_global_levels() - 1; + + MGLevelObject> constraints( + min_level, max_level); + + if (mg_constrained_dofs) + for (unsigned int l = min_level; l <= max_level; ++l) + mg_constrained_dofs->merge_constraints( + constraints[l], + l, + /*add_boundary_indices*/ true, + /*add_refinement_edge_indices*/ false, + /*add_level_constraints*/ true, + /*add_user_constraints*/ true); + + this->internal_transfer.resize(min_level, max_level); + + for (unsigned int l = min_level; l < max_level; ++l) + internal_transfer[l + 1].reinit_geometric_transfer( + dof_handler, dof_handler, constraints[l + 1], constraints[l], l + 1, l); +} + + + +template +template +void +MGTransferGlobalCoarsening::intitialize_transfer_references( + const MGLevelObject &transfer) { const unsigned int min_level = transfer.min_level(); const unsigned int max_level = transfer.max_level(); @@ -1031,8 +1231,6 @@ MGTransferGlobalCoarsening::MGTransferGlobalCoarsening( this->transfer[l] = &const_cast &>( static_cast &>( Utilities::get_underlying_value(transfer[l]))); - - this->build(initialize_dof_vector); } @@ -1094,6 +1292,34 @@ MGTransferGlobalCoarsening::build( +template +void +MGTransferGlobalCoarsening::build( + const DoFHandler &dof_handler, + const std::vector> + &external_partitioners) +{ + this->intitialize_internal_transfer(dof_handler, mg_constrained_dofs); + this->intitialize_transfer_references(internal_transfer); + this->build(external_partitioners); +} + + + +template +void +MGTransferGlobalCoarsening::build( + const DoFHandler &dof_handler, + const std::function + &initialize_dof_vector) +{ + this->intitialize_internal_transfer(dof_handler, mg_constrained_dofs); + this->intitialize_transfer_references(internal_transfer); + this->build(initialize_dof_vector); +} + + + template void MGTransferGlobalCoarsening::initialize_dof_vector( diff --git a/tests/multigrid-global-coarsening/multigrid_a_03.cc b/tests/multigrid-global-coarsening/multigrid_a_03.cc new file mode 100644 index 0000000000..0f4be4b029 --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_a_03.cc @@ -0,0 +1,153 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2022 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 global-coarsening multigrid so that it works also on refinement levels. + */ + +#include "multigrid_util.h" + +template +void +test(const unsigned int n_refinements, 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 unsigned int min_level = 0; + const unsigned int max_level = n_refinements; + + MGLevelObject> constraints(min_level, max_level); + MGLevelObject> operators(min_level, max_level); + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + fe = std::make_unique>(fe_degree_fine); + quad = std::make_unique>(fe_degree_fine + 1); + mapping = std::make_unique>(FE_Q(1)); + + // set up dofhandler + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(*fe); + dof_handler.distribute_mg_dofs(); + + const std::set dirichlet_boundary = {0}; + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, + dirichlet_boundary); + + // set up levels + for (auto l = min_level; l <= max_level; ++l) + { + auto &constraint = constraints[l]; + auto &op = operators[l]; + + // set up constraints + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + l, + relevant_dofs); + constraint.reinit(relevant_dofs); + constraint.add_lines(mg_constrained_dofs.get_boundary_indices(l)); + constraint.close(); + + constraint.print(std::cout); + + // set up operator + op.reinit(*mapping, dof_handler, *quad, constraint, l); + } + + // set up transfer operator + MGTransferGlobalCoarsening transfer(mg_constrained_dofs); + + transfer.build(dof_handler, [&](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_handler, + operators[max_level], + operators, + transfer); + + constraints[max_level].distribute(dst); + + deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' ' + << solver_control.last_step() << std::endl; + + return; + + static unsigned int counter = 0; + + MGLevelObject results(min_level, max_level); + + transfer.interpolate_to_mg(dof_handler, results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + 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) + test<2>(n_refinements, degree); +} diff --git a/tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output new file mode 100644 index 0000000000..e460d32df9 --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output @@ -0,0 +1,10 @@ + +DEAL:0::2 2 2 3 +DEAL:0::2 3 2 3 +DEAL:0::2 4 2 3 +DEAL:0::2 2 3 3 +DEAL:0::2 3 3 3 +DEAL:0::2 4 3 3 +DEAL:0::2 2 4 3 +DEAL:0::2 3 4 3 +DEAL:0::2 4 4 3 -- 2.39.5