AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
IndexDistribution.inc.hpp
1#pragma once
2
3#include <utility>
4
5#include <dune/common/timer.hh>
6
7#if HAVE_MPI
8#include <amdis/common/parallel/Collective.hpp>
9#include <amdis/common/parallel/RequestOperations.hpp>
10#endif
11
12namespace AMDiS {
13
14#if HAVE_MPI
15
16template <class GlobalId, class LI>
17 template <class Basis>
18void PetscParallelIndexDistribution<GlobalId,LI>::
19update(Basis const& basis)
20{
21 Dune::Timer t;
22
23 // clear all vectors and reset sizes
24 reset();
25 buildParallelIndexSet(basis, *indexSet_);
26
27 remoteIndices_->setIndexSets(*indexSet_, *indexSet_, world_);
28 remoteIndices_->template rebuild<true>();
29
30 // 1. insert and number owner DOFs
31 globalIndices_.resize(indexSet_->size(), 0);
32 owner_.resize(indexSet_->size(), false);
33 localSize_ = 0;
34 ghostSize_ = 0;
35 for (auto const& ip : *indexSet_) {
36 if (ip.local().attribute() == Attribute::owner) {
37 size_type idx = ip.local();
38 globalIndices_[idx] = localSize_++;
39 owner_[idx] = true;
40 } else {
41 ghostSize_++;
42 }
43 }
44
45 // communicate the local sizes from all processors
46 Mpi::all_gather(world_, localSize_, sizes_);
47
48 // at which global index do the local partitions start
49 starts_.resize(world_.size() + 1);
50 starts_[0] = 0;
51 std::partial_sum(sizes_.begin(), sizes_.end(), starts_.begin()+1);
52 globalSize_ = starts_.back();
53
54 // update the global index for all local indices by shifting by global start position
55 for (auto& i : globalIndices_)
56 i += starts_[world_.rank()];
57
58 // build up the communication of overlap DOFs that do not yet have a global index
59 // assigned. Therefore send the global index for all already computed owner DOFs
60 // to the neighboring remote processors. And receive from those their owner DOFs
61 // global indices.
62 using GlobalAssoc = std::pair<typename IndexSet::GlobalIndex, size_type>; // {globalId, globalIndex}
63 std::vector<std::vector<GlobalAssoc>> sendList(world_.size());
64 std::vector<std::size_t> receiveList(world_.size(), 0);
65
66 // Communicate attributes at the interface
67 for (auto const& rim : *remoteIndices_) {
68 int p = rim.first;
69
70 auto* sourceRemoteIndexList = rim.second.first;
71 auto* targetRemoteIndexList = rim.second.second;
72
73 // send to overlap
74 for (auto const& ri : *sourceRemoteIndexList) {
75 auto const& lip = ri.localIndexPair();
76 Attribute remoteAttr = ri.attribute();
77 Attribute myAttr = lip.local().attribute();
78 if (myAttr == Attribute::owner && remoteAttr != Attribute::owner) {
79 size_type globalIndex = globalIndices_[size_type(lip.local())];
80 sendList[p].push_back({lip.global(), globalIndex});
81 }
82 }
83
84 // receive from owner
85 for (auto const& ri : *targetRemoteIndexList) {
86 auto const& lip = ri.localIndexPair();
87 Attribute remoteAttr = ri.attribute();
88 Attribute myAttr = lip.local().attribute();
89 if (myAttr != Attribute::owner && remoteAttr == Attribute::owner) {
90 receiveList[p]++;
91 }
92 }
93 }
94 // all ghostDOFs must be communicated!
95 assert(ghostSize_ == std::accumulate(receiveList.begin(), receiveList.end(), 0u));
96
97 // send {globalId, globalIndex} to remote processors
98 std::vector<Mpi::Request> sendRequests;
99 for (int p = 0; p < world_.size(); ++p) {
100 if (!sendList[p].empty()) {
101 sendRequests.emplace_back( world_.isend(sendList[p], p, tag_) );
102 }
103 }
104
105 // receive {globalID, globalIndex} from remote processors
106 std::vector<Mpi::Request> recvRequests;
107 std::vector<std::vector<GlobalAssoc>> recvData(world_.size());
108 for (int p = 0; p < world_.size(); ++p) {
109 if (receiveList[p] > 0)
110 recvRequests.emplace_back( world_.irecv(recvData[p], p, tag_) );
111 }
112
113 Mpi::wait_all(recvRequests.begin(), recvRequests.end());
114
115 ghostIndices_.reserve(ghostSize_);
116 ghostLocalIndices_.resize(indexSet_->size(), LocalIndex(-1));
117
118 // insert all remote global indices into the map
119 std::size_t counter = 0;
120 for (int p = 0; p < world_.size(); ++p) {
121 auto const& data = recvData[p];
122 assert(data.size() == receiveList[p]);
123 for (auto const& d : data) {
124 typename IndexSet::IndexPair const& l = indexSet_->at(d.first);
125 assert(!owner_[size_type(l.local())]);
126
127 globalIndices_[size_type(l.local())] = d.second;
128 ghostIndices_.push_back(d.second);
129 ghostLocalIndices_[size_type(l.local())] = counter++;
130 }
131 }
132 assert(counter == ghostSize_);
133 assert(ghostSize_ + localSize_ == indexSet_->size());
134
135 Mpi::wait_all(sendRequests.begin(), sendRequests.end());
136 msg("update IndexDistribution need {} sec", t.elapsed());
137}
138
139
140template <class GlobalId, class LI>
141void PetscParallelIndexDistribution<GlobalId,LI>::
142debug() const
143{
144 int p = world_.rank();
145 std::cout << "[" << p << "] sizes_.size()=" << sizes_.size() << ", my_size=" << sizes_[p] << std::endl;
146 std::cout << "[" << p << "] starts_.size()=" << starts_.size() << ", my_start=" << starts_[p] << std::endl;
147 std::cout << "[" << p << "] localSize_=" << localSize_ << ", globalSize_=" << globalSize_ << ", ghostSize_=" << ghostSize_ << std::endl;
148 std::cout << "[" << p << "] globalIndices_.size()=" << globalIndices_.size() << std::endl;
149 std::cout << "[" << p << "] ghostIndices_.size()=" << ghostIndices_.size() << std::endl;
150 std::cout << "[" << p << "] ghostLocalIndices_.size()=" << ghostLocalIndices_.size() << std::endl;
151 std::cout << "[" << p << "] owner_.size()=" << owner_.size() << std::endl;
152}
153#endif
154
155} // end namespace AMDiS