observables: removed last MPI callbacks from observables. To capture …

…the runtime errors through parallel_try_catches of the operator() of all the observables, we are not allowed to use any MPI-callbacks anymore, so no MPI-deadlock is happing, because of a conflict of the order of communcations.

src/core/grid_based_algorithms/lb_interface.cpp

@@ -150,17 +150,18 @@ static Utils::VectorXd<9> get_pressure_tensor() {
  return lb_walberla()->get_pressure_tensor();
 }
 
-REGISTER_CALLBACK_REDUCTION(get_pressure_tensor, std::plus<>())
-
 } // namespace Walberla
 #endif // WALBERLA
 
-Utils::VectorXd<9> const get_pressure_tensor() {
+Utils::VectorXd<9> const
+get_pressure_tensor(boost::mpi::communicator const &comm) {
  if (lattice_switch == ActiveLB::WALBERLA_LB) {
 #ifdef WALBERLA
- return ::Communication::mpiCallbacks().call(
- ::Communication::Result::reduction, std::plus<>(),
- Walberla::get_pressure_tensor);
+ auto const local_pressure_tensor = Walberla::get_pressure_tensor();
+ std::remove_const_t<decltype(local_pressure_tensor)> pressure_tensor;
+ boost::mpi::reduce(comm, local_pressure_tensor, pressure_tensor,
+ std::plus<>(), 0);
+ return pressure_tensor;
 #endif
  }
  throw NoLBActive();

src/core/grid_based_algorithms/lb_interface.hpp

@@ -109,7 +109,8 @@ double get_lattice_speed();
  * over all nodes and dividing by the number of nodes.
  * Returns the lower triangle of the LB pressure tensor.
  */
-Utils::VectorXd<9> const get_pressure_tensor();
+Utils::VectorXd<9> const
+get_pressure_tensor(boost::mpi::communicator const &comm);
 
 Utils::Vector3d calc_fluid_momentum();
 

src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp

@@ -26,7 +26,7 @@
 #include <utils/Span.hpp>
 #include <utils/math/coordinate_transformation.hpp>
 
-#include <boost/mpi/collectives/gather.hpp>
+#include <boost/mpi/collectives/all_gather.hpp>
 
 #include <vector>
 
@@ -47,11 +47,7 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate(
  }
 
  std::vector<std::vector<pos_type>> global_folded_positions;
- boost::mpi::gather(comm, local_folded_positions, global_folded_positions, 0);
-
- if (comm.rank() != 0) {
- return {};
- }
+ boost::mpi::all_gather(comm, local_folded_positions, global_folded_positions);
 
  Utils::CylindricalHistogram<double, 3> histogram(n_bins(), limits());
  // First collect all positions (since we want to call the LB function to
@@ -60,18 +56,23 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate(
  for (auto const &pos_vec : global_folded_positions) {
  for (auto const &pos : pos_vec) {
  auto const v =
- LB::get_interpolated_velocity(pos) * LB::get_lattice_speed();
- auto const flux_dens = LB::get_interpolated_density(pos) * v;
+ LB::get_interpolated_velocity(comm, pos) * LB::get_lattice_speed();
+ auto const flux_dens = LB::get_interpolated_density(comm, pos) * v;
 
- histogram.update(Utils::transform_coordinate_cartesian_to_cylinder(
- pos - transform_params->center(),
- transform_params->axis(),
- transform_params->orientation()),
- Utils::transform_vector_cartesian_to_cylinder(
- flux_dens, transform_params->axis(),
- pos - transform_params->center()));
+ if (comm.rank() == 0) {
+ histogram.update(Utils::transform_coordinate_cartesian_to_cylinder(
+ pos - transform_params->center(),
+ transform_params->axis(),
+ transform_params->orientation()),
+ Utils::transform_vector_cartesian_to_cylinder(
+ flux_dens, transform_params->axis(),
+ pos - transform_params->center()));
+ }
  }
  }
+ if (comm.rank() != 0) {
+ return {};
+ }
 
  // normalize by number of hits per bin
  auto hist_tmp = histogram.get_histogram();

src/core/observables/CylindricalLBVelocityProfile.cpp

@@ -32,20 +32,22 @@ namespace Observables {
 
 std::vector<double> CylindricalLBVelocityProfile::operator()(
  boost::mpi::communicator const &comm) const {
- if (comm.rank() != 0) {
- return {};
- }
-
  Utils::CylindricalHistogram<double, 3> histogram(n_bins(), limits());
  for (auto const &p : sampling_positions) {
  auto const velocity =
- LB::get_interpolated_velocity(p) * LB::get_lattice_speed();
- auto const pos_shifted = p - transform_params->center();
- auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder(
- pos_shifted, transform_params->axis(), transform_params->orientation());
- histogram.update(pos_cyl,
- Utils::transform_vector_cartesian_to_cylinder(
- velocity, transform_params->axis(), pos_shifted));
+ LB::get_interpolated_velocity(comm, p) * LB::get_lattice_speed();
+
+ if (comm.rank() == 0) {
+ auto const pos_shifted = p - transform_params->center();
+ auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder(
+ pos_shifted, transform_params->axis(), transform_params->orientation());
+ histogram.update(pos_cyl,
+ Utils::transform_vector_cartesian_to_cylinder(
+ velocity, transform_params->axis(), pos_shifted));
+ }
+ }
+ if (comm.rank() != 0) {
+ return {};
  }
  auto hist_data = histogram.get_histogram();
  auto const tot_count = histogram.get_tot_count();

src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp

@@ -26,6 +26,8 @@
 #include <utils/Span.hpp>
 #include <utils/math/coordinate_transformation.hpp>
 
+#include <boost/mpi/collectives/all_gather.hpp>
+
 #include <cstddef>
 #include <vector>
 
@@ -45,27 +47,28 @@ std::vector<double> CylindricalLBVelocityProfileAtParticlePositions::evaluate(
  }
 
  std::vector<std::vector<pos_type>> global_folded_positions;
- boost::mpi::gather(comm, local_folded_positions, global_folded_positions, 0);
-
- if (comm.rank() != 0) {
- return {};
- }
+ boost::mpi::all_gather(comm, local_folded_positions, global_folded_positions);
 
  Utils::CylindricalHistogram<double, 3> histogram(n_bins(), limits());
 
  for (auto const &vec : global_folded_positions) {
  for (auto const &pos : vec) {
  auto const v =
- LB::get_interpolated_velocity(pos) * LB::get_lattice_speed();
+ LB::get_interpolated_velocity(comm, pos) * LB::get_lattice_speed();
 
- histogram.update(
- Utils::transform_coordinate_cartesian_to_cylinder(
- pos - transform_params->center(), transform_params->axis(),
- transform_params->orientation()),
- Utils::transform_vector_cartesian_to_cylinder(
- v, transform_params->axis(), pos - transform_params->center()));
+ if (comm.rank() == 0) {
+ histogram.update(
+ Utils::transform_coordinate_cartesian_to_cylinder(
+ pos - transform_params->center(), transform_params->axis(),
+ transform_params->orientation()),
+ Utils::transform_vector_cartesian_to_cylinder(
+ v, transform_params->axis(), pos - transform_params->center()));
+ }
  }
  }
+ if (comm.rank() != 0) {
+ return {};
+ }
 
  // normalize by number of hits per bin
  auto hist_tmp = histogram.get_histogram();

src/core/observables/LBFluidPressureTensor.hpp

@@ -33,12 +33,9 @@ class LBFluidPressureTensor : public Observable {
  std::vector<std::size_t> shape() const override { return {3, 3}; }
  std::vector<double>
  operator()(boost::mpi::communicator const &comm) const override {
- if (comm.rank() != 0) {
- return {};
- }
  auto const unit_conversion =
  1. / (LB::get_agrid() * Utils::sqr(LB::get_tau()));
- auto const tensor = LB::get_pressure_tensor() * unit_conversion;
+ auto const tensor = LB::get_pressure_tensor(comm) * unit_conversion;
  return tensor.as_vector();
  }
 };

src/core/observables/LBVelocityProfile.cpp

@@ -31,16 +31,16 @@ namespace Observables {
 
 std::vector<double>
 LBVelocityProfile::operator()(boost::mpi::communicator const &comm) const {
- // TODO: make a copy of LB::get_interpolated_velocity without MPI-callback and
- // replace this one
- if (comm.rank() != 0) {
- return {};
- }
-
+ // TODO: split this histogram according to domain-decomposition
  Utils::Histogram<double, 3> histogram(n_bins(), limits());
  for (auto const &p : sampling_positions) {
- const auto v = LB::get_interpolated_velocity(p) * LB::get_lattice_speed();
- histogram.update(p, v);
+ const auto v = LB::get_interpolated_velocity(comm, p) * LB::get_lattice_speed();
+ if (comm.rank() == 0) {
+ histogram.update(p, v);
+ }
+ }
+ if (comm.rank() != 0) {
+ return {};
  }
  auto hist_tmp = histogram.get_histogram();
  auto const tot_count = histogram.get_tot_count();

src/core/observables/RDF.cpp

@@ -29,58 +29,78 @@
 #include <utils/math/int_pow.hpp>
 
 #include <boost/range/algorithm/transform.hpp>
+#include <boost/range/combine.hpp>
 
 #include <cmath>
 #include <functional>
 #include <memory>
+#include <tuple>
 #include <vector>
 
 namespace Observables {
 std::vector<double>
 RDF::operator()(boost::mpi::communicator const &comm) const {
- if (comm.rank() != 0) {
- return {};
- }
-
- auto particles1 = old_fetch_particles(ids1());
- std::vector<const Particle *> particles_ptrs1(particles1.size());
- boost::transform(particles1, particles_ptrs1.begin(),
- [](auto const &p) { return std::addressof(p); });
+ auto const &local_particles_1 = fetch_particles(ids1());
 
  if (ids2().empty()) {
- return this->evaluate(comm, particles_ptrs1, {});
+ return this->evaluate(comm, local_particles_1, {}, {});
  }
 
- auto particles2 = old_fetch_particles(ids2());
- std::vector<const Particle *> particles_ptrs2(particles2.size());
- boost::transform(particles2, particles_ptrs2.begin(),
- [](auto const &p) { return std::addressof(p); });
- return this->evaluate(comm, particles_ptrs1, particles_ptrs2);
+ auto const &local_particles_2 = fetch_particles(ids2());
+
+ return this->evaluate(comm, local_particles_1, local_particles_2, {});
 }
 
 std::vector<double>
 RDF::evaluate(boost::mpi::communicator const &comm,
- Utils::Span<const Particle *const> particles1,
- Utils::Span<const Particle *const> particles2) const {
+ ParticleReferenceRange const &local_particles_1,
+ ParticleReferenceRange const &local_particles_2,
+ const ParticleObservables::traits<Particle> &traits) const {
+ auto const positions_1 = detail::get_all_particle_positions(
+ comm, local_particles_1, ids1(), traits, false);
+ auto const positions_2 = detail::get_all_particle_positions(
+ comm, local_particles_2, ids2(), traits, false);
+
+ if (comm.rank() != 0) {
+ return {};
+ }
+
  auto const bin_width = (max_r - min_r) / static_cast<double>(n_r_bins);
  auto const inv_bin_width = 1.0 / bin_width;
  std::vector<double> res(n_values(), 0.0);
  long int cnt = 0;
- auto op = [this, inv_bin_width, &cnt, &res](const Particle *const p1,
- const Particle *const p2) {
- auto const dist = box_geo.get_mi_vector(p1->pos(), p2->pos()).norm();
+ auto op = [this, inv_bin_width, &cnt, &res](auto const &pos1,
+ auto const &pos2) {
+ auto const dist = box_geo.get_mi_vector(pos1, pos2).norm();
  if (dist > min_r && dist < max_r) {
  auto const ind =
  static_cast<int>(std::floor((dist - min_r) * inv_bin_width));
  res[ind]++;
  }
  cnt++;
  };
- if (particles2.empty()) {
- Utils::for_each_pair(particles1, op);
+
+ if (local_particles_2.empty()) {
+ Utils::for_each_pair(positions_1, op);
  } else {
- auto cmp = std::not_equal_to<const Particle *const>();
- Utils::for_each_cartesian_pair_if(particles1, particles2, op, cmp);
+ auto const combine_1 = boost::combine(ids1(), positions_1);
+ auto const combine_2 = boost::combine(ids2(), positions_2);
+
+ auto op2 = [this, inv_bin_width, &cnt, &res, &op](auto const &it1,
+ auto const &it2) {
+ auto const &[id1, pos1] = it1;
+ auto const &[id2, pos2] = it2;
+
+ op(pos1, pos2);
+ };
+
+ auto cmp = [](auto const &it1, auto const &it2) {
+ auto const &[id1, pos1] = it1;
+ auto const &[id2, pos2] = it2;
+
+ return id1 != id2;
+ };
+ Utils::for_each_cartesian_pair_if(combine_1, combine_2, op2, cmp);
  }
  if (cnt == 0)
  return res;
@@ -97,4 +117,44 @@ RDF::evaluate(boost::mpi::communicator const &comm,
 
  return res;
 }
+
+// std::vector<double>
+// RDF::evaluate(boost::mpi::communicator const &comm,
+// Utils::Span<const Particle *const> particles1,
+// Utils::Span<const Particle *const> particles2) const {
+// auto const bin_width = (max_r - min_r) / static_cast<double>(n_r_bins);
+// auto const inv_bin_width = 1.0 / bin_width;
+// std::vector<double> res(n_values(), 0.0);
+// long int cnt = 0;
+// auto op = [this, inv_bin_width, &cnt, &res](const Particle *const p1,
+// const Particle *const p2) {
+// auto const dist = box_geo.get_mi_vector(p1->pos(), p2->pos()).norm();
+// if (dist > min_r && dist < max_r) {
+// auto const ind =
+// static_cast<int>(std::floor((dist - min_r) * inv_bin_width));
+// res[ind]++;
+// }
+// cnt++;
+// };
+// if (particles2.empty()) {
+// Utils::for_each_pair(particles1, op);
+// } else {
+// auto cmp = std::not_equal_to<const Particle *const>();
+// Utils::for_each_cartesian_pair_if(particles1, particles2, op, cmp);
+// }
+// if (cnt == 0)
+// return res;
+// // normalization
+// auto const volume = box_geo.volume();
+// for (int i = 0; i < n_r_bins; ++i) {
+// auto const r_in = i * bin_width + min_r;
+// auto const r_out = r_in + bin_width;
+// auto const bin_volume =
+// (4.0 / 3.0) * Utils::pi() *
+// (Utils::int_pow<3>(r_out) - Utils::int_pow<3>(r_in));
+// res[i] *= volume / (bin_volume * static_cast<double>(cnt));
+// }
+
+// return res;
+// }
 } // namespace Observables

src/core/observables/RDF.hpp

@@ -21,7 +21,8 @@
 #define OBSERVABLES_RDF_HPP
 
 #include "Observable.hpp"
-#include "Particle.hpp"
+
+#include "PidObservable.hpp"
 
 #include <utils/Span.hpp>
 
@@ -40,10 +41,15 @@ class RDF : public Observable {
  /** Identifiers of the distant particles */
  std::vector<int> m_ids2;
 
- virtual std::vector<double>
+ // virtual std::vector<double>
+ // evaluate(boost::mpi::communicator const &comm,
+ // Utils::Span<const Particle *const> particles1,
+ // Utils::Span<const Particle *const> particles2) const;
+ std::vector<double>
  evaluate(boost::mpi::communicator const &comm,
- Utils::Span<const Particle *const> particles1,
- Utils::Span<const Particle *const> particles2) const;
+ ParticleReferenceRange const &local_particles_1,
+ ParticleReferenceRange const &local_particles_2,
+ const ParticleObservables::traits<Particle> &traits) const;
 
 public:
  // Range of the profile.