From cad9769a35bc3bb4c4d643a0a9989010d6d331dd Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 25 Jun 2021 18:17:44 +0200 Subject: [PATCH] Make mapping in GridTools::distributed_compute_point_locations() unique --- source/grid/grid_tools.cc | 182 +++--------------- tests/remote_point_evaluation/mapping_02.cc | 106 ++++++++++ .../mapping_02.mpirun=1.output | 2 + .../mapping_02.mpirun=4.output | 2 + 4 files changed, 139 insertions(+), 153 deletions(-) create mode 100644 tests/remote_point_evaluation/mapping_02.cc create mode 100644 tests/remote_point_evaluation/mapping_02.mpirun=1.output create mode 100644 tests/remote_point_evaluation/mapping_02.mpirun=4.output diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 7b973c2f36..240c194d18 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -5640,7 +5640,7 @@ namespace GridTools { // run internal function ... const auto all = internal::distributed_compute_point_locations( - cache, points, global_bboxes, tolerance, false) + cache, points, global_bboxes, tolerance, false, true) .send_components; // ... and reshuffle the data @@ -5749,7 +5749,8 @@ namespace GridTools const Point & point, typename Triangulation::active_cell_iterator &cell_hint, const std::vector &marked_vertices, - const double tolerance) + const double tolerance, + const bool enforce_unique_mapping) { std::vector< std::pair::active_cell_iterator, @@ -5768,6 +5769,9 @@ namespace GridTools tolerance, &cache.get_locally_owned_cell_bounding_boxes_rtree()); + const unsigned int my_rank = Utilities::MPI::this_mpi_process( + cache.get_triangulation().get_communicator()); + cell_hint = first_cell.first; if (cell_hint.state() == IteratorState::valid) { @@ -5779,6 +5783,20 @@ namespace GridTools tolerance, first_cell); + if (enforce_unique_mapping) + { + // check if the rank of this process is the lowest of all cells + // if not, the other process will handle this cell and we don't + // have to do here anything in the case of unique mapping + unsigned int lowes_rank = numbers::invalid_unsigned_int; + + for (const auto &cell : active_cells_around_point) + lowes_rank = std::min(lowes_rank, cell.first->subdomain_id()); + + if (lowes_rank != my_rank) + return {}; + } + locally_owned_active_cells_around_point.reserve( active_cells_around_point.size()); @@ -5791,7 +5809,12 @@ namespace GridTools locally_owned_active_cells_around_point.end(), [](const auto &a, const auto &b) { return a.first < b.first; }); - return locally_owned_active_cells_around_point; + if (enforce_unique_mapping && + locally_owned_active_cells_around_point.size() > 1) + // in the case of unique mapping, we only need a single cell + return {locally_owned_active_cells_around_point.front()}; + else + return locally_owned_active_cells_around_point; } @@ -5806,8 +5829,6 @@ namespace GridTools const bool perform_handshake, const bool enforce_unique_mapping) { - Assert(!enforce_unique_mapping || perform_handshake, ExcInternalError()); - DistributedComputePointLocationsInternal result; auto &send_components = result.send_components; @@ -5879,7 +5900,8 @@ namespace GridTools index_and_point.second, cell_hint, marked_vertices, - tolerance); + tolerance, + enforce_unique_mapping); for (const auto &cell_and_reference_position : cells_and_reference_positions) @@ -5893,18 +5915,9 @@ namespace GridTools cell_and_reference_position.second, index_and_point.second, numbers::invalid_unsigned_int); - - if (enforce_unique_mapping) - break; // in the case of unique mapping, we only need a - // single cell (we take the first) } - if (perform_handshake) - request_buffer_temp[i] = - enforce_unique_mapping ? - std::min( - 1, cells_and_reference_positions.size()) : - cells_and_reference_positions.size(); + request_buffer_temp[i] = cells_and_reference_positions.size(); } if (perform_handshake) @@ -5946,143 +5959,6 @@ namespace GridTools process, cache.get_triangulation().get_communicator()) .run(); - // for unique mapping, we need to modify recv_components and - // send_components consistently - if (enforce_unique_mapping) - { - std::vector mask_recv(recv_components.size()); - std::vector mask_send(send_components.size()); - - // set up new recv_components and monitor which entries (we keep - // the entry of the lowest rank) have been eliminated so that we can - // communicate it - auto recv_components_copy = recv_components; - recv_components.clear(); - - for (unsigned int i = 0; i < recv_components_copy.size(); ++i) - std::get<2>(recv_components_copy[i]) = i; - - std::sort(recv_components_copy.begin(), - recv_components_copy.end(), - [&](const auto &a, const auto &b) { - if (std::get<0>(a) != std::get<0>(b)) // rank - return std::get<0>(a) < std::get<0>(b); - - return std::get<2>(a) < std::get<2>(b); // enumeration - }); - - std::vector unique(points.size(), false); - - std::vector recv_ranks; - std::vector recv_ptrs; - - for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int; - i < recv_components_copy.size(); - ++i) - { - if (dummy != std::get<0>(recv_components_copy[i])) - { - dummy = std::get<0>(recv_components_copy[i]); - recv_ranks.push_back(dummy); - recv_ptrs.push_back(i); - } - - if (unique[std::get<1>(recv_components_copy[i])] == false) - { - recv_components.emplace_back(recv_components_copy[i]); - mask_recv[i] = 1; - unique[std::get<1>(recv_components_copy[i])] = true; - } - else - { - mask_recv[i] = 0; - } - } - recv_ptrs.push_back(recv_components_copy.size()); - - Assert(std::all_of(unique.begin(), - unique.end(), - [](const auto &v) { return v; }), - ExcInternalError()); - - - // prepare send_components so that not needed entries can be - // eliminated later on - auto send_components_copy = send_components; - send_components.clear(); - - for (unsigned int i = 0; i < send_components_copy.size(); ++i) - std::get<5>(send_components_copy[i]) = i; - - std::sort(send_components_copy.begin(), - send_components_copy.end(), - [&](const auto &a, const auto &b) { - if (std::get<1>(a) != std::get<1>(b)) // rank - return std::get<1>(a) < std::get<1>(b); - - return std::get<5>(a) < std::get<5>(b); // enumeration - }); - - std::vector send_ranks; - std::vector send_ptrs; - - for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int; - i < send_components_copy.size(); - ++i) - { - if (dummy != std::get<1>(send_components_copy[i])) - { - dummy = std::get<1>(send_components_copy[i]); - send_ranks.push_back(dummy); - send_ptrs.push_back(i); - } - } - send_ptrs.push_back(send_components_copy.size()); - - // perform communication -#ifdef DEAL_II_WITH_MPI - std::vector req(send_ranks.size() + recv_ranks.size()); - - for (unsigned int i = 0; i < send_ranks.size(); ++i) - { - const auto ierr = - MPI_Irecv(mask_send.data() + send_ptrs[i], - send_ptrs[i + 1] - send_ptrs[i], - MPI_UNSIGNED, - send_ranks[i], - Utilities::MPI::internal::Tags:: - distributed_compute_point_locations, - cache.get_triangulation().get_communicator(), - &req[i]); - AssertThrowMPI(ierr); - } - - for (unsigned int i = 0; i < recv_ranks.size(); ++i) - { - const auto ierr = - MPI_Isend(mask_recv.data() + recv_ptrs[i], - recv_ptrs[i + 1] - recv_ptrs[i], - MPI_UNSIGNED, - recv_ranks[i], - Utilities::MPI::internal::Tags:: - distributed_compute_point_locations, - cache.get_triangulation().get_communicator(), - &req[i] + send_ranks.size()); - AssertThrowMPI(ierr); - } - - auto ierr = MPI_Waitall(req.size(), req.data(), MPI_STATUSES_IGNORE); - AssertThrowMPI(ierr); -#else - mask_send = mask_recv; -#endif - - // eliminate not needed entries - for (unsigned int i = 0; i < send_components_copy.size(); ++i) - if (mask_send[i] == 1) - send_components.emplace_back(send_components_copy[i]); - } - if (true) { // sort according to rank (and point index and cell) -> make diff --git a/tests/remote_point_evaluation/mapping_02.cc b/tests/remote_point_evaluation/mapping_02.cc new file mode 100644 index 0000000000..29f4e9dbc2 --- /dev/null +++ b/tests/remote_point_evaluation/mapping_02.cc @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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 GridTools::distributed_compute_point_locations() that it works with +// unique mapping. + +#include + +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::subdivided_hyper_cube(tria, 2); + tria.refine_global(1); + + FE_DGQ fe(0); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vec(dof_handler.n_dofs()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + vec[cell->global_active_cell_index()] = cell->global_active_cell_index(); + + MappingQ1 mapping; + + std::vector> local_boxes; + for (const auto &cell : tria.cell_iterators_on_level(0)) + local_boxes.push_back(mapping.get_bounding_box(cell)); + + // create r-tree of bounding boxes + 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, 0); + + // gather bounding boxes of other processes + const auto global_bboxes = + Utilities::MPI::all_gather(tria.get_communicator(), local_reduced_box); + + const GridTools::Cache cache(tria, mapping); + + std::vector> local_points; + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + for (unsigned int j = 0; j <= 4; ++j) + for (unsigned int i = 0; i <= 4; ++i) + local_points.emplace_back(i * 0.25, j * 0.25); + + const auto temp = + GridTools::distributed_compute_point_locations(cache, + local_points, + global_bboxes); + + unsigned int count = 0; + + for (const auto &i : std::get<1>(temp)) + count += i.size(); + + const auto global_count = Utilities::MPI::sum(count, MPI_COMM_WORLD); + const auto expected_global_count = + Utilities::MPI::sum(local_points.size(), MPI_COMM_WORLD); + + deallog << global_count << " " << expected_global_count << std::endl; +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + initlog(); + + test<2>(); +} diff --git a/tests/remote_point_evaluation/mapping_02.mpirun=1.output b/tests/remote_point_evaluation/mapping_02.mpirun=1.output new file mode 100644 index 0000000000..62654b94f7 --- /dev/null +++ b/tests/remote_point_evaluation/mapping_02.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL::25 25 diff --git a/tests/remote_point_evaluation/mapping_02.mpirun=4.output b/tests/remote_point_evaluation/mapping_02.mpirun=4.output new file mode 100644 index 0000000000..62654b94f7 --- /dev/null +++ b/tests/remote_point_evaluation/mapping_02.mpirun=4.output @@ -0,0 +1,2 @@ + +DEAL::25 25 -- 2.39.5