mpi_remote_point_evaluation.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2021 - 2023 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.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_mpi_mpi_remote_point_evaluation_h
#define dealii_mpi_mpi_remote_point_evaluation_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/mpi.h>
#include <deal.II/base/mpi_tags.h>
 
#include <deal.II/dofs/dof_handler.h>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
namespace GridTools
{
  namespace internal
  {
    template <int dim, int spacedim>
    struct DistributedComputePointLocationsInternal;
  }
} // namespace GridTools
 
namespace Utilities
{
  namespace MPI
  {
    template <int dim, int spacedim = dim>
    class RemotePointEvaluation
    {
    public:
      RemotePointEvaluation(
        const double       tolerance                              = 1e-6,
        const bool         enforce_unique_mapping                 = false,
        const unsigned int rtree_level                            = 0,
        const std::function<std::vector<bool>()> &marked_vertices = {});
 
      ~RemotePointEvaluation();
 
      void
      reinit(const std::vector<Point<spacedim>> &points,
             const Triangulation<dim, spacedim> &tria,
             const Mapping<dim, spacedim> &      mapping);
 
      void
      reinit(const GridTools::internal::
               DistributedComputePointLocationsInternal<dim, spacedim> &data,
             const Triangulation<dim, spacedim> &                       tria,
             const Mapping<dim, spacedim> &mapping);
 
      struct CellData
      {
        std::vector<std::pair<int, int>> cells;
 
        std::vector<unsigned int> reference_point_ptrs;
 
        std::vector<Point<dim>> reference_point_values;
      };
 
      const CellData &
      get_cell_data() const;
 
      template <typename T>
      void
      evaluate_and_process(
        std::vector<T> &output,
        std::vector<T> &buffer,
        const std::function<void(const ArrayView<T> &, const CellData &)>
          &evaluation_function) const;
 
      template <typename T>
      void
      process_and_evaluate(
        const std::vector<T> &input,
        std::vector<T> &      buffer,
        const std::function<void(const ArrayView<const T> &, const CellData &)>
          &evaluation_function) const;
 
      const std::vector<unsigned int> &
      get_point_ptrs() const;
 
      bool
      is_map_unique() const;
 
      bool
      all_points_found() const;
 
      bool
      point_found(const unsigned int i) const;
 
      const Triangulation<dim, spacedim> &
      get_triangulation() const;
 
      const Mapping<dim, spacedim> &
      get_mapping() const;
 
      bool
      is_ready() const;
 
    private:
      const double tolerance;
 
      const bool enforce_unique_mapping;
 
      const unsigned int rtree_level;
 
      const std::function<std::vector<bool>()> marked_vertices;
 
      boost::signals2::connection tria_signal;
 
      bool ready_flag;
 
      SmartPointer<const Triangulation<dim, spacedim>> tria;
 
      SmartPointer<const Mapping<dim, spacedim>> mapping;
 
      bool unique_mapping;
 
      bool all_points_found_flag;
 
      std::vector<unsigned int> point_ptrs;
 
      std::vector<unsigned int> recv_permutation;
 
      std::vector<unsigned int> recv_ptrs;
 
      std::vector<unsigned int> recv_ranks;
 
      CellData cell_data;
 
      std::vector<unsigned int> send_permutation;
 
      std::vector<unsigned int> send_ranks;
 
      std::vector<unsigned int> send_ptrs;
    };
 
 
    template <int dim, int spacedim>
    template <typename T>
    void
    RemotePointEvaluation<dim, spacedim>::evaluate_and_process(
      std::vector<T> &output,
      std::vector<T> &buffer,
      const std::function<void(const ArrayView<T> &, const CellData &)>
        &evaluation_function) const
    {
#ifndef DEAL_II_WITH_MPI
      Assert(false, ExcNeedsMPI());
      (void)output;
      (void)buffer;
      (void)evaluation_function;
#else
      static CollectiveMutex      mutex;
      CollectiveMutex::ScopedLock lock(mutex, tria->get_communicator());
 
      const unsigned int my_rank =
        Utilities::MPI::this_mpi_process(tria->get_communicator());
 
      // allocate memory for output and buffer
      output.resize(point_ptrs.back());
      buffer.resize(std::max(send_permutation.size() * 2,
                             point_ptrs.back() + send_permutation.size()));
 
      // ... for evaluation
      ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
 
      // ... for communication
      ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
                               send_permutation.size());
 
      // evaluate functions at points
      evaluation_function(buffer_eval, cell_data);
 
      // sort for communication
      unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
      unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
 
      const auto my_rank_local_recv_ptr =
        std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
 
