Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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 - 2022 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 
29 namespace Utilities
30 {
31  namespace MPI
32  {
43  template <int dim, int spacedim = dim>
45  {
46  public:
63  const double tolerance = 1e-6,
64  const bool enforce_unique_mapping = false,
65  const unsigned int rtree_level = 0,
66  const std::function<std::vector<bool>()> &marked_vertices = {});
67 
72 
84  void
85  reinit(const std::vector<Point<spacedim>> &points,
88 
92  struct CellData
93  {
97  std::vector<std::pair<int, int>> cells;
98 
103  std::vector<unsigned int> reference_point_ptrs;
104 
108  std::vector<Point<dim>> reference_point_values;
109  };
110 
125  template <typename T>
126  void
128  std::vector<T> &output,
129  std::vector<T> &buffer,
130  const std::function<void(const ArrayView<T> &, const CellData &)>
131  &evaluation_function) const;
132 
141  template <typename T>
142  void
144  const std::vector<T> &input,
145  std::vector<T> & buffer,
146  const std::function<void(const ArrayView<const T> &, const CellData &)>
147  &evaluation_function) const;
148 
153  const std::vector<unsigned int> &
154  get_point_ptrs() const;
155 
162  bool
163  is_map_unique() const;
164 
168  bool
169  all_points_found() const;
170 
174  bool
175  point_found(const unsigned int i) const;
176 
181  get_triangulation() const;
182 
186  const Mapping<dim, spacedim> &
187  get_mapping() const;
188 
194  bool
195  is_ready() const;
196 
197  private:
202  const double tolerance;
203 
209 
213  const unsigned int rtree_level;
214 
219  const std::function<std::vector<bool>()> marked_vertices;
220 
224  boost::signals2::connection tria_signal;
225 
232 
237 
242 
247 
252 
258  std::vector<unsigned int> point_ptrs;
259 
263  std::vector<unsigned int> recv_permutation;
264 
269  std::vector<unsigned int> recv_ptrs;
270 
274  std::vector<unsigned int> recv_ranks;
275 
281 
285  std::vector<unsigned int> send_permutation;
286 
290  std::vector<unsigned int> send_ranks;
291 
296  std::vector<unsigned int> send_ptrs;
297  };
298 
299 
300  template <int dim, int spacedim>
301  template <typename T>
302  void
304  std::vector<T> &output,
305  std::vector<T> &buffer,
306  const std::function<void(const ArrayView<T> &, const CellData &)>
307  &evaluation_function) const
308  {
309 #ifndef DEAL_II_WITH_MPI
310  Assert(false, ExcNeedsMPI());
311  (void)output;
312  (void)buffer;
313  (void)evaluation_function;
314 #else
315  static CollectiveMutex mutex;
317 
318  output.resize(point_ptrs.back());
319  buffer.resize(send_permutation.size() * 2);
320  ArrayView<T> buffer_1(buffer.data(), buffer.size() / 2);
321  ArrayView<T> buffer_2(buffer.data() + buffer.size() / 2,
322  buffer.size() / 2);
323 
324  // evaluate functions at points
325  evaluation_function(buffer_1, cell_data);
326 
327  // sort for communication
328  for (unsigned int i = 0; i < send_permutation.size(); ++i)
329  buffer_2[send_permutation[i]] = buffer_1[i];
330 
331  // process remote quadrature points and send them away
332  std::map<unsigned int, std::vector<char>> temp_map;
333 
334  std::vector<MPI_Request> requests;
335  requests.reserve(send_ranks.size());
336 
337  const unsigned int my_rank =
339 
340  std::map<unsigned int, std::vector<T>> temp_recv_map;
341 
342  for (unsigned int i = 0; i < send_ranks.size(); ++i)
343  {
344  if (send_ranks[i] == my_rank)
345  {
346  // process locally-owned values
347  temp_recv_map[my_rank] =
348  std::vector<T>(buffer_2.begin() + send_ptrs[i],
349  buffer_2.begin() + send_ptrs[i + 1]);
350  continue;
351  }
352 
353  temp_map[send_ranks[i]] =
354  Utilities::pack(std::vector<T>(buffer_2.begin() + send_ptrs[i],
355  buffer_2.begin() + send_ptrs[i + 1]),
356  false);
357 
358  auto &buffer = temp_map[send_ranks[i]];
359 
360  requests.push_back(MPI_Request());
361 
362  const int ierr = MPI_Isend(buffer.data(),
363  buffer.size(),
364  MPI_CHAR,
365  send_ranks[i],
368  &requests.back());
369  AssertThrowMPI(ierr);
370  }
371 
372  for (const auto recv_rank : recv_ranks)
373  {
374  if (recv_rank == my_rank)
375  continue;
376 
377  MPI_Status status;
378  int ierr = MPI_Probe(MPI_ANY_SOURCE,
381  &status);
382  AssertThrowMPI(ierr);
383 
384  int message_length;
385  ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
386  AssertThrowMPI(ierr);
387 
388  std::vector<char> buffer(message_length);
389 
390  ierr = MPI_Recv(buffer.data(),
391  buffer.size(),
392  MPI_CHAR,
393  status.MPI_SOURCE,
396  MPI_STATUS_IGNORE);
397  AssertThrowMPI(ierr);
398 
399  temp_recv_map[status.MPI_SOURCE] =
400  Utilities::unpack<std::vector<T>>(buffer, false);
401  }
402 
403  // make sure all messages have been sent
404  const int ierr =
405  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
406  AssertThrowMPI(ierr);
407 
408  // copy received data into output vector
409  auto it = recv_permutation.begin();
410  for (const auto &j : temp_recv_map)
411  for (const auto &i : j.second)
412  {
413  output[*it] = i;
414  it++;
415  }
416 #endif
417  }
418 
419 
420  template <int dim, int spacedim>
421  template <typename T>
422  void
424  const std::vector<T> &input,
425  std::vector<T> & buffer,
426  const std::function<void(const ArrayView<const T> &, const CellData &)>
427  &evaluation_function) const
428  {
429 #ifndef DEAL_II_WITH_MPI
430  Assert(false, ExcNeedsMPI());
431  (void)input;
432  (void)buffer;
433  (void)evaluation_function;
434 #else
435  static CollectiveMutex mutex;
437 
438  const auto &ptr = this->get_point_ptrs();
439 
440  std::map<unsigned int, std::vector<T>> temp_recv_map;
441 
442  for (unsigned int i = 0; i < recv_ranks.size(); ++i)
443  temp_recv_map[recv_ranks[i]].resize(recv_ptrs[i + 1] - recv_ptrs[i]);
444 
445  const unsigned int my_rank =
447 
448 # ifdef DEBUG
449  {
450  unsigned int i = 0;
451 
452  for (auto &j : temp_recv_map)
453  i += j.second.size();
454 
455  AssertDimension(recv_permutation.size(), i);
456  }
457 # endif
458 
459  {
460  // duplicate data to be able to sort it more easily in the next step
461  std::vector<T> buffer_(ptr.back());
462  for (unsigned int i = 0, c = 0; i < ptr.size() - 1; ++i)
463  {
464  const auto n_entries = ptr[i + 1] - ptr[i];
465 
466  for (unsigned int j = 0; j < n_entries; ++j, ++c)
467  buffer_[c] = input[i];
468  }
469 
470  // sort data according to the ranks
471  auto it = recv_permutation.begin();
472  for (auto &j : temp_recv_map)
473  for (auto &i : j.second)
474  {
475  i = buffer_[*it];
476  it++;
477  }
478  }
479 
480  // buffer.resize(point_ptrs.back());
481  buffer.resize(send_permutation.size() * 2);
482  ArrayView<T> buffer_1(buffer.data(), buffer.size() / 2);
483  ArrayView<T> buffer_2(buffer.data() + buffer.size() / 2,
484  buffer.size() / 2);
485 
486  // process remote quadrature points and send them away
487  std::map<unsigned int, std::vector<char>> temp_map;
488 
489  std::vector<MPI_Request> requests;
490  requests.reserve(recv_ranks.size());
491 
492  for (const auto recv_rank : recv_ranks)
493  {
494  if (recv_rank == my_rank)
495  continue;
496 
497  temp_map[recv_rank] =
498  Utilities::pack(temp_recv_map[recv_rank], false);
499 
500  auto &buffer_send = temp_map[recv_rank];
501 
502  requests.push_back(MPI_Request());
503 
504  const int ierr = MPI_Isend(buffer_send.data(),
505  buffer_send.size(),
506  MPI_CHAR,
507  recv_rank,
510  &requests.back());
511  AssertThrowMPI(ierr);
512  }
513 
514  for (unsigned int i = 0; i < send_ranks.size(); ++i)
515  {
516  if (send_ranks[i] == my_rank)
517  {
518  const auto &buffer_send = temp_recv_map[send_ranks[i]];
519  // process locally-owned values
520  const unsigned int j = std::distance(send_ranks.begin(),
521  std::find(send_ranks.begin(),
522  send_ranks.end(),
523  my_rank));
524 
525  AssertDimension(buffer_send.size(),
526  send_ptrs[j + 1] - send_ptrs[j]);
527 
528  for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
529  ++i, ++c)
530  buffer_1[i] = buffer_send[c];
531 
532  continue;
533  }
534 
535  MPI_Status status;
536  int ierr = MPI_Probe(MPI_ANY_SOURCE,
539  &status);
540  AssertThrowMPI(ierr);
541 
542  int message_length;
543  ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
544  AssertThrowMPI(ierr);
545 
546  std::vector<char> recv_buffer(message_length);
547 
548  ierr = MPI_Recv(recv_buffer.data(),
549  recv_buffer.size(),
550  MPI_CHAR,
551  status.MPI_SOURCE,
554  MPI_STATUS_IGNORE);
555  AssertThrowMPI(ierr);
556 
557  const auto recv_buffer_unpacked =
558  Utilities::unpack<std::vector<T>>(recv_buffer, false);
559 
560  auto ptr =
561  std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
562 
563  Assert(ptr != send_ranks.end(), ExcNotImplemented());
564 
565  const unsigned int j = std::distance(send_ranks.begin(), ptr);
566 
567  AssertDimension(recv_buffer_unpacked.size(),
568  send_ptrs[j + 1] - send_ptrs[j]);
569 
570  for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
571  ++i, ++c)
572  {
573  AssertIndexRange(i, buffer_1.size());
574  AssertIndexRange(c, recv_buffer_unpacked.size());
575  buffer_1[i] = recv_buffer_unpacked[c];
576  }
577  }
578 
579  const int ierr =
580  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
581  AssertThrowMPI(ierr);
582 
583  // sort for easy access during function call
584  for (unsigned int i = 0; i < send_permutation.size(); ++i)
585  buffer_2[i] = buffer_1[send_permutation[i]];
586 
587  // evaluate function at points
588  evaluation_function(buffer_2, cell_data);
589 #endif
590  }
591 
592  } // end of namespace MPI
593 } // end of namespace Utilities
594 
595 
597 
598 #endif
iterator begin() const
Definition: array_view.h:594
std::size_t size() const
Definition: array_view.h:576
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
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
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1886
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:161
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1334
const ::Triangulation< dim, spacedim > & tria