Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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>
23
25
27
28// Forward declarations
29namespace GridTools
30{
31 namespace internal
32 {
33 template <int dim, int spacedim>
35 }
36} // namespace GridTools
37
38namespace 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
172 template <typename T>
173 void
175 const std::vector<T> &input,
176 std::vector<T> & buffer,
177 const std::function<void(const ArrayView<const T> &, const CellData &)>
178 &evaluation_function) const;
179
184 const std::vector<unsigned int> &
185 get_point_ptrs() const;
186
193 bool
194 is_map_unique() const;
195
199 bool
200 all_points_found() const;
201
205 bool
206 point_found(const unsigned int i) const;
207
212 get_triangulation() const;
213
218 get_mapping() const;
219
225 bool
226 is_ready() const;
227
228 private:
233 const double tolerance;
234
240
244 const unsigned int rtree_level;
245
250 const std::function<std::vector<bool>()> marked_vertices;
251
255 boost::signals2::connection tria_signal;
256
263
268
273
278
283
289 std::vector<unsigned int> point_ptrs;
290
294 std::vector<unsigned int> recv_permutation;
295
300 std::vector<unsigned int> recv_ptrs;
301
305 std::vector<unsigned int> recv_ranks;
306
312
316 std::vector<unsigned int> send_permutation;
317
321 std::vector<unsigned int> send_ranks;
322
327 std::vector<unsigned int> send_ptrs;
328 };
329
330
331 template <int dim, int spacedim>
332 template <typename T>
333 void
335 std::vector<T> &output,
336 std::vector<T> &buffer,
337 const std::function<void(const ArrayView<T> &, const CellData &)>
338 &evaluation_function) const
339 {
340#ifndef DEAL_II_WITH_MPI
341 Assert(false, ExcNeedsMPI());
342 (void)output;
343 (void)buffer;
344 (void)evaluation_function;
345#else
346 static CollectiveMutex mutex;
348
349 const unsigned int my_rank =
351
352 // allocate memory for output and buffer
353 output.resize(point_ptrs.back());
354 buffer.resize(std::max(send_permutation.size() * 2,
355 point_ptrs.back() + send_permutation.size()));
356
357 // ... for evaluation
358 ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
359
360 // ... for communication
361 ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
362 send_permutation.size());
363
364 // evaluate functions at points
365 evaluation_function(buffer_eval, cell_data);
366
367 // sort for communication
368 unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
369 unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
370
371 const auto my_rank_local_recv_ptr =
372 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
373
374 if (my_rank_local_recv_ptr != recv_ranks.end())
375 {
376 my_rank_local_recv =
377 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
378 my_rank_local_send = std::distance(send_ranks.begin(),
379 std::find(send_ranks.begin(),
380 send_ranks.end(),
381 my_rank));
382 }
383
384 for (unsigned int i = 0; i < send_permutation.size(); ++i)
385 {
386 const unsigned int send_index = send_permutation[i];
387
388 // local data -> can be copied to output directly
389 if (my_rank_local_send != numbers::invalid_unsigned_int &&
390 (send_ptrs[my_rank_local_send] <= send_index &&
391 send_index < send_ptrs[my_rank_local_send + 1]))
392 output[recv_permutation[send_index - send_ptrs[my_rank_local_send] +
393 recv_ptrs[my_rank_local_recv]]] =
394 buffer_eval[i];
395 else // data to be sent
396 buffer_comm[send_index] = buffer_eval[i];
397 }
398
399 // send data
400 std::vector<std::vector<char>> send_buffer;
401 send_buffer.reserve(send_ranks.size());
402
403 std::vector<MPI_Request> send_requests;
404 send_requests.reserve(send_ranks.size());
405
406 for (unsigned int i = 0; i < send_ranks.size(); ++i)
407 {
408 if (send_ranks[i] == my_rank)
409 continue;
410
411 send_requests.emplace_back(MPI_Request());
412
413 send_buffer.emplace_back(Utilities::pack(
414 std::vector<T>(buffer_comm.begin() + send_ptrs[i],
415 buffer_comm.begin() + send_ptrs[i + 1]),
416 false));
417
418 const int ierr = MPI_Isend(send_buffer.back().data(),
419 send_buffer.back().size(),
420 MPI_CHAR,
421 send_ranks[i],
424 &send_requests.back());
425 AssertThrowMPI(ierr);
426 }
427
428 // receive data
429 std::vector<char> buffer_char;
430
431 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
432 {
433 if (recv_ranks[i] == my_rank)
434 continue;
435
436 // receive remote data
437 MPI_Status status;
438 int ierr = MPI_Probe(MPI_ANY_SOURCE,
441 &status);
442 AssertThrowMPI(ierr);
443
444 int message_length;
445 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
446 AssertThrowMPI(ierr);
447
448 buffer_char.resize(message_length);
449
450 ierr = MPI_Recv(buffer_char.data(),
451 buffer_char.size(),
452 MPI_CHAR,
453 status.MPI_SOURCE,
456 MPI_STATUS_IGNORE);
457 AssertThrowMPI(ierr);
458
459 // unpack data
460 const auto buffer =
461 Utilities::unpack<std::vector<T>>(buffer_char, false);
462
463 // write data into output vector
464 const auto ptr =
465 std::find(recv_ranks.begin(), recv_ranks.end(), status.MPI_SOURCE);
466
467 Assert(ptr != recv_ranks.end(), ExcNotImplemented());
468
469 const unsigned int j = std::distance(recv_ranks.begin(), ptr);
470
471 AssertDimension(buffer.size(), recv_ptrs[j + 1] - recv_ptrs[j]);
472
473 for (unsigned int i = recv_ptrs[j], c = 0; i < recv_ptrs[j + 1];
474 ++i, ++c)
475 output[recv_permutation[i]] = buffer[c];
476 }
477
478 // make sure all messages have been sent
479 const int ierr = MPI_Waitall(send_requests.size(),
480 send_requests.data(),
481 MPI_STATUSES_IGNORE);
482 AssertThrowMPI(ierr);
483#endif
484 }
485
486
487 template <int dim, int spacedim>
488 template <typename T>
489 void
491 const std::vector<T> &input,
492 std::vector<T> & buffer,
493 const std::function<void(const ArrayView<const T> &, const CellData &)>
494 &evaluation_function) const
495 {
496#ifndef DEAL_II_WITH_MPI
497 Assert(false, ExcNeedsMPI());
498 (void)input;
499 (void)buffer;
500 (void)evaluation_function;
501#else
502 static CollectiveMutex mutex;
504
505 const unsigned int my_rank =
507
508 // invert permutation matrices (TODO: precompute)
509 std::vector<unsigned int> recv_permutation_inv(recv_permutation.size());
510 for (unsigned int c = 0; c < recv_permutation.size(); ++c)
511 recv_permutation_inv[recv_permutation[c]] = c;
512
513 std::vector<unsigned int> send_permutation_inv(send_permutation.size());
514 for (unsigned int c = 0; c < send_permutation.size(); ++c)
515 send_permutation_inv[send_permutation[c]] = c;
516
517 // allocate memory for buffer
518 const auto &point_ptrs = this->get_point_ptrs();
519 AssertDimension(input.size(), point_ptrs.size() - 1);
520 buffer.resize(std::max(send_permutation.size() * 2,
521 point_ptrs.back() + send_permutation.size()));
522
523 // ... for evaluation
524 ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
525
526 // ... for communication
527 ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
528 point_ptrs.back());
529
530 // sort for communication (and duplicate data if necessary)
531 unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
532 unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
533
534 const auto my_rank_local_recv_ptr =
535 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
536
537 if (my_rank_local_recv_ptr != recv_ranks.end())
538 {
539 my_rank_local_recv =
540 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
541 my_rank_local_send = std::distance(send_ranks.begin(),
542 std::find(send_ranks.begin(),
543 send_ranks.end(),
544 my_rank));
545 }
546
547 for (unsigned int i = 0, c = 0; i < point_ptrs.size() - 1; ++i)
548 {
549 const auto n_entries = point_ptrs[i + 1] - point_ptrs[i];
550
551 for (unsigned int j = 0; j < n_entries; ++j, ++c)
552 {
553 const unsigned int recv_index = recv_permutation_inv[c];
554
555 // local data -> can be copied to final buffer directly
556 if (my_rank_local_recv != numbers::invalid_unsigned_int &&
557 (recv_ptrs[my_rank_local_recv] <= recv_index &&
558 recv_index < recv_ptrs[my_rank_local_recv + 1]))
559 buffer_eval[send_permutation_inv
560 [recv_index - recv_ptrs[my_rank_local_recv] +
561 send_ptrs[my_rank_local_send]]] = input[i];
562 else // data to be sent
563 buffer_comm[recv_index] = input[i];
564 }
565 }
566
567 // send data
568 std::vector<std::vector<char>> send_buffer;
569 send_buffer.reserve(recv_ranks.size());
570
571 std::vector<MPI_Request> send_requests;
572 send_requests.reserve(recv_ranks.size());
573
574 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
575 {
576 if (recv_ranks[i] == my_rank)
577 continue;
578
579 send_requests.push_back(MPI_Request());
580
581 send_buffer.emplace_back(Utilities::pack(
582 std::vector<T>(buffer_comm.begin() + recv_ptrs[i],
583 buffer_comm.begin() + recv_ptrs[i + 1]),
584 false));
585
586 const int ierr = MPI_Isend(send_buffer.back().data(),
587 send_buffer.back().size(),
588 MPI_CHAR,
589 recv_ranks[i],
592 &send_requests.back());
593 AssertThrowMPI(ierr);
594 }
595
596 // receive data
597 std::vector<char> recv_buffer;
598
599 for (unsigned int i = 0; i < send_ranks.size(); ++i)
600 {
601 if (send_ranks[i] == my_rank)
602 continue;
603
604 // receive remote data
605 MPI_Status status;
606 int ierr = MPI_Probe(MPI_ANY_SOURCE,
609 &status);
610 AssertThrowMPI(ierr);
611
612 int message_length;
613 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
614 AssertThrowMPI(ierr);
615
616 recv_buffer.resize(message_length);
617
618 ierr = MPI_Recv(recv_buffer.data(),
619 recv_buffer.size(),
620 MPI_CHAR,
621 status.MPI_SOURCE,
624 MPI_STATUS_IGNORE);
625 AssertThrowMPI(ierr);
626
627 // unpack data
628 const auto recv_buffer_unpacked =
629 Utilities::unpack<std::vector<T>>(recv_buffer, false);
630
631 // write data into buffer vector
632 const auto ptr =
633 std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
634
635 Assert(ptr != send_ranks.end(), ExcNotImplemented());
636
637 const unsigned int j = std::distance(send_ranks.begin(), ptr);
638
639 AssertDimension(recv_buffer_unpacked.size(),
640 send_ptrs[j + 1] - send_ptrs[j]);
641
642 for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
643 ++i, ++c)
644 buffer_eval[send_permutation_inv[i]] = recv_buffer_unpacked[c];
645 }
646
647 const int ierr = MPI_Waitall(send_requests.size(),
648 send_requests.data(),
649 MPI_STATUSES_IGNORE);
650 AssertThrowMPI(ierr);
651
652 // evaluate function at points
653 evaluation_function(buffer_eval, cell_data);
654#endif
655 }
656
657 } // end of namespace MPI
658} // end of namespace Utilities
659
660
662
663#endif
iterator begin() const
Definition array_view.h:594
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
virtual MPI_Comm get_communicator() const
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
const std::function< std::vector< bool >()> marked_vertices
void reinit(const std::vector< Point< spacedim > > &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
SmartPointer< const Triangulation< dim, spacedim > > tria
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
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:1351
static const unsigned int invalid_unsigned_int
Definition types.h:213
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const ::Triangulation< dim, spacedim > & tria