Remove more global variables (#4783)

Partial fix for #2628

Description of changes:
- encapsulate cell structure, box geometry and local box in the `System` class
- gather MPI-related global variables in a global struct
- rewrite position folding/unfolding free functions as member functions of `BoxGeometry`

I took this opportunity to cleanup and modernize code, by e.g. replacing hard-coded loops by `Utils::Vector` operations, improving type safety (no mixed signed and unsigned int in arithmetic expression), using more expressive variable names, passing global variables by function argument, adding a collision detection test case, etc.

samples/electrophoresis.py

@@ -33,9 +33,6 @@
 
 logging.basicConfig(level=logging.INFO)
 
-# Use a fixed int for a deterministic behavior
-np.random.seed()
-
 required_features = ["P3M", "EXTERNAL_FORCES", "WCA"]
 espressomd.assert_features(required_features)
 
@@ -133,7 +130,8 @@
 
 # activate electrostatics
 #############################################################
-p3m = espressomd.electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
+p3m_params = {"prefactor": 1., "accuracy": 1e-2}
+p3m = espressomd.electrostatics.P3M(**p3m_params)
 system.electrostatics.solver = p3m
 
 # apply external force (external electric field)
@@ -160,17 +158,15 @@
     obs=obs_bond_length, delta_N=1)
 system.auto_update_accumulators.add(acc_bond_length)
 
-# data storage for python analysis
-#############################################################
-pos = np.full((N_SAMPLES, N_MONOMERS, 3), np.nan)
+obs_pos = espressomd.observables.ParticlePositions(ids=range(N_MONOMERS))
+acc_pos = espressomd.accumulators.TimeSeries(obs=obs_pos, delta_N=100)
+system.auto_update_accumulators.add(acc_pos)
 
 # sampling Loop
 #############################################################
