ALL 0.9.3
A Loadbalacing Library
Loading...
Searching...
No Matches
ALL_ForceBased.hpp
Go to the documentation of this file.
1/*
2Copyright 2018-2020 Rene Halver, Forschungszentrum Juelich GmbH, Germany
3Copyright 2018-2020 Godehard Sutmann, Forschungszentrum Juelich GmbH, Germany
4
5Redistribution and use in source and binary forms, with or without modification,
6are permitted provided that the following conditions are met:
7
81. Redistributions of source code must retain the above copyright notice, this
9 list of conditions and the following disclaimer.
10
112. Redistributions in binary form must reproduce the above copyright notice,
12this list of conditions and the following disclaimer in the documentation and/or
13 other materials provided with the distribution.
14
153. Neither the name of the copyright holder nor the names of its contributors
16may be used to endorse or promote products derived from this software without
17 specific prior written permission.
18
19THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
20ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
23ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
26ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29*/
30
31#ifndef ALL_FORCEBASED_HEADER_INCLUDED
32#define ALL_FORCEBASED_HEADER_INCLUDED
33
34#define ALL_FORCE_OLD
35
36/*
37
38 ForceBased load-balancing scheme:
39
40 Requirements:
41 a) cartesian communicator
42 b) 2^(dim) vertices describing the domain (in theory could be more,
43 we assume a orthogonal grid as start point)
44 c) MPI >= 3.0 due to necessary non-blocking collectives
45
46
47 Method:
48
49 For each of the vertices of a domain a communicator is created, which
50 contains the surrounding processes. For one of the vertices, a domain
51 collects the centres of gravitiy and the work of the neighboring domains. In
52 a next step a force acting on the vertex is computed as: F = 1.0 / ( gamma *
53 n * sum(W_i) ) * sum( W_i * x_i ) with:
54
55 n: number of neighboring domains
56 W_i: work on domain i
57 x_i: vector pointing from the vertex to
58 the center of gravity of domain i
59 gamma: correction factor to control the
60 speed of the vertex shift
61
62 */
63
65#include "ALL_LB.hpp"
66#include "ALL_Point.hpp"
67#include "ALL_Defines.h"
68#include <algorithm>
69#include <cmath>
70#include <exception>
71#include <mpi.h>
72#include <vector>
73
74namespace ALL {
75
85template <class T, class W> class ForceBased_LB : public LB<T, W> {
86public:
93 ForceBased_LB(int d, W w, T g) : LB<T, W>(d, g) {
94 this->setWork(w);
95 }
96
98 ~ForceBased_LB() override;
99
101 void setup() override;
102
105 virtual void balance(int step) override;
106
107 // getter for variables (passed by reference to avoid
108 // direct access to private members of the object)
109
113 virtual std::vector<int> &getNeighbors() override;
114
117 virtual void setAdditionalData(known_unused const void *data) override {}
118
122 virtual W getEstimatedEfficiency() override {return (W)-1;}
123
124private:
125 // type for MPI communication
126 MPI_Datatype MPIDataTypeT;
127 MPI_Datatype MPIDataTypeW;
128
129 // array of MPI communicators for each vertex
130 std::vector<T> vertex_neighbors;
131
132 // MPI communicator to collect the work in the
133 // main dimension
134 MPI_Comm main_communicator;
135
136 // list of neighbors
137 std::vector<int> neighbors;
138
139 // find neighbors (more sophisticated than for tensor-product
140 // variant of LB)
141 void find_neighbors();
142
143 // number of vertices (since it is not equal for all domains)
144 int n_vertices;
145
146 // main dimension (correction in staggered style)
147 int main_dim;
148
149 // secondary dimensions (2D force shift)
150 int sec_dim[2];
151};
152
153template <class T, class W> ForceBased_LB<T, W>::~ForceBased_LB() {}
154
155// provide list of neighbors
156template <class T, class W>
158 return neighbors;
159}
160
161template <class T, class W> void ForceBased_LB<T, W>::setup() {
162 n_vertices = this->vertices.size();
163 vertex_neighbors.resize(n_vertices * 8);
164
165 // determine correct MPI data type for template T
166 if (std::is_same<T, double>::value)
167 MPIDataTypeT = MPI_DOUBLE;
168 else if (std::is_same<T, float>::value)
169 MPIDataTypeT = MPI_FLOAT;
170 else if (std::is_same<T, int>::value)
171 MPIDataTypeT = MPI_INT;
172 else if (std::is_same<T, long>::value)
173 MPIDataTypeT = MPI_LONG;
174 else {
176 __FILE__, __func__, __LINE__,
177 "Invalid data type for boundaries given (T)");
178 }
179
180 // determine correct MPI data type for template W
181 if (std::is_same<W, double>::value)
182 MPIDataTypeW = MPI_DOUBLE;
183 else if (std::is_same<W, float>::value)
184 MPIDataTypeW = MPI_FLOAT;
185 else if (std::is_same<W, int>::value)
186 MPIDataTypeW = MPI_INT;
187 else if (std::is_same<W, long>::value)
188 MPIDataTypeW = MPI_LONG;
189 else {
190 throw InvalidCommTypeException(__FILE__, __func__, __LINE__,
191 "Invalid data type for work given (W)");
192 }
193
194 MPI_Group all;
195 MPI_Comm_group(this->globalComm, &all);
196
197 // check if communicator is cartesian
198 int status;
199 MPI_Topo_test(this->globalComm, &status);
200
201 if (status != MPI_CART) {
203 __FILE__, __func__, __LINE__,
204 "Cartesian MPI communicator required, passed communicator is not "
205 "cartesian");
206 }
207
208 // get the local coordinates, this->periodicity and global size from the MPI
209 // communicator
210 MPI_Cart_get(this->globalComm, this->dimension, this->global_dims.data(),
211 this->periodicity.data(), this->local_coords.data());
212
213 // get the local rank from the MPI communicator
214 MPI_Cart_rank(this->globalComm, this->local_coords.data(), &this->localRank);
215
216 // groups required for new communicators
217 MPI_Group known_unused groups[n_vertices];
218 // arrays of processes belonging to group
219 int known_unused processes[n_vertices][n_vertices];
220
221 // shifted local coordinates to find neighboring processes
222 int known_unused shifted_coords[this->dimension];
223
224 // get main and secondary dimensions
225 if (this->global_dims.at(2) >= this->global_dims.at(1)) {
226 if (this->global_dims.at(2) >= this->global_dims.at(0)) {
227 main_dim = 2;
228 sec_dim[0] = 0;
229 sec_dim[1] = 1;
230 } else {
231 main_dim = 0;
232 sec_dim[0] = 1;
233 sec_dim[1] = 2;
234 }
235 } else {
236 if (this->global_dims.at(1) >= this->global_dims.at(0)) {
237 main_dim = 1;
238 sec_dim[0] = 0;
239 sec_dim[1] = 2;
240 } else {
241 main_dim = 0;
242 sec_dim[0] = 1;
243 sec_dim[1] = 2;
244 }
245 }
246
247 if (this->localRank == 0)
248 std::cout << "DEBUG: main_dim: " << main_dim << std::endl;
249
250 // create main communicator
251 MPI_Comm_split(this->globalComm, this->local_coords.at(main_dim),
252 this->local_coords.at(sec_dim[0]) +
253 this->local_coords.at(sec_dim[1]) *
254 this->global_dims.at(sec_dim[0]),
255 &main_communicator);
256
257 /*
258
259 // sequence of vertices in 1d
260 //
261 // 0 - - - 1
262
263 // sequence of vertices in 2d
264 //
265 // 2 - - - 3
266 // | |
267 // | |
268 // 0 - - - 1
269
270 // sequence of vertices in 3d
271 //
272 // 6 - - - 7
273 // /| /|
274 // / | / |
275 // 4 -2- - 5 3
276 // | / | /
277 // |/ |/
278 // 0 - - - 1
279
280 */
281
282#ifdef ALL_DEBUG_ENABLED
283 MPI_Barrier(this->globalComm);
284 if (this->localRank == 0)
285 std::cout << "ALL::ForceBased_LB<T,W>::setup() preparing communicators..."
286 << std::endl;
287#endif
288 std::vector<int> dim_vert(this->global_dims);
289
290 MPI_Comm known_unused tmp_comm;
291 int known_unused own_vertex;
292
293#ifdef ALL_DEBUG_ENABLED
294 MPI_Barrier(this->globalComm);
295 if (this->localRank == 0)
296 std::cout << "ALL::ForceBased_LB<T,W>::setup() computing communicators..."
297 << std::endl;
298 std::cout << "DEBUG: "
299 << " rank: " << this->localRank << " dim_vert: " << dim_vert.at(0)
300 << " " << dim_vert.at(1) << " " << dim_vert.at(2)
301 << " size(local_coords): " << this->local_coords.size() << " "
302 << " size(global_dims): " << this->global_dims.size() << std::endl;
303#endif
304 for (int iz = 0; iz < dim_vert.at(2); ++iz) {
305 for (int iy = 0; iy < dim_vert.at(1); ++iy) {
306 for (int ix = 0; ix < dim_vert.at(0); ++ix) {
307 bool affected[8];
308 for (auto &a : affected)
309 a = false;
310 int v_neighbors[8];
311 for (auto &vn : vertex_neighbors)
312 vn = -1;
313 if (ix == ((this->local_coords.at(0) + 0) % dim_vert.at(0)) &&
314 iy == ((this->local_coords.at(1) + 0) % dim_vert.at(1)) &&
315 iz == ((this->local_coords.at(2) + 0) % dim_vert.at(2))) {
316 affected[0] = true;
317 v_neighbors[0] = this->localRank;
318 }
319 if (ix == ((this->local_coords.at(0) + 1) % dim_vert.at(0)) &&
320 iy == ((this->local_coords.at(1) + 0) % dim_vert.at(1)) &&
321 iz == ((this->local_coords.at(2) + 0) % dim_vert.at(2))) {
322 affected[1] = true;
323 v_neighbors[1] = this->localRank;
324 }
325 if (ix == ((this->local_coords.at(0) + 0) % dim_vert.at(0)) &&
326 iy == ((this->local_coords.at(1) + 1) % dim_vert.at(1)) &&
327 iz == ((this->local_coords.at(2) + 0) % dim_vert.at(2))) {
328 affected[2] = true;
329 v_neighbors[2] = this->localRank;
330 }
331 if (ix == ((this->local_coords.at(0) + 1) % dim_vert.at(0)) &&
332 iy == ((this->local_coords.at(1) + 1) % dim_vert.at(1)) &&
333 iz == ((this->local_coords.at(2) + 0) % dim_vert.at(2))) {
334 affected[3] = true;
335 v_neighbors[3] = this->localRank;
336 }
337 if (ix == ((this->local_coords.at(0) + 0) % dim_vert.at(0)) &&
338 iy == ((this->local_coords.at(1) + 0) % dim_vert.at(1)) &&
339 iz == ((this->local_coords.at(2) + 1) % dim_vert.at(2))) {
340 affected[4] = true;
341 v_neighbors[4] = this->localRank;
342 }
343 if (ix == ((this->local_coords.at(0) + 1) % dim_vert.at(0)) &&
344 iy == ((this->local_coords.at(1) + 0) % dim_vert.at(1)) &&
345 iz == ((this->local_coords.at(2) + 1) % dim_vert.at(2))) {
346 affected[5] = true;
347 v_neighbors[5] = this->localRank;
348 }
349 if (ix == ((this->local_coords.at(0) + 0) % dim_vert.at(0)) &&
350 iy == ((this->local_coords.at(1) + 1) % dim_vert.at(1)) &&
351 iz == ((this->local_coords.at(2) + 1) % dim_vert.at(2))) {
352 affected[6] = true;
353 v_neighbors[6] = this->localRank;
354 }
355 if (ix == ((this->local_coords.at(0) + 1) % dim_vert.at(0)) &&
356 iy == ((this->local_coords.at(1) + 1) % dim_vert.at(1)) &&
357 iz == ((this->local_coords.at(2) + 1) % dim_vert.at(2))) {
358 affected[7] = true;
359 v_neighbors[7] = this->localRank;
360 }
361
362 MPI_Allreduce(MPI_IN_PLACE, v_neighbors, 8, MPI_INT, MPI_MAX,
363 this->globalComm);
364
365 for (int v = 0; v < n_vertices; ++v) {
366 if (affected[v]) {
367 for (int n = 0; n < 8; ++n) {
368 vertex_neighbors.at(8 * v + n) = v_neighbors[n];
369 }
370 }
371 }
372 }
373 }
374// ToDo: check if vertices correct
375#ifdef ALL_DEBUG_ENABLED
376 MPI_Barrier(this->globalComm);
377 if (this->localRank == 0)
378 std::cout << "ALL::ForceBased_LB<T,W>::setup() finished computing "
379 "communicators..."
380 << std::endl;
381#endif
382 }
383}
384
385// TODO: periodic boundary conditions (would require size of the system)
386// for now: do not change the vertices along the edge of the system
387
388template <class T, class W>
390
391 this->prevVertices = this->vertices;
392
393 // geometrical center and work of each neighbor for each vertex
394 // work is cast to vertex data type, since in the end a shift
395 // of the vertex is computed
396 T vertex_info[n_vertices * n_vertices * (this->dimension + 2)];
397 T center[this->dimension];
398
399 // local geometric center
400 T known_unused local_info[(this->dimension + 2) * n_vertices];
401
402 for (int i = 0; i < n_vertices * n_vertices * (this->dimension + 2); ++i)
403 vertex_info[i] = -1;
404
405 for (int d = 0; d < (this->dimension + 2) * n_vertices; ++d)
406 local_info[d] = (T)0;
407
408 for (int d = 0; d < this->dimension; ++d)
409 center[d] = (T)0;
410
411 // compute geometric center
412 for (int v = 0; v < n_vertices; ++v) {
413 for (int d = 0; d < this->dimension; ++d) {
414 center[d] += ((this->prevVertices.at(v))[d] / (T)n_vertices);
415 }
416 local_info[v * (this->dimension + 2) + this->dimension] =
417 (T)this->work.at(0);
418 }
419 for (int v = 0; v < n_vertices; ++v) {
420 for (int d = 0; d < this->dimension; ++d) {
421 // vector pointing to center
422 local_info[v * (this->dimension + 2) + d] =
423 center[d] - (this->prevVertices.at(v))[d];
424 }
425 }
426
427 // compute for each vertex the maximum movement
428 switch (main_dim) {
429 case 0:
430 local_info[0 * (this->dimension + 2) + this->dimension + 1] =
431 0.5 * std::min(this->prevVertices.at(0).d(this->prevVertices.at(2)),
432 this->prevVertices.at(0).d(this->prevVertices.at(4)));
433 local_info[1 * (this->dimension + 2) + this->dimension + 1] =
434 0.5 * std::min(this->prevVertices.at(1).d(this->prevVertices.at(3)),
435 this->prevVertices.at(1).d(this->prevVertices.at(5)));
436 local_info[2 * (this->dimension + 2) + this->dimension + 1] =
437 0.5 * std::min(this->prevVertices.at(2).d(this->prevVertices.at(0)),
438 this->prevVertices.at(2).d(this->prevVertices.at(6)));
439 local_info[3 * (this->dimension + 2) + this->dimension + 1] =
440 0.5 * std::min(this->prevVertices.at(3).d(this->prevVertices.at(1)),
441 this->prevVertices.at(3).d(this->prevVertices.at(7)));
442 local_info[4 * (this->dimension + 2) + this->dimension + 1] =
443 0.5 * std::min(this->prevVertices.at(4).d(this->prevVertices.at(0)),
444 this->prevVertices.at(4).d(this->prevVertices.at(6)));
445 local_info[5 * (this->dimension + 2) + this->dimension + 1] =
446 0.5 * std::min(this->prevVertices.at(5).d(this->prevVertices.at(1)),
447 this->prevVertices.at(5).d(this->prevVertices.at(7)));
448 local_info[6 * (this->dimension + 2) + this->dimension + 1] =
449 0.5 * std::min(this->prevVertices.at(6).d(this->prevVertices.at(2)),
450 this->prevVertices.at(6).d(this->prevVertices.at(4)));
451 local_info[7 * (this->dimension + 2) + this->dimension + 1] =
452 0.5 * std::min(this->prevVertices.at(7).d(this->prevVertices.at(3)),
453 this->prevVertices.at(7).d(this->prevVertices.at(5)));
454 break;
455 case 1:
456 local_info[0 * (this->dimension + 2) + this->dimension + 1] =
457 0.5 * std::min(this->prevVertices.at(0).d(this->prevVertices.at(1)),
458 this->prevVertices.at(0).d(this->prevVertices.at(4)));
459 local_info[1 * (this->dimension + 2) + this->dimension + 1] =
460 0.5 * std::min(this->prevVertices.at(1).d(this->prevVertices.at(0)),
461 this->prevVertices.at(1).d(this->prevVertices.at(5)));
462 local_info[2 * (this->dimension + 2) + this->dimension + 1] =
463 0.5 * std::min(this->prevVertices.at(2).d(this->prevVertices.at(3)),
464 this->prevVertices.at(2).d(this->prevVertices.at(6)));
465 local_info[3 * (this->dimension + 2) + this->dimension + 1] =
466 0.5 * std::min(this->prevVertices.at(3).d(this->prevVertices.at(2)),
467 this->prevVertices.at(3).d(this->prevVertices.at(7)));
468 local_info[4 * (this->dimension + 2) + this->dimension + 1] =
469 0.5 * std::min(this->prevVertices.at(4).d(this->prevVertices.at(0)),
470 this->prevVertices.at(4).d(this->prevVertices.at(5)));
471 local_info[5 * (this->dimension + 2) + this->dimension + 1] =
472 0.5 * std::min(this->prevVertices.at(5).d(this->prevVertices.at(1)),
473 this->prevVertices.at(5).d(this->prevVertices.at(4)));
474 local_info[6 * (this->dimension + 2) + this->dimension + 1] =
475 0.5 * std::min(this->prevVertices.at(6).d(this->prevVertices.at(2)),
476 this->prevVertices.at(6).d(this->prevVertices.at(7)));
477 local_info[7 * (this->dimension + 2) + this->dimension + 1] =
478 0.5 * std::min(this->prevVertices.at(7).d(this->prevVertices.at(3)),
479 this->prevVertices.at(7).d(this->prevVertices.at(6)));
480 break;
481 case 2:
482 local_info[0 * (this->dimension + 2) + this->dimension + 1] =
483 0.5 * std::min(this->prevVertices.at(0).d(this->prevVertices.at(1)),
484 this->prevVertices.at(0).d(this->prevVertices.at(2)));
485 local_info[1 * (this->dimension + 2) + this->dimension + 1] =
486 0.5 * std::min(this->prevVertices.at(1).d(this->prevVertices.at(0)),
487 this->prevVertices.at(1).d(this->prevVertices.at(3)));
488 local_info[2 * (this->dimension + 2) + this->dimension + 1] =
489 0.5 * std::min(this->prevVertices.at(2).d(this->prevVertices.at(0)),
490 this->prevVertices.at(2).d(this->prevVertices.at(3)));
491 local_info[3 * (this->dimension + 2) + this->dimension + 1] =
492 0.5 * std::min(this->prevVertices.at(3).d(this->prevVertices.at(1)),
493 this->prevVertices.at(3).d(this->prevVertices.at(2)));
494 local_info[4 * (this->dimension + 2) + this->dimension + 1] =
495 0.5 * std::min(this->prevVertices.at(4).d(this->prevVertices.at(5)),
496 this->prevVertices.at(4).d(this->prevVertices.at(6)));
497 local_info[5 * (this->dimension + 2) + this->dimension + 1] =
498 0.5 * std::min(this->prevVertices.at(5).d(this->prevVertices.at(4)),
499 this->prevVertices.at(5).d(this->prevVertices.at(7)));
500 local_info[6 * (this->dimension + 2) + this->dimension + 1] =
501 0.5 * std::min(this->prevVertices.at(6).d(this->prevVertices.at(4)),
502 this->prevVertices.at(6).d(this->prevVertices.at(7)));
503 local_info[7 * (this->dimension + 2) + this->dimension + 1] =
504 0.5 * std::min(this->prevVertices.at(7).d(this->prevVertices.at(5)),
505 this->prevVertices.at(7).d(this->prevVertices.at(6)));
506 break;
507 default:
509 __FILE__, __func__, __LINE__,
510 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
511 break;
512 }
513 // exchange information with all vertex neighbors
514 MPI_Request known_unused request[n_vertices];
515 MPI_Status known_unused status[n_vertices];
516
517 // compute new position for vertex 7 (if not periodic)
518 T known_unused total_work = (T)0;
519 T known_unused shift_vectors[this->dimension * n_vertices];
520
521 for (int v = 0; v < n_vertices; ++v) {
522 for (int d = 0; d < this->dimension; ++d) {
523 shift_vectors[v * this->dimension + d] = (T)0;
524 }
525 }
526
527 // TODO:
528 // 1.) collect work in main_dim communicators
529 // 2.) exchange work and extension in main_dim to neighbors in main_dim
530 // 3.) correct extension in main_dim
531 // 4.) do force-based correction in secondary dimensions
532 // 5.) find new neighbors in main_dim (hardest part, probably cell based)
533
534 int main_up;
535 int main_down;
536
537 // as each process computes the correction of the upper layer of vertices
538 // the source is 'main_up' and the target 'main_down'
539 MPI_Cart_shift(this->globalComm, main_dim, 1, &main_down, &main_up);
540
541 W local_work = this->work.at(0);
542 W remote_work;
543
544 // reduce work from local layer
545 MPI_Allreduce(MPI_IN_PLACE, &local_work, 1, MPIDataTypeW, MPI_SUM,
546 main_communicator);
547
548 // exchange local work with neighbor in main direction
549 MPI_Status state;
550 MPI_Sendrecv(&local_work, 1, MPIDataTypeW, main_down, 1010, &remote_work, 1,
551 MPIDataTypeW, main_up, 1010, this->globalComm, &state);
552
553 // transfer width of domains as well
554 T remote_width;
555 T local_width;
556 switch (main_dim) {
557 case 0:
558 local_width = this->prevVertices.at(1)[0] - this->prevVertices.at(0)[0];
559 break;
560 case 1:
561 local_width = this->prevVertices.at(2)[1] - this->prevVertices.at(0)[1];
562 break;
563 case 2:
564 local_width = this->prevVertices.at(4)[2] - this->prevVertices.at(0)[2];
565 break;
566 default:
568 __FILE__, __func__, __LINE__,
569 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
570 break;
571 }
572 MPI_Sendrecv(&local_width, 1, MPIDataTypeT, main_down, 1010, &remote_width, 1,
573 MPIDataTypeT, main_up, 1010, this->globalComm, &state);
574 T total_width = local_width + remote_width;
575 T max_main_shift = std::min(0.49 * std::min(local_width, remote_width),
576 std::min(local_width, remote_width) - 1.0);
577
578 // compute shift
579 T local_shift = (remote_work - local_work) / (remote_work + local_work) /
580 this->gamma / 2.0 * total_width;
581 local_shift = (std::abs(local_shift) > max_main_shift)
582 ? std::abs(local_shift) * max_main_shift / local_shift
583 : local_shift;
584 T remote_shift = 0;
585
586 // TODO: fix -> send to upper neighbor, receive from lower neighbor!!
587 // sendrecv not correct!
588
589 // transfer shift back
590 MPI_Sendrecv(&local_shift, 1, MPIDataTypeT, main_up, 2020, &remote_shift, 1,
591 MPIDataTypeT, main_down, 2020, this->globalComm, &state);
592
593 // apply shift in main direction
594 switch (main_dim) {
595 case 0:
596 if (this->local_coords.at(0) > 0) {
597 this->vertices.at(0)[0] = this->prevVertices.at(0)[0] + remote_shift;
598 this->vertices.at(2)[0] = this->prevVertices.at(2)[0] + remote_shift;
599 this->vertices.at(4)[0] = this->prevVertices.at(4)[0] + remote_shift;
600 this->vertices.at(6)[0] = this->prevVertices.at(6)[0] + remote_shift;
601 }
602 if (this->local_coords.at(0) < this->global_dims.at(0) - 1) {
603 this->vertices.at(1)[0] = this->prevVertices.at(1)[0] + local_shift;
604 this->vertices.at(3)[0] = this->prevVertices.at(3)[0] + local_shift;
605 this->vertices.at(5)[0] = this->prevVertices.at(5)[0] + local_shift;
606 this->vertices.at(7)[0] = this->prevVertices.at(7)[0] + local_shift;
607 }
608 break;
609 case 1:
610 if (this->local_coords.at(1) > 0) {
611 this->vertices.at(0)[1] = this->prevVertices.at(0)[1] + remote_shift;
612 this->vertices.at(1)[1] = this->prevVertices.at(1)[1] + remote_shift;
613 this->vertices.at(4)[1] = this->prevVertices.at(4)[1] + remote_shift;
614 this->vertices.at(5)[1] = this->prevVertices.at(5)[1] + remote_shift;
615 }
616 if (this->local_coords.at(1) < this->global_dims.at(1) - 1) {
617 this->vertices.at(2)[1] = this->prevVertices.at(2)[1] + local_shift;
618 this->vertices.at(3)[1] = this->prevVertices.at(3)[1] + local_shift;
619 this->vertices.at(6)[1] = this->prevVertices.at(6)[1] + local_shift;
620 this->vertices.at(7)[1] = this->prevVertices.at(7)[1] + local_shift;
621 }
622 break;
623 case 2:
624 if (this->local_coords.at(2) > 0) {
625 this->vertices.at(0)[2] = this->prevVertices.at(0)[2] + remote_shift;
626 this->vertices.at(1)[2] = this->prevVertices.at(1)[2] + remote_shift;
627 this->vertices.at(2)[2] = this->prevVertices.at(2)[2] + remote_shift;
628 this->vertices.at(3)[2] = this->prevVertices.at(3)[2] + remote_shift;
629 }
630 if (this->local_coords.at(2) < this->global_dims.at(2) - 1) {
631 this->vertices.at(4)[2] = this->prevVertices.at(4)[2] + local_shift;
632 this->vertices.at(5)[2] = this->prevVertices.at(5)[2] + local_shift;
633 this->vertices.at(6)[2] = this->prevVertices.at(6)[2] + local_shift;
634 this->vertices.at(7)[2] = this->prevVertices.at(7)[2] + local_shift;
635 }
636 break;
637 default:
639 __FILE__, __func__, __LINE__,
640 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
641 break;
642 }
643
644 if (this->localRank <= 1) {
645 std::cout << "DEBUG (main-dim shift): " << this->localRank << " "
646 << remote_work << " " << local_work << " " << total_width << " "
647 << local_shift << " " << remote_shift
648 << " 0: " << this->prevVertices.at(0)[2]
649 << " 0: " << this->vertices.at(0)[2]
650 << " 4: " << this->prevVertices.at(4)[2]
651 << " 4: " << this->vertices.at(4)[2] << " " << std::endl;
652 }
653
654 int last_vertex = (n_vertices - 1) * n_vertices * (this->dimension + 2);
655 int vertex_offset = this->dimension + 2;
656
657 // compute max shift
658 T max_shift = std::max(
659 0.49 * vertex_info[last_vertex + 0 * vertex_offset + this->dimension + 1],
660 1e-6);
661 for (int v = 1; v < n_vertices; ++v) {
662 max_shift = std::min(
663 max_shift,
664 0.49 *
665 vertex_info[last_vertex + v * vertex_offset + this->dimension + 1]);
666 }
667
668 // average work of neighboring processes for last vertex
669 T avg_work = 0.0;
670 for (int v = 0; v < n_vertices; ++v) {
671 avg_work += vertex_info[last_vertex + v * vertex_offset + this->dimension];
672 }
673 avg_work /= (T)n_vertices;
674
675 // compute shift vector
676 std::vector<T> vertex_shift(n_vertices * this->dimension, (T)0.0);
677 for (auto &sv : vertex_shift)
678 sv = (T)0.0;
679 int shift_offset = (n_vertices - 1) * this->dimension;
680
681 for (int v = 0; v < n_vertices; ++v) {
682 for (int d = 0; d < this->dimension; ++d) {
683 if (this->localRank == 0)
684 std::cout
685 << "DEBUG (shift x): " << v << ", " << d << " "
686 << (vertex_info[last_vertex + v * vertex_offset + this->dimension] -
687 avg_work) *
688 ((vertex_info[last_vertex + v * vertex_offset +
689 this->dimension] > avg_work)
690 ? 1.0
691 : -1.0) /
692 this->gamma *
693 (vertex_info[last_vertex + v * vertex_offset + d] -
694 this->prevVertices.at(v)[d])
695 << " => " << vertex_shift.at(shift_offset + d) <<
696
697 " max: " << max_shift
698 << " "
699 " avg_work: "
700 << avg_work
701 << " "
702 " work: "
703 << vertex_info[last_vertex + v * vertex_offset + this->dimension]
704 << " "
705 " 1/0: "
706 << ((vertex_info[last_vertex + v * vertex_offset +
707 this->dimension] > avg_work)
708 ? 1.0
709 : -1.0)
710
711 << std::endl;
712 vertex_shift.at(shift_offset + d) +=
713 (vertex_info[last_vertex + v * vertex_offset + this->dimension] -
714 avg_work) *
715 ((vertex_info[last_vertex + v * vertex_offset + this->dimension] >
716 avg_work)
717 ? 1.0
718 : -1.0) /
719 this->gamma *
720 (vertex_info[last_vertex + v * vertex_offset + d] -
721 this->prevVertices.at(v)[d]);
722 }
723 }
724
725 Point<T> shift_vector(3);
726 shift_vector[0] = vertex_shift.at(shift_offset);
727 shift_vector[1] = vertex_shift.at(shift_offset + 1);
728 shift_vector[2] = vertex_shift.at(shift_offset + 2);
729 T shift_length = shift_vector.norm();
730 if (this->localRank == 0)
731 std::cout << "DEBUG (shift length): " << shift_length
732 << " shift vector: " << shift_vector[0] << " " << shift_vector[1]
733 << " " << shift_vector[2] << " " << std::endl;
734
735 // apply correct length, if the shift vector is too large
736 if (shift_length > max_shift) {
737 for (int d = 0; d < this->dimension; ++d) {
738 vertex_shift.at(shift_offset + d) *= max_shift / shift_length;
739 }
740 }
741
742 if (this->localRank == 0) {
743 for (int v = 0; v < n_vertices; ++v) {
744 std::cout << "DEBUG shift vector vertex (before): " << v << ": ";
745 for (int d = 0; d < this->dimension; ++d) {
746 std::cout << vertex_shift.at(v * this->dimension + d) << " ";
747 }
748 std::cout << std::endl;
749 }
750 }
751
752 if (this->localRank == 0 || this->localRank == 1) {
753 for (int v = 0; v < n_vertices; ++v) {
754 std::cout << "DEBUG shift vector vertex " << v << ": ";
755 for (int d = 0; d < this->dimension; ++d) {
756 std::cout << vertex_shift.at(v * this->dimension + d) << " ";
757 }
758 std::cout << std::endl;
759 }
760 }
761
762 bool shift[3];
763 for (int z = 0; z < 2; ++z) {
764 if ((z == 0 && this->local_coords.at(2) == 0) ||
765 (z == 1 && this->local_coords.at(2) == this->global_dims.at(2) - 1))
766 shift[2] = false;
767 else
768 shift[2] = true;
769 for (int y = 0; y < 2; ++y) {
770 if ((y == 0 && this->local_coords.at(1) == 0) ||
771 (y == 1 && this->local_coords.at(1) == this->global_dims.at(1) - 1))
772 shift[1] = false;
773 else
774 shift[1] = true;
775 for (int x = 0; x < 2; ++x) {
776 if ((x == 0 && this->local_coords.at(0) == 0) ||
777 (x == 1 && this->local_coords.at(0) == this->global_dims.at(0) - 1))
778 shift[0] = false;
779 else
780 shift[0] = true;
781 for (int d = 0; d < this->dimension; ++d) {
782 if (d == main_dim)
783 continue;
784 int v = 4 * z + 2 * y + x;
785 if (shift[d])
786 this->vertices.at(v)[d] = this->prevVertices.at(v)[d] +
787 vertex_shift.at(v * this->dimension + d);
788 else
789 this->vertices.at(v)[d] = this->prevVertices.at(v)[d];
790 }
791 }
792 }
793 }
794}
795
796}//namespace ALL
797
798#endif
#define known_unused
Definition ALL_Defines.h:49
virtual W getEstimatedEfficiency() override
virtual void balance(int step) override
virtual std::vector< int > & getNeighbors() override
ForceBased_LB()
default constructor
void setup() override
setup internal data structures and parameters
ForceBased_LB(int d, W w, T g)
virtual void setAdditionalData(known_unused const void *data) override
~ForceBased_LB() override
default destructor
LB(const int dim, const T g)
Definition ALL_LB.hpp:49
int localRank
local rank in the used MPI communicator
Definition ALL_LB.hpp:212
std::vector< W > work
local work
Definition ALL_LB.hpp:219
MPI_Comm globalComm
used MPI communicator
Definition ALL_LB.hpp:210
int dimension
dimension of the used vertices
Definition ALL_LB.hpp:208
std::vector< int > global_dims
dimensions of the global process grid
Definition ALL_LB.hpp:225
virtual void setWork(const std::vector< W > &w)
Definition ALL_LB.hpp:112
T gamma
correction factor
Definition ALL_LB.hpp:206
std::vector< int > local_coords
cartesian coordinates of the local domain in the process grid
Definition ALL_LB.hpp:227
std::vector< Point< T > > prevVertices
original vertices before previous balancing step
Definition ALL_LB.hpp:223
std::vector< Point< T > > vertices
local vertices after previous balancing step
Definition ALL_LB.hpp:221
T norm(T nd=2)
Definition ALL.hpp:75