dune-grid 2.9.0
torus.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_GRID_YASPGRID_TORUS_HH
6#define DUNE_GRID_YASPGRID_TORUS_HH
7
8#include <array>
9#include <bitset>
10#include <cmath>
11#include <deque>
12#include <iostream>
13#include <vector>
14
15#if HAVE_MPI
16#include <mpi.h>
17#endif
18
19#include <dune/common/binaryfunctions.hh>
20#include <dune/common/streamoperators.hh>
22
23#include "partitioning.hh"
24
29namespace Dune
30{
31
45 template<class Communication, int d>
46 class Torus {
47 public:
49 typedef std::array<int, d> iTupel;
50
51
52 private:
53 struct CommPartner {
54 int rank;
55 iTupel delta;
56 int index;
57 };
58
59 struct CommTask {
60 int rank; // process to send to / receive from
61 int size; // size of buffer
62 void *buffer; // buffer to send / receive
63 };
64
65 public:
68 {}
69
71 Torus (Communication comm, int tag, iTupel size, int overlap, const Yasp::Partitioning<d>* partitioner)
72 : _comm(comm), _tag(tag)
73 {
74 // determine dimensions
75 partitioner->partition(size, _comm.size(), _dims, overlap);
76
77 // compute increments for lexicographic ordering
78 int inc = 1;
79 for (int i=0; i<d; i++)
80 {
81 _increment[i] = inc;
82 inc *= _dims[i];
83 }
84
85 // check whether the load balancing matches the size of the communicator
86 if (inc != _comm.size())
87 DUNE_THROW(Dune::Exception, "Communicator size and result of the given load balancer do not match!");
88
89 // make full schedule
90 proclists();
91 }
92
94 int rank () const
95 {
96 return _comm.rank();
97 }
98
100 iTupel coord () const
101 {
102 return rank_to_coord(_comm.rank());
103 }
104
106 int procs () const
107 {
108 return _comm.size();
109 }
110
112 const iTupel & dims () const
113 {
114 return _dims;
115 }
116
118 int dims (int i) const
119 {
120 return _dims[i];
121 }
122
124 Communication comm () const
125 {
126 return _comm;
127 }
128
130 int tag () const
131 {
132 return _tag;
133 }
134
136 bool inside (iTupel c) const
137 {
138 for (int i=d-1; i>=0; i--)
139 if (c[i]<0 || c[i]>=_dims[i]) return false;
140 return true;
141 }
142
145 {
147 rank = rank%_comm.size();
148 for (int i=d-1; i>=0; i--)
149 {
150 coord[i] = rank/_increment[i];
151 rank = rank%_increment[i];
152 }
153 return coord;
154 }
155
158 {
159 for (int i=0; i<d; i++) coord[i] = coord[i]%_dims[i];
160 int rank = 0;
161 for (int i=0; i<d; i++) rank += coord[i]*_increment[i];
162 return rank;
163 }
164
166 int rank_relative (int rank, int dir, int cnt) const
167 {
169 coord[dir] = (coord[dir]+_dims[dir]+cnt)%_dims[dir];
170 return coord_to_rank(coord);
171 }
172
174 int color (const iTupel & coord) const
175 {
176 int c = 0;
177 int power = 1;
178
179 // interior coloring
180 for (int i=0; i<d; i++)
181 {
182 if (coord[i]%2==1) c += power;
183 power *= 2;
184 }
185
186 // extra colors for boundary processes
187 for (int i=0; i<d; i++)
188 {
189 if (_dims[i]>1 && coord[i]==_dims[i]-1) c += power;
190 power *= 2;
191 }
192
193 return c;
194 }
195
197 int color (int rank) const
198 {
199 return color(rank_to_coord(rank));
200 }
201
203 int neighbors () const
204 {
205 int n=1;
206 for (int i=0; i<d; ++i)
207 n *= 3;
208 return n-1;
209 }
210
212 bool is_neighbor (iTupel delta, std::bitset<d> periodic) const
213 {
214 iTupel coord = rank_to_coord(_comm.rank()); // my own coordinate with 0 <= c_i < dims_i
215
216 for (int i=0; i<d; i++)
217 {
218 if (delta[i]<0)
219 {
220 // if I am on the boundary and domain is not periodic => no neighbor
221 if (coord[i]==0 && periodic[i]==false) return false;
222 }
223 if (delta[i]>0)
224 {
225 // if I am on the boundary and domain is not periodic => no neighbor
226 if (coord[i]==_dims[i]-1 && periodic[i]==false) return false;
227 }
228 }
229 return true;
230 }
231
239 double partition (int rank, iTupel origin_in, iTupel size_in, iTupel& origin_out, iTupel& size_out) const
240 {
242 double maxsize = 1;
243 double sz = 1;
244
245 // make a tensor product partition
246 for (int i=0; i<d; i++)
247 {
248 // determine
249 int m = size_in[i]/_dims[i];
250 int r = size_in[i]%_dims[i];
251
252 sz *= size_in[i];
253
254 if (coord[i]<_dims[i]-r)
255 {
256 origin_out[i] = origin_in[i] + coord[i]*m;
257 size_out[i] = m;
258 maxsize *= m;
259 }
260 else
261 {
262 origin_out[i] = origin_in[i] + (_dims[i]-r)*m + (coord[i]-(_dims[i]-r))*(m+1);
263 size_out[i] = m+1;
264 maxsize *= m+1;
265 }
266 }
267 return maxsize/(sz/_comm.size());
268 }
269
277 public:
279 ProcListIterator (typename std::deque<CommPartner>::const_iterator iter)
280 {
281 i = iter;
282 }
283
285 int rank () const
286 {
287 return i->rank;
288 }
289
291 iTupel delta () const
292 {
293 return i->delta;
294 }
295
297 int index () const
298 {
299 return i->index;
300 }
301
303 int distance () const
304 {
305 int dist = 0;
306 iTupel delta=i->delta;
307 for (int j=0; j<d; ++j)
308 dist += std::abs(delta[j]);
309 return dist;
310 }
311
313 bool operator== (const ProcListIterator& iter) const
314 {
315 return i == iter.i;
316 }
317
318
320 bool operator!= (const ProcListIterator& iter) const
321 {
322 return i != iter.i;
323 }
324
327 {
328 ++i;
329 return *this;
330 }
331
332 private:
333 typename std::deque<CommPartner>::const_iterator i;
334 };
335
338 {
339 return ProcListIterator(_sendlist.begin());
340 }
341
344 {
345 return ProcListIterator(_sendlist.end());
346 }
347
350 {
351 return ProcListIterator(_recvlist.begin());
352 }
353
356 {
357 return ProcListIterator(_recvlist.end());
358 }
359
361 void send (int rank, void* buffer, int size) const
362 {
363 CommTask task;
364 task.rank = rank;
365 task.buffer = buffer;
366 task.size = size;
367 if (rank!=_comm.rank())
368 _sendrequests.push_back(task);
369 else
370 _localsendrequests.push_back(task);
371 }
372
374 void recv (int rank, void* buffer, int size) const
375 {
376 CommTask task;
377 task.rank = rank;
378 task.buffer = buffer;
379 task.size = size;
380 if (rank!=_comm.rank())
381 _recvrequests.push_back(task);
382 else
383 _localrecvrequests.push_back(task);
384 }
385
387 void exchange () const
388 {
389 // handle local requests first
390 if (_localsendrequests.size()!=_localrecvrequests.size())
391 {
392 std::cout << "[" << rank() << "]: ERROR: local sends/receives do not match in exchange!" << std::endl;
393 return;
394 }
395 for (unsigned int i=0; i<_localsendrequests.size(); i++)
396 {
397 if (_localsendrequests[i].size!=_localrecvrequests[i].size)
398 {
399 std::cout << "[" << rank() << "]: ERROR: size in local sends/receive does not match in exchange!" << std::endl;
400 return;
401 }
402 memcpy(_localrecvrequests[i].buffer,_localsendrequests[i].buffer,_localsendrequests[i].size);
403 }
404 _localsendrequests.clear();
405 _localrecvrequests.clear();
406
407#if HAVE_MPI
408 // handle foreign requests
409
410 std::vector<MPI_Request> requests(_sendrequests.size() + _recvrequests.size());
411 MPI_Request* req = requests.data();
412
413 // issue sends to foreign processes
414 for (unsigned int i=0; i<_sendrequests.size(); i++)
415 if (_sendrequests[i].rank!=rank())
416 {
417 // std::cout << "[" << rank() << "]" << " send " << _sendrequests[i].size << " bytes "
418 // << "to " << _sendrequests[i].rank << " p=" << _sendrequests[i].buffer << std::endl;
419 MPI_Isend(_sendrequests[i].buffer, _sendrequests[i].size, MPI_BYTE,
420 _sendrequests[i].rank, _tag, _comm, req++);
421 }
422
423 // issue receives from foreign processes
424 for (unsigned int i=0; i<_recvrequests.size(); i++)
425 if (_recvrequests[i].rank!=rank())
426 {
427 // std::cout << "[" << rank() << "]" << " recv " << _recvrequests[i].size << " bytes "
428 // << "fm " << _recvrequests[i].rank << " p=" << _recvrequests[i].buffer << std::endl;
429 MPI_Irecv(_recvrequests[i].buffer, _recvrequests[i].size, MPI_BYTE,
430 _recvrequests[i].rank, _tag, _comm, req++);
431 }
432
433 // Wait for communication to complete
434 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
435
436 // clear request buffers
437 _sendrequests.clear();
438 _recvrequests.clear();
439#endif
440 }
441
443 double global_max (double x) const
444 {
445 double res = 0.0;
446 _comm.template allreduce<Dune::Max<double>,double>(&x, &res, 1);
447 return res;
448 }
449
451 void print (std::ostream& s) const
452 {
453 s << "[" << rank() << "]: Torus " << procs() << " processor(s) arranged as " << dims() << std::endl;
454 for (ProcListIterator i=sendbegin(); i!=sendend(); ++i)
455 {
456 s << "[" << rank() << "]: send to "
457 << "rank=" << i.rank()
458 << " index=" << i.index()
459 << " delta=" << i.delta() << " dist=" << i.distance() << std::endl;
460 }
461 for (ProcListIterator i=recvbegin(); i!=recvend(); ++i)
462 {
463 s << "[" << rank() << "]: recv from "
464 << "rank=" << i.rank()
465 << " index=" << i.index()
466 << " delta=" << i.delta() << " dist=" << i.distance() << std::endl;
467 }
468 }
469
470 private:
471
472 void proclists ()
473 {
474 // compile the full neighbor list
475 CommPartner cp;
476 iTupel delta;
477
478 std::fill(delta.begin(), delta.end(), -1);
479 bool ready = false;
480 iTupel me, nb;
481 me = rank_to_coord(_comm.rank());
482 int index = 0;
483 int last = neighbors()-1;
484 while (!ready)
485 {
486 // find neighbors coordinates
487 for (int i=0; i<d; i++)
488 nb[i] = ( me[i]+_dims[i]+delta[i] ) % _dims[i];
489
490 // find neighbors rank
491 int nbrank = coord_to_rank(nb);
492
493 // check if delta is not zero
494 for (int i=0; i<d; i++)
495 if (delta[i]!=0)
496 {
497 cp.rank = nbrank;
498 cp.delta = delta;
499 cp.index = index;
500 _recvlist.push_back(cp);
501 cp.index = last-index;
502 _sendlist.push_front(cp);
503 index++;
504 break;
505 }
506
507 // next neighbor
508 ready = true;
509 for (int i=0; i<d; i++)
510 if (delta[i]<1)
511 {
512 (delta[i])++;
513 ready=false;
514 break;
515 }
516 else
517 {
518 delta[i] = -1;
519 }
520 }
521
522 }
523
524 Communication _comm;
525
526 iTupel _dims;
527 iTupel _increment;
528 int _tag;
529 std::deque<CommPartner> _sendlist;
530 std::deque<CommPartner> _recvlist;
531
532 mutable std::vector<CommTask> _sendrequests;
533 mutable std::vector<CommTask> _recvrequests;
534 mutable std::vector<CommTask> _localsendrequests;
535 mutable std::vector<CommTask> _localrecvrequests;
536
537 };
538
540 template <class Communication, int d>
541 inline std::ostream& operator<< (std::ostream& s, const Torus<Communication, d> & t)
542 {
543 t.print(s);
544 return s;
545 }
546}
547
548#endif
This file provides tools to partition YaspGrids. If you want to write your own partitioner,...
std::ostream & operator<<(std::ostream &out, const PartitionType &type)
write a PartitionType to a stream
Definition: gridenums.hh:72
Include standard header files.
Definition: agrid.hh:60
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:328
constexpr Overlap overlap
PartitionSet for the overlap partition.
Definition: partitionset.hh:278
a base class for the yaspgrid partitioning strategy
Definition: partitioning.hh:39
virtual void partition(const iTupel &, int, iTupel &, int) const =0
Definition: torus.hh:46
Torus()
constructor making uninitialized object
Definition: torus.hh:67
int color(int rank) const
assign color to given rank
Definition: torus.hh:197
double partition(int rank, iTupel origin_in, iTupel size_in, iTupel &origin_out, iTupel &size_out) const
partition the given grid onto the torus and return the piece of the process with given rank; returns ...
Definition: torus.hh:239
int dims(int i) const
return dimensions of torus in direction i
Definition: torus.hh:118
iTupel coord() const
return own coordinates
Definition: torus.hh:100
int rank() const
return own rank
Definition: torus.hh:94
const iTupel & dims() const
return dimensions of torus
Definition: torus.hh:112
int rank_relative(int rank, int dir, int cnt) const
return rank of process where its coordinate in direction dir has offset cnt (handles periodic case)
Definition: torus.hh:166
void recv(int rank, void *buffer, int size) const
store a receive request; buffers are received in order; handles also local requests with memcpy
Definition: torus.hh:374
void send(int rank, void *buffer, int size) const
store a send request; buffers are sent in order; handles also local requests with memcpy
Definition: torus.hh:361
Torus(Communication comm, int tag, iTupel size, int overlap, const Yasp::Partitioning< d > *partitioner)
make partitioner from communicator and coarse mesh size
Definition: torus.hh:71
int neighbors() const
return the number of neighbors, which is
Definition: torus.hh:203
void print(std::ostream &s) const
print contents of torus object
Definition: torus.hh:451
double global_max(double x) const
global max
Definition: torus.hh:443
Communication comm() const
return communicator
Definition: torus.hh:124
int color(const iTupel &coord) const
assign color to given coordinate
Definition: torus.hh:174
ProcListIterator recvend() const
last process in receive list
Definition: torus.hh:355
bool is_neighbor(iTupel delta, std::bitset< d > periodic) const
return true if neighbor with given delta is a neighbor under the given periodicity
Definition: torus.hh:212
ProcListIterator sendend() const
end of send list
Definition: torus.hh:343
int procs() const
return number of processes
Definition: torus.hh:106
iTupel rank_to_coord(int rank) const
map rank to coordinate in torus using lexicographic ordering
Definition: torus.hh:144
std::array< int, d > iTupel
type used to pass tupels in and out
Definition: torus.hh:49
ProcListIterator sendbegin() const
first process in send list
Definition: torus.hh:337
void exchange() const
exchange messages stored in request buffers; clear request buffers afterwards
Definition: torus.hh:387
ProcListIterator recvbegin() const
first process in receive list
Definition: torus.hh:349
int coord_to_rank(iTupel coord) const
map coordinate in torus to rank using lexicographic ordering
Definition: torus.hh:157
bool inside(iTupel c) const
return true if coordinate is inside torus
Definition: torus.hh:136
int tag() const
return tag used by torus
Definition: torus.hh:130
Definition: torus.hh:276
iTupel delta() const
return distance vector
Definition: torus.hh:291
bool operator==(const ProcListIterator &iter) const
Return true when two iterators point to same member.
Definition: torus.hh:313
bool operator!=(const ProcListIterator &iter) const
Return true when two iterators do not point to same member.
Definition: torus.hh:320
int rank() const
return rank of neighboring process
Definition: torus.hh:285
ProcListIterator(typename std::deque< CommPartner >::const_iterator iter)
make an iterator
Definition: torus.hh:279
ProcListIterator & operator++()
Increment iterator to next cell.
Definition: torus.hh:326
int index() const
return index in proclist
Definition: torus.hh:297
int distance() const
return 1-norm of distance vector
Definition: torus.hh:303