-for i in range(N_SAMPLES):
-    if i % 100 == 0:
-        logging.info(f"\rsampling: {i:4d}")
-    system.integrator.run(N_INT_STEPS)
-    pos[i] = system.part.by_ids(range(N_MONOMERS)).pos
+for i in range(N_SAMPLES // 100):
+    logging.info(f"\rsampling: {100 * i:4d}")
+    system.integrator.run(100 * N_INT_STEPS)
 
 logging.info("\nsampling finished!\n")
 
@@ -179,6 +175,7 @@
 
 # calculate center of mass (COM) and its velocity
 #############################################################
+pos = acc_pos.time_series()
 COM = pos.sum(axis=1) / N_MONOMERS
 COM_v = (COM[1:] - COM[:-1]) / (N_INT_STEPS * system.time_step)
 
@@ -220,7 +217,7 @@ def exponential(x, lp):
 # ...second by using observables
 
 
-def persistence_length_obs(
+def persistence_length_from_obs(
         acc_bond_length, acc_persistence_angles, exponential):
     bond_lengths_obs = np.array(acc_bond_length.mean())
     sampling_positions_obs = np.insert(
@@ -233,7 +230,7 @@ def persistence_length_obs(
     return sampling_positions_obs, cos_thetas_obs, opt_obs[0]
 
 
-sampling_positions_obs, cos_thetas_obs, persistence_length_obs = persistence_length_obs(
+sampling_positions_obs, cos_thetas_obs, persistence_length_obs = persistence_length_from_obs(
     acc_bond_length, acc_persistence_angles, exponential)
 logging.info(f"persistence length (observables): {persistence_length_obs}")
 

src/core/BoxGeometry.hpp

@@ -16,8 +16,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef ESPRESSO_SRC_CORE_BOX_GEOMETRY_HPP
-#define ESPRESSO_SRC_CORE_BOX_GEOMETRY_HPP
+
+#pragma once
 
 #include "algorithm/periodic_fold.hpp"
 #include "lees_edwards/LeesEdwardsBC.hpp"
@@ -29,6 +29,7 @@
 #include <cassert>
 #include <cmath>
 #include <limits>
+#include <stdexcept>
 #include <utility>
 
 namespace detail {
@@ -68,6 +69,28 @@ template <typename T> T get_mi_coord(T a, T b, T box_length, bool periodic) {
   return get_mi_coord(a, b, box_length, 1. / box_length, 0.5 * box_length,
                       periodic);
 }
+
+/** @brief Calculate image box shift vector.
+ *  @param image_box  image box offset
+ *  @param box        box length
+ *  @return Image box coordinates.
+ */
+inline auto image_shift(Utils::Vector3i const &image_box,
+                        Utils::Vector3d const &box) {
+  return hadamard_product(image_box, box);
+}
+
+/** @brief Unfold particle coordinates to image box.
+ *  @param pos        coordinate to unfold
+ *  @param image_box  image box offset
+ *  @param box        box length
+ *  @return Unfolded coordinates.
+ */
+inline auto unfolded_position(Utils::Vector3d const &pos,
+                              Utils::Vector3i const &image_box,
+                              Utils::Vector3d const &box) {
+  return pos + image_shift(image_box, box);
+}
 } // namespace detail
 
 enum class BoxType { CUBOID = 0, LEES_EDWARDS = 1 };
@@ -120,7 +143,7 @@ class BoxGeometry {
    * @return true iff periodic in direction.
    */
   constexpr bool periodic(unsigned coord) const {
-    assert(coord <= 2);
+    assert(coord <= 2u);
     return m_periodic[coord];
   }
 
@@ -192,11 +215,11 @@ class BoxGeometry {
       auto a_tmp = a;
       auto b_tmp = b;
       a_tmp[shear_plane_normal] = Algorithm::periodic_fold(
-          a_tmp[shear_plane_normal], length()[shear_plane_normal]);
+          a_tmp[shear_plane_normal], m_length[shear_plane_normal]);
       b_tmp[shear_plane_normal] = Algorithm::periodic_fold(
-          b_tmp[shear_plane_normal], length()[shear_plane_normal]);
-      return lees_edwards_bc().distance(a_tmp - b_tmp, length(), length_half(),
-                                        length_inv(), m_periodic);
+          b_tmp[shear_plane_normal], m_length[shear_plane_normal]);
+      return lees_edwards_bc().distance(a_tmp - b_tmp, m_length, m_length_half,
+                                        m_length_inv, m_periodic);
     }
     assert(type() == BoxType::CUBOID);
     return {get_mi_coord(a[0], b[0], 0), get_mi_coord(a[1], b[1], 1),
@@ -238,61 +261,54 @@ class BoxGeometry {
     }
     return ret;
   }
-};
-
-/** @brief Fold a coordinate to primary simulation box.
- *  @param pos        coordinate to fold
- *  @param image_box  image box offset
- *  @param length     box length
- */
-inline std::pair<double, int> fold_coordinate(double pos, int image_box,
-                                              double const &length) {
-  std::tie(pos, image_box) = Algorithm::periodic_fold(pos, image_box, length);
-
-  if ((image_box == std::numeric_limits<int>::min()) ||
-      (image_box == std::numeric_limits<int>::max())) {
-    throw std::runtime_error(
-        "Overflow in the image box count while folding a particle coordinate "
-        "into the primary simulation box. Maybe a particle experienced a "
-        "huge force.");
-  }
 
-  return {pos, image_box};
-}
-
-/** @brief Fold particle coordinates to primary simulation box.
- *  Lees-Edwards offset is ignored.
- *  @param[in,out] pos        coordinate to fold
- *  @param[in,out] image_box  image box offset
- *  @param[in] box            box parameters (side lengths, periodicity)
- */
-inline void fold_position(Utils::Vector3d &pos, Utils::Vector3i &image_box,
-                          const BoxGeometry &box) {
-  for (unsigned int i = 0; i < 3; i++) {
-    if (box.periodic(i)) {
-      std::tie(pos[i], image_box[i]) =
-          fold_coordinate(pos[i], image_box[i], box.length()[i]);
+  /** @brief Fold coordinates to primary simulation box in-place.
+   *  Lees-Edwards offset is ignored.
+   *  @param[in,out] pos        coordinate to fold
+   *  @param[in,out] image_box  image box offset
+   */
+  void fold_position(Utils::Vector3d &pos, Utils::Vector3i &image_box) const {
+    for (unsigned int i = 0u; i < 3u; i++) {
+      if (m_periodic[i]) {
+        auto const result =
+            Algorithm::periodic_fold(pos[i], image_box[i], m_length[i]);
+        if (result.second == std::numeric_limits<int>::min() or
+            result.second == std::numeric_limits<int>::max()) {
+          throw std::runtime_error(
+              "Overflow in the image box count while folding a particle "
+              "coordinate into the primary simulation box. Maybe a particle "
+              "experienced a huge force.");
+        }
+        std::tie(pos[i], image_box[i]) = result;
+      }
     }
   }
-}
 
-/** @brief Fold particle coordinates to primary simulation box.
- *  @param p    coordinate to fold
- *  @param box  box parameters (side lengths, periodicity)
- *  @return Folded coordinates.
- */
-inline Utils::Vector3d folded_position(const Utils::Vector3d &p,
-                                       const BoxGeometry &box) {
-  Utils::Vector3d p_folded;
-  for (unsigned int i = 0; i < 3; i++) {
-    if (box.periodic(i)) {
-      p_folded[i] = Algorithm::periodic_fold(p[i], box.length()[i]);
-    } else {
-      p_folded[i] = p[i];
+  /** @brief Calculate coordinates folded to primary simulation box.
+   *  @param p    coordinate to fold
+   *  @return Folded coordinates.
+   */
+  auto folded_position(Utils::Vector3d const &p) const {
+    Utils::Vector3d p_folded;
+    for (unsigned int i = 0u; i < 3u; i++) {
+      if (m_periodic[i]) {
+        p_folded[i] = Algorithm::periodic_fold(p[i], m_length[i]);
+      } else {
+        p_folded[i] = p[i];
+      }
     }
+
+    return p_folded;
   }
 
-  return p_folded;
-}
+  /** @brief Calculate image box shift vector */
+  auto image_shift(Utils::Vector3i const &image_box) const {
+    return detail::image_shift(image_box, m_length);
+  }
 
-#endif
+  /** @brief Unfold particle coordinates to image box. */
+  auto unfolded_position(Utils::Vector3d const &pos,
+                         Utils::Vector3i const &image_box) const {
+    return detail::unfolded_position(pos, image_box, m_length);
+  }
+};

src/core/CMakeLists.txt

@@ -33,7 +33,6 @@ add_library(
   forcecap.cpp
   forces.cpp
   ghosts.cpp
-  grid.cpp
   immersed_boundaries.cpp
   interactions.cpp
   event.cpp

src/core/EspressoSystemStandAlone.cpp

@@ -19,10 +19,13 @@
 
 #include "config/config.hpp"
 
+#include "BoxGeometry.hpp"
 #include "EspressoSystemStandAlone.hpp"
+#include "cell_system/CellStructure.hpp"
+#include "cell_system/CellStructureType.hpp"
+#include "cells.hpp"
 #include "communication.hpp"
 #include "event.hpp"
-#include "grid.hpp"
 #include "integrate.hpp"
 #include "system/System.hpp"
 #include "system/System.impl.hpp"
@@ -48,16 +51,20 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) {
 #ifdef VIRTUAL_SITES
   set_virtual_sites(std::make_shared<VirtualSitesOff>());
 #endif
-  ::System::set_system(std::make_shared<::System::System>());
+  auto system = std::make_shared<::System::System>();
+  ::System::set_system(system);
+  auto &cell_structure = *system->cell_structure;
+  cells_re_init(cell_structure, CellStructureType::CELL_STRUCTURE_REGULAR);
 }
 
 void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const {
-  set_box_length(box_l);
+  System::get_system().box_geo->set_length(box_l);
+  on_boxl_change();
 }
 
 void EspressoSystemStandAlone::set_node_grid(
     Utils::Vector3i const &node_grid) const {
-  ::set_node_grid(node_grid);
+  ::communicator.set_node_grid(node_grid);
 }
 
 void EspressoSystemStandAlone::set_time_step(double time_step) const {

src/core/LocalBox.hpp

@@ -16,37 +16,37 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef ESPRESSO_SRC_CORE_LOCALBOX_HPP
-#define ESPRESSO_SRC_CORE_LOCALBOX_HPP
+
+#pragma once
 
 #include "cell_system/CellStructureType.hpp"
 
 #include <utils/Array.hpp>
 #include <utils/Vector.hpp>
 
-template <class T> class LocalBox {
-  Utils::Vector<T, 3> m_local_box_l = {1, 1, 1};
-  Utils::Vector<T, 3> m_lower_corner = {0, 0, 0};
-  Utils::Vector<T, 3> m_upper_corner = {1, 1, 1};
+class LocalBox {
+  Utils::Vector3d m_local_box_l = {1., 1., 1.};
+  Utils::Vector3d m_lower_corner = {0., 0., 0.};
+  Utils::Vector3d m_upper_corner = {1., 1., 1.};
   Utils::Array<int, 6> m_boundaries = {};
   CellStructureType m_cell_structure_type;
 
 public:
   LocalBox() = default;
-  LocalBox(Utils::Vector<T, 3> const &lower_corner,
-           Utils::Vector<T, 3> const &local_box_length,
+  LocalBox(Utils::Vector3d const &lower_corner,
+           Utils::Vector3d const &local_box_length,
            Utils::Array<int, 6> const &boundaries,
            CellStructureType const cell_structure_type)
       : m_local_box_l(local_box_length), m_lower_corner(lower_corner),
         m_upper_corner(lower_corner + local_box_length),
         m_boundaries(boundaries), m_cell_structure_type(cell_structure_type) {}
 
   /** Left (bottom, front) corner of this nodes local box. */
-  Utils::Vector<T, 3> const &my_left() const { return m_lower_corner; }
+  auto const &my_left() const { return m_lower_corner; }
   /** Right (top, back) corner of this nodes local box. */
-  Utils::Vector<T, 3> const &my_right() const { return m_upper_corner; }
+  auto const &my_right() const { return m_upper_corner; }
   /** Dimensions of the box a single node is responsible for. */
-  Utils::Vector<T, 3> const &length() const { return m_local_box_l; }
+  auto const &length() const { return m_local_box_l; }
   /** @brief Boundary information for the local box.
    *
    * This returns for each of the faces of the local box if
@@ -56,17 +56,32 @@ template <class T> class LocalBox {
    *
    * @return Array with boundary information.
    */
-  Utils::Array<int, 6> const &boundary() const { return m_boundaries; }
+  auto const &boundary() const { return m_boundaries; }
 
   /** Return cell structure type. */
-  CellStructureType const &cell_structure_type() const {
-    return m_cell_structure_type;
-  }
+  auto const &cell_structure_type() const { return m_cell_structure_type; }
 
   /** Set cell structure type. */
   void set_cell_structure_type(CellStructureType cell_structure_type) {
     m_cell_structure_type = cell_structure_type;
   }
-};
 
-#endif
+  static LocalBox make_regular_decomposition(Utils::Vector3d const &box_l,
+                                             Utils::Vector3i const &node_index,
+                                             Utils::Vector3i const &node_grid) {
+
+    auto const local_length = Utils::hadamard_division(box_l, node_grid);
+    auto const my_left = Utils::hadamard_product(node_index, local_length);
+
+    decltype(LocalBox::m_boundaries) boundaries;
+    for (unsigned int dir = 0u; dir < 3u; dir++) {
+      /* left boundary ? */
+      boundaries[2u * dir] = (node_index[dir] == 0);
+      /* right boundary ? */
+      boundaries[2u * dir + 1u] = -(node_index[dir] + 1 == node_grid[dir]);
+    }
+
+    return {my_left, local_length, boundaries,
+            CellStructureType::CELL_STRUCTURE_REGULAR};
+  }
+};