162 n_vertices = this->vertices.size();
163 vertex_neighbors.resize(n_vertices * 8);
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;
176 __FILE__, __func__, __LINE__,
177 "Invalid data type for boundaries given (T)");
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;
191 "Invalid data type for work given (W)");
195 MPI_Comm_group(this->globalComm, &all);
199 MPI_Topo_test(this->globalComm, &status);
201 if (status != MPI_CART) {
203 __FILE__, __func__, __LINE__,
204 "Cartesian MPI communicator required, passed communicator is not "
210 MPI_Cart_get(this->globalComm, this->dimension, this->global_dims.data(),
211 this->periodicity.data(), this->local_coords.data());
214 MPI_Cart_rank(this->globalComm, this->local_coords.data(), &this->localRank);
225 if (this->global_dims.at(2) >= this->global_dims.at(1)) {
226 if (this->global_dims.at(2) >= this->global_dims.at(0)) {
236 if (this->global_dims.at(1) >= this->global_dims.at(0)) {
247 if (this->localRank == 0)
248 std::cout <<
"DEBUG: main_dim: " << main_dim << std::endl;
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]),
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..."
288 std::vector<int> dim_vert(this->global_dims);
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..."
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;
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) {
308 for (
auto &a : affected)
311 for (
auto &vn : vertex_neighbors)
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))) {
317 v_neighbors[0] = this->localRank;
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))) {
323 v_neighbors[1] = this->localRank;
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))) {
329 v_neighbors[2] = this->localRank;
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))) {
335 v_neighbors[3] = this->localRank;
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))) {
341 v_neighbors[4] = this->localRank;
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))) {
347 v_neighbors[5] = this->localRank;
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))) {
353 v_neighbors[6] = this->localRank;
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))) {
359 v_neighbors[7] = this->localRank;
362 MPI_Allreduce(MPI_IN_PLACE, v_neighbors, 8, MPI_INT, MPI_MAX,
365 for (
int v = 0; v < n_vertices; ++v) {
367 for (
int n = 0; n < 8; ++n) {
368 vertex_neighbors.at(8 * v + n) = v_neighbors[n];
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 "
391 this->prevVertices = this->vertices;
396 T vertex_info[n_vertices * n_vertices * (this->dimension + 2)];
397 T center[this->dimension];
400 T
known_unused local_info[(this->dimension + 2) * n_vertices];
402 for (
int i = 0; i < n_vertices * n_vertices * (this->dimension + 2); ++i)
405 for (
int d = 0; d < (this->dimension + 2) * n_vertices; ++d)
406 local_info[d] = (T)0;
408 for (
int d = 0; d < this->dimension; ++d)
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);
416 local_info[v * (this->dimension + 2) + this->dimension] =
419 for (
int v = 0; v < n_vertices; ++v) {
420 for (
int d = 0; d < this->dimension; ++d) {
422 local_info[v * (this->dimension + 2) + d] =
423 center[d] - (this->prevVertices.at(v))[d];
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)));
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)));
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)));
509 __FILE__, __func__, __LINE__,
510 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
519 T
known_unused shift_vectors[this->dimension * n_vertices];
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;
539 MPI_Cart_shift(this->globalComm, main_dim, 1, &main_down, &main_up);
541 W local_work = this->work.at(0);
545 MPI_Allreduce(MPI_IN_PLACE, &local_work, 1, MPIDataTypeW, MPI_SUM,
550 MPI_Sendrecv(&local_work, 1, MPIDataTypeW, main_down, 1010, &remote_work, 1,
551 MPIDataTypeW, main_up, 1010, this->globalComm, &state);
558 local_width = this->prevVertices.at(1)[0] - this->prevVertices.at(0)[0];
561 local_width = this->prevVertices.at(2)[1] - this->prevVertices.at(0)[1];
564 local_width = this->prevVertices.at(4)[2] - this->prevVertices.at(0)[2];
568 __FILE__, __func__, __LINE__,
569 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
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);
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
590 MPI_Sendrecv(&local_shift, 1, MPIDataTypeT, main_up, 2020, &remote_shift, 1,
591 MPIDataTypeT, main_down, 2020, this->globalComm, &state);
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;
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;
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;
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;
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;
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;
639 __FILE__, __func__, __LINE__,
640 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
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;
654 int last_vertex = (n_vertices - 1) * n_vertices * (this->dimension + 2);
655 int vertex_offset = this->dimension + 2;
658 T max_shift = std::max(
659 0.49 * vertex_info[last_vertex + 0 * vertex_offset + this->dimension + 1],
661 for (
int v = 1; v < n_vertices; ++v) {
662 max_shift = std::min(
665 vertex_info[last_vertex + v * vertex_offset + this->dimension + 1]);
670 for (
int v = 0; v < n_vertices; ++v) {
671 avg_work += vertex_info[last_vertex + v * vertex_offset + this->dimension];
673 avg_work /= (T)n_vertices;
676 std::vector<T> vertex_shift(n_vertices * this->dimension, (T)0.0);
677 for (
auto &sv : vertex_shift)
679 int shift_offset = (n_vertices - 1) * this->dimension;
681 for (
int v = 0; v < n_vertices; ++v) {
682 for (
int d = 0; d < this->dimension; ++d) {
683 if (this->localRank == 0)
685 <<
"DEBUG (shift x): " << v <<
", " << d <<
" "
686 << (vertex_info[last_vertex + v * vertex_offset + this->dimension] -
688 ((vertex_info[last_vertex + v * vertex_offset +
689 this->dimension] > avg_work)
693 (vertex_info[last_vertex + v * vertex_offset + d] -
694 this->prevVertices.at(v)[d])
695 <<
" => " << vertex_shift.at(shift_offset + d) <<
697 " max: " << max_shift
703 << vertex_info[last_vertex + v * vertex_offset + this->dimension]
706 << ((vertex_info[last_vertex + v * vertex_offset +
707 this->dimension] > avg_work)
712 vertex_shift.at(shift_offset + d) +=
713 (vertex_info[last_vertex + v * vertex_offset + this->dimension] -
715 ((vertex_info[last_vertex + v * vertex_offset + this->dimension] >
720 (vertex_info[last_vertex + v * vertex_offset + d] -
721 this->prevVertices.at(v)[d]);
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;
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;
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) <<
" ";
748 std::cout << std::endl;
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) <<
" ";
758 std::cout << std::endl;
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))
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))
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))
781 for (
int d = 0; d < this->dimension; ++d) {
784 int v = 4 * z + 2 * y + x;
786 this->vertices.at(v)[d] = this->prevVertices.at(v)[d] +
787 vertex_shift.at(v * this->dimension + d);
789 this->vertices.at(v)[d] = this->prevVertices.at(v)[d];