Remove Observable_stat global variables

src/core/CMakeLists.txt

@@ -41,7 +41,6 @@ set(EspressoCore_SRC
     CellStructure.cpp
     PartCfg.cpp
     AtomDecomposition.cpp
-    reduce_observable_stat.cpp
     EspressoSystemStandAlone.cpp
     DomainDecomposition.cpp)
 

src/core/Observable_stat.cpp

@@ -24,32 +24,42 @@
 #include "config.hpp"
 
 #include "bonded_interactions/bonded_interaction_data.hpp"
+#include "communication.hpp"
+#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include <utils/Span.hpp>
+#include <utils/index.hpp>
+
+#include <boost/mpi/collectives/reduce.hpp>
 
 #include <cassert>
 #include <cstddef>
+#include <functional>
 #include <vector>
 
+inline std::size_t max_non_bonded_pairs() {
+  auto const n = static_cast<std::size_t>(max_seen_particle_type);
+  return (n * (n + 1)) / 2;
+}
+
 Observable_stat::Observable_stat(std::size_t chunk_size)
     : m_chunk_size(chunk_size) {
   // number of chunks for different interaction types
-  auto constexpr n_coulomb = 2;
-  auto constexpr n_dipolar = 2;
+  constexpr std::size_t n_coulomb = 2;
+  constexpr std::size_t n_dipolar = 2;
 #ifdef VIRTUAL_SITES
-  auto constexpr n_vs = 1;
+  constexpr std::size_t n_vs = 1;
 #else
-  auto constexpr n_vs = 0;
+  constexpr std::size_t n_vs = 0;
 #endif
   auto const n_bonded = bonded_ia_params.size();
   auto const n_non_bonded = max_non_bonded_pairs();
   constexpr std::size_t n_ext_fields = 1; // reduction over all fields
   constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions
 
-  // resize vector
-  auto const total = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
-                     n_dipolar + n_vs + n_ext_fields;
-  m_data = std::vector<double>(m_chunk_size * total);
+  auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
+                          n_dipolar + n_vs + n_ext_fields;
+  m_data = std::vector<double>(m_chunk_size * n_elements);
 
   // spans for the different contributions
   kinetic = Utils::Span<double>(m_data.data(), m_chunk_size);
@@ -65,3 +75,20 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
       Utils::Span<double>(non_bonded_intra.end(), n_non_bonded * m_chunk_size);
   assert(non_bonded_inter.end() == (m_data.data() + m_data.size()));
 }
+
+Utils::Span<double>
+Observable_stat::non_bonded_contribution(Utils::Span<double> base_pointer,
+                                         int type1, int type2) const {
+  auto const offset = static_cast<std::size_t>(Utils::upper_triangular(
+      std::min(type1, type2), std::max(type1, type2), max_seen_particle_type));
+  return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size};
+}
+
+void Observable_stat::mpi_reduce() {
+  if (comm_cart.rank() == 0) {
+    std::vector<double> temp(m_data);
+    boost::mpi::reduce(comm_cart, temp, m_data, std::plus<>{}, 0);
+  } else {
+    boost::mpi::reduce(comm_cart, m_data, std::plus<>{}, 0);
+  }
+}

src/core/Observable_stat.hpp

@@ -23,7 +23,6 @@
 #include <boost/range/numeric.hpp>
 
 #include <utils/Span.hpp>
-#include <utils/index.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -38,32 +37,14 @@ class Observable_stat {
   /** Number of doubles per data item */
   std::size_t m_chunk_size;
 
-  /** Calculate the maximal number of non-bonded interaction pairs in the
-   *  system.
-   */
-  static std::size_t max_non_bonded_pairs() {
-    extern int max_seen_particle_type;
-    return static_cast<std::size_t>(
-        (max_seen_particle_type * (max_seen_particle_type + 1)) / 2);
-  }
-
   /** Get contribution from a non-bonded interaction */
   Utils::Span<double> non_bonded_contribution(Utils::Span<double> base_pointer,
-                                              int type1, int type2) const {
-    extern int max_seen_particle_type;
-    int offset = Utils::upper_triangular(
-        std::min(type1, type2), std::max(type1, type2), max_seen_particle_type);
-    return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size};
-  }
+                                              int type1, int type2) const;
 
 public:
   explicit Observable_stat(std::size_t chunk_size);
 
