From e3d094b3a4acdcd12d9ec344e29125e4d8462ea0 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 22 Mar 2023 20:34:12 +0100 Subject: [PATCH] Communicate within guess_point_owner --- source/base/mpi_remote_point_evaluation.cc | 8 +- source/grid/grid_tools.cc | 102 ++++++++++++++++----- 2 files changed, 80 insertions(+), 30 deletions(-) diff --git a/source/base/mpi_remote_point_evaluation.cc b/source/base/mpi_remote_point_evaluation.cc index c5c7a5b759..debb3128f6 100644 --- a/source/base/mpi_remote_point_evaluation.cc +++ b/source/base/mpi_remote_point_evaluation.cc @@ -90,12 +90,8 @@ namespace Utilities const auto local_tree = pack_rtree(local_boxes); // compress r-tree to a minimal set of bounding boxes - const auto local_reduced_box = - extract_rtree_level(local_tree, rtree_level); - - // gather bounding boxes of other processes - const auto global_bboxes = - Utilities::MPI::all_gather(tria.get_communicator(), local_reduced_box); + std::vector>> global_bboxes(1); + global_bboxes[0] = extract_rtree_level(local_tree, rtree_level); const GridTools::Cache cache(tria, mapping); diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 78b73d0097..c6854b732d 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -5976,28 +5976,82 @@ namespace GridTools namespace internal { - template + /** + * Determine for each rank which entry of @p entities it + * might own. The first entry of the returned tuple is a list of + * ranks and the second and third entry give CRS data + * structure (pointers within a list of indices). + */ + template std::tuple, std::vector, std::vector> - guess_point_owner( + guess_owners_of_entities( + const MPI_Comm & comm, const std::vector>> &global_bboxes, - const std::vector> & points, + const std::vector & entities, const double tolerance) { + std::vector>> global_bboxes_temp; + auto *global_bboxes_to_be_used = &global_bboxes; + + if (global_bboxes.size() == 1) // TODO: and not ArborX installed + { + global_bboxes_temp = + Utilities::MPI::all_gather(comm, global_bboxes[0]); + global_bboxes_to_be_used = &global_bboxes_temp; + } + std::vector> ranks_and_indices; - ranks_and_indices.reserve(points.size()); + ranks_and_indices.reserve(entities.size()); - for (unsigned int i = 0; i < points.size(); ++i) + if (true) { - const auto &point = points[i]; - for (unsigned rank = 0; rank < global_bboxes.size(); ++rank) - for (const auto &box : global_bboxes[rank]) - if (box.point_inside(point, tolerance)) - { - ranks_and_indices.emplace_back(rank, i); - break; - } + // helper function to determine if a bounding box is valid + const auto is_valid = [](const auto &bb) { + for (unsigned int i = 0; i < spacedim; ++i) + if (bb.get_boundary_points().first[i] > + bb.get_boundary_points().second[i]) + return false; + + return true; + }; + + // linearize vector of vectors + std::vector, unsigned int>> + boxes_and_ranks; + + for (unsigned rank = 0; rank < global_bboxes_to_be_used->size(); + ++rank) + for (const auto &box : (*global_bboxes_to_be_used)[rank]) + if (is_valid(box)) + boxes_and_ranks.emplace_back(box, rank); + + // pack boxes into r-tree + const auto tree = pack_rtree(boxes_and_ranks); + + // loop over all entities + for (unsigned int i = 0; i < entities.size(); ++i) + { + // create a bounding box with tolerance + const auto bb = + BoundingBox(entities[i]).create_extended(tolerance); + + // determine ranks potentially owning point/bounding box + std::set my_ranks; + + for (const auto &box_and_rank : + tree | boost::geometry::index::adaptors::queried( + boost::geometry::index::intersects(bb))) + my_ranks.insert(box_and_rank.second); + + for (const auto rank : my_ranks) + ranks_and_indices.emplace_back(rank, i); + } + } + else + { + // TODO: use ArborX } // convert to CRS @@ -6007,14 +6061,14 @@ namespace GridTools std::vector ptr; std::vector indices; - unsigned int dummy_rank = numbers::invalid_unsigned_int; + unsigned int current_rank = numbers::invalid_unsigned_int; - for (const auto &i : ranks_and_indices) + for (const std::pair &i : ranks_and_indices) { - if (dummy_rank != i.first) + if (current_rank != i.first) { - dummy_rank = i.first; - ranks.push_back(dummy_rank); + current_rank = i.first; + ranks.push_back(current_rank); ptr.push_back(indices.size()); } @@ -6022,9 +6076,7 @@ namespace GridTools } ptr.push_back(indices.size()); - return std::make_tuple(std::move(ranks), - std::move(ptr), - std::move(indices)); + return {std::move(ranks), std::move(ptr), std::move(indices)}; } @@ -6128,8 +6180,10 @@ namespace GridTools auto &recv_ranks = result.recv_ranks; auto &recv_ptrs = result.recv_ptrs; - const auto potential_owners = - internal::guess_point_owner(global_bboxes, points, tolerance); + const auto comm = cache.get_triangulation().get_communicator(); + + const auto potential_owners = internal::guess_owners_of_entities( + comm, global_bboxes, points, tolerance); const auto &potential_owners_ranks = std::get<0>(potential_owners); const auto &potential_owners_ptrs = std::get<1>(potential_owners); @@ -6264,7 +6318,7 @@ namespace GridTools create_request, answer_request, process_answer, - cache.get_triangulation().get_communicator()); + comm); if (true) { -- 2.39.5