      if (my_rank_local_recv_ptr != recv_ranks.end())
        {
          my_rank_local_recv =
            std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
          my_rank_local_send = std::distance(send_ranks.begin(),
                                             std::find(send_ranks.begin(),
                                                       send_ranks.end(),
                                                       my_rank));
        }
 
      for (unsigned int i = 0; i < send_permutation.size(); ++i)
        {
          const unsigned int send_index = send_permutation[i];
 
          // local data -> can be copied to output directly
          if (my_rank_local_send != numbers::invalid_unsigned_int &&
              (send_ptrs[my_rank_local_send] <= send_index &&
               send_index < send_ptrs[my_rank_local_send + 1]))
            output[recv_permutation[send_index - send_ptrs[my_rank_local_send] +
                                    recv_ptrs[my_rank_local_recv]]] =
              buffer_eval[i];
          else // data to be sent
            buffer_comm[send_index] = buffer_eval[i];
        }
 
      // send data
      std::vector<std::vector<char>> send_buffer;
      send_buffer.reserve(send_ranks.size());
 
      std::vector<MPI_Request> send_requests;
      send_requests.reserve(send_ranks.size());
 
      for (unsigned int i = 0; i < send_ranks.size(); ++i)
        {
          if (send_ranks[i] == my_rank)
            continue;
 
          send_requests.emplace_back(MPI_Request());
 
          send_buffer.emplace_back(Utilities::pack(
            std::vector<T>(buffer_comm.begin() + send_ptrs[i],
                           buffer_comm.begin() + send_ptrs[i + 1]),
            false));
 
          const int ierr = MPI_Isend(send_buffer.back().data(),
                                     send_buffer.back().size(),
                                     MPI_CHAR,
                                     send_ranks[i],
                                     internal::Tags::remote_point_evaluation,
                                     tria->get_communicator(),
                                     &send_requests.back());
          AssertThrowMPI(ierr);
        }
 
      // receive data
      std::vector<char> buffer_char;
 
      for (unsigned int i = 0; i < recv_ranks.size(); ++i)
        {
          if (recv_ranks[i] == my_rank)
            continue;
 
          // receive remote data
          MPI_Status status;
          int        ierr = MPI_Probe(MPI_ANY_SOURCE,
                               internal::Tags::remote_point_evaluation,
                               tria->get_communicator(),
                               &status);
          AssertThrowMPI(ierr);
 
          int message_length;
          ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
          AssertThrowMPI(ierr);
 
          buffer_char.resize(message_length);
 
          ierr = MPI_Recv(buffer_char.data(),
                          buffer_char.size(),
                          MPI_CHAR,
                          status.MPI_SOURCE,
                          internal::Tags::remote_point_evaluation,
                          tria->get_communicator(),
                          MPI_STATUS_IGNORE);
          AssertThrowMPI(ierr);
 
          // unpack data
          const auto buffer =
            Utilities::unpack<std::vector<T>>(buffer_char, false);
 
          // write data into output vector
          const auto ptr =
            std::find(recv_ranks.begin(), recv_ranks.end(), status.MPI_SOURCE);
 
          Assert(ptr != recv_ranks.end(), ExcNotImplemented());
 
          const unsigned int j = std::distance(recv_ranks.begin(), ptr);
 
          AssertDimension(buffer.size(), recv_ptrs[j + 1] - recv_ptrs[j]);
 
          for (unsigned int i = recv_ptrs[j], c = 0; i < recv_ptrs[j + 1];
               ++i, ++c)
            output[recv_permutation[i]] = buffer[c];
        }
 
      // make sure all messages have been sent
      const int ierr = MPI_Waitall(send_requests.size(),
                                   send_requests.data(),
                                   MPI_STATUSES_IGNORE);
      AssertThrowMPI(ierr);
#endif
    }
 
 
    template <int dim, int spacedim>
    template <typename T>
    void
    RemotePointEvaluation<dim, spacedim>::process_and_evaluate(
      const std::vector<T> &input,
      std::vector<T> &      buffer,
      const std::function<void(const ArrayView<const T> &, const CellData &)>
        &evaluation_function) const
    {
#ifndef DEAL_II_WITH_MPI
      Assert(false, ExcNeedsMPI());
      (void)input;
      (void)buffer;
      (void)evaluation_function;
#else
      static CollectiveMutex      mutex;
      CollectiveMutex::ScopedLock lock(mutex, tria->get_communicator());
 
      const unsigned int my_rank =
        Utilities::MPI::this_mpi_process(tria->get_communicator());
 
      // invert permutation matrices (TODO: precompute)
      std::vector<unsigned int> recv_permutation_inv(recv_permutation.size());
      for (unsigned int c = 0; c < recv_permutation.size(); ++c)
        recv_permutation_inv[recv_permutation[c]] = c;
 
      std::vector<unsigned int> send_permutation_inv(send_permutation.size());
      for (unsigned int c = 0; c < send_permutation.size(); ++c)
        send_permutation_inv[send_permutation[c]] = c;
 
      // allocate memory for buffer
      const auto &point_ptrs = this->get_point_ptrs();
      AssertDimension(input.size(), point_ptrs.size() - 1);
      buffer.resize(std::max(send_permutation.size() * 2,
                             point_ptrs.back() + send_permutation.size()));
 
      // ... for evaluation
      ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
 
      // ... for communication
      ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
                               point_ptrs.back());
 
