improve testing

src/core/cell_system/CellStructure.cpp

@@ -44,6 +44,7 @@
 #include <cstddef>
 #include <iterator>
 #include <memory>
+#include <optional>
 #include <set>
 #include <stdexcept>
 #include <string>
@@ -249,8 +250,7 @@ void CellStructure::set_atom_decomposition() {
 }
 
 void CellStructure::set_regular_decomposition(
-    double range,
-    boost::optional<std::pair<int, int>> fully_connected_boundary) {
+    double range, std::optional<std::pair<int, int>> fully_connected_boundary) {
   auto &system = get_system();
   auto &local_geo = *system.local_geo;
   auto const &box_geo = *system.box_geo;

src/core/cell_system/CellStructure.hpp

@@ -46,6 +46,7 @@
 #include <cassert>
 #include <iterator>
 #include <memory>
+#include <optional>
 #include <set>
 #include <span>
 #include <stdexcept>
@@ -559,7 +560,7 @@ struct CellStructure : public System::Leaf<CellStructure> {
    */
   void set_regular_decomposition(
       double range,
-      boost::optional<std::pair<int, int>> fully_connected_boundary);
+      std::optional<std::pair<int, int>> fully_connected_boundary);
 
   /**
    * @brief Set the particle decomposition to @ref HybridDecomposition.

src/core/cell_system/HybridDecomposition.cpp

@@ -37,6 +37,7 @@
 #include <cstddef>
 #include <functional>
 #include <iterator>
+#include <optional>
 #include <set>
 #include <utility>
 
@@ -48,7 +49,7 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
                                          std::set<int> n_square_types)
     : m_comm(std::move(comm)), m_box(box_geo), m_cutoff_regular(cutoff_regular),
       m_regular_decomposition(RegularDecomposition(
-          m_comm, cutoff_regular + skin, m_box, local_box, {})),
+          m_comm, cutoff_regular + skin, m_box, local_box, std::nullopt)),
       m_n_square(AtomDecomposition(m_comm, m_box)),
       m_n_square_types(std::move(n_square_types)),
       m_get_global_ghost_flags(std::move(get_ghost_flags)) {

src/core/cell_system/RegularDecomposition.cpp

@@ -397,30 +397,29 @@ void RegularDecomposition::init_cell_interactions() {
   auto const &node_grid = ::communicator.node_grid;
   auto const global_halo_offset = hadamard_product(node_pos, cell_grid) - halo;
   auto const global_size = hadamard_product(node_grid, cell_grid);
-  auto at_boundary = [&global_size](int coord, Utils::Vector3i cell_idx) {
+  auto const at_boundary = [&global_size](int coord, Utils::Vector3i cell_idx) {
     return (cell_idx[coord] == 0 or cell_idx[coord] == global_size[coord]);
   };
 
   // For the fully connected feature (cells that don't share at least a corner)
   // only apply if one cell is a ghost cell (i.e. connections across the
-  // periodic boudnary.
-  auto fcb_is_inner_connection = [&global_size, this](Utils::Vector3i a,
-                                                      Utils::Vector3i b) {
-    if (not fully_connected_boundary())
-      return false;
-
-    const int fc_normal = (*fully_connected_boundary()).first;
-    const int fc_dir = (*fully_connected_boundary()).second;
-
-    bool involves_ghost_cell =
-        (a[fc_normal] == -1 or a[fc_normal] == global_size[fc_normal] or
-         b[fc_normal] == -1 or b[fc_normal] == global_size[fc_normal]);
-    if (involves_ghost_cell)
-      return false;
-    return (abs((a - b)[fc_dir]) > 1); // cells do not share at least a corner
+  // periodic boundary.
+  auto const fcb_is_inner_connection = [&global_size, this](Utils::Vector3i a,
+                                                            Utils::Vector3i b) {
+    if (fully_connected_boundary()) {
+      auto const [fc_normal, fc_dir] = *fully_connected_boundary();
+      auto const involves_ghost_cell =
+          (a[fc_normal] == -1 or a[fc_normal] == global_size[fc_normal] or
+           b[fc_normal] == -1 or b[fc_normal] == global_size[fc_normal]);
+      if (not involves_ghost_cell) {
+        // check if cells do not share at least a corner
+        return std::abs((a - b)[fc_dir]) > 1;
+      }
+    }
+    return false;
   };
 
-  /* Tanslate a node local index (relative to the origin of the local grid)
+  /* Translate a node local index (relative to the origin of the local grid)
    * to a global index. */
   auto global_index =
       [&](Utils::Vector3i const &local_index) -> Utils::Vector3i {
@@ -439,18 +438,19 @@ void RegularDecomposition::init_cell_interactions() {
       [&](Utils::Vector3i const &global_index) -> Utils::Vector3i {
     return (global_index - global_halo_offset);
   };
-  // Validation
+
+  // sanity checks
   if (fully_connected_boundary()) {
-    if ((*fully_connected_boundary()).first ==
-        (*fully_connected_boundary()).second) {
-      throw std::runtime_error("fully_connectd_boundary normal and connection "
-                               "coordinates need to differ.");
+    auto const [fc_normal, fc_dir] = *fully_connected_boundary();
+    if (fc_normal == fc_dir) {
+      throw std::domain_error("fully_connected_boundary normal and connection "
+                              "coordinates need to differ.");
     }
-    if (node_grid[(*fully_connected_boundary()).second] != 1) {
+    if (node_grid[fc_dir] != 1) {
       throw std::runtime_error(
           "The MPI nodegrid must be 1 in the fully connected direction.");
-    };
-  };
+    }
+  }
 
   /* We only consider local cells (e.g. not halo cells), which
    * span the range [(1,1,1), cell_grid) in local coordinates. */
@@ -469,8 +469,7 @@ void RegularDecomposition::init_cell_interactions() {
         * in the direction as neighbors, not only the nearest ones.
         //         */
         if (fully_connected_boundary()) {
-          int fc_boundary = (*fully_connected_boundary()).first;
-          int fc_direction = (*fully_connected_boundary()).second;
+          auto const [fc_boundary, fc_direction] = *fully_connected_boundary();
 
           // Fully connected is only needed at the box surface
           if (at_boundary(fc_boundary, {m, n, o})) {
@@ -671,7 +670,7 @@ GhostCommunicator RegularDecomposition::prepare_comm() {
 RegularDecomposition::RegularDecomposition(
     boost::mpi::communicator comm, double range, BoxGeometry const &box_geo,
     LocalBox const &local_geo,
-    boost::optional<std::pair<int, int>> fully_connected)
+    std::optional<std::pair<int, int>> fully_connected)
     : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo),
       m_fully_connected_boundary(std::move(fully_connected)) {
 

src/core/cell_system/RegularDecomposition.hpp

@@ -79,7 +79,7 @@ struct RegularDecomposition : public ParticleDecomposition {
   boost::mpi::communicator m_comm;
   BoxGeometry const &m_box;
   LocalBox m_local_box;
-  boost::optional<std::pair<int, int>> m_fully_connected_boundary = {};
+  std::optional<std::pair<int, int>> m_fully_connected_boundary = {};
   std::vector<Cell> cells;
   std::vector<Cell *> m_local_cells;
   std::vector<Cell *> m_ghost_cells;
@@ -89,11 +89,8 @@ struct RegularDecomposition : public ParticleDecomposition {
 public:
   RegularDecomposition(boost::mpi::communicator comm, double range,
                        BoxGeometry const &box_geo, LocalBox const &local_geo,
-                       boost::optional<std::pair<int, int>> fully_connected);
+                       std::optional<std::pair<int, int>> fully_connected);
 
-  boost::optional<std::pair<int, int>> fully_connected_boundary() const {
-    return m_fully_connected_boundary;
-  };
   GhostCommunicator const &exchange_ghosts_comm() const override {
     return m_exchange_ghosts_comm;
   }
@@ -120,6 +117,8 @@ struct RegularDecomposition : public ParticleDecomposition {
   Utils::Vector3d max_cutoff() const override;
   Utils::Vector3d max_range() const override;
 
+  auto fully_connected_boundary() const { return m_fully_connected_boundary; }
+
   std::optional<BoxGeometry> minimum_image_distance() const override {
     return {m_box};
   }

src/core/forces.cpp

@@ -221,7 +221,6 @@ void System::System::calculate_forces() {
     lb_couple_particles();
   }
 
-  particles = cell_structure->local_particles();
 #ifdef CUDA
 #ifdef CALIPER
   CALI_MARK_BEGIN("copy_forces_from_GPU");

src/python/espressomd/cell_system.py

@@ -115,14 +115,14 @@ def set_regular_decomposition(self, **kwargs):
         use_verlet_lists : :obj:`bool`, optional
             Activates or deactivates the usage of Verlet lists.
             Defaults to ``True``.
-        fully_connected_boundary : dict, optional
-            If set, connects all cells on a given boundary along the given direciton.
-            Example: {"boundary":"z", "direction":"x"} connects all cells on the boundary normal to 
-            the z-direcion along the x axis. This corresponds to
-            z as shear plane normal and x as shear direction in Lees-Edwards
-            boundary conditions.
-        """
+        fully_connected_boundary : :obj:`dict`, optional
+            If set, connects all cells on a given boundary along the given direction.
+            Example: ``{"boundary": "z", "direction": "x"}`` connects all
+            cells on the boundary normal to the z-direction along the x-axis.
+            This corresponds to z-axis as shear plane normal and x-axis as
+            shear direction in Lees-Edwards boundary conditions.
 
+        """
         self.call_method("initialize", name="regular_decomposition", **kwargs)
 
     def set_n_square(self, **kwargs):

src/script_interface/cell_system/CellSystem.cpp

@@ -40,7 +40,9 @@
 #include <boost/variant.hpp>
 
 #include <algorithm>
+#include <cassert>
 #include <iterator>
+#include <optional>
 #include <set>
 #include <sstream>
 #include <stdexcept>
@@ -49,29 +51,25 @@
 #include <utility>
 #include <vector>
 
-namespace {
-
-int coord(const std::string &s) {
+static int coord(std::string const &s) {
   if (s == "x")
     return 0;
   if (s == "y")
     return 1;
   if (s == "z")
     return 2;
-  throw std::runtime_error("Invalid argument for Cartesian coordinate: " + s);
+  throw std::invalid_argument("Invalid Cartesian coordinate: '" + s + "'");
 }
 
-std::string coord_letter(int c) {
+static std::string coord_letter(int c) {
   if (c == 0)
     return "x";
   if (c == 1)
     return "y";
-  if (c == 2)
-    return "z";
-  throw std::runtime_error("Invalid numeric coordinate");
+  assert(c == 2);
+  return "z";
 }
 
-} // namespace
 namespace ScriptInterface {
 namespace CellSystem {
 
@@ -299,19 +297,14 @@ void CellSystem::initialize(CellStructureType const &cs_type,
     auto n_square_types = std::set<int>{ns_types.begin(), ns_types.end()};
     m_cell_structure->set_hybrid_decomposition(cutoff_regular, n_square_types);
   } else if (cs_type == CellStructureType::REGULAR) {
-    boost::optional<std::pair<int, int>> fcb_pair = {};
-    if (params.count("fully_connected_boundary") == 0 ||
-        is_type<None>(params.at("fully_connected_boundary"))) {
-      // pass
-    } else {
-      auto const fully_connected_boundary =
-          get_value<std::unordered_map<std::string, Variant>>(
-              params, "fully_connected_boundary");
-      const int boundary = coord(
-          boost::get<std::string>(fully_connected_boundary.at("boundary")));
-      const int direction = coord(
-          boost::get<std::string>(fully_connected_boundary.at("direction")));
-      fcb_pair = {{boundary, direction}};
+    std::optional<std::pair<int, int>> fcb_pair = std::nullopt;
+    if (params.contains("fully_connected_boundary") and
+        not is_none(params.at("fully_connected_boundary"))) {
+      auto const variant = get_value<VariantMap>(params, "fully_connected_boundary");
+      context()->parallel_try_catch([&fcb_pair, &variant]() {
+        fcb_pair = {{coord(boost::get<std::string>(variant.at("boundary"))),
+                     coord(boost::get<std::string>(variant.at("direction")))}};
+      });
     }
     context()->parallel_try_catch([this, &fcb_pair]() {
       m_cell_structure->set_regular_decomposition(