HybridDecomposition: interface; coupling of child decompositions.

src/core/CellStructure.cpp

@@ -257,8 +257,8 @@ void CellStructure::set_domain_decomposition(
 
 void CellStructure::set_hybrid_decomposition(
     boost::mpi::communicator const &comm, double range, BoxGeometry const &box,
-    LocalBox<double> const &local_geo) {
-  set_particle_decomposition(
-      std::make_unique<HybridDecomposition>(comm, range, box, local_geo));
+    LocalBox<double> const &local_geo, std::set<int> n_square_types) {
+  set_particle_decomposition(std::make_unique<HybridDecomposition>(
+      comm, range, box, local_geo, n_square_types));
   m_type = CELL_STRUCTURE_HYBRID;
 }

src/core/CellStructure.hpp

@@ -506,7 +506,8 @@ struct CellStructure {
    */
   void set_hybrid_decomposition(boost::mpi::communicator const &comm,
                                 double range, BoxGeometry const &box,
-                                LocalBox<double> const &local_geo);
+                                LocalBox<double> const &local_geo,
+                                std::set<int> n_square_types);
 
 public:
   template <class BondKernel> void bond_loop(BondKernel const &bond_kernel) {

src/core/DomainDecomposition.cpp

@@ -571,7 +571,7 @@ DomainDecomposition::DomainDecomposition(boost::mpi::communicator comm,
                                          double range,
                                          const BoxGeometry &box_geo,
                                          const LocalBox<double> &local_geo)
-    : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo) {
+    : m_comm(comm), m_box(box_geo), m_local_box(local_geo) {
   /* set up new domain decomposition cell structure */
   create_cell_grid(range);
 

src/core/HybridDecomposition.cpp

@@ -27,11 +27,12 @@
 HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
                                          double range,
                                          const BoxGeometry &box_geo,
-                                         const LocalBox<double> &local_box)
-    : m_comm(std::move(comm)), m_box(box_geo),
-      m_domain_decomposition(
-          DomainDecomposition(m_comm, range, m_box, local_box)),
-      m_n_square(AtomDecomposition(m_comm, m_box)) {
+                                         const LocalBox<double> &local_box,
+                                         std::set<int> n_square_types)
+    : m_comm(comm), m_box(box_geo), m_domain_decomposition(DomainDecomposition(
+                                        m_comm, range, m_box, local_box)),
+      m_n_square(AtomDecomposition(m_comm, m_box)),
+      m_n_square_types(std::move(n_square_types)) {
 
   /* Vector containing cells of both child decompositions */
   m_local_cells = m_domain_decomposition.get_local_cells();
@@ -59,4 +60,18 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
   std::copy(collect_ghost_force_comm_n_square.communications.begin(),
             collect_ghost_force_comm_n_square.communications.end(),
             std::back_inserter(m_collect_ghost_force_comm.communications));
+
+  /* coupling between the child decompositions via neighborship relation */
+  std::vector<Cell *> additional_reds = m_n_square.get_local_cells();
+  std::copy(ghost_cells_n_square.begin(), ghost_cells_n_square.end(),
+            std::back_inserter(additional_reds));
+  for (auto &local_cell : m_domain_decomposition.local_cells()) {
+    std::vector<Cell *> red_neighbors(local_cell->m_neighbors.red().begin(),
+                                      local_cell->m_neighbors.red().end());
+    std::vector<Cell *> black_neighbors(local_cell->m_neighbors.black().begin(),
+                                        local_cell->m_neighbors.black().end());
+    std::copy(additional_reds.begin(), additional_reds.end(),
+              std::back_inserter(red_neighbors));
+    local_cell->m_neighbors = Neighbors<Cell *>(red_neighbors, black_neighbors);
+  }
 }

src/core/HybridDecomposition.hpp

@@ -23,6 +23,7 @@
 #define ESPRESSO_HYBRID_DECOMPOSITION_HPP
 
 #include "AtomDecomposition.hpp"
+#include "Cell.hpp"
 #include "DomainDecomposition.hpp"
 #include "Particle.hpp"
 #include "ParticleDecomposition.hpp"
@@ -52,12 +53,13 @@ class HybridDecomposition : public ParticleDecomposition {
   /** N-Square Decomposition to hold large particles */
   AtomDecomposition m_n_square;
   /** Set containing the types that should be handled using n_quare */
-  std::set<int> m_n_square_types = {1};
+  std::set<int> const m_n_square_types;
 
 public:
   HybridDecomposition(boost::mpi::communicator comm, double range,
                       const BoxGeometry &box_geo,
-                      const LocalBox<double> &local_box);
+                      const LocalBox<double> &local_box,
+                      std::set<int> n_square_types);
 
   Utils::Vector3i get_cell_grid() const {
     return m_domain_decomposition.get_cell_grid();
@@ -91,13 +93,10 @@ class HybridDecomposition : public ParticleDecomposition {
   }
 
   Cell *particle_to_cell(Particle const &p) override {
-    Cell *id_to_cell;
-    if (m_n_square_types.find(p.identity()) != m_n_square_types.end()) {
-      id_to_cell = m_domain_decomposition.particle_to_cell(p);
-    } else {
-      id_to_cell = m_n_square.particle_to_cell(p);
+    if (m_n_square_types.find(p.p.type) != m_n_square_types.end()) {
+      return m_domain_decomposition.particle_to_cell(p);
     }
-    return id_to_cell;
+    return m_n_square.particle_to_cell(p);
   }
 
   Utils::Vector3d max_cutoff() const override {

src/core/cells.cpp

@@ -43,6 +43,7 @@
 
 #include <algorithm>
 #include <boost/range/algorithm/min_element.hpp>
+#include <boost/serialization/set.hpp>
 #include <functional>
 #include <stdexcept>
 #include <utility>
@@ -206,30 +207,59 @@ std::vector<int> mpi_resort_particles(int global_flag) {
   return n_parts;
 }
 
-void cells_re_init(int new_cs) {
-  switch (new_cs) {
+void set_domain_decomposition() {
+  cell_structure.set_domain_decomposition(comm_cart, interaction_range(),
+                                          box_geo, local_geo);
+  on_cell_structure_change();
+}
+
+void set_atom_decomposition() {
+  cell_structure.set_atom_decomposition(comm_cart, box_geo);
+  on_cell_structure_change();
+}
+
+void set_hybrid_decomposition(std::set<int> n_square_types) {
+  cell_structure.set_hybrid_decomposition(comm_cart, interaction_range(),
+                                          box_geo, local_geo, n_square_types);
+  on_cell_structure_change();
+}
+
+REGISTER_CALLBACK(set_domain_decomposition)
+REGISTER_CALLBACK(set_atom_decomposition)
+REGISTER_CALLBACK(set_hybrid_decomposition)
+
+void mpi_set_domain_decomposition() { mpi_call_all(set_domain_decomposition); }
+void mpi_set_atom_decomposition() { mpi_call_all(set_atom_decomposition); }
+void mpi_set_hybrid_decomposition(std::set<int> n_square_types) {
+  mpi_call_all(set_hybrid_decomposition, n_square_types);
+}
+
+void cells_re_init() {
+  switch (cell_structure.decomposition_type()) {
   case CELL_STRUCTURE_DOMDEC:
     cell_structure.set_domain_decomposition(comm_cart, interaction_range(),
                                             box_geo, local_geo);
     break;
   case CELL_STRUCTURE_NSQUARE:
     cell_structure.set_atom_decomposition(comm_cart, box_geo);
     break;
-  case CELL_STRUCTURE_HYBRID:
-    cell_structure.set_hybrid_decomposition(comm_cart, interaction_range(),
-                                            box_geo, local_geo);
+  case CELL_STRUCTURE_HYBRID: {
+    /* Get current HybridDecomposition to extract n_square_types */
+    auto &current_hybrid_decomposition =
+        dynamic_cast<const HybridDecomposition &>(
+            Utils::as_const(cell_structure).decomposition());
+    cell_structure.set_hybrid_decomposition(
+        comm_cart, interaction_range(), box_geo, local_geo,
+        current_hybrid_decomposition.get_n_square_types());
     break;
+  }
   default:
     throw std::runtime_error("Unknown cell system type");
   }
 
   on_cell_structure_change();
 }
 
-REGISTER_CALLBACK(cells_re_init)
-
-void mpi_bcast_cell_structure(int cs) { mpi_call_all(cells_re_init, cs); }
-
 void check_resort_particles() {
   const double skin2 = Utils::sqr(skin / 2.0);
 

src/core/cells.hpp

@@ -36,6 +36,14 @@
  *    regardless their spatial position (see \ref AtomDecomposition.hpp).
  *    This is suitable for long range interactions that cannot be treated by a
  *    special method like P3M.
+ *  - hybrid decomposition: Initializes both domain decomposition and nsquare
+ *    at the same time and and has is given a set of particle types
+ *    n_square_types (see \ref HybridDecomposition.hpp). By default, particles
+ *    will be distributed using the domain decomposition. For particles of the
+ *    types defined as n_square_types the nsquare method is used. This is
+ *    suitable for systems containing lots of small particles with short
+ *    interaction range mixed with a few large particles with long interaction
+ *    range. There, the large particles should be treated using nsquare.
  */
 
 #include "Cell.hpp"
@@ -60,13 +68,30 @@ enum {
 /** Type of cell structure in use. */
 extern CellStructure cell_structure;
 
-/** Reinitialize the cell structures.
- *  @param new_cs The new topology to use afterwards.
+/** Initialize cell sturcture DomainDecomposition */
+void set_domain_decomposition();
+
+/** Initialize cell structure AtomDecomposition */
+void set_atom_decomposition();
+
+/** Initialize cell structure HybridDecomposition
+ *  @param n_square_types Set of particle types that will use n_square.
+ */
+void set_hybrid_decomposition(std::set<int> n_square_types);
+
+/** Change cell structure to DomainDecomposition on all nodes. */
+void mpi_set_domain_decomposition();
+
+/** Change cell structure to AtomDecomposition on all nodes. */
+void mpi_set_atom_decomposition();
+
+/** Change cell structure to HybridDecomposition on all nodes.
+ *  @param n_square_types Set of particle types that will use n_square.
  */
-void cells_re_init(int new_cs);
+void mpi_set_hybrid_decomposition(std::set<int> n_square_types);
 
-/** Change the cell structure on all nodes. */
-void mpi_bcast_cell_structure(int cs);
+/** Reinitialize the cell structures (using same topology). */
+void cells_re_init();
 
 /**
  * @brief Set @ref CellStructure::use_verlet_list

src/core/event.cpp

@@ -91,7 +91,7 @@ void on_program_start() {
   init_node_grid();
 
   /* initially go for domain decomposition */
-  cells_re_init(CELL_STRUCTURE_DOMDEC);
+  set_domain_decomposition();
 
   if (this_node == 0) {
     /* make sure interaction 0<->0 always exists */
@@ -222,8 +222,7 @@ void on_coulomb_change() {
 }
 
 void on_short_range_ia_change() {
-  cells_re_init(cell_structure.decomposition_type());
-
+  cells_re_init();
   recalc_forces = true;
 }
 
@@ -241,7 +240,7 @@ void on_boxl_change(bool skip_method_adaption) {
   grid_changed_box_l(box_geo);
   /* Electrostatics cutoffs mostly depend on the system size,
    * therefore recalculate them. */
-  cells_re_init(cell_structure.decomposition_type());
+  cells_re_init();
 
   if (not skip_method_adaption) {
     /* Now give methods a chance to react to the change in box length */
@@ -303,7 +302,7 @@ void on_periodicity_change() {
 }
 
 void on_skin_change() {
-  cells_re_init(cell_structure.decomposition_type());
+  cells_re_init();
   on_coulomb_change();
 }
 
@@ -320,7 +319,7 @@ void on_forcecap_change() { recalc_forces = true; }
 
 void on_nodegrid_change() {
   grid_changed_n_nodes();
-  cells_re_init(cell_structure.decomposition_type());
+  cells_re_init();
 }
 
 /**

src/python/espressomd/cellsystem.pxd

@@ -18,7 +18,7 @@
 #
 
 from libcpp cimport bool
-from libcpp.set cimport set
+from libcpp.set cimport set as cpp_set
 from libcpp.pair cimport pair, tuple
 from libcpp.vector cimport vector
 from .utils cimport Vector3i, Vector3d
@@ -54,7 +54,9 @@ cdef extern from "cells.hpp":
     vector[pair[int, int]] mpi_get_pairs_of_types(double distance, vector[int] types) except +
     vector[PairInfo] mpi_non_bonded_loop_trace() 
     vector[int] mpi_resort_particles(int global_flag)
-    void mpi_bcast_cell_structure(int cs)
+    void mpi_set_domain_decomposition()
+    void mpi_set_atom_decomposition()
+    void mpi_set_hybrid_decomposition(cpp_set[int] n_square_types)
     void mpi_set_use_verlet_lists(bool use_verlet_lists)
 
 cdef extern from "tuning.hpp":
@@ -74,7 +76,7 @@ cdef extern from "HybridDecomposition.hpp":
     cppclass HybridDecomposition:
         Vector3i get_cell_grid()
         Vector3d get_cell_size()
-        set[int] get_n_square_types()
+        cpp_set[int] get_n_square_types()
 
 cdef extern from "grid.hpp":
     void mpi_set_node_grid(const Vector3i & node_grid)

src/python/espressomd/cellsystem.pyx

@@ -42,7 +42,7 @@ cdef class CellSystem:
 
         """
         mpi_set_use_verlet_lists(use_verlet_lists)
-        mpi_bcast_cell_structure(CELL_STRUCTURE_DOMDEC)
+        mpi_set_domain_decomposition()
 
         handle_errors("Error while initializing the cell system.")
         return True
@@ -59,12 +59,12 @@ cdef class CellSystem:
 
         """
         mpi_set_use_verlet_lists(use_verlet_lists)
-        mpi_bcast_cell_structure(CELL_STRUCTURE_NSQUARE)
+        mpi_set_atom_decomposition()
 
         return True
 
     def set_hybrid_decomposition(
-            self, n_square_types=[], use_verlet_lists=True):
+            self, n_square_types=None, use_verlet_lists=True):
         """
         Activates the hybrid domain decomposition.
 
@@ -78,7 +78,9 @@ cdef class CellSystem:
 
         """
         mpi_set_use_verlet_lists(use_verlet_lists)
-        mpi_bcast_cell_structure(CELL_STRUCTURE_HYBRID)
+        if n_square_types is None:
+            n_square_types = set()
+        mpi_set_hybrid_decomposition(n_square_types)
 
         return True