hyper.deal
Loading...
Searching...
No Matches
vector_partitioner.h
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 by the hyper.deal authors
4//
5// This file is part of the hyper.deal library.
6//
7// The hyper.deal 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 3.0 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 hyper.deal.
13//
14// ---------------------------------------------------------------------
15
16#ifndef DEALII_LINEARALGEBRA_SHAREDMPI_PARTITIONER
17#define DEALII_LINEARALGEBRA_SHAREDMPI_PARTITIONER
18
19#include <hyper.deal/base/config.h>
20
21#include <deal.II/base/conditional_ostream.h>
22#include <deal.II/base/mpi.h>
23#include <deal.II/base/mpi.templates.h>
24#include <deal.II/base/mpi_compute_index_owner_internal.h>
25#include <deal.II/base/mpi_consensus_algorithms.h>
26#include <deal.II/base/timer.h>
27
28#ifdef DEAL_II_WITH_64BIT_INDICES
29# include <deal.II/base/mpi_consensus_algorithms.templates.h>
30#endif
31
32#include <deal.II/matrix_free/vector_data_exchange.h>
33
34#include <hyper.deal/base/mpi.h>
35#include <hyper.deal/base/mpi_tags.h>
36#include <hyper.deal/matrix_free/shape_info.h>
37
38#include <map>
39#include <vector>
40
41
42DEAL_II_NAMESPACE_OPEN
43
44namespace internal
45{
46 namespace MatrixFreeFunctions
47 {
48 namespace VectorDataExchange
49 {
54 class Contiguous
55 : public dealii::internal::MatrixFreeFunctions::VectorDataExchange::Base
56 {
57 const dealii::types::global_dof_index dofs_per_cell;
58 const dealii::types::global_dof_index dofs_per_face;
59 const std::vector<std::vector<unsigned int>> &face_to_cell_index_nodal;
60
61 public:
62 using RankType = unsigned int;
63 using LocalDoFType = unsigned int;
64 using CellIdType = dealii::types::global_dof_index;
65 using FaceIdType = std::pair<CellIdType, unsigned int>;
66
70 template <typename ShapeInfo>
71 Contiguous(const ShapeInfo &shape_info);
72
77 inline void
78 reinit(const std::vector<dealii::types::global_dof_index> local_cells,
79 const std::vector<std::pair<dealii::types::global_dof_index,
80 std::vector<unsigned int>>>
81 local_ghost_faces,
82 const MPI_Comm comm,
83 const MPI_Comm comm_sm,
84 const bool do_buffering);
85
86 unsigned int
87 locally_owned_size() const override
88 {
89 return n_local_elements;
90 }
91
92 unsigned int
93 n_ghost_indices() const override
94 {
95 return n_ghost_elements;
96 }
97
98 unsigned int
99 n_import_indices() const override
100 {
101 return send_ptr.back() * dofs_per_ghost;
102 }
103
104 void
105 reset_ghost_values(
106 const dealii::ArrayView<double> &ghost_array) const override
107 {
108 (void)ghost_array;
109 // TODO
110 }
111
112 void
113 reset_ghost_values(
114 const dealii::ArrayView<float> &ghost_array) const override
115 {
116 (void)ghost_array;
117 // TODO
118 }
119
120 unsigned int
121 n_import_sm_procs() const override
122 {
123 AssertThrow(false, ExcNotImplemented());
124 return 0; // TODO
125 }
126
127 types::global_dof_index
128 size() const override
129 {
130 AssertThrow(false, ExcNotImplemented());
131 return 0; // TODO
132 }
133
137 inline void
138 export_to_ghosted_array_start(
139 const unsigned int communication_channel,
140 const dealii::ArrayView<const double> &locally_owned_array,
141 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
142 const dealii::ArrayView<double> & ghost_array,
143 const dealii::ArrayView<double> & temporary_storage,
144 std::vector<MPI_Request> &requests) const override;
145
149 inline void
150 export_to_ghosted_array_finish(
151 const dealii::ArrayView<const double> &locally_owned_array,
152 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
153 const dealii::ArrayView<double> & ghost_array,
154 std::vector<MPI_Request> &requests) const override;
155
159 inline void
160 import_from_ghosted_array_start(
161 const dealii::VectorOperation::values vector_operation,
162 const unsigned int communication_channel,
163 const dealii::ArrayView<const double> &locally_owned_array,
164 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
165 const dealii::ArrayView<double> & ghost_array,
166 const dealii::ArrayView<double> & temporary_storage,
167 std::vector<MPI_Request> &requests) const override;
168
172 inline void
173 import_from_ghosted_array_finish(
174 const dealii::VectorOperation::values vector_operation,
175 const dealii::ArrayView<double> & locally_owned_storage,
176 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
177 const dealii::ArrayView<double> & ghost_array,
178 const dealii::ArrayView<const double> & temporary_storage,
179 std::vector<MPI_Request> &requests) const override;
180
184 inline void
185 export_to_ghosted_array_start(
186 const unsigned int communication_channel,
187 const dealii::ArrayView<const float> &locally_owned_array,
188 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
189 const dealii::ArrayView<float> & ghost_array,
190 const dealii::ArrayView<float> & temporary_storage,
191 std::vector<MPI_Request> &requests) const override;
192
196 inline void
197 export_to_ghosted_array_finish(
198 const dealii::ArrayView<const float> &locally_owned_array,
199 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
200 const dealii::ArrayView<float> & ghost_array,
201 std::vector<MPI_Request> &requests) const override;
202
206 inline void
207 import_from_ghosted_array_start(
208 const dealii::VectorOperation::values vector_operation,
209 const unsigned int communication_channel,
210 const dealii::ArrayView<const float> &locally_owned_array,
211 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
212 const dealii::ArrayView<float> & ghost_array,
213 const dealii::ArrayView<float> & temporary_storage,
214 std::vector<MPI_Request> &requests) const override;
215
219 inline void
220 import_from_ghosted_array_finish(
221 const dealii::VectorOperation::values vector_operation,
222 const dealii::ArrayView<float> & locally_owned_storage,
223 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
224 const dealii::ArrayView<float> & ghost_array,
225 const dealii::ArrayView<const float> & temporary_storage,
226 std::vector<MPI_Request> &requests) const override;
227
231 template <typename Number>
232 void
233 export_to_ghosted_array_finish_0(
234 const dealii::ArrayView<const Number> &locally_owned_array,
235 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
236 const dealii::ArrayView<Number> & ghost_array,
237 std::vector<MPI_Request> & requests) const;
238
242 template <typename Number>
243 void
244 export_to_ghosted_array_finish_1(
245 const dealii::ArrayView<const Number> &locally_owned_array,
246 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
247 const dealii::ArrayView<Number> & ghost_array,
248 std::vector<MPI_Request> & requests) const;
249
250 private:
255 template <typename Number>
256 void
257 export_to_ghosted_array_start_impl(
258 const unsigned int communication_channel,
259 const dealii::ArrayView<const Number> &locally_owned_array,
260 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
261 const dealii::ArrayView<Number> & ghost_array,
262 const dealii::ArrayView<Number> & temporary_storage,
263 std::vector<MPI_Request> & requests) const;
264
269 template <typename Number>
270 void
271 export_to_ghosted_array_finish_impl(
272 const dealii::ArrayView<const Number> &locally_owned_array,
273 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
274 const dealii::ArrayView<Number> & ghost_array,
275 std::vector<MPI_Request> & requests) const;
276
281 template <typename Number>
282 void
283 import_from_ghosted_array_start_impl(
284 const dealii::VectorOperation::values vector_operation,
285 const unsigned int communication_channel,
286 const dealii::ArrayView<const Number> &locally_owned_array,
287 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
288 const dealii::ArrayView<Number> & ghost_array,
289 const dealii::ArrayView<Number> & temporary_storage,
290 std::vector<MPI_Request> & requests) const;
291
296 template <typename Number>
297 void
298 import_from_ghosted_array_finish_impl(
299 const dealii::VectorOperation::values vector_operation,
300 const dealii::ArrayView<Number> & locally_owned_storage,
301 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
302 const dealii::ArrayView<Number> & ghost_array,
303 const dealii::ArrayView<const Number> & temporary_storage,
304 std::vector<MPI_Request> & requests) const;
305
306 public:
310 inline const std::map<dealii::types::global_dof_index,
311 std::pair<unsigned int, unsigned int>> &
312 get_maps() const;
313
314
318 inline const std::map<
319 std::pair<dealii::types::global_dof_index, unsigned int>,
320 std::pair<unsigned int, unsigned int>> &
321 get_maps_ghost() const;
322
328 inline std::size_t
329 memory_consumption() const;
330
334 inline void
335 sync(const unsigned int tag = 0) const;
336
337 private:
341 MPI_Comm comm;
342
346 MPI_Comm comm_sm;
347
351 unsigned int n_mpi_processes_;
352
356 unsigned int n_local_elements;
357
361 unsigned int n_ghost_elements;
362
363 // I) configuration parameters
364 bool do_buffering; // buffering vs. non-buffering modus
365 unsigned int dofs_per_ghost; // ghost face or ghost cell
366
367 // II) MPI-communicator related stuff
368 unsigned int sm_size;
369 unsigned int sm_rank;
370
371 // III) access cells and ghost faces
372 std::map<CellIdType, std::pair<RankType, LocalDoFType>> maps;
373 std::map<FaceIdType, std::pair<RankType, LocalDoFType>> maps_ghost;
374
375 // III) information to pack/unpack buffers
376 std::vector<unsigned int> send_ranks;
377 std::vector<dealii::types::global_dof_index> send_ptr;
378 std::vector<dealii::types::global_dof_index> send_data_id;
379 std::vector<unsigned int> send_data_face_no;
380
381 std::vector<unsigned int> recv_ranks;
382 std::vector<dealii::types::global_dof_index> recv_ptr;
383 std::vector<dealii::types::global_dof_index> recv_size;
384
385 std::vector<unsigned int> sm_targets;
386 std::vector<unsigned int> sm_sources;
387
388
389 std::vector<dealii::types::global_dof_index> sm_send_ptr;
390 std::vector<unsigned int> sm_send_rank;
391 std::vector<unsigned int> sm_send_offset_1;
392 std::vector<unsigned int> sm_send_offset_2;
393 std::vector<unsigned int> sm_send_no;
394
395 std::vector<dealii::types::global_dof_index> sm_recv_ptr;
396 std::vector<unsigned int> sm_recv_rank;
397 std::vector<unsigned int> sm_recv_offset_1;
398 std::vector<unsigned int> sm_recv_offset_2;
399 std::vector<unsigned int> sm_recv_no;
400 };
401
402
403 inline void
404 Contiguous::export_to_ghosted_array_start(
405 const unsigned int communication_channel,
406 const dealii::ArrayView<const double> &locally_owned_array,
407 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
408 const dealii::ArrayView<double> & ghost_array,
409 const dealii::ArrayView<double> & temporary_storage,
410 std::vector<MPI_Request> & requests) const
411 {
412 export_to_ghosted_array_start_impl(communication_channel,
413 locally_owned_array,
414 shared_arrays,
415 ghost_array,
416 temporary_storage,
417 requests);
418 }
419
420
421
422 inline void
423 Contiguous::export_to_ghosted_array_finish(
424 const dealii::ArrayView<const double> & locally_owned_array,
425 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
426 const dealii::ArrayView<double> & ghost_array,
427 std::vector<MPI_Request> & requests) const
428 {
429 export_to_ghosted_array_finish_impl(locally_owned_array,
430 shared_arrays,
431 ghost_array,
432 requests);
433 }
434
435
436
437 inline void
438 Contiguous::import_from_ghosted_array_start(
439 const dealii::VectorOperation::values vector_operation,
440 const unsigned int communication_channel,
441 const dealii::ArrayView<const double> &locally_owned_array,
442 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
443 const dealii::ArrayView<double> & ghost_array,
444 const dealii::ArrayView<double> & temporary_storage,
445 std::vector<MPI_Request> & requests) const
446 {
447 import_from_ghosted_array_start_impl(vector_operation,
448 communication_channel,
449 locally_owned_array,
450 shared_arrays,
451 ghost_array,
452 temporary_storage,
453 requests);
454 }
455
456
457
458 inline void
459 Contiguous::import_from_ghosted_array_finish(
460 const dealii::VectorOperation::values vector_operation,
461 const dealii::ArrayView<double> & locally_owned_storage,
462 const std::vector<dealii::ArrayView<const double>> &shared_arrays,
463 const dealii::ArrayView<double> & ghost_array,
464 const dealii::ArrayView<const double> & temporary_storage,
465 std::vector<MPI_Request> & requests) const
466 {
467 import_from_ghosted_array_finish_impl(vector_operation,
468 locally_owned_storage,
469 shared_arrays,
470 ghost_array,
471 temporary_storage,
472 requests);
473 }
474
475
476
477 inline void
478 Contiguous::export_to_ghosted_array_start(
479 const unsigned int communication_channel,
480 const dealii::ArrayView<const float> &locally_owned_array,
481 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
482 const dealii::ArrayView<float> & ghost_array,
483 const dealii::ArrayView<float> & temporary_storage,
484 std::vector<MPI_Request> & requests) const
485 {
486 export_to_ghosted_array_start_impl(communication_channel,
487 locally_owned_array,
488 shared_arrays,
489 ghost_array,
490 temporary_storage,
491 requests);
492 }
493
494
495
496 inline void
497 Contiguous::export_to_ghosted_array_finish(
498 const dealii::ArrayView<const float> & locally_owned_array,
499 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
500 const dealii::ArrayView<float> & ghost_array,
501 std::vector<MPI_Request> & requests) const
502 {
503 export_to_ghosted_array_finish_impl(locally_owned_array,
504 shared_arrays,
505 ghost_array,
506 requests);
507 }
508
509
510
511 inline void
512 Contiguous::import_from_ghosted_array_start(
513 const dealii::VectorOperation::values vector_operation,
514 const unsigned int communication_channel,
515 const dealii::ArrayView<const float> &locally_owned_array,
516 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
517 const dealii::ArrayView<float> & ghost_array,
518 const dealii::ArrayView<float> & temporary_storage,
519 std::vector<MPI_Request> & requests) const
520 {
521 import_from_ghosted_array_start_impl(vector_operation,
522 communication_channel,
523 locally_owned_array,
524 shared_arrays,
525 ghost_array,
526 temporary_storage,
527 requests);
528 }
529
530
531
532 inline void
533 Contiguous::import_from_ghosted_array_finish(
534 const dealii::VectorOperation::values vector_operation,
535 const dealii::ArrayView<float> & locally_owned_storage,
536 const std::vector<dealii::ArrayView<const float>> &shared_arrays,
537 const dealii::ArrayView<float> & ghost_array,
538 const dealii::ArrayView<const float> & temporary_storage,
539 std::vector<MPI_Request> & requests) const
540 {
541 import_from_ghosted_array_finish_impl(vector_operation,
542 locally_owned_storage,
543 shared_arrays,
544 ghost_array,
545 temporary_storage,
546 requests);
547 }
548
549
550
551 inline void
552 Contiguous::sync(const unsigned int tag) const
553 {
554 std::vector<MPI_Request> req(sm_targets.size() + sm_sources.size());
555
556 for (unsigned int i = 0; i < sm_targets.size(); i++)
557 {
558 int dummy;
559 MPI_Isend(&dummy,
560 0,
561 MPI_INT,
562 sm_targets[i],
563 tag,
564 this->comm_sm,
565 req.data() + i);
566 }
567
568 for (unsigned int i = 0; i < sm_sources.size(); i++)
569 {
570 int dummy;
571 MPI_Irecv(&dummy,
572 0,
573 MPI_INT,
574 sm_sources[i],
575 tag,
576 this->comm_sm,
577 req.data() + i + sm_targets.size());
578 }
579
580 MPI_Waitall(req.size(), req.data(), MPI_STATUSES_IGNORE);
581 }
582
583
584
585 namespace internal
586 {
587 template <typename T, typename U>
588 std::vector<std::pair<T, U>>
589 MPI_Allgather_Pairs(const std::vector<std::pair<T, U>> &src,
590 const MPI_Comm & comm)
591 {
592 int size;
593 MPI_Comm_size(comm, &size);
594
595 std::vector<T> src_1;
596 std::vector<U> src_2;
597
598 for (auto i : src)
599 {
600 src_1.push_back(i.first);
601 src_2.push_back(i.second);
602 }
603
604
605 unsigned int len_local = src_1.size();
606 std::vector<int> len_global(
607 size); // actually unsigned int but MPI wants int
608 MPI_Allgather(
609 &len_local,
610 1,
611 MPI_INT,
612 &len_global[0],
613 1,
614 dealii::Utilities::MPI::mpi_type_id_for_type<decltype(len_local)>,
615 comm);
616
617
618 std::vector<int> displs; // actually unsigned int but MPI wants int
619 displs.push_back(0);
620
621 int total_size = 0;
622
623 for (auto i : len_global)
624 {
625 displs.push_back(i + displs.back());
626 total_size += i;
627 }
628
629 std::vector<T> dst_1(total_size);
630 std::vector<U> dst_2(total_size);
631 MPI_Allgatherv(
632 src_1.data(),
633 len_local,
634 dealii::Utilities::MPI::mpi_type_id_for_type<decltype(src_1[0])>,
635 dst_1.data(),
636 len_global.data(),
637 displs.data(),
638 dealii::Utilities::MPI::mpi_type_id_for_type<decltype(dst_1[0])>,
639 comm);
640 MPI_Allgatherv(
641 src_2.data(),
642 len_local,
643 dealii::Utilities::MPI::mpi_type_id_for_type<decltype(src_2[0])>,
644 dst_2.data(),
645 len_global.data(),
646 displs.data(),
647 dealii::Utilities::MPI::mpi_type_id_for_type<decltype(dst_2[0])>,
648 comm);
649
650 std::vector<std::pair<T, U>> dst(total_size);
651
652 for (unsigned int i = 0; i < dst_1.size(); i++)
653 dst[i] = {dst_1[i], dst_2[i]};
654
655 return dst;
656 }
657 } // namespace internal
658
659
660
661 template <typename ShapeInfo>
662 Contiguous::Contiguous(const ShapeInfo &shape_info)
663 : dofs_per_cell(shape_info.dofs_per_cell)
664 , dofs_per_face(shape_info.dofs_per_face)
665 , face_to_cell_index_nodal(shape_info.face_to_cell_index_nodal)
666 {}
667
668
669
670 inline void
671 Contiguous::reinit(
672 const std::vector<dealii::types::global_dof_index> local_cells,
673 const std::vector<
674 std::pair<dealii::types::global_dof_index, std::vector<unsigned int>>>
675 local_ghost_faces,
676 const MPI_Comm comm,
677 const MPI_Comm comm_sm,
678 const bool do_buffering)
679 {
680 // fill some information needed by PartitionerBase
681 this->comm = comm;
682 this->comm_sm = comm_sm;
683 this->n_mpi_processes_ = dealii::Utilities::MPI::n_mpi_processes(comm);
684
685 this->do_buffering = do_buffering;
686
687 AssertThrow(local_cells.size() > 0,
688 dealii::ExcMessage("No local cells!"));
689
690 this->n_local_elements = local_cells.size() * dofs_per_cell;
691
692 // 1) determine if ghost faces or ghost cells are needed
693 const dealii::types::global_dof_index dofs_per_ghost = [&]() {
694 unsigned int result = dofs_per_face;
695
696 for (const auto &ghost_faces : local_ghost_faces)
697 for (const auto ghost_face : ghost_faces.second)
698 if (ghost_face == dealii::numbers::invalid_unsigned_int)
699 result = dofs_per_cell;
700 return dealii::Utilities::MPI::max(result, comm);
701 }();
702
703 this->dofs_per_ghost = dofs_per_ghost;
704
705 // const auto this->comm_sm = create_sm(comm);
706 const auto sm_procs = hyperdeal::mpi::procs_of_sm(comm, this->comm_sm);
707 const auto sm_rank = [&]() {
708 const auto ptr =
709 std::find(sm_procs.begin(),
710 sm_procs.end(),
711 dealii::Utilities::MPI::this_mpi_process(comm));
712
713 AssertThrow(ptr != sm_procs.end(),
714 dealii::ExcMessage("Proc not found!"));
715
716 return std::distance(sm_procs.begin(), ptr);
717 }();
718
719 this->sm_rank = sm_rank;
720 this->sm_size = sm_procs.size();
721
722
723 for (unsigned int i = 0; i < local_cells.size(); i++)
724 this->maps[local_cells[i]] = {sm_rank, i * dofs_per_cell};
725
726
727 // 2) determine which ghost face is shared or remote
728 std::vector<
729 std::pair<dealii::types::global_dof_index, std::vector<unsigned int>>>
730 local_ghost_faces_remote, local_ghost_faces_shared;
731 {
732 const auto n_total_cells = dealii::Utilities::MPI::sum(
733 static_cast<dealii::types::global_dof_index>(local_cells.size()),
734 comm);
735
736 dealii::IndexSet is_local_cells(n_total_cells);
737 is_local_cells.add_indices(local_cells.begin(), local_cells.end());
738
739 dealii::IndexSet is_ghost_cells(n_total_cells);
740 for (const auto &ghost_faces : local_ghost_faces)
741 is_ghost_cells.add_index(ghost_faces.first);
742
743 for (unsigned int i = 0; i < local_ghost_faces.size(); i++)
744 AssertThrow(local_ghost_faces[i].first ==
745 is_ghost_cells.nth_index_in_set(i),
746 dealii::ExcMessage("PROBLEM!"));
747
748 AssertThrow(
749 local_ghost_faces.size() == is_ghost_cells.n_elements(),
750 dealii::ExcMessage(
751 "Dimensions " + std::to_string(local_ghost_faces.size()) + " " +
752 std::to_string(is_ghost_cells.n_elements()) + " do not match!"));
753
754 std::vector<unsigned int> owning_ranks_of_ghosts(
755 is_ghost_cells.n_elements());
756
757 // set up dictionary
758 dealii::Utilities::MPI::internal::ComputeIndexOwner::
759 ConsensusAlgorithmsPayload process(is_local_cells,
760 is_ghost_cells,
761 comm,
762 owning_ranks_of_ghosts,
763 false);
764
765 dealii::Utilities::MPI::ConsensusAlgorithms::Selector<
766 std::vector<std::pair<dealii::types::global_dof_index,
767 dealii::types::global_dof_index>>,
768 std::vector<unsigned int>>
769 consensus_algorithm;
770 consensus_algorithm.run(process, comm);
771
772 std::map<unsigned int, std::vector<dealii::types::global_dof_index>>
773 shared_procs_to_cells;
774
775 for (unsigned int i = 0; i < owning_ranks_of_ghosts.size(); i++)
776 {
777 AssertThrow(dealii::Utilities::MPI::this_mpi_process(comm) !=
778 owning_ranks_of_ghosts[i],
779 dealii::ExcMessage(
780 "Locally owned cells should be not ghosted!"));
781
782 const auto ptr = std::find(sm_procs.begin(),
783 sm_procs.end(),
784 owning_ranks_of_ghosts[i]);
785
786 if (ptr == sm_procs.end())
787 local_ghost_faces_remote.push_back(local_ghost_faces[i]);
788 else
789 {
790 local_ghost_faces_shared.push_back(local_ghost_faces[i]);
791 shared_procs_to_cells[std::distance(sm_procs.begin(), ptr)]
792 .emplace_back(local_ghost_faces[i].first);
793 }
794 }
795
796 // determine (ghost sm cell) -> (sm rank, offset)
797 dealii::Utilities::MPI::ConsensusAlgorithms::selector<
798 std::vector<dealii::types::global_dof_index>,
799 std::vector<unsigned int>>(
800 [&]() {
801 std::vector<unsigned int> result;
802 for (auto &i : shared_procs_to_cells)
803 result.push_back(i.first);
804 return result;
805 }(),
806 [&](const auto other_rank) {
807 return shared_procs_to_cells[other_rank];
808 },
809 [&](const auto &other_rank, const auto &buffer_recv) {
810 std::vector<unsigned int> request_buffer(buffer_recv.size());
811
812 for (unsigned int i = 0; i < buffer_recv.size(); i++)
813 {
814 const auto value = buffer_recv[i];
815 const auto ptr =
816 std::find(local_cells.begin(), local_cells.end(), value);
817
818 AssertThrow(ptr != local_cells.end(),
819 dealii::ExcMessage(
820 "Cell " + std::to_string(value) + " at index " +
821 std::to_string(i) + " on rank " +
822 std::to_string(
823 dealii::Utilities::MPI::this_mpi_process(
824 comm)) +
825 " requested by rank " +
826 std::to_string(other_rank) + " not found!"));
827
828 request_buffer[i] = std::distance(local_cells.begin(), ptr);
829 }
830 return request_buffer;
831 },
832 [&](const auto other_rank, const auto &recv_buffer) {
833 for (unsigned int i = 0; i < recv_buffer.size(); i++)
834 {
835 const dealii::types::global_dof_index cell =
836 shared_procs_to_cells[other_rank][i];
837 const unsigned int offset = recv_buffer[i];
838
839 Assert(maps.find(cell) == maps.end(),
840 dealii::ExcMessage("Cell " + std::to_string(cell) +
841 " is already in maps!"));
842
843 this->maps[cell] = {other_rank, offset * dofs_per_cell};
844 }
845 },
846 this->comm_sm);
847 }
848
849
850 // 3) merge local_ghost_faces_remote and sort -> ghost_faces_remote
851 const auto local_ghost_faces_remote_pairs_global =
852 [&local_ghost_faces_remote, this]() {
853 std::vector<
854 std::pair<dealii::types::global_dof_index, unsigned int>>
855 local_ghost_faces_remote_pairs_local;
856
857 // convert vector<pair<U, std::vector<V>>> ->.vector<std::pair<U,
858 // V>>>
859 for (const auto &ghost_faces : local_ghost_faces_remote)
860 for (const auto ghost_face : ghost_faces.second)
861 local_ghost_faces_remote_pairs_local.emplace_back(
862 ghost_faces.first, ghost_face);
863
864 // collect all on which are shared
865 std::vector<
866 std::pair<dealii::types::global_dof_index, unsigned int>>
867 local_ghost_faces_remote_pairs_global =
868 internal::MPI_Allgather_Pairs(
869 local_ghost_faces_remote_pairs_local, this->comm_sm);
870
871 // sort
872 std::sort(local_ghost_faces_remote_pairs_global.begin(),
873 local_ghost_faces_remote_pairs_global.end());
874
875 return local_ghost_faces_remote_pairs_global;
876 }();
877
878
879 // 4) distributed ghost_faces_remote,
880 auto distributed_local_ghost_faces_remote_pairs_global =
881 [&local_ghost_faces_remote_pairs_global, &sm_procs]() {
882 std::vector<std::vector<
883 std::pair<dealii::types::global_dof_index, unsigned int>>>
884 result(sm_procs.size());
885
886 unsigned int counter = 0;
887 const unsigned int faces_per_process =
888 (local_ghost_faces_remote_pairs_global.size() + sm_procs.size() -
889 1) /
890 sm_procs.size();
891 for (auto p : local_ghost_faces_remote_pairs_global)
892 result[(counter++) / faces_per_process].push_back(p);
893
894 return result;
895 }();
896
897 // revert partitioning of ghost faces (TODO)
898 {
899 std::vector<std::pair<dealii::types::global_dof_index, unsigned int>>
900 local_ghost_faces_remote_pairs_local;
901
902 for (const auto &ghost_faces : local_ghost_faces_remote)
903 for (const auto ghost_face : ghost_faces.second)
904 local_ghost_faces_remote_pairs_local.emplace_back(
905 ghost_faces.first, ghost_face);
906
907 distributed_local_ghost_faces_remote_pairs_global.clear();
908 distributed_local_ghost_faces_remote_pairs_global.resize(
909 sm_procs.size());
910 distributed_local_ghost_faces_remote_pairs_global[sm_rank] =
911 local_ghost_faces_remote_pairs_local;
912 }
913
914
915 // ... update ghost size, and
916 this->n_ghost_elements =
917 (distributed_local_ghost_faces_remote_pairs_global[sm_rank].size() +
918 (do_buffering ?
919 std::accumulate(local_ghost_faces_shared.begin(),
920 local_ghost_faces_shared.end(),
921 std::size_t(0),
922 [](const std::size_t &a, const auto &b) {
923 return a + b.second.size();
924 }) :
925 std::size_t(0))) *
926 dofs_per_ghost;
927
928
929 // ... update ghost map
930 this->maps_ghost = [&]() {
931 std::map<std::pair<dealii::types::global_dof_index, unsigned int>,
932 std::pair<unsigned int, unsigned int>>
933 maps_ghost;
934
935 // counter for offset of ghost faces
936 unsigned int my_offset = local_cells.size() * dofs_per_cell;
937
938 // buffering-mode: insert shared faces into maps_ghost
939 if (do_buffering)
940 for (const auto &pair : local_ghost_faces_shared)
941 for (const auto &face : pair.second)
942 {
943 maps_ghost[{pair.first, face}] = {sm_rank, my_offset};
944 my_offset += dofs_per_ghost;
945 }
946
947 // create map (cell, face_no) -> (sm rank, offset)
948 const auto maps_ghost_inverse =
949 [&distributed_local_ghost_faces_remote_pairs_global,
950 &dofs_per_ghost,
951 this,
952 &sm_procs,
953 &my_offset]() {
954 std::vector<unsigned int> offsets(sm_procs.size());
955
956 MPI_Allgather(&my_offset,
957 1,
958 dealii::Utilities::MPI::mpi_type_id_for_type<
959 decltype(my_offset)>,
960 offsets.data(),
961 1,
962 dealii::Utilities::MPI::mpi_type_id_for_type<
963 decltype(my_offset)>,
964 this->comm_sm);
965
966 std::map<std::pair<unsigned int, unsigned int>,
967 std::pair<dealii::types::global_dof_index, unsigned int>>
968 maps_ghost_inverse;
969
970 for (unsigned int i = 0; i < sm_procs.size(); i++)
971 for (unsigned int j = 0;
972 j < distributed_local_ghost_faces_remote_pairs_global[i]
973 .size();
974 j++)
975 maps_ghost_inverse
976 [distributed_local_ghost_faces_remote_pairs_global[i][j]] =
977 {i, offsets[i] + j * dofs_per_ghost};
978
979 return maps_ghost_inverse;
980 }();
981
982 // create map (cell, face_no) -> (sm rank, offset) with only
983 // ghost faces needed for evaluation
984 for (const auto &i : local_ghost_faces_remote)
985 for (const auto &j : i.second)
986 maps_ghost[{i.first, j}] = maps_ghost_inverse.at({i.first, j});
987
988 for (const auto &i :
989 distributed_local_ghost_faces_remote_pairs_global[sm_rank])
990 maps_ghost[i] = maps_ghost_inverse.at(i);
991
992 return maps_ghost;
993 }();
994
995 // buffering-mode: pre-compute information for memcopy of ghost faces
996 std::vector<std::array<LocalDoFType, 5>> ghost_list_shared_precomp;
997 std::vector<std::array<LocalDoFType, 5>> maps_ghost_inverse_precomp;
998
999 if (do_buffering)
1000 {
1001 std::map<unsigned int, std::vector<std::array<LocalDoFType, 5>>>
1002 send_data;
1003
1004 for (const auto &pair : local_ghost_faces_shared)
1005 for (const auto &face : pair.second)
1006 {
1007 const auto ptr1 = maps.find(pair.first);
1008 AssertThrow(ptr1 != maps.end(),
1009 dealii::ExcMessage("Entry not found!"));
1010
1011 const auto ptr2 = maps_ghost.find({pair.first, face});
1012 AssertThrow(ptr2 != maps_ghost.end(),
1013 dealii::ExcMessage("Entry not found!"));
1014
1015 std::array<LocalDoFType, 5> v{{ptr1->second.first,
1016 ptr1->second.second,
1017 face,
1018 ptr2->second.first,
1019 ptr2->second.second}};
1020
1021 ghost_list_shared_precomp.push_back(v);
1022
1023 send_data[ptr1->second.first].push_back(v);
1024 }
1025
1026 dealii::Utilities::MPI::ConsensusAlgorithms::selector<
1027 std::vector<LocalDoFType>>(
1028 [&]() {
1029 std::vector<unsigned int> result;
1030 for (auto &i : send_data)
1031 result.push_back(i.first);
1032 return result;
1033 }(),
1034 [&](const auto other_rank) {
1035 std::vector<LocalDoFType> send_buffer;
1036
1037 for (const auto &i : send_data[other_rank])
1038 for (const auto &j : i)
1039 send_buffer.push_back(j);
1040
1041 return send_buffer;
1042 },
1043 [&](const auto & /*other_rank*/, const auto &buffer_recv) {
1044 for (unsigned int i = 0; i < buffer_recv.size(); i += 5)
1045 maps_ghost_inverse_precomp.push_back({{buffer_recv[i],
1046 buffer_recv[i + 1],
1047 buffer_recv[i + 2],
1048 buffer_recv[i + 3],
1049 buffer_recv[i + 4]}});
1050 },
1051 this->comm_sm);
1052
1053 std::sort(maps_ghost_inverse_precomp.begin(),
1054 maps_ghost_inverse_precomp.end());
1055 }
1056
1057 // rank -> pair(offset, size)
1058 std::map<unsigned int, std::pair<unsigned int, unsigned int>>
1059 receive_info;
1060 std::map<unsigned int, std::vector<std::array<unsigned int, 3>>>
1061 requests_from_relevant_precomp;
1062
1063 // 5) setup communication patterns (during update_ghost_values &
1064 // compress)
1065 [&local_cells,
1066 &distributed_local_ghost_faces_remote_pairs_global,
1067 &sm_rank,
1068 &comm,
1069 &dofs_per_ghost](auto & requests_from_relevant_precomp,
1070 auto & receive_info,
1071 const auto &maps,
1072 const auto &maps_ghost) {
1073 // determine of the owner of cells of remote ghost faces
1074 const auto n_total_cells = dealii::Utilities::MPI::sum(
1075 static_cast<dealii::types::global_dof_index>(local_cells.size()),
1076 comm);
1077
1078 // owned cells (TODO: generalize so that local_cells is also
1079 // partitioned)
1080 dealii::IndexSet is_local_cells(n_total_cells);
1081 is_local_cells.add_indices(local_cells.begin(), local_cells.end());
1082
1083 // needed (ghost) cell
1084 dealii::IndexSet is_ghost_cells(n_total_cells);
1085 for (const auto &ghost_faces :
1086 distributed_local_ghost_faces_remote_pairs_global[sm_rank])
1087 is_ghost_cells.add_index(ghost_faces.first);
1088
1089 // determine rank of (ghost) cells
1090 const auto owning_ranks_of_ghosts = [&]() {
1091 std::vector<unsigned int> owning_ranks_of_ghosts(
1092 is_ghost_cells.n_elements());
1093
1094 dealii::Utilities::MPI::internal::ComputeIndexOwner::
1095 ConsensusAlgorithmsPayload process(is_local_cells,
1096 is_ghost_cells,
1097 comm,
1098 owning_ranks_of_ghosts,
1099 false);
1100
1101 dealii::Utilities::MPI::ConsensusAlgorithms::Selector<
1102 std::vector<std::pair<dealii::types::global_dof_index,
1103 dealii::types::global_dof_index>>,
1104 std::vector<unsigned int>>
1105 consensus_algorithm;
1106 consensus_algorithm.run(process, comm);
1107
1108 return owning_ranks_of_ghosts;
1109 }();
1110
1111 // determine targets
1112 const auto send_ranks = [&]() {
1113 std::set<unsigned int> send_ranks_set;
1114
1115 for (const auto &i : owning_ranks_of_ghosts)
1116 send_ranks_set.insert(i);
1117
1118 const std::vector<unsigned int> send_ranks(send_ranks_set.begin(),
1119 send_ranks_set.end());
1120
1121 return send_ranks;
1122 }();
1123
1124 // collect ghost faces (separated for each target)
1125 const auto send_data = [&]() {
1126 std::vector<std::vector<std::pair<dealii::types::global_dof_index,
1127 dealii::types::global_dof_index>>>
1128 send_data(send_ranks.size());
1129
1130 unsigned int index = 0;
1131 unsigned int index_cell = dealii::numbers::invalid_unsigned_int;
1132
1133 for (const auto &ghost_faces :
1134 distributed_local_ghost_faces_remote_pairs_global[sm_rank])
1135 {
1136 if (index_cell != ghost_faces.first)
1137 {
1138 index_cell = ghost_faces.first;
1139 const unsigned int index_rank =
1140 owning_ranks_of_ghosts[is_ghost_cells.index_within_set(
1141 ghost_faces.first)];
1142 index = std::distance(send_ranks.begin(),
1143 std::find(send_ranks.begin(),
1144 send_ranks.end(),
1145 index_rank));
1146 }
1147 send_data[index].emplace_back(ghost_faces.first,
1148 ghost_faces.second);
1149 }
1150
1151 return send_data;
1152 }();
1153
1154 // send ghost faces to the owners
1155 std::vector<MPI_Request> send_requests(send_ranks.size());
1156
1157 for (unsigned int i = 0; i < send_ranks.size(); i++)
1158 {
1159 MPI_Isend(send_data[i].data(),
1160 2 * send_data[i].size(),
1161 dealii::Utilities::MPI::mpi_type_id_for_type<
1162 dealii::types::global_dof_index>,
1163 send_ranks[i],
1164 105,
1165 comm,
1166 send_requests.data() + i);
1167
1168 receive_info[send_ranks[i]] = {
1169 send_data[i].size() * dofs_per_ghost,
1170 maps_ghost.at(send_data[i][0]).second};
1171 }
1172
1173 MPI_Barrier(comm);
1174
1175 const auto targets = dealii::Utilities::MPI::
1176 compute_point_to_point_communication_pattern(comm, send_ranks);
1177
1178 // process requests
1179 for (unsigned int i = 0; i < targets.size(); i++)
1180 {
1181 // wait for any request
1182 MPI_Status status;
1183 auto ierr = MPI_Probe(MPI_ANY_SOURCE, 105, comm, &status);
1184 AssertThrowMPI(ierr);
1185
1186 // determine number of ghost faces * 2 (since we are considering
1187 // pairs)
1188 int len;
1189 MPI_Get_count(&status,
1190 dealii::Utilities::MPI::mpi_type_id_for_type<
1191 dealii::types::global_dof_index>,
1192 &len);
1193
1194 AssertThrow(len % 2 == 0,
1195 dealii::ExcMessage("Length " + std::to_string(len) +
1196 " is not a multiple of two!"));
1197
1198 // allocate memory for the incoming vector
1199 std::vector<std::pair<dealii::types::global_dof_index,
1200 dealii::types::global_dof_index>>
1201 recv_data(len / 2);
1202
1203 // receive data
1204 ierr = MPI_Recv(recv_data.data(),
1205 len,
1206 dealii::Utilities::MPI::mpi_type_id_for_type<
1207 dealii::types::global_dof_index>,
1208 status.MPI_SOURCE,
1209 status.MPI_TAG,
1210 comm,
1211 &status);
1212 AssertThrowMPI(ierr);
1213
1214 // setup pack and unpack info
1215 requests_from_relevant_precomp[status.MPI_SOURCE] = [&]() {
1216 std::vector<std::array<unsigned int, 3>> temp(len / 2);
1217 for (unsigned int i = 0; i < static_cast<unsigned int>(len) / 2;
1218 i++)
1219 {
1220 const CellIdType cell = recv_data[i].first;
1221 const unsigned int face_no = recv_data[i].second;
1222
1223 const auto ptr = maps.find(cell);
1224 AssertThrow(ptr != maps.end(),
1225 dealii::ExcMessage("Entry " +
1226 std::to_string(cell) +
1227 " not found!"));
1228
1229 temp[i] = std::array<unsigned int, 3>{
1230 {ptr->second.first,
1231 (unsigned int)ptr->second.second,
1232 face_no}};
1233 }
1234 return temp;
1235 }();
1236 }
1237
1238 // make sure requests have been sent away
1239 MPI_Waitall(send_requests.size(),
1240 send_requests.data(),
1241 MPI_STATUSES_IGNORE);
1242 }(requests_from_relevant_precomp,
1243 receive_info,
1244 this->maps,
1245 this->maps_ghost);
1246
1247
1248 {
1249 recv_ptr.clear();
1250 recv_size.clear();
1251 recv_ranks.clear();
1252
1253 for (const auto &i : receive_info)
1254 {
1255 recv_ranks.push_back(i.first);
1256 recv_size.push_back(i.second.first);
1257 recv_ptr.push_back(i.second.second -
1258 local_cells.size() * dofs_per_cell);
1259 }
1260 }
1261
1262 {
1263 // TODO: clear
1264 send_ptr.push_back(0);
1265
1266 for (const auto &i : requests_from_relevant_precomp)
1267 {
1268 send_ranks.push_back(i.first);
1269
1270 for (const auto &j : i.second)
1271 {
1272 AssertThrow(j[0] == sm_rank,
1273 dealii::StandardExceptions::ExcNotImplemented());
1274 send_data_id.push_back(j[1]);
1275 send_data_face_no.push_back(j[2]);
1276 }
1277
1278 send_ptr.push_back(send_data_id.size());
1279 }
1280 }
1281
1282 {
1283 std::set<unsigned int> temp;
1284
1285 for (const auto &i : maps)
1286 if (i.second.first != sm_rank)
1287 temp.insert(i.second.first);
1288
1289 for (const auto &i : temp)
1290 this->sm_sources.push_back(i);
1291
1292 this->sm_targets = dealii::Utilities::MPI::
1293 compute_point_to_point_communication_pattern(this->comm_sm,
1294 sm_sources);
1295 }
1296
1297 if (do_buffering)
1298 {
1299 auto temp = ghost_list_shared_precomp;
1300 std::sort(temp.begin(), temp.end());
1301
1302 std::vector<
1303 std::vector<std::array<dealii::types::global_dof_index, 3>>>
1304 temp_(sm_size);
1305
1306 for (auto t : temp)
1307 {
1308 AssertThrow(sm_rank == t[3],
1309 dealii::StandardExceptions::ExcNotImplemented());
1310 temp_[t[0]].emplace_back(
1311 std::array<dealii::types::global_dof_index, 3>{
1312 {t[1], t[2], t[4]}});
1313 }
1314
1315
1316 sm_send_ptr.push_back(0);
1317
1318 for (unsigned int i = 0; i < temp_.size(); i++)
1319 {
1320 if (temp_[i].size() == 0)
1321 continue;
1322
1323 sm_send_rank.push_back(i);
1324
1325 for (const auto &v : temp_[i])
1326 {
1327 sm_send_offset_1.push_back(v[2] - local_cells.size() *
1328 dofs_per_cell);
1329 sm_send_offset_2.push_back(v[0]);
1330 sm_send_no.push_back(v[1]);
1331 }
1332
1333 sm_send_ptr.push_back(sm_send_no.size());
1334 }
1335
1336 AssertThrow(sm_send_rank.size() == sm_sources.size(),
1337 dealii::StandardExceptions::ExcNotImplemented());
1338 }
1339
1340 if (do_buffering)
1341 {
1342 auto temp = maps_ghost_inverse_precomp;
1343 std::sort(temp.begin(), temp.end());
1344
1345 std::vector<
1346 std::vector<std::array<dealii::types::global_dof_index, 3>>>
1347 temp_(sm_size);
1348
1349 for (auto t : temp)
1350 {
1351 AssertThrow(sm_rank == t[0],
1352 dealii::StandardExceptions::ExcNotImplemented());
1353 temp_[t[3]].emplace_back(
1354 std::array<dealii::types::global_dof_index, 3>{
1355 {t[4], t[2], t[1]}});
1356 }
1357
1358
1359 sm_recv_ptr.push_back(0);
1360
1361 for (unsigned int i = 0; i < temp_.size(); i++)
1362 {
1363 if (temp_[i].size() == 0)
1364 continue;
1365
1366 sm_recv_rank.push_back(i);
1367
1368 for (const auto &v : temp_[i])
1369 {
1370 sm_recv_offset_1.push_back(v[2]);
1371 sm_recv_offset_2.push_back(v[0]);
1372 sm_recv_no.push_back(v[1]);
1373 }
1374
1375 sm_recv_ptr.push_back(sm_recv_no.size());
1376 }
1377
1378 AssertThrow(sm_recv_rank.size() == sm_targets.size(),
1379 dealii::StandardExceptions::ExcNotImplemented());
1380 }
1381 }
1382
1383
1384
1385 template <typename Number>
1386 void
1387 Contiguous::export_to_ghosted_array_start_impl(
1388 const unsigned int communication_channel,
1389 const dealii::ArrayView<const Number> &locally_owned_array,
1390 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
1391 const dealii::ArrayView<Number> & ghost_array,
1392 const dealii::ArrayView<Number> & temporary_storage,
1393 std::vector<MPI_Request> & requests) const
1394 {
1395 (void)shared_arrays;
1396
1397 AssertThrow(temporary_storage.size() ==
1398 send_ptr.back() * dofs_per_ghost,
1399 dealii::StandardExceptions::ExcDimensionMismatch(
1400 temporary_storage.size(),
1401 send_ptr.back() * dofs_per_ghost));
1402
1403 requests.resize(sm_sources.size() + sm_targets.size() +
1404 recv_ranks.size() + send_ranks.size());
1405
1406 // 1) notify relevant shared processes that local data is available
1407 if (sm_size > 1)
1408 {
1409 int dummy;
1410 for (unsigned int i = 0; i < sm_targets.size(); i++)
1411 MPI_Isend(&dummy,
1412 0,
1413 MPI_INT,
1414 sm_targets[i],
1415 communication_channel + 21,
1416 this->comm_sm,
1417 requests.data() + i + sm_sources.size());
1418
1419 for (unsigned int i = 0; i < sm_sources.size(); i++)
1420 MPI_Irecv(&dummy,
1421 0,
1422 MPI_INT,
1423 sm_sources[i],
1424 communication_channel + 21,
1425 this->comm_sm,
1426 requests.data() + i);
1427 }
1428
1429 // 2) start receiving form (remote) processes
1430 {
1431 for (unsigned int i = 0; i < recv_ranks.size(); i++)
1432 MPI_Irecv(ghost_array.data() + recv_ptr[i],
1433 recv_size[i],
1434 MPI_DOUBLE,
1435 recv_ranks[i],
1436 communication_channel + 22,
1437 comm,
1438 requests.data() + i + sm_sources.size() +
1439 sm_targets.size());
1440 }
1441
1442 // 3) fill buffers and start sending to (remote) processes
1443 for (unsigned int c = 0; c < send_ranks.size(); c++)
1444 {
1445 auto buffer =
1446 temporary_storage.data() + send_ptr[c] * dofs_per_ghost;
1447
1448 for (unsigned int i = send_ptr[c]; i < send_ptr[c + 1];
1449 i++, buffer += dofs_per_ghost)
1450 if (dofs_per_ghost == dofs_per_face)
1451 {
1452 auto *__restrict dst = buffer;
1453 const auto *__restrict src =
1454 locally_owned_array.data() + send_data_id[i];
1455 const auto *__restrict idx =
1456 face_to_cell_index_nodal[send_data_face_no[i]].data();
1457
1458 for (unsigned int i = 0; i < dofs_per_face; i++)
1459 dst[i] = src[idx[i]];
1460 }
1461 else if (dofs_per_ghost == dofs_per_cell)
1462 {
1463 AssertThrow(false,
1464 dealii::StandardExceptions::ExcNotImplemented());
1465 }
1466
1467 MPI_Issend(temporary_storage.data() + send_ptr[c] * dofs_per_ghost,
1468 (send_ptr[c + 1] - send_ptr[c]) * dofs_per_ghost,
1469 MPI_DOUBLE,
1470 send_ranks[c],
1471 communication_channel + 22,
1472 comm,
1473 requests.data() + c + sm_sources.size() +
1474 sm_targets.size() + recv_ranks.size());
1475 }
1476 }
1477
1478
1479
1480 template <typename Number>
1481 void
1482 Contiguous::export_to_ghosted_array_finish_impl(
1483 const dealii::ArrayView<const Number> & locally_owned_array,
1484 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
1485 const dealii::ArrayView<Number> & ghost_array,
1486 std::vector<MPI_Request> & requests) const
1487 {
1488 (void)locally_owned_array;
1489
1490 AssertDimension(requests.size(),
1491 sm_sources.size() + sm_targets.size() +
1492 recv_ranks.size() + send_ranks.size());
1493
1494 if (do_buffering) // deal with shared faces if buffering is requested
1495 {
1496 // update ghost values of shared cells (if requested)
1497 for (unsigned int c = 0; c < sm_sources.size(); c++)
1498 {
1499 int i;
1500 MPI_Status status;
1501 const auto ierr =
1502 MPI_Waitany(sm_sources.size(), requests.data(), &i, &status);
1503 AssertThrowMPI(ierr);
1504
1505 for (unsigned int j = sm_send_ptr[i]; j < sm_send_ptr[i + 1];
1506 j++)
1507 if (dofs_per_ghost == dofs_per_face)
1508 {
1509 auto *__restrict dst =
1510 ghost_array.data() + sm_send_offset_1[j];
1511 const auto *__restrict src =
1512 shared_arrays[sm_send_rank[i]].data() +
1513 sm_send_offset_2[j];
1514 const auto *__restrict idx =
1515 face_to_cell_index_nodal[sm_send_no[j]].data();
1516
1517 for (unsigned int i = 0; i < dofs_per_face; i++)
1518 dst[i] = src[idx[i]];
1519 }
1520 else if (dofs_per_ghost == dofs_per_cell)
1521 {
1522 AssertThrow(
1523 false, dealii::StandardExceptions::ExcNotImplemented());
1524 }
1525 }
1526 }
1527
1528 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1529 }
1530
1531
1532
1533 template <typename Number>
1534 void
1535 Contiguous::export_to_ghosted_array_finish_0(
1536 const dealii::ArrayView<const Number> & locally_owned_array,
1537 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
1538 const dealii::ArrayView<Number> & ghost_array,
1539 std::vector<MPI_Request> & requests) const
1540 {
1541 (void)locally_owned_array;
1542
1543 AssertDimension(requests.size(),
1544 sm_sources.size() + sm_targets.size() +
1545 recv_ranks.size() + send_ranks.size());
1546
1547 if (do_buffering) // deal with shared faces if buffering is requested
1548 {
1549 // update ghost values of shared cells (if requested)
1550 for (unsigned int c = 0; c < sm_sources.size(); c++)
1551 {
1552 int i;
1553 MPI_Status status;
1554 const auto ierr =
1555 MPI_Waitany(sm_sources.size(), requests.data(), &i, &status);
1556 AssertThrowMPI(ierr);
1557
1558 for (unsigned int j = sm_send_ptr[i]; j < sm_send_ptr[i + 1];
1559 j++)
1560 if (dofs_per_ghost == dofs_per_face)
1561 {
1562 auto *__restrict dst =
1563 ghost_array.data() + sm_send_offset_1[j];
1564 const auto *__restrict src =
1565 shared_arrays[sm_send_rank[i]].data() +
1566 sm_send_offset_2[j];
1567 const auto *__restrict idx =
1568 face_to_cell_index_nodal[sm_send_no[j]].data();
1569
1570 for (unsigned int i = 0; i < dofs_per_face; i++)
1571 dst[i] = src[idx[i]];
1572 }
1573 else if (dofs_per_ghost == dofs_per_cell)
1574 {
1575 AssertThrow(
1576 false, dealii::StandardExceptions::ExcNotImplemented());
1577 }
1578 }
1579 }
1580 else
1581 {
1582 MPI_Waitall(sm_sources.size(),
1583 requests.data(),
1584 MPI_STATUSES_IGNORE);
1585 }
1586 }
1587
1588
1589
1590 template <typename Number>
1591 void
1592 Contiguous::export_to_ghosted_array_finish_1(
1593 const dealii::ArrayView<const Number> & locally_owned_array,
1594 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
1595 const dealii::ArrayView<Number> & ghost_array,
1596 std::vector<MPI_Request> & requests) const
1597 {
1598 (void)locally_owned_array;
1599 (void)shared_arrays;
1600 (void)ghost_array;
1601
1602 AssertDimension(requests.size(),
1603 sm_sources.size() + sm_targets.size() +
1604 recv_ranks.size() + send_ranks.size());
1605
1606 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1607 }
1608
1609
1610
1611 template <typename Number>
1612 void
1613 Contiguous::import_from_ghosted_array_start_impl(
1614 const dealii::VectorOperation::values operation,
1615 const unsigned int communication_channel,
1616 const dealii::ArrayView<const Number> &locally_owned_array,
1617 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
1618 const dealii::ArrayView<Number> & ghost_array,
1619 const dealii::ArrayView<Number> & temporary_storage,
1620 std::vector<MPI_Request> & requests) const
1621 {
1622 (void)locally_owned_array;
1623 (void)shared_arrays;
1624 (void)communication_channel;
1625 (void)ghost_array;
1626
1627 AssertThrow(operation == dealii::VectorOperation::add,
1628 dealii::ExcMessage("Not yet implemented."));
1629
1630 AssertThrow(temporary_storage.size() ==
1631 send_ptr.back() * dofs_per_ghost,
1632 dealii::StandardExceptions::ExcDimensionMismatch(
1633 temporary_storage.size(),
1634 send_ptr.back() * dofs_per_ghost));
1635
1636 requests.resize(sm_sources.size() + sm_targets.size() +
1637 recv_ranks.size() + send_ranks.size());
1638
1639 // 1) notify relevant shared processes that data is available
1640 if (sm_size > 1)
1641 {
1642 int dummy;
1643
1644 for (unsigned int i = 0; i < sm_sources.size(); i++)
1645 MPI_Isend(&dummy,
1646 0,
1647 MPI_INT,
1648 sm_sources[i],
1649 communication_channel + 21,
1650 this->comm_sm,
1651 requests.data() + i);
1652
1653 for (unsigned int i = 0; i < sm_targets.size(); i++)
1654 MPI_Irecv(&dummy,
1655 0,
1656 MPI_INT,
1657 sm_targets[i],
1658 communication_channel + 21,
1659 this->comm_sm,
1660 requests.data() + i + sm_sources.size());
1661 }
1662
1663 // request receive
1664 {
1665 for (unsigned int i = 0; i < recv_ranks.size(); i++)
1666 MPI_Isend(ghost_array.data() + recv_ptr[i],
1667 recv_size[i],
1668 MPI_DOUBLE,
1669 recv_ranks[i],
1670 0,
1671 comm,
1672 requests.data() + i + sm_sources.size() +
1673 sm_targets.size());
1674 }
1675
1676 // fill buffers and request send
1677 for (unsigned int i = 0; i < send_ranks.size(); i++)
1678 MPI_Irecv(temporary_storage.data() + send_ptr[i] * dofs_per_ghost,
1679 (send_ptr[i + 1] - send_ptr[i]) * dofs_per_ghost,
1680 MPI_DOUBLE,
1681 send_ranks[i],
1682 0,
1683 comm,
1684 requests.data() + i + sm_sources.size() +
1685 sm_targets.size() + recv_ranks.size());
1686 }
1687
1688
1689
1690 template <typename Number>
1691 void
1692 Contiguous::import_from_ghosted_array_finish_impl(
1693 const dealii::VectorOperation::values operation,
1694 const dealii::ArrayView<Number> & locally_owned_array,
1695 const std::vector<dealii::ArrayView<const Number>> &shared_arrays,
1696 const dealii::ArrayView<Number> & ghost_array,
1697 const dealii::ArrayView<const Number> & temporary_storage,
1698 std::vector<MPI_Request> & requests) const
1699 {
1700 (void)ghost_array;
1701
1702 AssertThrow(operation == dealii::VectorOperation::add,
1703 dealii::ExcMessage("Not yet implemented."));
1704
1705 AssertThrow(temporary_storage.size() ==
1706 send_ptr.back() * dofs_per_ghost,
1707 dealii::StandardExceptions::ExcDimensionMismatch(
1708 temporary_storage.size(),
1709 send_ptr.back() * dofs_per_ghost));
1710
1711 AssertDimension(requests.size(),
1712 sm_sources.size() + sm_targets.size() +
1713 recv_ranks.size() + send_ranks.size());
1714
1715 // 1) compress for shared faces
1716 if (do_buffering)
1717 {
1718 for (unsigned int c = 0; c < sm_targets.size(); c++)
1719 {
1720 int i;
1721 MPI_Status status;
1722 const auto ierr =
1723 MPI_Waitany(sm_targets.size(),
1724 requests.data() + sm_sources.size(),
1725 &i,
1726 &status);
1727 AssertThrowMPI(ierr);
1728
1729 for (unsigned int j = sm_recv_ptr[i]; j < sm_recv_ptr[i + 1];
1730 j++)
1731 if (dofs_per_ghost == dofs_per_face)
1732 {
1733 auto *__restrict dst =
1734 locally_owned_array.data() + sm_recv_offset_1[j];
1735 const auto *__restrict src =
1736 shared_arrays[sm_recv_rank[i]].data() +
1737 sm_recv_offset_2[j];
1738 const auto *__restrict idx =
1739 face_to_cell_index_nodal[sm_recv_no[j]].data();
1740
1741 for (unsigned int i = 0; i < dofs_per_face; i++)
1742 dst[idx[i]] += src[i];
1743 }
1744 else if (dofs_per_ghost == dofs_per_cell)
1745 {
1746 AssertThrow(
1747 false, dealii::StandardExceptions::ExcNotImplemented());
1748 }
1749 }
1750 }
1751
1752 // 2) receive data and compress for remote faces
1753 for (unsigned int c = 0; c < send_ranks.size(); c++)
1754 {
1755 int r;
1756 MPI_Status status;
1757 const auto ierr =
1758 MPI_Waitany(send_ranks.size(),
1759 requests.data() + sm_sources.size() +
1760 sm_targets.size() + recv_ranks.size(),
1761 &r,
1762 &status);
1763 AssertThrowMPI(ierr);
1764
1765 auto buffer =
1766 temporary_storage.data() + send_ptr[r] * dofs_per_ghost;
1767
1768 for (unsigned int i = send_ptr[r]; i < send_ptr[r + 1];
1769 i++, buffer += dofs_per_ghost)
1770 if (dofs_per_ghost == dofs_per_face)
1771 {
1772 auto *__restrict dst =
1773 locally_owned_array.data() + send_data_id[i];
1774 const auto *__restrict src = buffer;
1775 const auto *__restrict idx =
1776 face_to_cell_index_nodal[send_data_face_no[i]].data();
1777
1778 for (unsigned int i = 0; i < dofs_per_face; i++)
1779 dst[idx[i]] += src[i];
1780 }
1781 else if (dofs_per_ghost == dofs_per_cell)
1782 {
1783 AssertThrow(false,
1784 dealii::StandardExceptions::ExcNotImplemented());
1785 }
1786 }
1787
1788 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1789 }
1790
1791
1792
1793 inline const std::map<dealii::types::global_dof_index,
1794 std::pair<unsigned int, unsigned int>> &
1795 Contiguous::get_maps() const
1796 {
1797 return maps;
1798 }
1799
1800
1801
1802 inline const std::map<
1803 std::pair<dealii::types::global_dof_index, unsigned int>,
1804 std::pair<unsigned int, unsigned int>> &
1805 Contiguous::get_maps_ghost() const
1806 {
1807 return maps_ghost;
1808 }
1809
1810
1811
1812 inline std::size_t
1813 Contiguous::memory_consumption() const
1814 {
1815 // [TODO] not counting maps and maps_ghost
1816
1817 return dealii::MemoryConsumption::memory_consumption(send_ranks) +
1818 dealii::MemoryConsumption::memory_consumption(send_ptr) +
1819 dealii::MemoryConsumption::memory_consumption(send_data_id) +
1820 dealii::MemoryConsumption::memory_consumption(
1821 send_data_face_no) +
1822 dealii::MemoryConsumption::memory_consumption(recv_ranks) +
1823 dealii::MemoryConsumption::memory_consumption(recv_ptr) +
1824 dealii::MemoryConsumption::memory_consumption(recv_size) +
1825 dealii::MemoryConsumption::memory_consumption(sm_targets) +
1826 dealii::MemoryConsumption::memory_consumption(sm_sources) +
1827 dealii::MemoryConsumption::memory_consumption(sm_send_ptr) +
1828 dealii::MemoryConsumption::memory_consumption(sm_send_rank) +
1829 dealii::MemoryConsumption::memory_consumption(sm_send_offset_1) +
1830 dealii::MemoryConsumption::memory_consumption(sm_send_offset_2) +
1831 dealii::MemoryConsumption::memory_consumption(sm_send_no) +
1832 dealii::MemoryConsumption::memory_consumption(sm_recv_ptr) +
1833 dealii::MemoryConsumption::memory_consumption(sm_recv_rank) +
1834 dealii::MemoryConsumption::memory_consumption(sm_recv_offset_1) +
1835 dealii::MemoryConsumption::memory_consumption(sm_recv_offset_2) +
1836 dealii::MemoryConsumption::memory_consumption(sm_recv_no);
1837 }
1838
1839 } // namespace VectorDataExchange
1840 } // namespace MatrixFreeFunctions
1841} // namespace internal
1842
1843DEAL_II_NAMESPACE_CLOSE
1844
1845#endif
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::Contiguous
Contiguous(const ShapeInfo &shape_info)
Definition vector_partitioner.h:662
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::reinit
void reinit(const std::vector< dealii::types::global_dof_index > local_cells, const std::vector< std::pair< dealii::types::global_dof_index, std::vector< unsigned int > > > local_ghost_faces, const MPI_Comm comm, const MPI_Comm comm_sm, const bool do_buffering)
Definition vector_partitioner.h:671
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::export_to_ghosted_array_start
void export_to_ghosted_array_start(const unsigned int communication_channel, const dealii::ArrayView< const double > &locally_owned_array, const std::vector< dealii::ArrayView< const double > > &shared_arrays, const dealii::ArrayView< double > &ghost_array, const dealii::ArrayView< double > &temporary_storage, std::vector< MPI_Request > &requests) const override
Definition vector_partitioner.h:404
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::sync
void sync(const unsigned int tag=0) const
Definition vector_partitioner.h:552
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::get_maps_ghost
const std::map< std::pair< dealii::types::global_dof_index, unsigned int >, std::pair< unsigned int, unsigned int > > & get_maps_ghost() const
Definition vector_partitioner.h:1805
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::import_from_ghosted_array_finish
void import_from_ghosted_array_finish(const dealii::VectorOperation::values vector_operation, const dealii::ArrayView< double > &locally_owned_storage, const std::vector< dealii::ArrayView< const double > > &shared_arrays, const dealii::ArrayView< double > &ghost_array, const dealii::ArrayView< const double > &temporary_storage, std::vector< MPI_Request > &requests) const override
Definition vector_partitioner.h:459
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::export_to_ghosted_array_finish_1
void export_to_ghosted_array_finish_1(const dealii::ArrayView< const Number > &locally_owned_array, const std::vector< dealii::ArrayView< const Number > > &shared_arrays, const dealii::ArrayView< Number > &ghost_array, std::vector< MPI_Request > &requests) const
Definition vector_partitioner.h:1592
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::export_to_ghosted_array_finish
void export_to_ghosted_array_finish(const dealii::ArrayView< const double > &locally_owned_array, const std::vector< dealii::ArrayView< const double > > &shared_arrays, const dealii::ArrayView< double > &ghost_array, std::vector< MPI_Request > &requests) const override
Definition vector_partitioner.h:423
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::memory_consumption
std::size_t memory_consumption() const
Definition vector_partitioner.h:1813
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::export_to_ghosted_array_finish_0
void export_to_ghosted_array_finish_0(const dealii::ArrayView< const Number > &locally_owned_array, const std::vector< dealii::ArrayView< const Number > > &shared_arrays, const dealii::ArrayView< Number > &ghost_array, std::vector< MPI_Request > &requests) const
Definition vector_partitioner.h:1535
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::import_from_ghosted_array_start
void import_from_ghosted_array_start(const dealii::VectorOperation::values vector_operation, const unsigned int communication_channel, const dealii::ArrayView< const double > &locally_owned_array, const std::vector< dealii::ArrayView< const double > > &shared_arrays, const dealii::ArrayView< double > &ghost_array, const dealii::ArrayView< double > &temporary_storage, std::vector< MPI_Request > &requests) const override
Definition vector_partitioner.h:438
internal::MatrixFreeFunctions::VectorDataExchange::Contiguous::get_maps
const std::map< dealii::types::global_dof_index, std::pair< unsigned int, unsigned int > > & get_maps() const
Definition vector_partitioner.h:1795
Generated by doxygen 1.12.0