-  auto chunk_size() const { return m_chunk_size; }
-  Utils::Span<double> data_() { return {m_data.data(), m_data.size()}; }
-  Utils::Span<const double> data_() const {
-    return {m_data.data(), m_data.size()};
-  }
+  auto get_chunk_size() const { return m_chunk_size; }
 
   /** Accumulate values.
    *  @param acc    Initial value for the accumulator.
@@ -97,7 +78,7 @@ class Observable_stat {
   Utils::Span<double> coulomb;
   /** Contribution(s) from dipolar interactions. */
   Utils::Span<double> dipolar;
-  /** Contribution(s) from virtual sites (accumulated). */
+  /** Contribution from virtual sites (accumulated). */
   Utils::Span<double> virtual_sites;
   /** Contribution from external fields (accumulated). */
   Utils::Span<double> external_fields;
@@ -107,17 +88,16 @@ class Observable_stat {
   Utils::Span<double> non_bonded_inter;
 
   /** Get contribution from a bonded interaction */
-  Utils::Span<double> bonded_contribution(int bond_id) {
-    return Utils::Span<double>(bonded.data() + m_chunk_size * bond_id,
-                               m_chunk_size);
+  Utils::Span<double> bonded_contribution(int bond_id) const {
+    auto const offset = m_chunk_size * static_cast<std::size_t>(bond_id);
+    return Utils::Span<double>(bonded.data() + offset, m_chunk_size);
   }
 
   void add_non_bonded_contribution(int type1, int type2,
                                    Utils::Span<const double> data) {
-    auto const dest =
-        (type1 == type2)
-            ? non_bonded_contribution(non_bonded_intra, type1, type2)
-            : non_bonded_contribution(non_bonded_inter, type1, type2);
+    assert(data.size() == m_chunk_size);
+    auto const source = (type1 == type2) ? non_bonded_intra : non_bonded_inter;
+    auto const dest = non_bonded_contribution(source, type1, type2);
 
     boost::transform(dest, data, dest.begin(), std::plus<>{});
   }
@@ -137,6 +117,9 @@ class Observable_stat {
                                                     int type2) const {
     return non_bonded_contribution(non_bonded_inter, type1, type2);
   }
+
+  /** MPI reduction. */
+  void mpi_reduce();
 };
 
 #endif // ESPRESSO_OBSERVABLE_STAT_HPP

src/core/constraints/ShapeBasedConstraint.cpp

@@ -20,6 +20,7 @@
 #include "ShapeBasedConstraint.hpp"
 
 #include "BoxGeometry.hpp"
+#include "Observable_stat.hpp"
 #include "communication.hpp"
 #include "config.hpp"
 #include "dpd.hpp"

src/core/energy.cpp

@@ -33,25 +33,24 @@
 #include "forces.hpp"
 #include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "reduce_observable_stat.hpp"
 
 #include "short_range_loop.hpp"
 
 #include "electrostatics_magnetostatics/coulomb.hpp"
 #include "electrostatics_magnetostatics/dipole.hpp"
 
+#include <memory>
+
 ActorList energyActors;
 
-/** Energy of the system */
-Observable_stat obs_energy{1};
+static std::shared_ptr<Observable_stat> calculate_energy_local() {
 
-Observable_stat const &get_obs_energy() { return obs_energy; }
+  auto obs_energy_ptr = std::make_shared<Observable_stat>(1);
 
-void energy_calc(const double time) {
   if (!interactions_sanity_checks())
-    return;
+    return obs_energy_ptr;
 
-  obs_energy = Observable_stat{1};
+  auto &obs_energy = *obs_energy_ptr;
 
 #ifdef CUDA
   clear_energy_on_GPU();
@@ -65,12 +64,15 @@ void energy_calc(const double time) {
 
   on_observable_calc();
 
-  for (auto const &p : cell_structure.local_particles()) {
+  auto const local_parts = cell_structure.local_particles();
+
+  for (auto const &p : local_parts) {
     obs_energy.kinetic[0] += calc_kinetic_energy(p);
   }
 
   short_range_loop(
-      [](Particle &p1, int bond_id, Utils::Span<Particle *> partners) {
+      [&obs_energy](Particle const &p1, int bond_id,
+                    Utils::Span<Particle *> partners) {
         auto const &iaparams = bonded_ia_params[bond_id];
         auto const result = calc_bonded_energy(iaparams, p1, partners);
         if (result) {
@@ -79,16 +81,23 @@ void energy_calc(const double time) {
         }
         return true;
       },
-      [](Particle const &p1, Particle const &p2, Distance const &d) {
+      [&obs_energy](Particle const &p1, Particle const &p2, Distance const &d) {
         add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2,
                                    obs_energy);
       },
       maximal_cutoff(), maximal_cutoff_bonded());
 
-  calc_long_range_energies(cell_structure.local_particles());
+#ifdef ELECTROSTATICS
+  /* calculate k-space part of electrostatic interaction. */
+  obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(local_parts);
+#endif
+
+#ifdef DIPOLES
+  /* calculate k-space part of magnetostatic interaction. */
+  obs_energy.dipolar[1] = Dipole::calc_energy_long_range(local_parts);
+#endif
 
-  auto local_parts = cell_structure.local_particles();
-  Constraints::constraints.add_energy(local_parts, time, obs_energy);
+  Constraints::constraints.add_energy(local_parts, get_sim_time(), obs_energy);
 
 #ifdef CUDA
   auto const energy_host = copy_energy_from_GPU();
@@ -98,37 +107,22 @@ void energy_calc(const double time) {
     obs_energy.dipolar[1] += energy_host.dipolar;
 #endif
 
-  /* gather data */
-  auto obs_energy_res = reduce(comm_cart, obs_energy);
-  if (obs_energy_res) {
-    std::swap(obs_energy, *obs_energy_res);
-  }
+  obs_energy.mpi_reduce();
+  return obs_energy_ptr;
 }
 
-void update_energy_local(int, int) { energy_calc(get_sim_time()); }
+REGISTER_CALLBACK_MASTER_RANK(calculate_energy_local)
 
-REGISTER_CALLBACK(update_energy_local)
-
-void update_energy() { mpi_call_all(update_energy_local, -1, -1); }
-
-void calc_long_range_energies(const ParticleRange &particles) {
-#ifdef ELECTROSTATICS
-  /* calculate k-space part of electrostatic interaction. */
-  obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(particles);
-#endif
-
-#ifdef DIPOLES
-  /* calculate k-space part of magnetostatic interaction. */
-  obs_energy.dipolar[1] = Dipole::calc_energy_long_range(particles);
-#endif
+std::shared_ptr<Observable_stat> calculate_energy() {
+  return mpi_call(Communication::Result::master_rank, calculate_energy_local);
 }
 
 double calculate_current_potential_energy_of_system() {
-  update_energy();
-  return obs_energy.accumulate(-obs_energy.kinetic[0]);
+  auto const obs_energy = calculate_energy();
+  return obs_energy->accumulate(-obs_energy->kinetic[0]);
 }
 
 double observable_compute_energy() {
-  update_energy();
-  return obs_energy.accumulate(0);
+  auto const obs_energy = calculate_energy();
+  return obs_energy->accumulate(0);
 }

src/core/energy.hpp

@@ -28,22 +28,14 @@
 #define _ENERGY_H
 
 #include "Observable_stat.hpp"
-#include "ParticleRange.hpp"
 #include "actor/ActorList.hpp"
 
+#include <memory>
+
 extern ActorList energyActors;
 
 /** Parallel energy calculation. */
-void energy_calc(double time);
-
-/** Run @ref energy_calc in parallel. */
-void update_energy();
-
-/** Return the energy observable. */
-Observable_stat const &get_obs_energy();
-
-/** Calculate long-range energies (P3M, ...). */
-void calc_long_range_energies(const ParticleRange &particles);
+std::shared_ptr<Observable_stat> calculate_energy();
 
 /** Calculate the total energy of the system. */
 double calculate_current_potential_energy_of_system();