Reference documentation for deal.II version GIT d6cf33b98c 2023-09-22 19:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mpi_remote_point_evaluation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_mpi_mpi_remote_point_evaluation_h
17 #define dealii_mpi_mpi_remote_point_evaluation_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/mpi.h>
22 #include <deal.II/base/mpi_tags.h>
23 
25 
27 
28 // Forward declarations
29 namespace GridTools
30 {
31  namespace internal
32  {
33  template <int dim, int spacedim>
35  }
36 } // namespace GridTools
37 
38 namespace Utilities
39 {
40  namespace MPI
41  {
52  template <int dim, int spacedim = dim>
54  {
55  public:
72  const double tolerance = 1e-6,
73  const bool enforce_unique_mapping = false,
74  const unsigned int rtree_level = 0,
75  const std::function<std::vector<bool>()> &marked_vertices = {});
76 
81 
93  void
94  reinit(const std::vector<Point<spacedim>> &points,
97 
107  void
109  DistributedComputePointLocationsInternal<dim, spacedim> &data,
112 
116  struct CellData
117  {
121  std::vector<std::pair<int, int>> cells;
122 
127  std::vector<unsigned int> reference_point_ptrs;
128 
132  std::vector<Point<dim>> reference_point_values;
133  };
134 
139  const CellData &
140  get_cell_data() const;
141 
156  template <typename T>
157  void
159  std::vector<T> &output,
160  std::vector<T> &buffer,
161  const std::function<void(const ArrayView<T> &, const CellData &)>
162  &evaluation_function) const;
163 
168  template <typename T>
169  std::vector<T>
171  const std::function<void(const ArrayView<T> &, const CellData &)>
172  &evaluation_function) const;
173 
182  template <typename T>
183  void
185  const std::vector<T> &input,
186  std::vector<T> &buffer,
187  const std::function<void(const ArrayView<const T> &, const CellData &)>
188  &evaluation_function) const;
189 
194  template <typename T>
195  void
197  const std::vector<T> &input,
198  const std::function<void(const ArrayView<const T> &, const CellData &)>
199  &evaluation_function) const;
200 
205  const std::vector<unsigned int> &
206  get_point_ptrs() const;
207 
214  bool
215  is_map_unique() const;
216 
220  bool
221  all_points_found() const;
222 
226  bool
227  point_found(const unsigned int i) const;
228 
233  get_triangulation() const;
234 
238  const Mapping<dim, spacedim> &
239  get_mapping() const;
240 
246  bool
247  is_ready() const;
248 
249  private:
254  const double tolerance;
255 
261 
265  const unsigned int rtree_level;
266 
271  const std::function<std::vector<bool>()> marked_vertices;
272 
276  boost::signals2::connection tria_signal;
277 
284 
289 
294 
299 
304 
310  std::vector<unsigned int> point_ptrs;
311 
315  std::vector<unsigned int> recv_permutation;
316 
321  std::vector<unsigned int> recv_ptrs;
322 
326  std::vector<unsigned int> recv_ranks;
327 
333 
337  std::vector<unsigned int> send_permutation;
338 
342  std::vector<unsigned int> send_ranks;
343 
348  std::vector<unsigned int> send_ptrs;
349  };
350 
351 
352  template <int dim, int spacedim>
353  template <typename T>
354  void
356  std::vector<T> &output,
357  std::vector<T> &buffer,
358  const std::function<void(const ArrayView<T> &, const CellData &)>
359  &evaluation_function) const
360  {
361 #ifndef DEAL_II_WITH_MPI
362  Assert(false, ExcNeedsMPI());
363  (void)output;
364  (void)buffer;
365  (void)evaluation_function;
366 #else
367  static CollectiveMutex mutex;
369 
370  const unsigned int my_rank =
372 
373  // allocate memory for output and buffer
374  output.resize(point_ptrs.back());
375  buffer.resize(std::max(send_permutation.size() * 2,
376  point_ptrs.back() + send_permutation.size()));
377 
378  // ... for evaluation
379  ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
380 
381  // ... for communication
382  ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
383  send_permutation.size());
384 
385  // evaluate functions at points
386  evaluation_function(buffer_eval, cell_data);
387 
388  // sort for communication
389  unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
390  unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
391 
392  const auto my_rank_local_recv_ptr =
393  std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
394 
395  if (my_rank_local_recv_ptr != recv_ranks.end())
396  {
397  my_rank_local_recv =
398  std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
399  my_rank_local_send = std::distance(send_ranks.begin(),
400  std::find(send_ranks.begin(),
401  send_ranks.end(),
402  my_rank));
403  }
404 
405  for (unsigned int i = 0; i < send_permutation.size(); ++i)
406  {
407  const unsigned int send_index = send_permutation[i];
408 
409  // local data -> can be copied to output directly
410  if (my_rank_local_send != numbers::invalid_unsigned_int &&
411  (send_ptrs[my_rank_local_send] <= send_index &&
412  send_index < send_ptrs[my_rank_local_send + 1]))
413  output[recv_permutation[send_index - send_ptrs[my_rank_local_send] +
414  recv_ptrs[my_rank_local_recv]]] =
415  buffer_eval[i];
416  else // data to be sent
417  buffer_comm[send_index] = buffer_eval[i];
418  }
419 
420  // send data
421  std::vector<std::vector<char>> send_buffer;
422  send_buffer.reserve(send_ranks.size());
423 
424  std::vector<MPI_Request> send_requests;
425  send_requests.reserve(send_ranks.size());
426 
427  for (unsigned int i = 0; i < send_ranks.size(); ++i)
428  {
429  if (send_ranks[i] == my_rank)
430  continue;
431 
432  send_requests.emplace_back(MPI_Request());
433 
434  send_buffer.emplace_back(Utilities::pack(
435  std::vector<T>(buffer_comm.begin() + send_ptrs[i],
436  buffer_comm.begin() + send_ptrs[i + 1]),
437  false));
438 
439  const int ierr = MPI_Isend(send_buffer.back().data(),
440  send_buffer.back().size(),
441  MPI_CHAR,
442  send_ranks[i],
445  &send_requests.back());
446  AssertThrowMPI(ierr);
447  }
448 
449  // receive data
450  std::vector<char> buffer_char;
451 
452  for (unsigned int i = 0; i < recv_ranks.size(); ++i)
453  {
454  if (recv_ranks[i] == my_rank)
455  continue;
456 
457  // receive remote data
458  MPI_Status status;
459  int ierr = MPI_Probe(MPI_ANY_SOURCE,
462  &status);
463  AssertThrowMPI(ierr);
464 
465  int message_length;
466  ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
467  AssertThrowMPI(ierr);
468 
469  buffer_char.resize(message_length);
470 
471  ierr = MPI_Recv(buffer_char.data(),
472  buffer_char.size(),
473  MPI_CHAR,
474  status.MPI_SOURCE,
477  MPI_STATUS_IGNORE);
478  AssertThrowMPI(ierr);
479 
480  // unpack data
481  const auto buffer =
482  Utilities::unpack<std::vector<T>>(buffer_char, false);
483 
484  // write data into output vector
485  const auto ptr =
486  std::find(recv_ranks.begin(), recv_ranks.end(), status.MPI_SOURCE);
487 
488  Assert(ptr != recv_ranks.end(), ExcNotImplemented());
489 
490  const unsigned int j = std::distance(recv_ranks.begin(), ptr);
491 
492  AssertDimension(buffer.size(), recv_ptrs[j + 1] - recv_ptrs[j]);
493 
494  for (unsigned int i = recv_ptrs[j], c = 0; i < recv_ptrs[j + 1];
495  ++i, ++c)
496  output[recv_permutation[i]] = buffer[c];
497  }
498 
499  // make sure all messages have been sent
500  if (!send_requests.empty())
501  {
502  const int ierr = MPI_Waitall(send_requests.size(),
503  send_requests.data(),
504  MPI_STATUSES_IGNORE);
505  AssertThrowMPI(ierr);
506  }
507 #endif
508  }
509 
510 
511  template <int dim, int spacedim>
512  template <typename T>
513  std::vector<T>
515  const std::function<void(const ArrayView<T> &, const CellData &)>
516  &evaluation_function) const
517  {
518  std::vector<T> output;
519  std::vector<T> buffer;
520 
521  this->evaluate_and_process(output, buffer, evaluation_function);
522 
523  return output;
524  }
525 
526 
527 
528  template <int dim, int spacedim>
529  template <typename T>
530  void
532  const std::vector<T> &input,
533  std::vector<T> &buffer,
534  const std::function<void(const ArrayView<const T> &, const CellData &)>
535  &evaluation_function) const
536  {
537 #ifndef DEAL_II_WITH_MPI
538  Assert(false, ExcNeedsMPI());
539  (void)input;
540  (void)buffer;
541  (void)evaluation_function;
542 #else
543  static CollectiveMutex mutex;
545 
546  const unsigned int my_rank =
548 
549  // invert permutation matrices (TODO: precompute)
550  std::vector<unsigned int> recv_permutation_inv(recv_permutation.size());
551  for (unsigned int c = 0; c < recv_permutation.size(); ++c)
552  recv_permutation_inv[recv_permutation[c]] = c;
553 
554  std::vector<unsigned int> send_permutation_inv(send_permutation.size());
555  for (unsigned int c = 0; c < send_permutation.size(); ++c)
556  send_permutation_inv[send_permutation[c]] = c;
557 
558  // allocate memory for buffer
559  const auto &point_ptrs = this->get_point_ptrs();
560  AssertDimension(input.size(), point_ptrs.size() - 1);
561  buffer.resize(std::max(send_permutation.size() * 2,
562  point_ptrs.back() + send_permutation.size()));
563 
564  // ... for evaluation
565  ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
566 
567  // ... for communication
568  ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
569  point_ptrs.back());
570 
571  // sort for communication (and duplicate data if necessary)
572  unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
573  unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
574 
575  const auto my_rank_local_recv_ptr =
576  std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
577 
578  if (my_rank_local_recv_ptr != recv_ranks.end())
579  {
580  my_rank_local_recv =
581  std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
582  my_rank_local_send = std::distance(send_ranks.begin(),
583  std::find(send_ranks.begin(),
584  send_ranks.end(),
585  my_rank));
586  }
587 
588  for (unsigned int i = 0, c = 0; i < point_ptrs.size() - 1; ++i)
589  {
590  const auto n_entries = point_ptrs[i + 1] - point_ptrs[i];
591 
592  for (unsigned int j = 0; j < n_entries; ++j, ++c)
593  {
594  const unsigned int recv_index = recv_permutation_inv[c];
595 
596  // local data -> can be copied to final buffer directly
597  if (my_rank_local_recv != numbers::invalid_unsigned_int &&
598  (recv_ptrs[my_rank_local_recv] <= recv_index &&
599  recv_index < recv_ptrs[my_rank_local_recv + 1]))
600  buffer_eval[send_permutation_inv
601  [recv_index - recv_ptrs[my_rank_local_recv] +
602  send_ptrs[my_rank_local_send]]] = input[i];
603  else // data to be sent
604  buffer_comm[recv_index] = input[i];
605  }
606  }
607 
608  // send data
609  std::vector<std::vector<char>> send_buffer;
610  send_buffer.reserve(recv_ranks.size());
611 
612  std::vector<MPI_Request> send_requests;
613  send_requests.reserve(recv_ranks.size());
614 
615  for (unsigned int i = 0; i < recv_ranks.size(); ++i)
616  {
617  if (recv_ranks[i] == my_rank)
618  continue;
619 
620  send_requests.push_back(MPI_Request());
621 
622  send_buffer.emplace_back(Utilities::pack(
623  std::vector<T>(buffer_comm.begin() + recv_ptrs[i],
624  buffer_comm.begin() + recv_ptrs[i + 1]),
625  false));
626 
627  const int ierr = MPI_Isend(send_buffer.back().data(),
628  send_buffer.back().size(),
629  MPI_CHAR,
630  recv_ranks[i],
633  &send_requests.back());
634  AssertThrowMPI(ierr);
635  }
636 
637  // receive data
638  std::vector<char> recv_buffer;
639 
640  for (unsigned int i = 0; i < send_ranks.size(); ++i)
641  {
642  if (send_ranks[i] == my_rank)
643  continue;
644 
645  // receive remote data
646  MPI_Status status;
647  int ierr = MPI_Probe(MPI_ANY_SOURCE,
650  &status);
651  AssertThrowMPI(ierr);
652 
653  int message_length;
654  ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
655  AssertThrowMPI(ierr);
656 
657  recv_buffer.resize(message_length);
658 
659  ierr = MPI_Recv(recv_buffer.data(),
660  recv_buffer.size(),
661  MPI_CHAR,
662  status.MPI_SOURCE,
665  MPI_STATUS_IGNORE);
666  AssertThrowMPI(ierr);
667 
668  // unpack data
669  const auto recv_buffer_unpacked =
670  Utilities::unpack<std::vector<T>>(recv_buffer, false);
671 
672  // write data into buffer vector
673  const auto ptr =
674  std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
675 
676  Assert(ptr != send_ranks.end(), ExcNotImplemented());
677 
678  const unsigned int j = std::distance(send_ranks.begin(), ptr);
679 
680  AssertDimension(recv_buffer_unpacked.size(),
681  send_ptrs[j + 1] - send_ptrs[j]);
682 
683  for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
684  ++i, ++c)
685  buffer_eval[send_permutation_inv[i]] = recv_buffer_unpacked[c];
686  }
687 
688  if (!send_requests.empty())
689  {
690  const int ierr = MPI_Waitall(send_requests.size(),
691  send_requests.data(),
692  MPI_STATUSES_IGNORE);
693  AssertThrowMPI(ierr);
694  }
695 
696  // evaluate function at points
697  evaluation_function(buffer_eval, cell_data);
698 #endif
699  }
700 
701 
702 
703  template <int dim, int spacedim>
704  template <typename T>
705  void
707  const std::vector<T> &input,
708  const std::function<void(const ArrayView<const T> &, const CellData &)>
709  &evaluation_function) const
710  {
711  std::vector<T> buffer;
712  this->process_and_evaluate(input, buffer, evaluation_function);
713  }
714 
715  } // end of namespace MPI
716 } // end of namespace Utilities
717 
718 
720 
721 #endif
iterator begin() const
Definition: array_view.h:591
virtual MPI_Comm get_communicator() const
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={})
const std::function< std::vector< bool >)> marked_vertices
void evaluate_and_process(std::vector< T > &output, std::vector< T > &buffer, const std::function< void(const ArrayView< T > &, const CellData &)> &evaluation_function) const
void process_and_evaluate(const std::vector< T > &input, const std::function< void(const ArrayView< const T > &, const CellData &)> &evaluation_function) const
const Triangulation< dim, spacedim > & get_triangulation() const
SmartPointer< const Mapping< dim, spacedim > > mapping
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
const Mapping< dim, spacedim > & get_mapping() const
SmartPointer< const Triangulation< dim, spacedim > > tria
std::vector< T > evaluate_and_process(const std::function< void(const ArrayView< T > &, const CellData &)> &evaluation_function) const
bool point_found(const unsigned int i) const
void reinit(const std::vector< Point< spacedim >> &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:164
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1343
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const ::Triangulation< dim, spacedim > & tria