      // sort for communication (and duplicate data if necessary)
      unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
      unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
 
      const auto my_rank_local_recv_ptr =
        std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
 
      if (my_rank_local_recv_ptr != recv_ranks.end())
        {
          my_rank_local_recv =
            std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
          my_rank_local_send = std::distance(send_ranks.begin(),
                                             std::find(send_ranks.begin(),
                                                       send_ranks.end(),
                                                       my_rank));
        }
 
      for (unsigned int i = 0, c = 0; i < point_ptrs.size() - 1; ++i)
        {
          const auto n_entries = point_ptrs[i + 1] - point_ptrs[i];
 
          for (unsigned int j = 0; j < n_entries; ++j, ++c)
            {
              const unsigned int recv_index = recv_permutation_inv[c];
 
              // local data -> can be copied to final buffer directly
              if (my_rank_local_recv != numbers::invalid_unsigned_int &&
                  (recv_ptrs[my_rank_local_recv] <= recv_index &&
                   recv_index < recv_ptrs[my_rank_local_recv + 1]))
                buffer_eval[send_permutation_inv
                              [recv_index - recv_ptrs[my_rank_local_recv] +
                               send_ptrs[my_rank_local_send]]] = input[i];
              else // data to be sent
                buffer_comm[recv_index] = input[i];
            }
        }
 
      // send data
      std::vector<std::vector<char>> send_buffer;
      send_buffer.reserve(recv_ranks.size());
 
      std::vector<MPI_Request> send_requests;
      send_requests.reserve(recv_ranks.size());
 
      for (unsigned int i = 0; i < recv_ranks.size(); ++i)
        {
          if (recv_ranks[i] == my_rank)
            continue;
 
          send_requests.push_back(MPI_Request());
 
          send_buffer.emplace_back(Utilities::pack(
            std::vector<T>(buffer_comm.begin() + recv_ptrs[i],
                           buffer_comm.begin() + recv_ptrs[i + 1]),
            false));
 
          const int ierr = MPI_Isend(send_buffer.back().data(),
                                     send_buffer.back().size(),
                                     MPI_CHAR,
                                     recv_ranks[i],
                                     internal::Tags::remote_point_evaluation,
                                     tria->get_communicator(),
                                     &send_requests.back());
          AssertThrowMPI(ierr);
        }
 
      // receive data
      std::vector<char> recv_buffer;
 
      for (unsigned int i = 0; i < send_ranks.size(); ++i)
        {
          if (send_ranks[i] == my_rank)
            continue;
 
          // receive remote data
          MPI_Status status;
          int        ierr = MPI_Probe(MPI_ANY_SOURCE,
                               internal::Tags::remote_point_evaluation,
                               tria->get_communicator(),
                               &status);
          AssertThrowMPI(ierr);
 
          int message_length;
          ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
          AssertThrowMPI(ierr);
 
          recv_buffer.resize(message_length);
 
          ierr = MPI_Recv(recv_buffer.data(),
                          recv_buffer.size(),
                          MPI_CHAR,
                          status.MPI_SOURCE,
                          internal::Tags::remote_point_evaluation,
                          tria->get_communicator(),
                          MPI_STATUS_IGNORE);
          AssertThrowMPI(ierr);
 
          // unpack data
          const auto recv_buffer_unpacked =
            Utilities::unpack<std::vector<T>>(recv_buffer, false);
 
          // write data into buffer vector
          const auto ptr =
            std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
 
          Assert(ptr != send_ranks.end(), ExcNotImplemented());
 
          const unsigned int j = std::distance(send_ranks.begin(), ptr);
 
          AssertDimension(recv_buffer_unpacked.size(),
                          send_ptrs[j + 1] - send_ptrs[j]);
 
          for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
               ++i, ++c)
            buffer_eval[send_permutation_inv[i]] = recv_buffer_unpacked[c];
        }
 
      const int ierr = MPI_Waitall(send_requests.size(),
                                   send_requests.data(),
                                   MPI_STATUSES_IGNORE);
      AssertThrowMPI(ierr);
 
      // evaluate function at points
      evaluation_function(buffer_eval, cell_data);
#endif
    }
 
  } // end of namespace MPI
} // end of namespace Utilities
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif