Reference documentation for deal.II version GIT acb749f4f5 2024-02-21 16:10: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// 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:
61 {
62 public:
67 const double tolerance = 1e-6,
68 const bool enforce_unique_mapping = false,
69 const unsigned int rtree_level = 0,
70 const std::function<std::vector<bool>()> &marked_vertices = {});
71
80 double tolerance;
81
87
91 unsigned int rtree_level;
92
97 std::function<std::vector<bool>()> marked_vertices;
98 };
99
107
126 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
127 "Use the constructor with AdditionalData struct.")
129 const double tolerance,
130 const bool enforce_unique_mapping = false,
131 const unsigned int rtree_level = 0,
132 const std::function<std::vector<bool>()> &marked_vertices = {});
133
138
150 void
151 reinit(const std::vector<Point<spacedim>> &points,
154
164 void
166 DistributedComputePointLocationsInternal<dim, spacedim> &data,
167 const Triangulation<dim, spacedim> &tria,
168 const Mapping<dim, spacedim> &mapping);
169
175 {
176 public:
183
210 cell_indices() const;
211
216 get_active_cell_iterator(const unsigned int cell) const;
217
222 get_unit_points(const unsigned int cell) const;
223
227 template <typename T>
229 get_data_view(const unsigned int cell,
230 const ArrayView<T> &values) const;
231
235 std::vector<std::pair<int, int>> cells;
236
241 std::vector<unsigned int> reference_point_ptrs;
242
246 std::vector<Point<dim>> reference_point_values;
247
248 private:
254 };
255
260 const CellData &
261 get_cell_data() const;
262
277 template <typename T>
278 void
280 std::vector<T> &output,
281 std::vector<T> &buffer,
282 const std::function<void(const ArrayView<T> &, const CellData &)>
283 &evaluation_function) const;
284
289 template <typename T>
290 std::vector<T>
292 const std::function<void(const ArrayView<T> &, const CellData &)>
293 &evaluation_function) const;
294
303 template <typename T>
304 void
306 const std::vector<T> &input,
307 std::vector<T> &buffer,
308 const std::function<void(const ArrayView<const T> &, const CellData &)>
309 &evaluation_function) const;
310
315 template <typename T>
316 void
318 const std::vector<T> &input,
319 const std::function<void(const ArrayView<const T> &, const CellData &)>
320 &evaluation_function) const;
321
326 const std::vector<unsigned int> &
327 get_point_ptrs() const;
328
335 bool
336 is_map_unique() const;
337
341 bool
342 all_points_found() const;
343
347 bool
348 point_found(const unsigned int i) const;
349
354 get_triangulation() const;
355
360 get_mapping() const;
361
367 bool
368 is_ready() const;
369
370 private:
375
379 boost::signals2::connection tria_signal;
380
387
392
397
402
407
413 std::vector<unsigned int> point_ptrs;
414
418 std::vector<unsigned int> recv_permutation;
419
424 std::vector<unsigned int> recv_ptrs;
425
429 std::vector<unsigned int> recv_ranks;
430
435 std::unique_ptr<CellData> cell_data;
436
440 std::vector<unsigned int> send_permutation;
441
445 std::vector<unsigned int> send_ranks;
446
451 std::vector<unsigned int> send_ptrs;
452 };
453
454
455
456 template <int dim, int spacedim>
457 template <typename T>
460 const unsigned int cell,
461 const ArrayView<T> &values) const
462 {
463 AssertIndexRange(cell, cells.size());
464 return {values.data() + reference_point_ptrs[cell],
466 }
467
468
469
470 template <int dim, int spacedim>
471 template <typename T>
472 void
474 std::vector<T> &output,
475 std::vector<T> &buffer,
476 const std::function<void(const ArrayView<T> &, const CellData &)>
477 &evaluation_function) const
478 {
479#ifndef DEAL_II_WITH_MPI
480 Assert(false, ExcNeedsMPI());
481 (void)output;
482 (void)buffer;
483 (void)evaluation_function;
484#else
485 static CollectiveMutex mutex;
487
488 const unsigned int my_rank =
490
491 // allocate memory for output and buffer
492 output.resize(point_ptrs.back());
493 buffer.resize(std::max(send_permutation.size() * 2,
494 point_ptrs.back() + send_permutation.size()));
495
496 // ... for evaluation
497 ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
498
499 // ... for communication
500 ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
501 send_permutation.size());
502
503 // evaluate functions at points
504 evaluation_function(buffer_eval, *cell_data);
505
506 // sort for communication
507 unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
508 unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
509
510 const auto my_rank_local_recv_ptr =
511 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
512
513 if (my_rank_local_recv_ptr != recv_ranks.end())
514 {
515 my_rank_local_recv =
516 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
517 my_rank_local_send = std::distance(send_ranks.begin(),
518 std::find(send_ranks.begin(),
519 send_ranks.end(),
520 my_rank));
521 }
522
523 for (unsigned int i = 0; i < send_permutation.size(); ++i)
524 {
525 const unsigned int send_index = send_permutation[i];
526
527 // local data -> can be copied to output directly
528 if (my_rank_local_send != numbers::invalid_unsigned_int &&
529 (send_ptrs[my_rank_local_send] <= send_index &&
530 send_index < send_ptrs[my_rank_local_send + 1]))
531 output[recv_permutation[send_index - send_ptrs[my_rank_local_send] +
532 recv_ptrs[my_rank_local_recv]]] =
533 buffer_eval[i];
534 else // data to be sent
535 buffer_comm[send_index] = buffer_eval[i];
536 }
537
538 // send data
539 std::vector<std::vector<char>> send_buffer;
540 send_buffer.reserve(send_ranks.size());
541
542 std::vector<MPI_Request> send_requests;
543 send_requests.reserve(send_ranks.size());
544
545 for (unsigned int i = 0; i < send_ranks.size(); ++i)
546 {
547 if (send_ranks[i] == my_rank)
548 continue;
549
550 send_requests.emplace_back(MPI_Request());
551
552 send_buffer.emplace_back(Utilities::pack(
553 std::vector<T>(buffer_comm.begin() + send_ptrs[i],
554 buffer_comm.begin() + send_ptrs[i + 1]),
555 false));
556
557 const int ierr = MPI_Isend(send_buffer.back().data(),
558 send_buffer.back().size(),
559 MPI_CHAR,
560 send_ranks[i],
563 &send_requests.back());
564 AssertThrowMPI(ierr);
565 }
566
567 // receive data
568 std::vector<char> buffer_char;
569
570 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
571 {
572 if (recv_ranks[i] == my_rank)
573 continue;
574
575 // receive remote data
576 MPI_Status status;
577 int ierr = MPI_Probe(MPI_ANY_SOURCE,
580 &status);
581 AssertThrowMPI(ierr);
582
583 int message_length;
584 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
585 AssertThrowMPI(ierr);
586
587 buffer_char.resize(message_length);
588
589 ierr = MPI_Recv(buffer_char.data(),
590 buffer_char.size(),
591 MPI_CHAR,
592 status.MPI_SOURCE,
595 MPI_STATUS_IGNORE);
596 AssertThrowMPI(ierr);
597
598 // unpack data
599 const auto buffer =
600 Utilities::unpack<std::vector<T>>(buffer_char, false);
601
602 // write data into output vector
603 const auto ptr =
604 std::find(recv_ranks.begin(), recv_ranks.end(), status.MPI_SOURCE);
605
606 Assert(ptr != recv_ranks.end(), ExcNotImplemented());
607
608 const unsigned int j = std::distance(recv_ranks.begin(), ptr);
609
610 AssertDimension(buffer.size(), recv_ptrs[j + 1] - recv_ptrs[j]);
611
612 for (unsigned int i = recv_ptrs[j], c = 0; i < recv_ptrs[j + 1];
613 ++i, ++c)
614 output[recv_permutation[i]] = buffer[c];
615 }
616
617 // make sure all messages have been sent
618 if (!send_requests.empty())
619 {
620 const int ierr = MPI_Waitall(send_requests.size(),
621 send_requests.data(),
622 MPI_STATUSES_IGNORE);
623 AssertThrowMPI(ierr);
624 }
625#endif
626 }
627
628
629 template <int dim, int spacedim>
630 template <typename T>
631 std::vector<T>
633 const std::function<void(const ArrayView<T> &, const CellData &)>
634 &evaluation_function) const
635 {
636 std::vector<T> output;
637 std::vector<T> buffer;
638
639 this->evaluate_and_process(output, buffer, evaluation_function);
640
641 return output;
642 }
643
644
645
646 template <int dim, int spacedim>
647 template <typename T>
648 void
650 const std::vector<T> &input,
651 std::vector<T> &buffer,
652 const std::function<void(const ArrayView<const T> &, const CellData &)>
653 &evaluation_function) const
654 {
655#ifndef DEAL_II_WITH_MPI
656 Assert(false, ExcNeedsMPI());
657 (void)input;
658 (void)buffer;
659 (void)evaluation_function;
660#else
661 static CollectiveMutex mutex;
663
664 const unsigned int my_rank =
666
667 // invert permutation matrices (TODO: precompute)
668 std::vector<unsigned int> recv_permutation_inv(recv_permutation.size());
669 for (unsigned int c = 0; c < recv_permutation.size(); ++c)
670 recv_permutation_inv[recv_permutation[c]] = c;
671
672 std::vector<unsigned int> send_permutation_inv(send_permutation.size());
673 for (unsigned int c = 0; c < send_permutation.size(); ++c)
674 send_permutation_inv[send_permutation[c]] = c;
675
676 // allocate memory for buffer
677 const auto &point_ptrs = this->get_point_ptrs();
678 AssertDimension(input.size(), point_ptrs.size() - 1);
679 buffer.resize(std::max(send_permutation.size() * 2,
680 point_ptrs.back() + send_permutation.size()));
681
682 // ... for evaluation
683 ArrayView<T> buffer_eval(buffer.data(), send_permutation.size());
684
685 // ... for communication
686 ArrayView<T> buffer_comm(buffer.data() + send_permutation.size(),
687 point_ptrs.back());
688
689 // sort for communication (and duplicate data if necessary)
690 unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
691 unsigned int my_rank_local_send = numbers::invalid_unsigned_int;
692
693 const auto my_rank_local_recv_ptr =
694 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
695
696 if (my_rank_local_recv_ptr != recv_ranks.end())
697 {
698 my_rank_local_recv =
699 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
700 my_rank_local_send = std::distance(send_ranks.begin(),
701 std::find(send_ranks.begin(),
702 send_ranks.end(),
703 my_rank));
704 }
705
706 for (unsigned int i = 0, c = 0; i < point_ptrs.size() - 1; ++i)
707 {
708 const auto n_entries = point_ptrs[i + 1] - point_ptrs[i];
709
710 for (unsigned int j = 0; j < n_entries; ++j, ++c)
711 {
712 const unsigned int recv_index = recv_permutation_inv[c];
713
714 // local data -> can be copied to final buffer directly
715 if (my_rank_local_recv != numbers::invalid_unsigned_int &&
716 (recv_ptrs[my_rank_local_recv] <= recv_index &&
717 recv_index < recv_ptrs[my_rank_local_recv + 1]))
718 buffer_eval[send_permutation_inv
719 [recv_index - recv_ptrs[my_rank_local_recv] +
720 send_ptrs[my_rank_local_send]]] = input[i];
721 else // data to be sent
722 buffer_comm[recv_index] = input[i];
723 }
724 }
725
726 // send data
727 std::vector<std::vector<char>> send_buffer;
728 send_buffer.reserve(recv_ranks.size());
729
730 std::vector<MPI_Request> send_requests;
731 send_requests.reserve(recv_ranks.size());
732
733 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
734 {
735 if (recv_ranks[i] == my_rank)
736 continue;
737
738 send_requests.push_back(MPI_Request());
739
740 send_buffer.emplace_back(Utilities::pack(
741 std::vector<T>(buffer_comm.begin() + recv_ptrs[i],
742 buffer_comm.begin() + recv_ptrs[i + 1]),
743 false));
744
745 const int ierr = MPI_Isend(send_buffer.back().data(),
746 send_buffer.back().size(),
747 MPI_CHAR,
748 recv_ranks[i],
751 &send_requests.back());
752 AssertThrowMPI(ierr);
753 }
754
755 // receive data
756 std::vector<char> recv_buffer;
757
758 for (unsigned int i = 0; i < send_ranks.size(); ++i)
759 {
760 if (send_ranks[i] == my_rank)
761 continue;
762
763 // receive remote data
764 MPI_Status status;
765 int ierr = MPI_Probe(MPI_ANY_SOURCE,
768 &status);
769 AssertThrowMPI(ierr);
770
771 int message_length;
772 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
773 AssertThrowMPI(ierr);
774
775 recv_buffer.resize(message_length);
776
777 ierr = MPI_Recv(recv_buffer.data(),
778 recv_buffer.size(),
779 MPI_CHAR,
780 status.MPI_SOURCE,
783 MPI_STATUS_IGNORE);
784 AssertThrowMPI(ierr);
785
786 // unpack data
787 const auto recv_buffer_unpacked =
788 Utilities::unpack<std::vector<T>>(recv_buffer, false);
789
790 // write data into buffer vector
791 const auto ptr =
792 std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
793
794 Assert(ptr != send_ranks.end(), ExcNotImplemented());
795
796 const unsigned int j = std::distance(send_ranks.begin(), ptr);
797
798 AssertDimension(recv_buffer_unpacked.size(),
799 send_ptrs[j + 1] - send_ptrs[j]);
800
801 for (unsigned int i = send_ptrs[j], c = 0; i < send_ptrs[j + 1];
802 ++i, ++c)
803 buffer_eval[send_permutation_inv[i]] = recv_buffer_unpacked[c];
804 }
805
806 if (!send_requests.empty())
807 {
808 const int ierr = MPI_Waitall(send_requests.size(),
809 send_requests.data(),
810 MPI_STATUSES_IGNORE);
811 AssertThrowMPI(ierr);
812 }
813
814 // evaluate function at points
815 evaluation_function(buffer_eval, *cell_data);
816#endif
817 }
818
819
820
821 template <int dim, int spacedim>
822 template <typename T>
823 void
825 const std::vector<T> &input,
826 const std::function<void(const ArrayView<const T> &, const CellData &)>
827 &evaluation_function) const
828 {
829 std::vector<T> buffer;
830 this->process_and_evaluate(input, buffer, evaluation_function);
831 }
832
833 } // end of namespace MPI
834} // end of namespace Utilities
835
836
838
839#endif
iterator begin() const
Definition array_view.h:703
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
Definition tria.h:1595
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:144
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1365
static const unsigned int invalid_unsigned_int
Definition types.h:221
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:46
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)