Line data Source code
1 : namespace ippl{
2 :
3 : template <typename T>
4 86 : FEMVector<T>::FEMVector(size_t n, std::vector<size_t> neighbors,
5 : std::vector< Kokkos::View<size_t*> > sendIdxs,
6 : std::vector< Kokkos::View<size_t*> > recvIdxs) :
7 258 : data_m("FEMVector::data", n), boundaryInfo_m(new BoundaryInfo(std::move(neighbors),
8 172 : std::move(sendIdxs), std::move(recvIdxs))) {
9 :
10 86 : }
11 :
12 :
13 : template <typename T>
14 22 : FEMVector<T>::FEMVector(size_t n) : data_m("FEMVector::data", n), boundaryInfo_m(nullptr){
15 :
16 22 : }
17 :
18 :
19 :
20 : template <typename T>
21 1308 : FEMVector<T>::FEMVector(const FEMVector<T>& other) : data_m(other.data_m),
22 1308 : boundaryInfo_m(other.boundaryInfo_m) {
23 :
24 1308 : }
25 :
26 :
27 : template <typename T>
28 86 : FEMVector<T>::BoundaryInfo::BoundaryInfo(std::vector<size_t> neighbors,
29 : std::vector< Kokkos::View<size_t*> > sendIdxs,
30 : std::vector< Kokkos::View<size_t*> > recvIdxs) :
31 86 : neighbors_m(neighbors), sendIdxs_m(sendIdxs), recvIdxs_m(recvIdxs) {
32 :
33 86 : }
34 :
35 :
36 : template <typename T>
37 10 : void FEMVector<T>::fillHalo() {
38 : // check that we have halo information
39 10 : if (!boundaryInfo_m) {
40 0 : throw IpplException(
41 : "FEMVector::fillHalo()",
42 : "Cannot do halo operations, as no MPI communication information is provided. "
43 : "Did you use the correct constructor to construct the FEMVector?");
44 : }
45 :
46 : using memory_space = typename Kokkos::View<size_t*>::memory_space;
47 : // List of MPI requests which we send
48 10 : std::vector<MPI_Request> requests(boundaryInfo_m->neighbors_m.size());
49 :
50 : // Send loop.
51 22 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
52 12 : size_t neighborRank = boundaryInfo_m->neighbors_m[i];
53 12 : size_t nsends = boundaryInfo_m->sendIdxs_m[i].extent(0);
54 12 : size_t tag = mpi::tag::FEMVECTOR + ippl::Comm->rank();
55 :
56 : // pack data, i.e., copy values from data_m to commBuffer_m
57 12 : pack(boundaryInfo_m->sendIdxs_m[i]);
58 :
59 :
60 : // ippl MPI communication which sends data to neighbors
61 12 : mpi::Communicator::buffer_type<memory_space> archive =
62 : ippl::Comm->getBuffer<memory_space, T>(nsends);
63 12 : ippl::Comm->isend(neighborRank, tag, boundaryInfo_m->commBuffer_m, *archive,
64 : requests[i], nsends);
65 12 : archive->resetWritePos();
66 : }
67 :
68 : // Recieve loop
69 22 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
70 12 : size_t neighborRank = boundaryInfo_m->neighbors_m[i];
71 12 : size_t nrecvs = boundaryInfo_m->recvIdxs_m[i].extent(0);
72 12 : size_t tag = mpi::tag::FEMVECTOR + boundaryInfo_m->neighbors_m[i];
73 :
74 : // ippl MPI communication which will recive data from neighbors,
75 : // will put data into commBuffer_m
76 12 : mpi::Communicator::buffer_type<memory_space> archive =
77 : ippl::Comm->getBuffer<memory_space, T>(nrecvs);
78 12 : ippl::Comm->recv(neighborRank, tag, boundaryInfo_m->commBuffer_m, *archive,
79 : nrecvs * sizeof(T), nrecvs);
80 12 : archive->resetReadPos();
81 :
82 : // unpack recieved data, i.e., copy from commBuffer_m to data_m
83 12 : unpack<Assign>(boundaryInfo_m->recvIdxs_m[i]);
84 : }
85 :
86 10 : if (requests.size() > 0) {
87 10 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
88 : }
89 10 : ippl::Comm->freeAllBuffers();
90 10 : }
91 :
92 :
93 : template <typename T>
94 10 : void FEMVector<T>::accumulateHalo() {
95 : // check that we have halo information
96 10 : if (!boundaryInfo_m) {
97 0 : throw IpplException(
98 : "FEMVector::accumulateHalo()",
99 : "Cannot do halo operations, as no MPI communication information is provided. "
100 : "Did you use the correct constructor to construct the FEMVector?");
101 : }
102 :
103 : using memory_space = typename Kokkos::View<size_t*>::memory_space;
104 : // List of MPI requests which we send
105 10 : std::vector<MPI_Request> requests(boundaryInfo_m->neighbors_m.size());
106 :
107 : // Send loop.
108 22 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
109 12 : size_t neighborRank = boundaryInfo_m->neighbors_m[i];
110 12 : size_t nsends = boundaryInfo_m->recvIdxs_m[i].extent(0);
111 12 : size_t tag = mpi::tag::FEMVECTOR + ippl::Comm->rank();
112 :
113 : // pack data, i.e., copy values from data_m to commBuffer_m
114 12 : pack(boundaryInfo_m->recvIdxs_m[i]);
115 :
116 :
117 : // ippl MPI communication which sends data to neighbors
118 12 : mpi::Communicator::buffer_type<memory_space> archive =
119 : ippl::Comm->getBuffer<memory_space, T>(nsends);
120 12 : ippl::Comm->isend(neighborRank, tag, boundaryInfo_m->commBuffer_m, *archive,
121 : requests[i], nsends);
122 12 : archive->resetWritePos();
123 : }
124 :
125 : // Receive loop
126 22 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
127 12 : size_t neighborRank = boundaryInfo_m->neighbors_m[i];
128 12 : size_t nrecvs = boundaryInfo_m->sendIdxs_m[i].extent(0);
129 12 : size_t tag = mpi::tag::FEMVECTOR + boundaryInfo_m->neighbors_m[i];
130 :
131 : // ippl MPI communication which will receive data from neighbors,
132 : // will put data into commBuffer_m
133 12 : mpi::Communicator::buffer_type<memory_space> archive =
134 : ippl::Comm->getBuffer<memory_space, T>(nrecvs);
135 12 : ippl::Comm->recv(neighborRank, tag, boundaryInfo_m->commBuffer_m, *archive,
136 : nrecvs * sizeof(T), nrecvs);
137 12 : archive->resetReadPos();
138 :
139 : // unpack recieved data, i.e., copy from commBuffer_m to data_m
140 12 : unpack<AssignAdd>(boundaryInfo_m->sendIdxs_m[i]);
141 : }
142 :
143 10 : if (requests.size() > 0) {
144 10 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
145 : }
146 10 : ippl::Comm->freeAllBuffers();
147 10 : }
148 :
149 :
150 : template <typename T>
151 106 : void FEMVector<T>::setHalo(T setValue) {
152 : // check that we have halo information
153 106 : if (!boundaryInfo_m) {
154 102 : return;
155 : throw IpplException(
156 : "FEMVector::setHalo()",
157 : "Cannot do halo operations, as no MPI communication information is provided. "
158 : "Did you use the correct constructor to construct the FEMVector?");
159 : }
160 :
161 12 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
162 8 : auto& view = boundaryInfo_m->recvIdxs_m[i];
163 8 : Kokkos::parallel_for("FEMVector::setHalo()",view.extent(0),
164 16 : KOKKOS_CLASS_LAMBDA(const size_t& j){
165 32 : data_m[view(j)] = setValue;
166 : }
167 : );
168 : }
169 : }
170 :
171 :
172 : template <typename T>
173 48 : FEMVector<T>& FEMVector<T>::operator= (T value) {
174 48 : Kokkos::parallel_for("FEMVector::operator=(T value)", data_m.extent(0),
175 2128 : KOKKOS_CLASS_LAMBDA(const size_t& i){
176 400 : data_m[i] = value;
177 : }
178 : );
179 48 : return *this;
180 : }
181 :
182 :
183 : template <typename T>
184 : template <typename E, size_t N>
185 168 : FEMVector<T>& FEMVector<T>::operator= (const detail::Expression<E, N>& expr) {
186 : using capture_type = detail::CapturedExpression<E, N>;
187 168 : capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
188 168 : Kokkos::parallel_for("FEMVector::operator=(Expression)", data_m.extent(0),
189 3696 : KOKKOS_CLASS_LAMBDA(const size_t& i){
190 1680 : data_m[i] = expr_(i);
191 : }
192 : );
193 168 : return *this;
194 : }
195 :
196 :
197 : template <typename T>
198 76 : FEMVector<T>& FEMVector<T>::operator= (const FEMVector<T>& v) {
199 76 : auto view = v.getView();
200 76 : Kokkos::parallel_for("FEMVector::operator=(FEMVector)", data_m.extent(0),
201 152 : KOKKOS_CLASS_LAMBDA(const size_t& i){
202 1520 : data_m[i] = view(i);
203 : }
204 : );
205 76 : return *this;
206 76 : }
207 :
208 :
209 : template <typename T>
210 3300 : T FEMVector<T>::operator[] (size_t i) const {
211 6600 : return data_m(i);
212 : }
213 :
214 :
215 : template <typename T>
216 3300 : T FEMVector<T>::operator() (size_t i) const {
217 3300 : return this->operator[](i);
218 : }
219 :
220 :
221 : template <typename T>
222 436 : const Kokkos::View<T*>& FEMVector<T>::getView() const {
223 436 : return data_m;
224 : }
225 :
226 :
227 : template <typename T>
228 78 : size_t FEMVector<T>::size() const {
229 78 : return data_m.extent(0);
230 : }
231 :
232 : template <typename T>
233 20 : FEMVector<T> FEMVector<T>::deepCopy() const {
234 : // We have to check if we have boundary information or not
235 20 : if (boundaryInfo_m) {
236 : // The neighbor_m can be simply passed to the new vector, the sendIdxs_m
237 : // and recvIdxs_m need to be explicitly copied.
238 2 : std::vector< Kokkos::View<size_t*> > newSendIdxs;
239 2 : std::vector< Kokkos::View<size_t*> > newRecvIdxs;
240 :
241 6 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
242 4 : newSendIdxs.emplace_back(Kokkos::View<size_t*>(boundaryInfo_m->sendIdxs_m[i].label(),
243 4 : boundaryInfo_m->sendIdxs_m[i].extent(0)));
244 4 : Kokkos::deep_copy(newSendIdxs[i], boundaryInfo_m->sendIdxs_m[i]);
245 :
246 4 : newRecvIdxs.emplace_back(Kokkos::View<size_t*>(boundaryInfo_m->recvIdxs_m[i].label(),
247 4 : boundaryInfo_m->recvIdxs_m[i].extent(0)));
248 :
249 4 : Kokkos::deep_copy(newRecvIdxs[i], boundaryInfo_m->recvIdxs_m[i]);
250 : }
251 :
252 : // create the new FEMVector
253 2 : FEMVector<T> newVector(size(), boundaryInfo_m->neighbors_m, newSendIdxs, newRecvIdxs);
254 : // copy over the values
255 2 : newVector = *this;
256 :
257 2 : return newVector;
258 0 : } else {
259 18 : FEMVector<T> newVector(size());
260 : // copy over the values
261 18 : newVector = *this;
262 :
263 18 : return newVector;
264 18 : }
265 : }
266 :
267 : template <typename T>
268 : template <typename K>
269 0 : FEMVector<K> FEMVector<T>::skeletonCopy() const {
270 : // We have to check if we have boundary information or not
271 0 : if (boundaryInfo_m) {
272 : // The neighbor_m can be simply passed to the new vector, the sendIdxs_m
273 : // and recvIdxs_m need to be explicitly copied.
274 0 : std::vector< Kokkos::View<size_t*> > newSendIdxs;
275 0 : std::vector< Kokkos::View<size_t*> > newRecvIdxs;
276 :
277 0 : for (size_t i = 0; i < boundaryInfo_m->neighbors_m.size(); ++i) {
278 0 : newSendIdxs.emplace_back(Kokkos::View<size_t*>(boundaryInfo_m->sendIdxs_m[i].label(),
279 0 : boundaryInfo_m->sendIdxs_m[i].extent(0)));
280 0 : Kokkos::deep_copy(newSendIdxs[i], boundaryInfo_m->sendIdxs_m[i]);
281 :
282 0 : newRecvIdxs.emplace_back(Kokkos::View<size_t*>(boundaryInfo_m->recvIdxs_m[i].label(),
283 0 : boundaryInfo_m->recvIdxs_m[i].extent(0)));
284 :
285 0 : Kokkos::deep_copy(newRecvIdxs[i], boundaryInfo_m->recvIdxs_m[i]);
286 : }
287 :
288 : // create the new FEMVector
289 0 : FEMVector<K> newVector(size(), boundaryInfo_m->neighbors_m, newSendIdxs, newRecvIdxs);
290 :
291 0 : return newVector;
292 0 : } else {
293 0 : FEMVector<K> newVector(size());
294 :
295 0 : return newVector;
296 0 : }
297 : }
298 :
299 :
300 : template <typename T>
301 24 : void FEMVector<T>::pack(const Kokkos::View<size_t*>& idxStore) {
302 : // check that we have halo information
303 24 : if (!boundaryInfo_m) {
304 0 : throw IpplException(
305 : "FEMVector::pack()",
306 : "Cannot do halo operations, as no MPI communication information is provided. "
307 : "Did you use the correct constructor to construct the FEMVector?");
308 : }
309 :
310 24 : size_t nIdxs = idxStore.extent(0);
311 24 : auto& bufferData = boundaryInfo_m->commBuffer_m.buffer;
312 :
313 24 : if (bufferData.size() < nIdxs) {
314 12 : int overalloc = Comm->getDefaultOverallocation();
315 12 : Kokkos::realloc(bufferData, nIdxs * overalloc);
316 : }
317 :
318 24 : Kokkos::parallel_for("FEMVector::pack()", nIdxs,
319 688 : KOKKOS_CLASS_LAMBDA(const size_t& i) {
320 1008 : bufferData(i) = data_m(idxStore(i));
321 : }
322 : );
323 24 : Kokkos::fence();
324 24 : }
325 :
326 :
327 : template <typename T>
328 : template <typename Op>
329 24 : void FEMVector<T>::unpack(const Kokkos::View<size_t*>& idxStore) {
330 : // check that we have halo information
331 24 : if (!boundaryInfo_m) {
332 0 : throw IpplException(
333 : "FEMVector::unpack()",
334 : "Cannot do halo operations, as no MPI communication information is provided. "
335 : "Did you use the correct constructor to construct the FEMVector?");
336 : }
337 :
338 24 : size_t nIdxs = idxStore.extent(0);
339 24 : auto& bufferData = boundaryInfo_m->commBuffer_m.buffer;
340 24 : if (bufferData.size() < nIdxs) {
341 0 : int overalloc = Comm->getDefaultOverallocation();
342 0 : Kokkos::realloc(bufferData, nIdxs * overalloc);
343 : }
344 :
345 : Op op;
346 24 : Kokkos::parallel_for("FEMVector::unpack()", nIdxs,
347 720 : KOKKOS_CLASS_LAMBDA(const size_t& i) {
348 1344 : op(data_m(idxStore(i)), bufferData(i));
349 : }
350 : );
351 24 : Kokkos::fence();
352 24 : }
353 :
354 :
355 :
356 : } // namespace ippl
|
