ALL 0.9.3
A Loadbalacing Library
Loading...
Searching...
No Matches
ALL.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_MAIN_HEADER_INC
32#define ALL_MAIN_HEADER_INC
33
35#include "ALL_LB.hpp"
36#include "ALL_Point.hpp"
37#include "ALL_Tensor.hpp"
38#ifdef ALL_VORONOI_ACTIVE
39#include "ALL_Voronoi.hpp"
40#endif
41#include "ALL_ForceBased.hpp"
42#include "ALL_Histogram.hpp"
43#include "ALL_Staggered.hpp"
44#include <algorithm>
45#include <iomanip>
46#include <memory>
47#include <numeric>
48#include <sstream>
49#include <string>
50#include <vector>
51
52#include <cerrno>
53#include <sys/stat.h>
54
55#ifdef ALL_VTK_OUTPUT
56#include <limits>
57#include <vtkCellArray.h>
58#include <vtkCellData.h>
59#include <vtkFloatArray.h>
60#include <vtkInformation.h>
61#include <vtkIntArray.h>
62#include <vtkMPIController.h>
63#include <vtkPolyhedron.h>
64#include <vtkProgrammableFilter.h>
65#include <vtkSmartPointer.h>
66#include <vtkUnstructuredGrid.h>
67#include <vtkVoxel.h>
68#include <vtkXMLPUnstructuredGridWriter.h>
69#include <vtkXMLUnstructuredGridWriter.h>
70#ifdef VTK_CELL_ARRAY_V2
71#include <vtkNew.h>
72#endif
73#endif
74
75namespace ALL {
76
77#define ALL_ESTIMATION_LIMIT_DEFAULT 0
78
80enum LB_t : int {
84 TENSOR = 1,
87#ifdef ALL_VORONOI_ACTIVE
90#endif
97};
98
101
102template <class T, class W> class ALL {
103public:
107 : loadbalancing_step(0)
108#ifdef ALL_VTK_OUTPUT
109 ,
110 vtk_init(false)
111#endif
112 {
113 // determine correct MPI data type for template T
114 if (std::is_same<T, double>::value)
115 MPIDataTypeT = MPI_DOUBLE;
116 else if (std::is_same<T, float>::value)
117 MPIDataTypeT = MPI_FLOAT;
118 else if (std::is_same<T, int>::value)
119 MPIDataTypeT = MPI_INT;
120 else if (std::is_same<T, long>::value)
121 MPIDataTypeT = MPI_LONG;
122
123 // determine correct MPI data type for template W
124 if (std::is_same<W, double>::value)
125 MPIDataTypeW = MPI_DOUBLE;
126 else if (std::is_same<W, float>::value)
127 MPIDataTypeW = MPI_FLOAT;
128 else if (std::is_same<W, int>::value)
129 MPIDataTypeW = MPI_INT;
130 else if (std::is_same<W, long>::value)
131 MPIDataTypeW = MPI_LONG;
132 }
133
142 ALL(const LB_t m, const int d, const T g) : ALL() {
143 method = m;
144 switch (method) {
145 case LB_t::TENSOR:
146 balancer.reset(new Tensor_LB<T, W>(d, (W)0, g));
147 break;
148 case LB_t::STAGGERED:
149 balancer.reset(new Staggered_LB<T, W>(d, (W)0, g));
150 break;
151 case LB_t::FORCEBASED:
152 balancer.reset(new ForceBased_LB<T, W>(d, (W)0, g));
153 break;
154#ifdef ALL_VORONOI_ACTIVE
155 case LB_t::VORONOI:
156 balancer.reset(new Voronoi_LB<T, W>(d, (W)0, g));
157 break;
158#endif
159 case LB_t::HISTOGRAM:
160 balancer.reset(new Histogram_LB<T, W>(d, std::vector<W>(10), g));
161 break;
162 case LB_t::TENSOR_MAX:
163 balancer.reset(new TensorMax_LB<T, W>(d, (W)0, g));
164 break;
165 default:
166 throw InvalidArgumentException(__FILE__, __func__, __LINE__,
167 "Unknown type of loadbalancing passed.");
168 }
169 balancer->setDimension(d);
170 balancer->setGamma(g);
171 estimation_limit = ALL_ESTIMATION_LIMIT_DEFAULT;
172 }
173
183 ALL(const LB_t m, const int d, const std::vector<Point<T>> &inp, const T g)
184 : ALL(m, d, g) {
185 balancer->setVertices(inp);
186 calculate_outline();
187 }
188
190 ~ALL() = default;
191
195 void setVertices(const std::vector<Point<T>> &inp) {
196 int dim = inp.at(0).getDimension();
197 if (dim != balancer->getDimension())
199 __FILE__, __func__, __LINE__,
200 "Dimension of ALL::Points in input vector do not match dimension of "
201 "ALL "
202 "object.");
203 balancer->setVertices(inp);
204 calculate_outline();
205 }
206
211 void setCommunicator(const MPI_Comm comm) {
212 int comm_type;
213 MPI_Topo_test(comm, &comm_type);
214 if (comm_type == MPI_CART) {
215 communicator = comm;
216 balancer->setCommunicator(communicator);
217 // set procGridLoc and procGridDim
218 const int dimension = balancer->getDimension();
219 int location[dimension];
220 int size[dimension];
221 int periodicity[dimension];
222 MPI_Cart_get(communicator, dimension, size, periodicity, location);
223 procGridLoc.assign(location, location + dimension);
224 procGridDim.assign(size, size + dimension);
225 procGridSet = true;
226 } else {
227 int rank;
228 MPI_Comm_rank(comm, &rank);
229 if (procGridSet) {
230 // compute 1D index from the location in the process grid (using z-y-x
231 // ordering, as MPI_Dims_create)
232 int idx;
233 switch (balancer->getDimension()) {
234 case 3:
235 idx = procGridLoc.at(2) + procGridLoc.at(1) * procGridDim.at(2) +
236 procGridLoc.at(0) * procGridDim.at(2) * procGridDim.at(1);
237 break;
238 case 2:
239 idx = procGridLoc.at(1) + procGridLoc.at(0) * procGridDim.at(1);
240 break;
241 case 1:
242 idx = procGridLoc.at(0);
243 break;
244 default:
246 __FILE__, __func__, __LINE__,
247 "ALL cannot deal with more then three dimensions (or less then "
248 "one).");
249 break;
250 }
251
252 // create new temporary communicator with correct ranks
253 MPI_Comm temp_comm;
254 MPI_Comm_split(comm, 0, idx, &temp_comm);
255
256 std::vector<int> periods(3, 1);
257
258 // transform temporary communicator to cartesian communicator
259 MPI_Cart_create(temp_comm, balancer->getDimension(), procGridDim.data(),
260 periods.data(), 0, &communicator);
261 balancer->setCommunicator(communicator);
262 } else {
264 __FILE__, __func__, __LINE__,
265 "When using a non-cartesian communicator process grid parameters "
266 "required (not set).");
267 }
268 }
269 }
270
273 T getGamma() { return balancer->getGamma(); };
274
277 void setGamma(const T g) { balancer->setGamma(g); };
278
281 void setWork(const W work) {
282 // set new value for work (single value for whole domain)
283 balancer->setWork(work);
284 }
285
286 // method to set a multi-dimensional work for the local domain
288 void setWork(const std::vector<W> &work) { balancer->setWork(work); }
289
292 void setup() { balancer->setup(); }
293
299 std::vector<Point<T>> &balance() {
300 nVertices = balancer->getNVertices();
301 switch (method) {
302 case LB_t::TENSOR:
303 balancer->balance(loadbalancing_step);
304 calculate_outline();
305 break;
306 case LB_t::STAGGERED:
307
308 /* ToDo: Check if repetition is useful at all and
309 * change it to be included for all methods,
310 * not only STAGGERED */
311 /*
312 // estimation to improve speed of convergence
313 // changing vertices during estimation
314 std::vector<Point<T>> tmp_vertices = balancer->getVertices();
315 // saved original vertices
316 std::vector<Point<T>> old_vertices = balancer->getVertices();
317 // old temporary vertices
318 std::vector<Point<T>> old_tmp_vertices = balancer->getVertices();
319 // result temporary vertices
320 std::vector<Point<T>> result_vertices = balancer->getVertices();
321 // temporary work
322 W tmp_work = balancer->getWork().at(0);
323 // old temporary work
324 W old_tmp_work = balancer->getWork().at(0);
325
326 W max_work;
327 W sum_work;
328 double d_max;
329 int n_ranks;
330
331 // compute work density on local domain
332 double V = 1.0;
333 for (int i = 0; i < balancer->getDimension(); ++i)
334 V *= ( outline.at(1)[i] - outline.at(0)[i] );
335 double rho = workArray.at(0) / V;
336
337 // collect maximum work in system
338 MPI_Allreduce(workArray.data(), &max_work, 1, MPIDataTypeW, MPI_MAX,
339 communicator); MPI_Allreduce(workArray.data(), &sum_work, 1,
340 MPIDataTypeW, MPI_SUM, communicator);
341 MPI_Comm_size(communicator,&n_ranks);
342 d_max = (double)(max_work) * (double)n_ranks / (double)sum_work - 1.0;
343
344
345 // internal needs to be readded to argument list
346 const bool internal = false
347 for (int i = 0; i < estimation_limit && d_max > 0.1 && !internal; ++i)
348 {
349 balancer->setWork(tmp_work);
350 balancer->setup();
351 balancer->balance(true);
352 bool sane = true;
353 // check if the estimated boundaries are not too deformed
354 for (int j = 0; j < balancer->getDimension(); ++j)
355 sane = sane && (result_vertices.at(0)[j] <=
356 old_vertices.at(1)[j]);
357 MPI_Allreduce(MPI_IN_PLACE,&sane,1,MPI_CXX_BOOL,MPI_LAND,communicator);
358 if (sane)
359 {
360 old_tmp_vertices = tmp_vertices;
361 tmp_vertices = result_vertices;
362 V = 1.0;
363 for (int i = 0; i < balancer->getDimension(); ++i)
364 V *= ( tmp_vertices.at(1)[i] - tmp_vertices.at(0)[i] );
365 old_tmp_work = tmp_work;
366 tmp_work = rho * V;
367 }
368 else if (!sane || i == estimation_limit - 1)
369 {
370 balancer->getVertices() = old_tmp_vertices;
371 workArray->at(0) = old_tmp_work;
372 calculate_outline();
373 i = estimation_limit;
374 }
375 }
376
377 ((Staggered_LB<T,W>*)balancer.get())->setVertices(outline->data());
378 */
379 balancer->balance(loadbalancing_step);
380 calculate_outline();
381 break;
382 case LB_t::FORCEBASED:
383 balancer->balance(loadbalancing_step);
384 calculate_outline();
385 break;
386#ifdef ALL_VORONOI_ACTIVE
387 case LB_t::VORONOI:
388 balancer->balance(loadbalancing_step);
389 break;
390#endif
391 case LB_t::HISTOGRAM:
392 balancer->balance(loadbalancing_step);
393 calculate_outline();
394 break;
395 case LB_t::TENSOR_MAX:
396 balancer->balance(loadbalancing_step);
397 calculate_outline();
398 break;
399
400 default:
401 throw InvalidArgumentException(__FILE__, __func__, __LINE__,
402 "Unknown type of loadbalancing passed.");
403 }
404 loadbalancing_step++;
405 return balancer->getVertices();
406 }
407
408 /**********************/
409 /* getter functions */
410 /**********************/
411
415 std::vector<Point<T>> &getPrevVertices() {
416 return balancer->getPrevVertices();
417 };
418
423 std::vector<Point<T>> &getVertices() { return balancer->getVertices(); };
424
428 int getDimension() { return balancer->getDimension(); }
429
433 void getWork(std::vector<W> &result) { result = balancer->getWork(); }
434
437 W getWork() { return balancer->getWork().at(0); }
438
442 std::vector<int> &getNeighbors() { return balancer->getNeighbors(); }
443
448 std::vector<T> &getNeighborVertices() {
449 return balancer->getNeighborVertices();
450 }
451
455 return balancer->getEfficiency();
456 }
457
461 // currently only implemented for HISTOGRAM
462 return balancer->getEstimatedEfficiency();
463 }
464
465 /**********************/
466 /* setter functions */
467 /**********************/
468
474 void setSysSize(const std::vector<T> &sysSize) {
475 balancer.get()->setSysSize(sysSize);
476 }
477
486 void setMethodData(const void *data) {
487 balancer.get()->setAdditionalData(data);
488 }
489
494 void setProcGridParams(const std::vector<int> &loc,
495 const std::vector<int> &size) {
496 // todo(s.schulz): We should probably throw if the dimension does not match
497 procGridLoc.insert(procGridLoc.begin(), loc.begin(), loc.end());
498 procGridDim.insert(procGridDim.begin(), size.begin(), size.end());
499 procGridSet = true;
500 }
501
504 void setProcTag(int tag) { procTag = tag; }
505
510 void setMinDomainSize(const std::vector<T> &minSize) {
511 (balancer.get())->setMinDomainSize(minSize);
512 }
513
514#ifdef ALL_VTK_OUTPUT
519 void printVTKoutlines(const int step) {
520 // define grid points, i.e. vertices of local domain
521 auto points = vtkSmartPointer<vtkPoints>::New();
522 // allocate space for eight points (describing the cuboid around the domain)
523 points->Allocate(8 * balancer->getDimension());
524
525 int n_ranks;
526 int local_rank;
527
528 if (!vtk_init) {
529 controller = vtkMPIController::New();
530 controller->Initialize();
531 controller->SetGlobalController(controller);
532 vtk_init = true;
533 }
534
535 MPI_Comm_rank(communicator, &local_rank);
536 MPI_Comm_size(communicator, &n_ranks);
537
538 std::vector<Point<W>> tmp_outline(2);
539 std::vector<W> tmp_0(
540 {outline.at(0)[0], outline.at(0)[1], outline.at(0)[2]});
541 std::vector<W> tmp_1(
542 {outline.at(1)[0], outline.at(1)[1], outline.at(1)[2]});
543 tmp_outline.at(0) = Point<W>(tmp_0);
544 tmp_outline.at(1) = Point<W>(tmp_1);
545
546 for (auto z = 0; z <= 1; ++z)
547 for (auto y = 0; y <= 1; ++y)
548 for (auto x = 0; x <= 1; ++x) {
549 points->InsertPoint(x + 2 * y + 4 * z, tmp_outline.at(x)[0],
550 tmp_outline.at(y)[1], tmp_outline.at(z)[2]);
551 }
552
553 auto hexa = vtkSmartPointer<vtkVoxel>::New();
554 for (int i = 0; i < 8; ++i)
555 hexa->GetPointIds()->SetId(i, i);
556
557 auto cellArray = vtkSmartPointer<vtkCellArray>::New();
558 cellArray->InsertNextCell(hexa);
559
560 // define work array (length = 1)
561 auto work = vtkSmartPointer<vtkFloatArray>::New();
562 work->SetNumberOfComponents(1);
563 work->SetNumberOfTuples(1);
564 work->SetName("work");
565 W total_work = std::accumulate(balancer->getWork().begin(),
566 balancer->getWork().end(), (W)0);
567 work->SetValue(0, total_work);
568
569 // define rank array (length = 4, x,y,z, rank)
570 int rank = 0;
571 MPI_Comm_rank(communicator, &rank);
572 auto coords = vtkSmartPointer<vtkFloatArray>::New();
573 coords->SetNumberOfComponents(3);
574 coords->SetNumberOfTuples(1);
575 coords->SetName("coords");
576 coords->SetValue(0, procGridLoc.at(0));
577 coords->SetValue(1, procGridLoc.at(1));
578 coords->SetValue(2, procGridLoc.at(2));
579
580 auto rnk = vtkSmartPointer<vtkFloatArray>::New();
581 rnk->SetNumberOfComponents(1);
582 rnk->SetNumberOfTuples(1);
583 rnk->SetName("MPI rank");
584 rnk->SetValue(0, rank);
585
586 // define tag array (length = 1)
587 auto tag = vtkSmartPointer<vtkIntArray>::New();
588 tag->SetNumberOfComponents(1);
589 tag->SetNumberOfTuples(1);
590 tag->SetName("tag");
591 tag->SetValue(0, procTag);
592
593 // determine extent of system
594 W known_unused global_extent[6];
595 W local_min[3];
596 W local_max[3];
597 W global_min[3];
598 W global_max[3];
599 W width[3];
600 W volume = (W)1.0;
601
602 for (int i = 0; i < 3; ++i) {
603 local_min[i] = outline.at(0)[i];
604 local_max[i] = outline.at(1)[i];
605 width[i] = local_max[i] - local_min[i];
606 volume *= width[i];
607 }
608
609 W surface = (W)2.0 * width[0] * width[1] + (W)2.0 * width[1] * width[2] +
610 (W)2.0 * width[0] * width[2];
611
612 auto expanse = vtkSmartPointer<vtkFloatArray>::New();
613 expanse->SetNumberOfComponents(6);
614 expanse->SetNumberOfTuples(1);
615 expanse->SetName("expanse");
616 expanse->SetValue(0, width[0]);
617 expanse->SetValue(1, width[0]);
618 expanse->SetValue(2, width[0]);
619 expanse->SetValue(3, volume);
620 expanse->SetValue(4, surface);
621 expanse->SetValue(5, surface / volume);
622
623 MPI_Allreduce(&local_min, &global_min, 3, MPIDataTypeW, MPI_MIN,
624 communicator);
625 MPI_Allreduce(&local_max, &global_max, 3, MPIDataTypeW, MPI_MAX,
626 communicator);
627
628 for (int i = 0; i < 3; ++i) {
629 global_extent[2 * i] = global_min[i];
630 global_extent[2 * i + 1] = global_max[i];
631 }
632
633 // create a structured grid and assign data to it
634 auto unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
635 unstructuredGrid->SetPoints(points);
636 unstructuredGrid->SetCells(VTK_VOXEL, cellArray);
637 unstructuredGrid->GetCellData()->AddArray(work);
638 unstructuredGrid->GetCellData()->AddArray(expanse);
639 unstructuredGrid->GetCellData()->AddArray(rnk);
640 unstructuredGrid->GetCellData()->AddArray(coords);
641 unstructuredGrid->GetCellData()->AddArray(tag);
642
643 createDirectory("vtk_outline");
644 std::ostringstream ss_local;
645 ss_local << "vtk_outline/ALL_vtk_outline_" << std::setw(7)
646 << std::setfill('0') << step << "_" << local_rank << ".vtu";
647
648 auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
649 writer->SetInputData(unstructuredGrid);
650 writer->SetFileName(ss_local.str().c_str());
651 writer->SetDataModeToAscii();
652 // writer->SetDataModeToBinary();
653 writer->Write();
654
655 // if (local_rank == 0)
656 //{
657 std::ostringstream ss_para;
658 ss_para << "vtk_outline/ALL_vtk_outline_" << std::setw(7)
659 << std::setfill('0') << step << ".pvtu";
660 // create the parallel writer
661 auto parallel_writer =
662 vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
663 parallel_writer->SetFileName(ss_para.str().c_str());
664 parallel_writer->SetNumberOfPieces(n_ranks);
665 parallel_writer->SetStartPiece(local_rank);
666 parallel_writer->SetEndPiece(local_rank);
667 parallel_writer->SetInputData(unstructuredGrid);
668 parallel_writer->SetDataModeToAscii();
669 // parallel_writer->SetDataModeToBinary();
670 parallel_writer->Write();
671 //}
672 }
673
678 void printVTKvertices(const int step) {
679 int n_ranks;
680 int local_rank;
681
682 MPI_Comm_rank(communicator, &local_rank);
683 MPI_Comm_size(communicator, &n_ranks);
684
685 // local vertices
686 // (vertices + work)
687 T local_vertices[nVertices * balancer->getDimension() + 1];
688
689 for (int v = 0; v < nVertices; ++v) {
690 for (int d = 0; d < balancer->getDimension(); ++d) {
691 local_vertices[v * balancer->getDimension() + d] =
692 balancer->getVertices().at(v)[d];
693 }
694 }
695 local_vertices[nVertices * balancer->getDimension()] =
696 (T)balancer->getWork().at(0);
697
698 /*
699 T *global_vertices;
700 if (local_rank == 0) {
701 global_vertices =
702 new T[n_ranks * (nVertices * balancer->getDimension() + 1)];
703 }
704 */
705 T global_vertices[n_ranks * (nVertices * balancer->getDimension() + 1)];
706
707 // collect all works and vertices on a single process
708 MPI_Gather(local_vertices, nVertices * balancer->getDimension() + 1,
709 MPIDataTypeT, global_vertices,
710 nVertices * balancer->getDimension() + 1, MPIDataTypeT, 0,
711 communicator);
712
713 if (local_rank == 0) {
714 auto vtkpoints = vtkSmartPointer<vtkPoints>::New();
715#ifdef VTK_CELL_ARRAY_V2
716 vtkNew<vtkUnstructuredGrid> unstructuredGrid;
717 unstructuredGrid->Allocate(n_ranks + 10);
718#else
719 auto unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
720#endif
721
722 // enter vertices into unstructured grid
723 for (int i = 0; i < n_ranks; ++i) {
724 for (int v = 0; v < nVertices; ++v) {
725 vtkpoints->InsertNextPoint(
726 global_vertices[i * (nVertices * balancer->getDimension() + 1) +
727 v * balancer->getDimension() + 0],
728 global_vertices[i * (nVertices * balancer->getDimension() + 1) +
729 v * balancer->getDimension() + 1],
730 global_vertices[i * (nVertices * balancer->getDimension() + 1) +
731 v * balancer->getDimension() + 2]);
732 }
733 }
734 unstructuredGrid->SetPoints(vtkpoints);
735
736 // data arrays for work and cell id
737 auto work = vtkSmartPointer<vtkFloatArray>::New();
738 auto cell = vtkSmartPointer<vtkFloatArray>::New();
739 work->SetNumberOfComponents(1);
740 work->SetNumberOfTuples(n_ranks);
741 work->SetName("work");
742 cell->SetNumberOfComponents(1);
743 cell->SetNumberOfTuples(n_ranks);
744 cell->SetName("cell id");
745
746 for (int n = 0; n < n_ranks; ++n) {
747
748#ifdef VTK_CELL_ARRAY_V2
749 // define grid points, i.e. vertices of local domain
750 vtkIdType pointIds[8] = {8 * n + 0, 8 * n + 1, 8 * n + 2, 8 * n + 3,
751 8 * n + 4, 8 * n + 5, 8 * n + 6, 8 * n + 7};
752
753 vtkIdType faces[48] = {
754 3, 8 * n + 0, 8 * n + 2, 8 * n + 1,
755 3, 8 * n + 1, 8 * n + 2, 8 * n + 3,
756 3, 8 * n + 0, 8 * n + 4, 8 * n + 2,
757 3, 8 * n + 2, 8 * n + 4, 8 * n + 6,
758 3, 8 * n + 2, 8 * n + 6, 8 * n + 3,
759 3, 8 * n + 3, 8 * n + 6, 8 * n + 7,
760 3, 8 * n + 1, 8 * n + 5, 8 * n + 3,
761 3, 8 * n + 3, 8 * n + 5, 8 * n + 7,
762 3, 8 * n + 0, 8 * n + 4, 8 * n + 1,
763 3, 8 * n + 1, 8 * n + 4, 8 * n + 5,
764 3, 8 * n + 4, 8 * n + 6, 8 * n + 5,
765 3, 8 * n + 5, 8 * n + 6, 8 * n + 7};
766
767 unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
768 faces);
769#else
770 // define grid points, i.e. vertices of local domain
771 vtkIdType pointIds[8] = {8 * n + 0, 8 * n + 1, 8 * n + 2, 8 * n + 3,
772 8 * n + 4, 8 * n + 5, 8 * n + 6, 8 * n + 7};
773
774 auto faces = vtkSmartPointer<vtkCellArray>::New();
775 // setup faces of polyhedron
776 vtkIdType f0[3] = {8 * n + 0, 8 * n + 2, 8 * n + 1};
777 vtkIdType f1[3] = {8 * n + 1, 8 * n + 2, 8 * n + 3};
778
779 vtkIdType f2[3] = {8 * n + 0, 8 * n + 4, 8 * n + 2};
780 vtkIdType f3[3] = {8 * n + 2, 8 * n + 4, 8 * n + 6};
781
782 vtkIdType f4[3] = {8 * n + 2, 8 * n + 6, 8 * n + 3};
783 vtkIdType f5[3] = {8 * n + 3, 8 * n + 6, 8 * n + 7};
784
785 vtkIdType f6[3] = {8 * n + 1, 8 * n + 5, 8 * n + 3};
786 vtkIdType f7[3] = {8 * n + 3, 8 * n + 5, 8 * n + 7};
787
788 vtkIdType f8[3] = {8 * n + 0, 8 * n + 4, 8 * n + 1};
789 vtkIdType f9[3] = {8 * n + 1, 8 * n + 4, 8 * n + 5};
790
791 vtkIdType fa[3] = {8 * n + 4, 8 * n + 6, 8 * n + 5};
792 vtkIdType fb[3] = {8 * n + 5, 8 * n + 6, 8 * n + 7};
793
794 faces->InsertNextCell(3, f0);
795 faces->InsertNextCell(3, f1);
796 faces->InsertNextCell(3, f2);
797 faces->InsertNextCell(3, f3);
798 faces->InsertNextCell(3, f4);
799 faces->InsertNextCell(3, f5);
800 faces->InsertNextCell(3, f6);
801 faces->InsertNextCell(3, f7);
802 faces->InsertNextCell(3, f8);
803 faces->InsertNextCell(3, f9);
804 faces->InsertNextCell(3, fa);
805 faces->InsertNextCell(3, fb);
806
807 unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
808 faces->GetPointer());
809#endif
810 work->SetValue(
811 n, global_vertices[n * (nVertices * balancer->getDimension() + 1) +
812 8 * balancer->getDimension()]);
813 cell->SetValue(n, (T)n);
814 }
815 unstructuredGrid->GetCellData()->AddArray(work);
816 unstructuredGrid->GetCellData()->AddArray(cell);
817
818 createDirectory("vtk_vertices");
819 std::ostringstream filename;
820 filename << "vtk_vertices/ALL_vtk_vertices_" << std::setw(7)
821 << std::setfill('0') << step << ".vtu";
822 auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
823 writer->SetInputData(unstructuredGrid);
824 writer->SetFileName(filename.str().c_str());
825 writer->SetDataModeToAscii();
826 // writer->SetDataModeToBinary();
827 writer->Write();
828
829 // delete[] global_vertices;
830 }
831 }
832#endif
833
834private:
836 LB_t method;
837
841 std::vector<Point<T>> outline;
842
844 std::unique_ptr<LB<T, W>> balancer;
845
848 std::vector<W> workArray;
849
851 MPI_Comm communicator;
852
854 int nVertices;
855
857 int nNeighbors;
858
861 void calculate_outline() {
862 // calculation only possible if there are vertices defining the domain
863 if (balancer->getVertices().size() > 0) {
864 outline.resize(2);
865 // setup the outline with the first point
866 for (int i = 0; i < balancer->getDimension(); ++i) {
867 outline.at(0) = balancer->getVertices().at(0);
868 outline.at(1) = balancer->getVertices().at(0);
869 }
870 // compare value of each outline point with all vertices to find the
871 // maximum extend of the domain in each dimension
872 for (int i = 1; i < (int)balancer->getVertices().size(); ++i) {
873 Point<T> p = balancer->getVertices().at(i);
874 for (int j = 0; j < balancer->getDimension(); ++j) {
875 outline.at(0)[j] = std::min(outline.at(0)[j], p[j]);
876 outline.at(1)[j] = std::max(outline.at(1)[j], p[j]);
877 }
878 }
879 }
880 }
881
885 void createDirectory(const char *filename) {
886 errno = 0;
887 mode_t mode = S_IRWXU | S_IRWXG; // rwx for user and group
888 if (mkdir(filename, mode)) {
889 if (errno != EEXIST) {
890 throw FilesystemErrorException(__FILE__, __func__, __LINE__,
891 "Could not create output directory.");
892 }
893 }
894 // check permission in case directory already existed, but had wrong
895 // permissions
896 struct stat attr;
897 stat(filename, &attr);
898 if ((attr.st_mode & S_IRWXU) != S_IRWXU) {
900 __FILE__, __func__, __LINE__,
901 "Possibly already existing directory does not have correct "
902 "permissions (rwx) for user");
903 }
904 }
905
909 int estimation_limit;
910
912 std::vector<int> procGridLoc;
914 std::vector<int> procGridDim;
916 bool procGridSet;
918 int procTag;
919
921 MPI_Datatype MPIDataTypeT;
922
924 MPI_Datatype MPIDataTypeW;
925
927 int loadbalancing_step;
928
929#ifdef ALL_VTK_OUTPUT
931 bool vtk_init;
932
934 vtkMPIController *controller;
935
937 vtkMPIController *controller_vertices;
938#endif
939};
940
941} // namespace ALL
942
943#endif
944
945// VIM modeline for indentation
946// vim: sw=2 et
#define ALL_ESTIMATION_LIMIT_DEFAULT
Definition ALL.hpp:77
#define known_unused
Definition ALL_Defines.h:49
std::vector< T > & getNeighborVertices()
Definition ALL.hpp:448
T getGamma()
Definition ALL.hpp:273
void setSysSize(const std::vector< T > &sysSize)
Definition ALL.hpp:474
ALL(const LB_t m, const int d, const std::vector< Point< T > > &inp, const T g)
Definition ALL.hpp:183
void setProcGridParams(const std::vector< int > &loc, const std::vector< int > &size)
Definition ALL.hpp:494
std::vector< Point< T > > & getPrevVertices()
Definition ALL.hpp:415
void setMethodData(const void *data)
Definition ALL.hpp:486
void setCommunicator(const MPI_Comm comm)
Definition ALL.hpp:211
std::vector< int > & getNeighbors()
Definition ALL.hpp:442
void printVTKvertices(const int step)
Definition ALL.hpp:678
void getWork(std::vector< W > &result)
Definition ALL.hpp:433
ALL()
Definition ALL.hpp:106
std::vector< Point< T > > & getVertices()
Definition ALL.hpp:423
W getEstimatedEfficiency()
Definition ALL.hpp:460
W getEfficiency()
Definition ALL.hpp:454
W getWork()
Definition ALL.hpp:437
void setMinDomainSize(const std::vector< T > &minSize)
Definition ALL.hpp:510
void setWork(const W work)
Definition ALL.hpp:281
void setup()
Definition ALL.hpp:292
std::vector< Point< T > > & balance()
Definition ALL.hpp:299
void setVertices(const std::vector< Point< T > > &inp)
Definition ALL.hpp:195
void printVTKoutlines(const int step)
Definition ALL.hpp:519
int getDimension()
Definition ALL.hpp:428
ALL(const LB_t m, const int d, const T g)
Definition ALL.hpp:142
void setWork(const std::vector< W > &work)
Definition ALL.hpp:288
void setGamma(const T g)
Definition ALL.hpp:277
void setProcTag(int tag)
Definition ALL.hpp:504
~ALL()=default
destructor
Definition ALL.hpp:75
LB_t
enum type to describe the different balancing methods
Definition ALL.hpp:80
@ FORCEBASED
unstructured-mesh load balancing
Definition ALL.hpp:86
@ STAGGERED
staggered grid load balancing
Definition ALL.hpp:82
@ VORONOI
voronoi cell based load balancing
Definition ALL.hpp:89
@ UNIMPLEMENTED
Unimplemented load balancing NEVER SUPPORTED.
Definition ALL.hpp:96
@ HISTOGRAM
histogram based load balancing
Definition ALL.hpp:92
@ TENSOR_MAX
tensor based load balancing using maximum values
Definition ALL.hpp:94
@ TENSOR
tensor based load balancing
Definition ALL.hpp:84
Exception to be used for missmatches in dimension for ALL::Point class.