ALL 0.9.3
A Loadbalacing Library
Loading...
Searching...
No Matches
ALL_Voronoi.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_VORONOI_HEADER_INCLUDED
32#define ALL_VORONOI_HEADER_INCLUDED
33
34// number of maximum neighboring Voronoi cells
35#define ALL_VORONOI_MAX_NEIGHBORS 32
36// number of subblock division for Voronoi cell creation
37#define ALL_VORONOI_SUBBLOCKS 20
38
39// depth used to find neighbor cells
40#define ALL_VORONOI_NEIGHBOR_SEARCH_DEPTH 2
41
42#ifdef ALL_VORONOI_ACTIVE
44#include "ALL_LB.hpp"
45#include "ALL_Point.hpp"
46#include <algorithm>
47#include <exception>
48#include <iomanip>
49#include <map>
50#include <mpi.h>
51#include <sstream>
52
53#pragma GCC diagnostic push
54#pragma GCC diagnostic ignored "-Wunused-parameter"
55#include <voro++.hh>
56#pragma GCC diagnostic pop
57
58namespace ALL {
59
69template <class T, class W> class Voronoi_LB : public LB<T, W> {
70public:
77 Voronoi_LB(int d, W w, T g) : LB<T, W>(d, g) { this->setWork(w); }
78
80 ~Voronoi_LB() override;
81
82 voro::voronoicell vc;
83
85 void setup() override;
86
89 virtual void balance(int step) override;
90
91 // getter for variables (passed by reference to avoid
92 // direct access to private members of the object)
93
97 virtual std::vector<int> &getNeighbors() override;
98
101 virtual void setAdditionalData(known_unused const void *data) override {}
102
105 std::vector<T> &getNeighborVertices();
106
110 virtual W getEstimatedEfficiency() override {return (W)-1;};
111
112private:
113 // MPI values
114
115 // type for MPI communication
116 MPI_Datatype MPIDataTypeT;
117 MPI_Datatype MPIDataTypeW;
118
119 // list of neighbors
120 std::vector<int> neighbors;
121 int nNeighbors[1];
122
123 // collection of generator points
124 std::vector<T> generator_points;
125
126 // number of ranks in global communicator
127 int nDomains;
128};
129
130template <class T, class W> Voronoi_LB<T, W>::~Voronoi_LB() {}
131
132// setup routine for the tensor-based load-balancing scheme
133// requires:
134// this->globalComm (int): cartesian MPI communicator, from
135// which separate sub communicators
136// are derived in order to represent
137// each plane of domains in the system
138template <class T, class W> void Voronoi_LB<T, W>::setup() {
139 // store global communicator
140 int dimension = this->getDimension();
141
142 // no special communicator required
143
144 // create array to store information about generator points
145 // TODO: hierarchical scheme or better preselection which
146 // generator points are required for correct creation
147 // of the local cell -> large-scale Voronoi grid creation
148 // is too expensive to be done
149 // every step
150 MPI_Comm_size(this->globalComm, &nDomains);
151 MPI_Comm_rank(this->globalComm, &this->localRank);
152 generator_points.resize(nDomains * (2 * dimension + 1));
153
154 // determine correct MPI data type for template T
155 if (std::is_same<T, double>::value)
156 MPIDataTypeT = MPI_DOUBLE;
157 else if (std::is_same<T, float>::value)
158 MPIDataTypeT = MPI_FLOAT;
159 else if (std::is_same<T, int>::value)
160 MPIDataTypeT = MPI_INT;
161 else if (std::is_same<T, long>::value)
162 MPIDataTypeT = MPI_LONG;
163
164 // determine correct MPI data type for template W
165 if (std::is_same<W, double>::value)
166 MPIDataTypeW = MPI_DOUBLE;
167 else if (std::is_same<W, float>::value)
168 MPIDataTypeW = MPI_FLOAT;
169 else if (std::is_same<W, int>::value)
170 MPIDataTypeW = MPI_INT;
171 else if (std::is_same<W, long>::value)
172 MPIDataTypeW = MPI_LONG;
173}
174
175template <class T, class W>
176void Voronoi_LB<T, W>::balance(int known_unused loadbalancing_step) {
177 int dimension = this->getDimension();
178 // collect local information in array
179 int info_length = 2 * dimension + 1;
180 std::vector<T> local_info(info_length);
181 // Todo(s.schulz): prevVertices is unused, might be a bug?
182 std::vector<Point<T>> &prevVertices = this->getVertices();
183 std::vector<Point<T>> &vertices = this->getVertices();
184 std::vector<W> &work = this->getWork();
185
186 for (auto i = 0; i < 2; ++i)
187 for (auto d = 0; d < dimension; ++d) {
188 local_info.at(d) = vertices.at(i)[d];
189 }
190 local_info.at(2 * dimension) = (T)work.at(0);
191
192 MPI_Allgather(local_info.data(), info_length, MPIDataTypeT,
193 generator_points.data(), info_length, MPIDataTypeT,
194 this->globalComm);
195
196 // create Voronoi-cells
197 // TODO: check if system is periodic or not -> true / false!
198 voro::container con_old(
199 this->sysSize.at(0), this->sysSize.at(1), this->sysSize.at(2),
200 this->sysSize.at(3), this->sysSize.at(4), this->sysSize.at(5),
202 false, false, false, nDomains + 10);
203
204 // add generator points to container
205 for (auto d = 0; d < nDomains; ++d) {
206 con_old.put(d, generator_points.at(info_length * d),
207 generator_points.at(info_length * d + 1),
208 generator_points.at(info_length * d + 2));
209 }
210
211#ifdef ALL_DEBUG_ENABLED
212 // print old voronoi cell structure
213 std::ostringstream ss_local_gp_pov;
214 ss_local_gp_pov << "voronoi/generator_points_" << std::setw(7)
215 << std::setfill('0') << loadbalancing_step << ".pov";
216 std::ostringstream ss_local_gp_gnu;
217 ss_local_gp_gnu << "voronoi/generator_points_" << std::setw(7)
218 << std::setfill('0') << loadbalancing_step << ".gnu";
219 std::ostringstream ss_local_vc_pov;
220 ss_local_vc_pov << "voronoi/voronoi_cells_" << std::setw(7)
221 << std::setfill('0') << loadbalancing_step << ".pov";
222 std::ostringstream ss_local_vc_gnu;
223 ss_local_vc_gnu << "voronoi/voronoi_cells_" << std::setw(7)
224 << std::setfill('0') << loadbalancing_step << ".gnu";
225
226 loadbalancing_step++;
227
228 if (this->localRank == 0) {
229 con_old.draw_particles(ss_local_gp_gnu.str().c_str());
230 con_old.draw_cells_gnuplot(ss_local_vc_gnu.str().c_str());
231 con_old.draw_particles_pov(ss_local_gp_pov.str().c_str());
232 con_old.draw_cells_pov(ss_local_vc_pov.str().c_str());
233 }
234#endif
235
236 // vector to store neighbor information
237 std::vector<double> cell_vertices;
238 voro::voronoicell_neighbor c;
239
240 // generate loop over all old generator points
241 voro::c_loop_all cl_old(con_old);
242
243 int pid;
244 double x, y, z, r;
245 // reset list of neighbors
246 neighbors.clear();
247
248 std::map<int, double> neig_area;
249
250 for (auto i = 0; i < 1; ++i) {
251 // find next neighbors
252 std::vector<std::vector<int>> next_neighbors(
253 std::max((int)neighbors.size(), 1));
254 std::vector<std::vector<double>> neighbors_area(
255 std::max((int)neighbors.size(), 1));
256 int idx = 0;
257 if (cl_old.start()) {
258 do {
259 // compute only local cell
260 cl_old.pos(pid, x, y, z, r);
261 if (i == 0) {
262 if (pid == this->localRank) {
263 con_old.compute_cell(c, cl_old);
264 c.neighbors(next_neighbors.at(idx));
265 c.face_areas(neighbors_area.at(idx));
266 for (int j = 0; j < (int)next_neighbors.at(idx).size(); ++j)
267 neig_area.insert(std::pair<int, double>(
268 next_neighbors.at(idx).at(j), neighbors_area.at(idx).at(j)));
269 idx++;
270 if (idx == (int)neighbors.size())
271 break;
272 }
273 } else {
274 if (std::count(neighbors.begin(), neighbors.end(), pid) == 1) {
275 con_old.compute_cell(c, cl_old);
276 c.neighbors(next_neighbors.at(idx));
277 idx++;
278 if (idx == (int)neighbors.size())
279 break;
280 }
281 }
282 } while (cl_old.inc());
283 }
284 for (auto it = next_neighbors.begin(); it != next_neighbors.end(); ++it) {
285 neighbors.insert(neighbors.begin(), it->begin(), it->end());
286 }
287
288 std::sort(neighbors.begin(), neighbors.end());
289 auto uniq_it = std::unique(neighbors.begin(), neighbors.end());
290 neighbors.resize(std::distance(neighbors.begin(), uniq_it));
291 neighbors.erase(
292 std::remove(neighbors.begin(), neighbors.end(), this->localRank),
293 neighbors.end());
294 neighbors.erase(std::remove_if(neighbors.begin(), neighbors.end(),
295 [](int x) { return x < 0; }),
296 neighbors.end());
297 }
298
299 std::vector<double> volumes(neighbors.size());
300 std::vector<double> surfaces(neighbors.size());
301
302 for (int i = 0; i < (int)neighbors.size(); ++i) {
303 auto it = neig_area.find(neighbors.at(i));
304 surfaces.at(i) = it->second;
305 }
306
307 // Todo(s.schulz): is local_volume needed for future work or can we remove it?
308 double local_volume;
309 // get volumes of each neighbor cell
310 if (cl_old.start()) {
311 do {
312 cl_old.pos(pid, x, y, z, r);
313 auto it = std::find(neighbors.begin(), neighbors.end(), pid);
314
315 if (it != neighbors.end()) {
316 int idx = std::distance(neighbors.begin(), it);
317 con_old.compute_cell(c, cl_old);
318 volumes.at(idx) = c.volume();
319 }
320 if (pid == this->localRank) {
321 con_old.compute_cell(c, cl_old);
322 local_volume = c.volume();
323 }
324 } while (cl_old.inc());
325 }
326#ifdef ALL_DEBUG_ENABLED
327 MPI_Barrier(this->globalComm);
328 if (this->localRank == 0)
329 std::cout << "ALL::Voronoi_LB<T,W>::balance begin checking..." << std::endl;
330#endif
331
332 // compute shift
333 std::vector<T> shift(dimension, 0.0);
334 T norm;
335
336 T work_load = (T)local_info.at(info_length - 1);
337 for (auto it = neighbors.begin(); it != neighbors.end(); ++it) {
338 int neighbor = *it;
339 work_load += generator_points.at(info_length * neighbor + info_length - 1);
340 }
341
342#ifdef ALL_DEBUG_ENABLED
343 MPI_Barrier(this->globalComm);
344 if (this->localRank == 0)
345 std::cout << "ALL::Voronoi_LB<T,W>::balance find neighbors..." << std::endl;
346#endif
347
348 T max_diff = 20.0;
349
350 for (auto it = neighbors.begin(); it != neighbors.end(); ++it) {
351 int neighbor = *it;
352 std::vector<T> diff(dimension);
353 bool correct = true;
354 for (int d = 0; d < dimension && correct; ++d) {
355 diff.at(d) = generator_points.at(info_length * neighbor + dimension + d) -
356 local_info.at(d);
357
358 if (diff.at(d) >
359 0.5 * (this->sysSize[2 * d + 1] - this->sysSize[2 * d])) {
360 diff.at(d) -= (this->sysSize[2 * d + 1] - this->sysSize[2 * d]);
361 correct = false;
362 } else if (diff.at(d) <
363 -0.5 * (this->sysSize[2 * d + 1] - this->sysSize[2 * d])) {
364 diff.at(d) += (this->sysSize[2 * d + 1] - this->sysSize[2 * d]);
365 correct = false;
366 }
367 }
368 max_diff = std::min(max_diff,
369 sqrt(diff.at(0) * diff.at(0) + diff.at(1) * diff.at(1) +
370 diff.at(2) * diff.at(2)));
371 if (correct) {
372 // compute difference in work load
373 T work_diff = (T)0.0;
374 work_diff =
375 (generator_points.at(info_length * neighbor + info_length - 1) -
376 local_info.at(info_length - 1));
377 T work_density =
378 (generator_points.at(info_length * neighbor + info_length - 1) +
379 local_info.at(info_length - 1));
380
381 if (work_density < 1e-6)
382 work_diff = (T)0;
383 else
384 work_diff /= work_density;
385 for (int d = 0; d < dimension; ++d) {
386 shift.at(d) += 0.5 * work_diff * diff.at(d) / this->gamma;
387 }
388 }
389 }
390
391#ifdef ALL_DEBUG_ENABLED
392 MPI_Barrier(this->globalComm);
393 if (this->localRank == 0)
394 std::cout << "ALL::Voronoi_LB<T,W>::balance after neighbors found..."
395 << std::endl;
396#endif
397
398 norm = sqrt(shift.at(0) * shift.at(0) + shift.at(1) * shift.at(1) +
399 shift.at(2) * shift.at(2));
400
401 T scale = 1.0;
402
403 if (norm > 0.45 * max_diff)
404 scale = 0.45 * max_diff / norm;
405
406#ifdef ALL_DEBUG_ENABLED
407 MPI_Barrier(this->globalComm);
408 if (this->localRank == 0)
409 std::cout << "ALL::Voronoi_LB<T,W>::balance find new neighbors..." << std::endl;
410#endif
411
412 // to find new neighbors
413 for (int d = 0; d < dimension; ++d) {
414 local_info.at(d) += scale * shift.at(d);
415
416 // periodic correction of points
417 local_info.at(d) = (local_info.at(d) < this->sysSize.at(2 * d))
418 ? this->sysSize.at(2 * d) + 1.0
419 : ((local_info.at(d) >= this->sysSize.at(2 * d + 1))
420 ? this->sysSize.at(2 * d + 1) - 1.0
421 : local_info.at(d));
422 }
423
424 MPI_Allgather(local_info.data(), info_length, MPIDataTypeT,
425 generator_points.data(), info_length, MPIDataTypeT,
426 this->globalComm);
427
428#ifdef ALL_DEBUG_ENABLED
429 MPI_Barrier(this->globalComm);
430 if (this->localRank == 0)
431 std::cout << "ALL::Voronoi_LB<T,W>::balance create new containers..."
432 << std::endl;
433#endif
434
435 // create Voronoi-cells
436 // TODO: check if system is periodic or not -> true / false!
437 voro::container con_new(
438 this->sysSize.at(0), this->sysSize.at(1), this->sysSize.at(2),
439 this->sysSize.at(3), this->sysSize.at(4), this->sysSize.at(5),
441 false, false, false, nDomains + 10);
442
443 // add generator points to container
444 for (int d = 0; d < nDomains; ++d) {
445 con_new.put(d, generator_points.at(info_length * d),
446 generator_points.at(info_length * d + 1),
447 generator_points.at(info_length * d + 2));
448 }
449
450#ifdef ALL_DEBUG_ENABLED
451 std::ostringstream ss_local_gp_pov2;
452 ss_local_gp_pov2 << "voronoi/generator_points_" << std::setw(7)
453 << std::setfill('0') << loadbalancing_step << ".pov";
454 std::ostringstream ss_local_gp_gnu2;
455 ss_local_gp_gnu2 << "voronoi/generator_points_" << std::setw(7)
456 << std::setfill('0') << loadbalancing_step << ".gnu";
457 std::ostringstream ss_local_vc_pov2;
458 ss_local_vc_pov2 << "voronoi/voronoi_cells_" << std::setw(7)
459 << std::setfill('0') << loadbalancing_step << ".pov";
460 std::ostringstream ss_local_vc_gnu2;
461 ss_local_vc_gnu2 << "voronoi/voronoi_cells_" << std::setw(7)
462 << std::setfill('0') << loadbalancing_step << ".gnu";
463
464 loadbalancing_step++;
465
466 if (this->localRank == 0) {
467 con_new.draw_particles(ss_local_gp_gnu2.str().c_str());
468 con_new.draw_cells_gnuplot(ss_local_vc_gnu2.str().c_str());
469 con_new.draw_particles_pov(ss_local_gp_pov2.str().c_str());
470 con_new.draw_cells_pov(ss_local_vc_pov2.str().c_str());
471 }
472#endif
473
474#ifdef ALL_DEBUG_ENABLED
475 MPI_Barrier(this->globalComm);
476 if (this->localRank == 0)
477 std::cout << "ALL::Voronoi_LB<T,W>::balance storing shifted vertices..."
478 << " " << vertices.size() << " " << dimension << std::endl;
479#endif
480
481 Point<T> tmp_pnt(dimension, local_info.data());
482
483 // compute new neighboring cells and generator points
484
485 // generate loop over all new generator points
486 voro::c_loop_all cl_new(con_new);
487
488 // reset list of neighbors
489 neighbors.clear();
490
491#ifdef ALL_DEBUG_ENABLED
492 MPI_Barrier(this->globalComm);
493 if (this->localRank == 0)
494 std::cout << "ALL::Voronoi_LB<T,W>::balance searching neighboring vertices..."
495 << std::endl;
496#endif
497
498 for (auto i = 0; i < ALL_VORONOI_NEIGHBOR_SEARCH_DEPTH; ++i) {
499 // find next neighbors
500 std::vector<std::vector<int>> next_neighbors(
501 std::max((int)neighbors.size(), 1));
502 int idx = 0;
503 if (cl_new.start()) {
504 do {
505 // compute next voronoi cell
506 cl_new.pos(pid, x, y, z, r);
507 if (i == 0) {
508 if (pid == this->localRank) {
509 con_new.compute_cell(c, cl_new);
510 c.neighbors(next_neighbors.at(idx));
511 idx++;
512 if (idx == (int)neighbors.size())
513 break;
514 }
515 } else {
516 if (std::count(neighbors.begin(), neighbors.end(), pid) == 1) {
517 con_new.compute_cell(c, cl_new);
518 c.neighbors(next_neighbors.at(idx));
519 idx++;
520 if (idx == (int)neighbors.size())
521 break;
522 }
523 }
524 } while (cl_new.inc());
525 }
526 for (auto it = next_neighbors.begin(); it != next_neighbors.end(); ++it) {
527 neighbors.insert(neighbors.begin(), it->begin(), it->end());
528 }
529
530 std::sort(neighbors.begin(), neighbors.end());
531 auto uniq_it = std::unique(neighbors.begin(), neighbors.end());
532 neighbors.resize(std::distance(neighbors.begin(), uniq_it));
533 neighbors.erase(
534 std::remove(neighbors.begin(), neighbors.end(), this->localRank),
535 neighbors.end());
536 neighbors.erase(std::remove_if(neighbors.begin(), neighbors.end(),
537 [](int x) { return x < 0; }),
538 neighbors.end());
539 }
540
541 // determine number of neighbors
542 nNeighbors[0] = neighbors.size();
543
544 // clear old neighbor vertices
545 this->neighborVertices.clear();
546
547 // find vertices of neighbors
548
549 for (auto n : neighbors) {
550 for (int i = 0; i < dimension; ++i)
551 this->neighborVertices.push_back(generator_points.at(info_length * n + i));
552 }
553}
554
555// provide list of neighbors
556template <class T, class W>
558 return neighbors;
559}
560
561
562}//namespace ALL
563
564#endif
565#endif
#define known_unused
Definition ALL_Defines.h:49
#define ALL_VORONOI_NEIGHBOR_SEARCH_DEPTH
#define ALL_VORONOI_SUBBLOCKS
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< T > sysSize
(orthogonal) system size
Definition ALL_LB.hpp:217
virtual void setWork(const std::vector< W > &w)
Definition ALL_LB.hpp:112
virtual std::vector< W > & getWork()
Definition ALL_LB.hpp:144
T gamma
correction factor
Definition ALL_LB.hpp:206
virtual int getDimension()
Definition ALL_LB.hpp:99
std::vector< T > neighborVertices
vertices describing neighboring domains
Definition ALL_LB.hpp:231
virtual std::vector< Point< T > > & getVertices()
Definition ALL_LB.hpp:148
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
voro::voronoicell vc
Voronoi_LB(int d, W w, T g)
std::vector< T > & getNeighborVertices()
void setup() override
setup internal data structures and parameters
virtual void setAdditionalData(known_unused const void *data) override
Voronoi_LB()
default constructor
~Voronoi_LB() override
default destructor
virtual void balance(int step) override
virtual std::vector< int > & getNeighbors() override
virtual W getEstimatedEfficiency() override
Definition ALL.hpp:75