Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20: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\}}\)
Loading...
Searching...
No Matches
mpi_remote_point_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_mpi_mpi_remote_point_evaluation_h
16#define dealii_mpi_mpi_remote_point_evaluation_h
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/base/mpi.h>
22
24
26
27// Forward declarations
28namespace GridTools
29{
30 namespace internal
31 {
32 template <int dim, int spacedim>
34 }
35} // namespace GridTools
36
37namespace Utilities
38{
39 namespace MPI
40 {
51 template <int dim, int spacedim = dim>
53 {
54 public:
60 {
61 public:
66 const double tolerance = 1e-6,
67 const bool enforce_unique_mapping = false,
68 const unsigned int rtree_level = 0,
69 const std::function<std::vector<bool>()> &marked_vertices = {});
70
79 double tolerance;
80
86
90 unsigned int rtree_level;
91
96 std::function<std::vector<bool>()> marked_vertices;
97 };
98
106
125 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
126 "Use the constructor with AdditionalData struct.")
128 const double tolerance,
129 const bool enforce_unique_mapping = false,
130 const unsigned int rtree_level = 0,
131 const std::function<std::vector<bool>()> &marked_vertices = {});
132
137
149 void
150 reinit(const std::vector<Point<spacedim>> &points,
153
163 void
165 DistributedComputePointLocationsInternal<dim, spacedim> &data,
166 const Triangulation<dim, spacedim> &tria,
167 const Mapping<dim, spacedim> &mapping);
168
174 {
175 public:
182
209 cell_indices() const;
210
215 get_active_cell_iterator(const unsigned int cell) const;
216
221 get_unit_points(const unsigned int cell) const;
222
226 template <typename T>
228 get_data_view(const unsigned int cell,
229 const ArrayView<T> &values) const;
230
234 std::vector<std::pair<int, int>> cells;
235
240 std::vector<unsigned int> reference_point_ptrs;
241
245 std::vector<Point<dim>> reference_point_values;
246
247 private:
253 };
254
259 const CellData &
260 get_cell_data() const;
261
276 template <typename T>
277 void
279 std::vector<T> &output,
280 std::vector<T> &buffer,
281 const std::function<void(const ArrayView<T> &, const CellData &)>
282 &evaluation_function) const;
283
288 template <typename T>
289 std::vector<T>
291 const std::function<void(const ArrayView<T> &, const CellData &)>
292 &evaluation_function) const;
293
302 template <typename T>
303 void
305 const std::vector<T> &input,
306 std::vector<T> &buffer,
307 const std::function<void(const ArrayView<const T> &, const CellData &)>
308 &evaluation_function) const;
309
314 template <typename T>
315 void
317 const std::vector<T> &input,
318 const std::function<void(const ArrayView<const T> &, const CellData &)>
319 &evaluation_function) const;
320
325 const std::vector<unsigned int> &
326 get_point_ptrs() const;
327
334 bool
335 is_map_unique() const;
336
340 bool
341 all_points_found() const;
342
346 bool
347 point_found(const unsigned int i) const;
348
353 get_triangulation() const;
354
359 get_mapping() const;
360
366 bool
367 is_ready() const;
368
369 private:
374
378 boost::signals2::connection tria_signal;
379
386
391
396
401
406
412 std::vector<unsigned int> point_ptrs;
413
417 std::vector<unsigned int> recv_permutation;
418
423 std::vector<unsigned int> recv_ptrs;
424
428 std::vector<unsigned int> recv_ranks;
429
434 std::unique_ptr<CellData> cell_data;
435
439 std::vector<unsigned int> send_permutation;
440
444 std::vector<unsigned int> send_ranks;
445
450 std::vector<unsigned int> send_ptrs;
451 };
452
453
454
455 template <int dim, int spacedim>
456 template <typename T>
459 const unsigned int cell,
460 const ArrayView<T> &values) const
461 {
462 AssertIndexRange(cell, cells.size());
463 return {values.data() + reference_point_ptrs[cell],
465 }
466
467
468
469 template <int dim, int spacedim>
470 template <typename T>
471 void
473 std::vector<T> &output,
474 std::vector<T> &buffer,
475 const std::function<void(const ArrayView<T> &, const CellData &)>
476 &evaluation_function) const
477 {
478#ifndef DEAL_II_WITH_MPI
479 Assert(false, ExcNeedsMPI());
480 (void)output;
481 (void)buffer;
482 (void)evaluation_function;
483#else
484 static CollectiveMutex mutex;
486
487 const unsigned int my_rank =
489
490 // allocate memory for output and buffer
491 output.resize(point_ptrs.back());
492 buffer.resize(std::max(send_permutation.size() * 2,
493 point_ptrs.back() + send_permutation.size()));
494
495 // ... for evaluation
496 ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
497
498 // ... for communication
499 ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
500 send_permutation.size());
501
502 // evaluate functions at points
503 evaluation_function(buffer_eval, *cell_data);
504
505 // sort for communication
506 const auto my_rank_local_recv_ptr =
507 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
508
509 if (my_rank_local_recv_ptr != recv_ranks.end())
510 {
511 const unsigned int my_rank_local_recv =
512 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
513 const unsigned int my_rank_local_send = std::distance(
514 send_ranks.begin(),
515 std::find(send_ranks.begin(), send_ranks.end(), my_rank));
516 const unsigned int start = send_ptrs[my_rank_local_send];
517 const unsigned int end = send_ptrs[my_rank_local_send + 1];
518 const unsigned int *recv_ptr =
519 recv_permutation.data() + recv_ptrs[my_rank_local_recv];
520 for (unsigned int i = 0; i < send_permutation.size(); ++i)
521 {
522 const unsigned int send_index = send_permutation[i];
523
524 // local data -> can be copied to output directly
525 if (start <= send_index && send_index < end)
526 output[recv_ptr[send_index - start]] = buffer_eval[i];
527 else // data to be sent
528 buffer_comm[send_index] = buffer_eval[i];
529 }
530 }
531 else
532 {
533 for (unsigned int i = 0; i < send_permutation.size(); ++i)
534 buffer_comm[send_permutation[i]] = buffer_eval[i];
535 }
536
537 // send data
538 std::vector<std::vector<char>> send_buffer;
539 send_buffer.reserve(send_ranks.size());
540
541 std::vector<MPI_Request> send_requests;
542 send_requests.reserve(send_ranks.size());
543
544 for (unsigned int i = 0; i < send_ranks.size(); ++i)
545 {
546 if (send_ranks[i] == my_rank)
547 continue;
548
549 send_requests.emplace_back(MPI_Request());
550
551 send_buffer.emplace_back(Utilities::pack(
552 std::vector<T>(buffer_comm.begin() + send_ptrs[i],
553 buffer_comm.begin() + send_ptrs[i + 1]),
554 false));
555
556 const int ierr = MPI_Isend(send_buffer.back().data(),
557 send_buffer.back().size(),
558 MPI_CHAR,
559 send_ranks[i],
562 &send_requests.back());
563 AssertThrowMPI(ierr);
564 }
565
566 // receive data
567 std::vector<char> buffer_char;
568
569 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
570 {
571 if (recv_ranks[i] == my_rank)
572 continue;
573
574 // receive remote data
575 MPI_Status status;
576 int ierr = MPI_Probe(MPI_ANY_SOURCE,
579 &status);
580 AssertThrowMPI(ierr);
581
582 int message_length;
583 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
584 AssertThrowMPI(ierr);
585
586 buffer_char.resize(message_length);
587
588 ierr = MPI_Recv(buffer_char.data(),
589 buffer_char.size(),
590 MPI_CHAR,
591 status.MPI_SOURCE,
594 MPI_STATUS_IGNORE);
595 AssertThrowMPI(ierr);
596
597 // unpack data
598 const auto buffer =
599 Utilities::unpack<std::vector<T>>(buffer_char, false);
600
601 // write data into output vector
602 const auto ptr =
603 std::find(recv_ranks.begin(), recv_ranks.end(), status.MPI_SOURCE);
604
605 Assert(ptr != recv_ranks.end(), ExcNotImplemented());
606
607 const unsigned int j = std::distance(recv_ranks.begin(), ptr);
608
609 AssertDimension(buffer.size(), recv_ptrs[j + 1] - recv_ptrs[j]);
610
611 for (unsigned int i = recv_ptrs[j], c = 0; i < recv_ptrs[j + 1];
612 ++i, ++c)
613 output[recv_permutation[i]] = buffer[c];
614 }
615
616 // make sure all messages have been sent
617 if (!send_requests.empty())
618 {
619 const int ierr = MPI_Waitall(send_requests.size(),
620 send_requests.data(),
621 MPI_STATUSES_IGNORE);
622 AssertThrowMPI(ierr);
623 }
624#endif
625 }
626
627
628 template <int dim, int spacedim>
629 template <typename T>
630 std::vector<T>
632 const std::function<void(const ArrayView<T> &, const CellData &)>
633 &evaluation_function) const
634 {
635 std::vector<T> output;
636 std::vector<T> buffer;
637
638 this->evaluate_and_process(output, buffer, evaluation_function);
639
640 return output;
641 }
642
643
644
645 template <int dim, int spacedim>
646 template <typename T>
647 void
649 const std::vector<T> &input,
650 std::vector<T> &buffer,
651 const std::function<void(const ArrayView<const T> &, const CellData &)>
652 &evaluation_function) const
653 {
654#ifndef DEAL_II_WITH_MPI
655 Assert(false, ExcNeedsMPI());
656 (void)input;
657 (void)buffer;
658 (void)evaluation_function;
659#else
660 static CollectiveMutex mutex;
662
663 const unsigned int my_rank =
665
666 // invert permutation matrices (TODO: precompute)
667 std::vector<unsigned int> recv_permutation_inv(recv_permutation.size());
668 for (unsigned int c = 0; c < recv_permutation.size(); ++c)
669 recv_permutation_inv[recv_permutation[c]] = c;
670
671 std::vector<unsigned int> send_permutation_inv(send_permutation.size());
672 for (unsigned int c = 0; c < send_permutation.size(); ++c)
673 send_permutation_inv[send_permutation[c]] = c;
674
675 // allocate memory for buffer
676 const auto &point_ptrs = this->get_point_ptrs();
677 AssertDimension(input.size(), point_ptrs.size() - 1);
678 buffer.resize(std::max(send_permutation.size() * 2,
679 point_ptrs.back() + send_permutation.size()));
680
681 // ... for evaluation
682 ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
683
684 // ... for communication
685 ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
686 point_ptrs.back());
687
688 // sort for communication (and duplicate data if necessary)
689
690 const auto my_rank_local_recv_ptr =
691 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
692
693 if (my_rank_local_recv_ptr != recv_ranks.end())
694 {
695 // optimize the case if we have our own rank among the possible list
696 const unsigned int my_rank_local_recv =
697 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
698 const unsigned int my_rank_local_send = std::distance(
699 send_ranks.begin(),
700 std::find(send_ranks.begin(), send_ranks.end(), my_rank));
701
702 const unsigned int start = recv_ptrs[my_rank_local_recv];
703 const unsigned int end = recv_ptrs[my_rank_local_recv + 1];
704 const unsigned int *send_ptr =
705 send_permutation_inv.data() + send_ptrs[my_rank_local_send];
706 for (unsigned int i = 0, c = 0; i < point_ptrs.size() - 1; ++i)
707 {
708 const unsigned int next = point_ptrs[i + 1];
709 for (unsigned int j = point_ptrs[i]; j < next; ++j, ++c)
710 {
711 const unsigned int recv_index = recv_permutation_inv[c];
712
713 // local data -> can be copied to final buffer directly
714 if (start <= recv_index && recv_index < end)
715 buffer_eval[send_ptr[recv_index - start]] = input[i];
716 else // data to be sent
717 buffer_comm[recv_index] = input[i];
718 }
719 }
720 }
721 else
722 {
723 for (unsigned int i = 0, c = 0; i < point_ptrs.size() - 1; ++i)
724 for (unsigned int j = point_ptrs[i]; j < point_ptrs[i + 1];
725 ++j, ++c)
726 buffer_comm[recv_permutation_inv[c]] = input[i];
727 }
728
729 // send data
730 std::vector<std::vector<char>> send_buffer;
731 send_buffer.reserve(recv_ranks.size());
732
733 std::vector<MPI_Request> send_requests;
734 send_requests.reserve(recv_ranks.size());
735
736 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
737 {
738 if (recv_ranks[i] == my_rank)
739 continue;
740
741 send_requests.push_back(MPI_Request());
742
743 send_buffer.emplace_back(Utilities::pack(
744 std::vector<T>(buffer_comm.begin() + recv_ptrs[i],
745 buffer_comm.begin() + recv_ptrs[i + 1]),
746 false));
747
748 const int ierr = MPI_Isend(send_buffer.back().data(),
749 send_buffer.back().size(),
750 MPI_CHAR,
751 recv_ranks[i],
754 &send_requests.back());
755 AssertThrowMPI(ierr);
756 }
757
758 // receive data
759 std::vector<char> recv_buffer;
760
761 for (unsigned int i = 0; i < send_ranks.size(); ++i)
762 {
763 if (send_ranks[i] == my_rank)
764 continue;
765
766 // receive remote data
767 MPI_Status status;
768 int ierr = MPI_Probe(MPI_ANY_SOURCE,
771 &status);
772 AssertThrowMPI(ierr);
773
774 int message_length;
775 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
776 AssertThrowMPI(ierr);
777
778 recv_buffer.resize(message_length);
779
780 ierr = MPI_Recv(recv_buffer.data(),
781 recv_buffer.size(),
782 MPI_CHAR,
783 status.MPI_SOURCE,
786 MPI_STATUS_IGNORE);
787 AssertThrowMPI(ierr);
788
789 // unpack data
790 const auto recv_buffer_unpacked =
791 Utilities::unpack<std::vector<T>>(recv_buffer, false);
792
793 // write data into buffer vector
794 const auto ptr =
795 std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
796
797 Assert(ptr != send_ranks.end(), ExcNotImplemented());
798
799 const unsigned int j = std::distance(send_ranks.begin(), ptr);
800
801 AssertDimension(recv_buffer_unpacked.size(),
802 send_ptrs[j + 1] - send_ptrs[j]);
803
804 for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
805 ++i, ++c)
806 buffer_eval[send_permutation_inv[i]] = recv_buffer_unpacked[c];
807 }
808
809 if (!send_requests.empty())
810 {
811 const int ierr = MPI_Waitall(send_requests.size(),
812 send_requests.data(),
813 MPI_STATUSES_IGNORE);
814 AssertThrowMPI(ierr);
815 }
816
817 // evaluate function at points
818 evaluation_function(buffer_eval, *cell_data);
819#endif
820 }
821
822
823
824 template <int dim, int spacedim>
825 template <typename T>
826 void
828 const std::vector<T> &input,
829 const std::function<void(const ArrayView<const T> &, const CellData &)>
830 &evaluation_function) const
831 {
832 std::vector<T> buffer;
833 this->process_and_evaluate(input, buffer, evaluation_function);
834 }
835
836 } // end of namespace MPI
837} // end of namespace Utilities
838
839
841
842#endif
iterator begin() const
Definition array_view.h:702
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
virtual MPI_Comm get_communicator() const
ArrayView< T > get_data_view(const unsigned int cell, const ArrayView< T > &values) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > cell_indices() const
ArrayView< const Point< dim > > get_unit_points(const unsigned int cell) const
Triangulation< dim, spacedim >::active_cell_iterator get_active_cell_iterator(const unsigned int cell) 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
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
void reinit(const std::vector< Point< spacedim > > &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
SmartPointer< const Triangulation< dim, spacedim > > tria
std::vector< T > evaluate_and_process(const std::function< void(const ArrayView< T > &, const CellData &)> &evaluation_function) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1381
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)