From ad09d6fdcc231255ea4bea6c2af240f278ec2fd1 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 8 Mar 2021 22:13:23 +0100 Subject: [PATCH] MGLevelGlobalTransfer/MGTransferMatrixFree: allow user to init vectors --- doc/news/changes/minor/20210309Munch | 7 +++ include/deal.II/multigrid/mg_transfer.h | 7 +++ .../deal.II/multigrid/mg_transfer.templates.h | 32 ++++++++----- .../multigrid/mg_transfer_matrix_free.h | 10 ++++ source/multigrid/mg_transfer_matrix_free.cc | 48 +++++++++++++++++++ tests/matrix_free/multigrid_dg_sip_01.cc | 10 ++-- 6 files changed, 99 insertions(+), 15 deletions(-) create mode 100644 doc/news/changes/minor/20210309Munch diff --git a/doc/news/changes/minor/20210309Munch b/doc/news/changes/minor/20210309Munch new file mode 100644 index 0000000000..e6b3b6019a --- /dev/null +++ b/doc/news/changes/minor/20210309Munch @@ -0,0 +1,7 @@ +New: The class MGTransferMatrixFree::build() now also +accepts an optional function for initializing the internal level vectors. +This is useful if one uses the the transfer operators in the context of +smoothers that are built around MatrixFree objects. +
+(Peter Munch, 2021/03/09) + diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h index b20acc46aa..021f53e262 100644 --- a/include/deal.II/multigrid/mg_transfer.h +++ b/include/deal.II/multigrid/mg_transfer.h @@ -611,6 +611,13 @@ protected: mutable MGLevelObject> solution_ghosted_level_vector; + /** + * Function to initialize internal level vectors. + */ + std::function &)> + initialize_dof_vector; + private: /** * This function is called to make sure that build() has been invoked. diff --git a/include/deal.II/multigrid/mg_transfer.templates.h b/include/deal.II/multigrid/mg_transfer.templates.h index 2b8c186bf4..4941974466 100644 --- a/include/deal.II/multigrid/mg_transfer.templates.h +++ b/include/deal.II/multigrid/mg_transfer.templates.h @@ -438,14 +438,15 @@ MGLevelGlobalTransfer>::copy_to_mg( // In case a ghosted level vector has been initialized, we can simply // use that as a template for the vector partitioning. If not, we // resort to the locally owned range of the dof handler. - if (level <= ghosted_level_vector.max_level() && - ghosted_level_vector[level].size() == dof_handler.n_dofs(level)) + if (initialize_dof_vector) + initialize_dof_vector(level, dst[level]); + else if (level <= ghosted_level_vector.max_level() && + ghosted_level_vector[level].size() == + dof_handler.n_dofs(level)) dst[level].reinit(ghosted_level_vector[level], false); else - { - dst[level].reinit(dof_handler.locally_owned_mg_dofs(level), - dof_handler.get_communicator()); - } + dst[level].reinit(dof_handler.locally_owned_mg_dofs(level), + dof_handler.get_communicator()); } else if ((perform_plain_copy == false && perform_renumbered_plain_copy == false) || @@ -561,20 +562,29 @@ MGLevelGlobalTransfer>::copy_from_mg( { // the ghosted vector should already have the correct local size (but // different parallel layout) - AssertDimension(ghosted_level_vector[level].locally_owned_size(), - src[level].locally_owned_size()); + if (ghosted_level_vector[level].size() > 0) + AssertDimension(ghosted_level_vector[level].locally_owned_size(), + src[level].locally_owned_size()); // the first time around, we copy the source vector to the temporary // vector that we hold for the purpose of data exchange LinearAlgebra::distributed::Vector &ghosted_vector = ghosted_level_vector[level]; - ghosted_vector = src[level]; - ghosted_vector.update_ghost_values(); + + if (ghosted_level_vector[level].size() > 0) + ghosted_vector = src[level]; + + const auto ghosted_vector_ptr = (ghosted_level_vector[level].size() > 0) ? + &ghosted_vector : + &src[level]; + + ghosted_vector_ptr->update_ghost_values(); + auto copy_unknowns = [&](const Table<2, unsigned int> &indices) { for (unsigned int i = 0; i < indices.n_cols(); ++i) dst.local_element(indices(0, i)) = - ghosted_vector.local_element(indices(1, i)); + ghosted_vector_ptr->local_element(indices(1, i)); }; // first copy local unknowns diff --git a/include/deal.II/multigrid/mg_transfer_matrix_free.h b/include/deal.II/multigrid/mg_transfer_matrix_free.h index f87fff7d48..3e558e0954 100644 --- a/include/deal.II/multigrid/mg_transfer_matrix_free.h +++ b/include/deal.II/multigrid/mg_transfer_matrix_free.h @@ -105,6 +105,16 @@ public: &external_partitioners = std::vector>()); + /** + * Same as above but taking a lambda for initializing vector instead of + * partitioners. + */ + void + build(const DoFHandler &dof_handler, + const std::function &)> + &initialize_dof_vector); + /** * Prolongate a vector from level to_level-1 to level * to_level using the embedding matrices of the underlying finite diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index 8010b0fa6f..6cb3169617 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -93,6 +93,7 @@ MGTransferMatrixFree::clear() } + template void MGTransferMatrixFree::build( @@ -106,6 +107,20 @@ MGTransferMatrixFree::build( "distribute_mg_dofs() function called, but this is a prerequisite " "for multigrid transfers. You will need to call this function, " "probably close to where you already call distribute_dofs().")); + if (external_partitioners.size() > 0) + { + Assert( + this->initialize_dof_vector == nullptr, + ExcMessage( + "A initialize_dof_vector function has already been registered in the constructor!")); + + this->initialize_dof_vector = + [external_partitioners]( + const unsigned int level, + LinearAlgebra::distributed::Vector &vec) { + vec.reinit(external_partitioners[level]); + }; + } this->fill_and_communicate_copy_indices(dof_handler); @@ -186,6 +201,39 @@ MGTransferMatrixFree::build( +template +void +MGTransferMatrixFree::build( + const DoFHandler &dof_handler, + const std::function &)> + &initialize_dof_vector) +{ + if (initialize_dof_vector) + { + const unsigned int n_levels = + dof_handler.get_triangulation().n_global_levels(); + + std::vector> + external_partitioners(n_levels); + + for (unsigned int level = 0; level < n_levels; ++level) + { + LinearAlgebra::distributed::Vector vector; + initialize_dof_vector(level, vector); + external_partitioners[level] = vector.get_partitioner(); + } + + build(dof_handler, external_partitioners); + } + else + { + build(dof_handler); + } +} + + + template void MGTransferMatrixFree::prolongate( diff --git a/tests/matrix_free/multigrid_dg_sip_01.cc b/tests/matrix_free/multigrid_dg_sip_01.cc index 381b615514..63a436459a 100644 --- a/tests/matrix_free/multigrid_dg_sip_01.cc +++ b/tests/matrix_free/multigrid_dg_sip_01.cc @@ -596,9 +596,11 @@ do_test(const DoFHandler &dof, const unsigned n_q_points_1d) mg_constrained_dofs.initialize(dof); mg_constrained_dofs.make_zero_boundary_constraints(dof, {0}); - MGTransferMF mg_transfer(mg_matrices, - mg_constrained_dofs); - mg_transfer.build(dof); + MGTransferMatrixFree mg_transfer(mg_constrained_dofs); + + mg_transfer.build(dof, [&](const unsigned int level, auto &vec) { + mg_matrices[level].initialize_dof_vector(vec); + }); mg::Matrix> mg_matrix(mg_matrices); @@ -606,7 +608,7 @@ do_test(const DoFHandler &dof, const unsigned n_q_points_1d) mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); PreconditionMG, - MGTransferMF> + MGTransferMatrixFree> preconditioner(dof, mg, mg_transfer); { -- 2.39.5