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)");
201 if (status != MPI_CART) {
203 __FILE__, __func__, __LINE__,
204 "Cartesian MPI communicator required, passed communicator is not "
211 this->periodicity.data(), this->local_coords.data());
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)) {
248 std::cout <<
"DEBUG: main_dim: " << main_dim << std::endl;
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
285 std::cout <<
"ALL::ForceBased_LB<T,W>::setup() preparing communicators..."
293#ifdef ALL_DEBUG_ENABLED
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))) {
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))) {
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))) {
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))) {
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))) {
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))) {
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))) {
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))) {
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
378 std::cout <<
"ALL::ForceBased_LB<T,W>::setup() finished computing "
396 T vertex_info[n_vertices * n_vertices * (this->
dimension + 2)];
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] =
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);
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);
602 if (this->
local_coords.at(0) < this->global_dims.at(0) - 1) {
616 if (this->
local_coords.at(1) < this->global_dims.at(1) - 1) {
630 if (this->
local_coords.at(2) < this->global_dims.at(2) - 1) {
639 __FILE__, __func__, __LINE__,
640 "Invalid main dimension provided (numerical value not 0, 1 or 2).");
645 std::cout <<
"DEBUG (main-dim shift): " << this->
localRank <<
" "
646 << remote_work <<
" " << local_work <<
" " << total_width <<
" "
647 << local_shift <<
" " << remote_shift
649 <<
" 0: " << this->
vertices.at(0)[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) {
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] -
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] -
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();
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;
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;
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) {
765 (z == 1 && this->local_coords.at(2) == this->global_dims.at(2) - 1))
769 for (
int y = 0; y < 2; ++y) {
771 (y == 1 && this->local_coords.at(1) == this->global_dims.at(1) - 1))
775 for (
int x = 0; x < 2; ++x) {
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;
787 vertex_shift.at(v * this->dimension + d);