From: Martin Kronbichler Date: Wed, 25 May 2022 15:21:00 +0000 (+0200) Subject: Use GridTools::exchange_cell_data_on_ghosts for DoF renumbering X-Git-Tag: v9.4.0-rc1~112^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=487c2a67e7968dc66570b00e56225c84371848b2;p=dealii.git Use GridTools::exchange_cell_data_on_ghosts for DoF renumbering --- diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index 96fc91bf24..c27693fd9e 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -2533,16 +2533,12 @@ namespace internal const std::vector &new_numbers, const IndexSet & indices_we_care_about, DoFHandler &dof_handler, - const unsigned int level, - const bool check_validity) + const unsigned int level) { - (void)check_validity; Assert(level < dof_handler.get_triangulation().n_levels(), ExcInternalError()); - for (typename std::vector< - typename DoFHandler::MGVertexDoFs>::iterator i = - dof_handler.mg_vertex_dofs.begin(); + for (auto i = dof_handler.mg_vertex_dofs.begin(); i != dof_handler.mg_vertex_dofs.end(); ++i) // if the present vertex lives on the current level @@ -2559,10 +2555,9 @@ namespace internal if (idx != numbers::invalid_dof_index) { - Assert(check_validity == false || - (indices_we_care_about.size() > 0 ? - indices_we_care_about.is_element(idx) : - (idx < new_numbers.size())), + Assert(indices_we_care_about.size() > 0 ? + indices_we_care_about.is_element(idx) : + (idx < new_numbers.size()), ExcInternalError()); i->set_index(level, d, @@ -2823,8 +2818,7 @@ namespace internal renumber_vertex_mg_dofs(new_numbers, indices_we_care_about, dof_handler, - level, - check_validity); + level); }); tasks += Threads::new_task([&]() { renumber_face_mg_dofs(new_numbers, @@ -3617,15 +3611,12 @@ namespace internal { template void - communicate_mg_ghost_cells( - const typename dealii::parallel:: - DistributedTriangulationBase &tria, - DoFHandler & dof_handler) + communicate_mg_ghost_cells(DoFHandler &dof_handler) { - (void)tria; const auto pack = [&](const auto &cell) { // why would somebody request a cell that is not ours? - Assert(cell->level_subdomain_id() == tria.locally_owned_subdomain(), + Assert(cell->level_subdomain_id() == + dof_handler.get_triangulation().locally_owned_subdomain(), ExcInternalError()); std::vector data( @@ -3686,9 +3677,7 @@ namespace internal template void - communicate_mg_ghost_cells(const typename dealii::parallel:: - distributed::Triangulation<1, spacedim> &, - DoFHandler<1, spacedim> &) + communicate_mg_ghost_cells(DoFHandler<1, spacedim> &) { Assert(false, ExcNotImplemented()); } @@ -4173,20 +4162,18 @@ namespace internal // mark all ghost cells for transfer { for (const auto &cell : dof_handler->cell_iterators()) - if (cell->level_subdomain_id() != - dealii::numbers::artificial_subdomain_id && - !cell->is_locally_owned_on_level()) + if (cell->is_ghost_on_level()) cell->set_user_flag(); } // Phase 1. Request all marked cells from corresponding owners. If we // managed to get every DoF, remove the user_flag, otherwise we // will request them again in the step below. - communicate_mg_ghost_cells(*triangulation, *dof_handler); + communicate_mg_ghost_cells(*dof_handler); // Phase 2, only request the cells that were not completed // in Phase 1. - communicate_mg_ghost_cells(*triangulation, *dof_handler); + communicate_mg_ghost_cells(*dof_handler); # ifdef DEBUG // make sure we have removed all flags: @@ -4256,229 +4243,119 @@ namespace internal // We start by checking whether only the numbering within the MPI - // ranks changed. In that case, we can apply the renumbering with some - // local renumbering only (this is similar to the renumber_mg_dofs() - // function below) - const bool locally_owned_set_changes = + // ranks changed, in which case we do not need to find a new index + // set. + const IndexSet &owned_dofs = dof_handler->locally_owned_dofs(); + const bool locally_owned_set_changes = std::any_of(new_numbers.cbegin(), new_numbers.cend(), - [this](const types::global_dof_index i) { - return dof_handler->locally_owned_dofs().is_element( - i) == false; + [&owned_dofs](const types::global_dof_index i) { + return owned_dofs.is_element(i) == false; }); - if (Utilities::MPI::sum(static_cast( - locally_owned_set_changes), - triangulation->get_communicator()) == 0) + IndexSet my_locally_owned_new_dof_indices = owned_dofs; + if (locally_owned_set_changes && owned_dofs.n_elements() > 0) { - // Since only the order within the local subdomains has changed, - // all we need to do is to propagate the knowledge about the - // numbers from the locally owned dofs (given by the new_numbers - // array) to all ghosted dofs on neighboring processors. We can do - // this by ghost layer exchange routines as in parallel vectors: - // We create an IndexSet for the relevant dofs and then export - // into an array of those values via Utilities::MPI::Partitioner. - IndexSet relevant_dofs; - DoFTools::extract_locally_relevant_dofs(*dof_handler, - relevant_dofs); - std::vector ghosted_new_numbers( - relevant_dofs.n_elements()); - { - Utilities::MPI::Partitioner partitioner( - dof_handler->locally_owned_dofs(), - relevant_dofs, - triangulation->get_communicator()); - - // choose some number that makes it unlikely to get conflicts - // with other ongoing non-blocking communication (there - // shouldn't be any at this place in most programs). - const unsigned int communication_channel = 19; - std::vector temp_array( - partitioner.n_import_indices()); - std::vector requests; - partitioner.export_to_ghosted_array_start( - communication_channel, - make_array_view(new_numbers), - make_array_view(temp_array), - ArrayView( - ghosted_new_numbers.data() + new_numbers.size(), - partitioner.n_ghost_indices()), - requests); - partitioner.export_to_ghosted_array_finish( - ArrayView( - ghosted_new_numbers.data() + new_numbers.size(), - partitioner.n_ghost_indices()), - requests); - - // we need to fill the indices of the locally owned part into - // the new numbers array, which is not provided by the parallel - // partitioner. their right position is somewhere in the middle - // of the array, so we first copy the ghosted part from smaller - // ranks to the front, then insert the data in the middle. - unsigned int n_ghosts_on_smaller_ranks = 0; - for (std::pair t : - partitioner.ghost_targets()) - { - if (t.first > partitioner.this_mpi_process()) - break; - n_ghosts_on_smaller_ranks += t.second; - } - if (n_ghosts_on_smaller_ranks > 0) - { - Assert(ghosted_new_numbers.data() != nullptr, - ExcInternalError()); - std::memmove(ghosted_new_numbers.data(), - ghosted_new_numbers.data() + new_numbers.size(), - sizeof(types::global_dof_index) * - n_ghosts_on_smaller_ranks); - } - if (new_numbers.size() > 0) - { - Assert(new_numbers.data() != nullptr, ExcInternalError()); - std::memcpy(ghosted_new_numbers.data() + - n_ghosts_on_smaller_ranks, - new_numbers.data(), - sizeof(types::global_dof_index) * - new_numbers.size()); - } - } + std::vector new_numbers_sorted = + new_numbers; + std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end()); - // In case we do not carry any relevant dof (but only some remote - // processor), we do not need to call the renumbering. We call the - // version without validity check because vertex dofs will be - // set already in the artificial region. - if (relevant_dofs.n_elements() > 0) - Implementation::renumber_dofs(ghosted_new_numbers, - relevant_dofs, - *dof_handler, - /*check_validity=*/false); - - NumberCache number_cache; - number_cache.locally_owned_dofs = dof_handler->locally_owned_dofs(); - number_cache.n_global_dofs = dof_handler->n_dofs(); - number_cache.n_locally_owned_dofs = - number_cache.locally_owned_dofs.n_elements(); - return number_cache; + my_locally_owned_new_dof_indices = IndexSet(dof_handler->n_dofs()); + my_locally_owned_new_dof_indices.add_indices( + new_numbers_sorted.begin(), new_numbers_sorted.end()); + my_locally_owned_new_dof_indices.compress(); + + Assert(my_locally_owned_new_dof_indices.n_elements() == + new_numbers.size(), + ExcInternalError()); } - else - { - // Now back to the more complicated case - // - // First figure out the new set of locally owned DoF indices. - // If we own no DoFs, we still need to go through this function, - // but we can skip this calculation. - // - // The IndexSet::add_indices() function is substantially more - // efficient if the set of indices is already sorted because - // it can then insert ranges instead of individual elements. - // consequently, pre-sort the array of new indices - IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs()); - if (dof_handler->n_locally_owned_dofs() > 0) + + // delete all knowledge of DoF indices that are not locally + // owned. we do so by getting DoF indices on cells, checking + // whether they are locally owned, if not, setting them to + // an invalid value, and then setting them again on the current + // cell + // + // DoFs we (i) know about, and (ii) don't own locally must be + // located either on ghost cells, or on the interface between a + // locally owned cell and a ghost cell. In any case, it is + // sufficient to kill them only from the ghost side cell, so loop + // only over ghost cells + { + std::vector local_dof_indices; + + for (auto cell : dof_handler->active_cell_iterators()) + if (cell->is_ghost()) { - std::vector - new_numbers_sorted = new_numbers; - std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end()); + local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(local_dof_indices); - my_locally_owned_new_dof_indices.add_indices( - new_numbers_sorted.begin(), new_numbers_sorted.end()); - my_locally_owned_new_dof_indices.compress(); + // delete a DoF index if it has not already been deleted + // (e.g., by visiting a neighboring cell, if it is on the + // boundary), and if we don't own it + for (types::global_dof_index &index : local_dof_indices) + if ((index != numbers::invalid_dof_index) && + (!owned_dofs.is_element(index))) + index = numbers::invalid_dof_index; - Assert(my_locally_owned_new_dof_indices.n_elements() == - new_numbers.size(), - ExcInternalError()); + cell->set_dof_indices(local_dof_indices); } + } - // delete all knowledge of DoF indices that are not locally - // owned. we do so by getting DoF indices on cells, checking - // whether they are locally owned, if not, setting them to - // an invalid value, and then setting them again on the current - // cell - // - // DoFs we (i) know about, and (ii) don't own locally must be - // located either on ghost cells, or on the interface between a - // locally owned cell and a ghost cell. In any case, it is - // sufficient to kill them only from the ghost side cell, so loop - // only over ghost cells - { - std::vector local_dof_indices; - for (auto cell : dof_handler->active_cell_iterators()) - if (cell->is_ghost()) - { - local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); - cell->get_dof_indices(local_dof_indices); - - for (unsigned int i = 0; - i < cell->get_fe().n_dofs_per_cell(); - ++i) - // delete a DoF index if it has not already been deleted - // (e.g., by visiting a neighboring cell, if it is on the - // boundary), and if we don't own it - if ((local_dof_indices[i] != - numbers::invalid_dof_index) && - (!dof_handler->locally_owned_dofs().is_element( - local_dof_indices[i]))) - local_dof_indices[i] = numbers::invalid_dof_index; - - cell->set_dof_indices(local_dof_indices); - } - } + // renumber. Skip when there is nothing to do because we own no DoF. + if (owned_dofs.n_elements() > 0) + Implementation::renumber_dofs(new_numbers, + owned_dofs, + *dof_handler, + /*check_validity=*/false, + /*update_cache=*/false); + // Communicate newly assigned DoF indices to other processors + // and get the same information for our own ghost cells. + // + // This is the same as phase 5+6 in the distribute_dofs() algorithm, + // taking into account that we have to unify a few DoFs in between + // then communication phases if we do hp-numbering + { + std::vector user_flags; + triangulation->save_user_flags(user_flags); + triangulation->clear_user_flags(); - // renumber. Skip when there is nothing to do because we own no DoF. - if (dof_handler->locally_owned_dofs().n_elements() > 0) - Implementation::renumber_dofs(new_numbers, - dof_handler->locally_owned_dofs(), - *dof_handler, - /*check_validity=*/false, - /*update_cache=*/false); + // mark all own cells for transfer + for (const auto &cell : dof_handler->active_cell_iterators()) + if (cell->is_ghost()) + cell->set_user_flag(); - // Communicate newly assigned DoF indices to other processors - // and get the same information for our own ghost cells. - // - // This is the same as phase 5+6 in the distribute_dofs() algorithm, - // taking into account that we have to unify a few DoFs in between - // then communication phases if we do hp-numbering - { - std::vector user_flags; - triangulation->save_user_flags(user_flags); - triangulation->clear_user_flags(); - - // mark all own cells for transfer - for (const auto &cell : dof_handler->active_cell_iterators()) - if (cell->is_ghost()) - cell->set_user_flag(); - - - // Send and receive cells. After this, only the local cells - // are marked, that received new data. This has to be - // communicated in a second communication step. - // - // as explained in the 'distributed' paper, this has to be - // done twice - communicate_dof_indices_on_marked_cells(*dof_handler); - - // if the DoFHandler has hp-capabilities then we may have - // received valid indices of degrees of freedom that are - // dominated by a FE object adjacent to a ghost interface. - // thus, we overwrite the remaining invalid indices with the - // valid ones in this step. - Implementation::merge_invalid_dof_indices_on_ghost_interfaces( - *dof_handler); - - communicate_dof_indices_on_marked_cells(*dof_handler); - - triangulation->load_user_flags(user_flags); - update_all_active_cell_dof_indices_caches(*this->dof_handler); - } - NumberCache number_cache; - number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices; - number_cache.n_global_dofs = dof_handler->n_dofs(); - number_cache.n_locally_owned_dofs = - number_cache.locally_owned_dofs.n_elements(); - return number_cache; - } + // Send and receive cells. After this, only the local cells + // are marked, that received new data. This has to be + // communicated in a second communication step. + // + // as explained in the 'distributed' paper, this has to be + // done twice + communicate_dof_indices_on_marked_cells(*dof_handler); + + // if the DoFHandler has hp-capabilities then we may have + // received valid indices of degrees of freedom that are + // dominated by a FE object adjacent to a ghost interface. + // thus, we overwrite the remaining invalid indices with the + // valid ones in this step. + Implementation::merge_invalid_dof_indices_on_ghost_interfaces( + *dof_handler); + + communicate_dof_indices_on_marked_cells(*dof_handler); + + triangulation->load_user_flags(user_flags); + update_all_active_cell_dof_indices_caches(*this->dof_handler); + } + + NumberCache number_cache; + number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices; + number_cache.n_global_dofs = dof_handler->n_dofs(); + number_cache.n_locally_owned_dofs = + number_cache.locally_owned_dofs.n_elements(); + return number_cache; #endif } @@ -4490,110 +4367,105 @@ namespace internal const unsigned int level, const std::vector &new_numbers) const { - // we only implement the case where the multigrid numbers are - // renumbered within the processor's partition, rather than the most - // general case - const IndexSet index_set = dof_handler->locally_owned_mg_dofs(level); +#ifndef DEAL_II_WITH_MPI -#ifdef DEAL_II_WITH_MPI + (void)level; + (void)new_numbers; - const dealii::parallel::TriangulationBase *tr = - (dynamic_cast - *>(&this->dof_handler->get_triangulation())); - Assert(tr != nullptr, ExcInternalError()); + Assert(false, ExcNotImplemented()); + return NumberCache(); +#else - const unsigned int my_rank = - Utilities::MPI::this_mpi_process(tr->get_communicator()); + dealii::parallel::DistributedTriangulationBase + *triangulation = + (dynamic_cast< + dealii::parallel::DistributedTriangulationBase *>( + const_cast *>( + &dof_handler->get_triangulation()))); + Assert(triangulation != nullptr, ExcInternalError()); -# ifdef DEBUG - for (types::global_dof_index i : new_numbers) + // This code is very close to the respective code in renumber_dofs, + // with the difference that we work on different entities with + // different objects. + const IndexSet &owned_dofs = dof_handler->locally_owned_mg_dofs(level); + AssertDimension(new_numbers.size(), owned_dofs.n_elements()); + + const bool locally_owned_set_changes = + std::any_of(new_numbers.cbegin(), + new_numbers.cend(), + [&owned_dofs](const types::global_dof_index i) { + return owned_dofs.is_element(i) == false; + }); + + IndexSet my_locally_owned_new_dof_indices = owned_dofs; + if (locally_owned_set_changes && owned_dofs.n_elements() > 0) { - Assert(index_set.is_element(i), - ExcNotImplemented( - "Renumberings that change the locally owned mg dofs " - "partitioning are currently not implemented for " - "the multigrid levels")); + std::vector new_numbers_sorted = + new_numbers; + std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end()); + + my_locally_owned_new_dof_indices = + IndexSet(dof_handler->n_dofs(level)); + my_locally_owned_new_dof_indices.add_indices( + new_numbers_sorted.begin(), new_numbers_sorted.end()); + my_locally_owned_new_dof_indices.compress(); + + Assert(my_locally_owned_new_dof_indices.n_elements() == + new_numbers.size(), + ExcInternalError()); } -# endif - // we need to access all locally relevant degrees of freedom. we - // use Utilities::MPI::Partitioner for handling the data exchange - // of the new numbers, which is simply the extraction of ghost data - IndexSet relevant_dofs; - DoFTools::extract_locally_relevant_level_dofs(*dof_handler, - level, - relevant_dofs); - std::vector ghosted_new_numbers( - relevant_dofs.n_elements()); + // delete all knowledge of DoF indices that are not locally + // owned { - Utilities::MPI::Partitioner partitioner(index_set, - relevant_dofs, - tr->get_communicator()); - std::vector temp_array( - partitioner.n_import_indices()); - const unsigned int communication_channel = 17; - std::vector requests; - partitioner.export_to_ghosted_array_start( - communication_channel, - make_array_view(new_numbers), - make_array_view(temp_array), - ArrayView(ghosted_new_numbers.data() + - new_numbers.size(), - partitioner.n_ghost_indices()), - requests); - partitioner.export_to_ghosted_array_finish( - ArrayView(ghosted_new_numbers.data() + - new_numbers.size(), - partitioner.n_ghost_indices()), - requests); - - // we need to fill the indices of the locally owned part into the - // new numbers array. their right position is somewhere in the - // middle of the array, so we first copy the ghosted part from - // smaller ranks to the front, then insert the data in the middle. - unsigned int n_ghosts_on_smaller_ranks = 0; - for (std::pair t : - partitioner.ghost_targets()) - { - if (t.first > my_rank) - break; - n_ghosts_on_smaller_ranks += t.second; - } - if (n_ghosts_on_smaller_ranks > 0) - { - Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError()); - std::memmove(ghosted_new_numbers.data(), - ghosted_new_numbers.data() + new_numbers.size(), - sizeof(types::global_dof_index) * - n_ghosts_on_smaller_ranks); - } - if (new_numbers.size() > 0) - { - Assert(new_numbers.data() != nullptr, ExcInternalError()); - std::memcpy(ghosted_new_numbers.data() + - n_ghosts_on_smaller_ranks, - new_numbers.data(), - sizeof(types::global_dof_index) * new_numbers.size()); - } + std::vector local_dof_indices; + + for (auto cell : dof_handler->cell_iterators_on_level(level)) + if (cell->is_ghost_on_level()) + { + local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_mg_dof_indices(local_dof_indices); + + for (types::global_dof_index &index : local_dof_indices) + if ((index != numbers::invalid_dof_index) && + (!owned_dofs.is_element(index))) + index = numbers::invalid_dof_index; + + cell->set_mg_dof_indices(local_dof_indices); + } } - // in case we do not own any of the given level (but only some remote - // processor), we do not need to call the renumbering - if (level < this->dof_handler->get_triangulation().n_levels() && - relevant_dofs.n_elements() > 0) + + // renumber. Skip when there is nothing to do because we own no DoF. + if (level < triangulation->n_levels() && owned_dofs.n_elements() > 0) Implementation::renumber_mg_dofs( - ghosted_new_numbers, relevant_dofs, *dof_handler, level, true); -#else - (void)new_numbers; - Assert(false, ExcNotImplemented()); -#endif + new_numbers, owned_dofs, *dof_handler, level, false); + + // communicate newly assigned DoF indices with other processors + { + std::vector user_flags; + triangulation->save_user_flags(user_flags); + triangulation->clear_user_flags(); + + // mark only cells on the chosen level + for (const auto &cell : dof_handler->cell_iterators_on_level(level)) + if (cell->is_ghost_on_level()) + cell->set_user_flag(); + + communicate_mg_ghost_cells(*dof_handler); + + communicate_mg_ghost_cells(*dof_handler); + + triangulation->load_user_flags(user_flags); + } NumberCache number_cache; - number_cache.locally_owned_dofs = index_set; + number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices; number_cache.n_global_dofs = dof_handler->n_dofs(level); number_cache.n_locally_owned_dofs = number_cache.locally_owned_dofs.n_elements(); return number_cache; +#endif } } // namespace Policy } // namespace DoFHandlerImplementation