Implementation of HybridDecomposition (#4373)

Fixes #4351 

Description of changes:
- Implement HybridDecomposition
- Add tests and doc
- Make DomainDecomposition `class` (instead of `struct`)

doc/sphinx/system_setup.rst

@@ -337,3 +337,53 @@ this node is twice as high. For 3 processors, the interactions are 0-0,
 
 Therefore it is highly recommended that you use N-squared only with an
 odd number of nodes, if with multiple processors at all.
+
+.. _Hybrid:
+
+Hybrid decomposition
+^^^^^^^^^^^^^^^^^^^^
+
+If for a simulation setup the interaction range is much smaller than the
+system size, use of a :ref:`Regular decomposition` leads to efficient
+scaling behavior (order :math:`N` instead of order :math:`N^2`).
+Consider a system with many small particles, e.g. a polymer solution.
+There, already the addition of one single large particle increases the maximum
+interaction range and thus the minimum cell size of the decomposition.
+Due to this larger cell size, throughout the simulation box a large number
+of non-interacting pairs of small particles is visited during the short
+range calculation. This can considerably increase the computational cost of
+the simulation.
+
+For such simulation setups, i.e. systems with a few large particles and much
+more small particles, the hybrid decomposition can be used. This hybrid
+decomposition is backed by two coupled particle decompositions which can
+be used to efficiently deal with the differently sized particles.
+Specifically that means putting the small particles into a
+:ref:`Regular decomposition`. There, the minimum cell size is limited only
+by the maximum interaction range of all particles within this decomposition.
+The few large particles are put into a :ref:`N-squared` cellsystem. Particles
+within this decomposition interact both, amongst each other and with all small
+particles in the :ref:`Regular decomposition`. The hybrid decomposition can therefore
+effectively recover the computational efficiency of the regular decomposition,
+given that only a few large particles have been added.
+
+Invoking :py:meth:`~espressomd.cellsystem.CellSystem.set_hybrid_decomposition`
+selects the hybrid decomposition. ::
+
+    system = espressomd.System(box_l=[10, 10, 10])
+    system.cell_system.set_hybrid_decomposition(n_square_types={1, 3}, cutoff_regular=1.2)
+
+Here, ``n_square_types`` is a python set containing the types of particles to
+put into the :ref:`N-squared` cellsystem, i.e. the particle types of the
+large particles. Particles with other types will by default be put into the
+:ref:`Regular decomposition`. Note that for now it is also necessary to manually set
+the maximum cutoff to consider for interactions within the
+:ref:`Regular decomposition`, i.e. the maximum interaction range among all
+small particle types. Set this via the ``cutoff_regular`` parameter.
+
+.. note::
+
+  The hybrid particle decomposition has been added to |es| only recently and
+  for now should be considered an experimental feature. If you notice some unexpected
+  behavior please let us know via github or the mailing list.
+

src/core/AtomDecomposition.hpp

@@ -80,6 +80,10 @@ class AtomDecomposition : public ParticleDecomposition {
     return Utils::make_span(m_ghost_cells);
   }
 
+  /* Getter needed for HybridDecomposition */
+  std::vector<Cell *> get_local_cells() const { return m_local_cells; }
+  std::vector<Cell *> get_ghost_cells() const { return m_ghost_cells; }
+
   /**
    * @brief Determine which cell a particle id belongs to.
    *

src/core/CMakeLists.txt

@@ -15,6 +15,7 @@ set(EspressoCore_SRC
     galilei.cpp
     ghosts.cpp
     grid.cpp
+    HybridDecomposition.cpp
     immersed_boundaries.cpp
     interactions.cpp
     event.cpp

src/core/CellStructure.cpp

@@ -23,6 +23,7 @@
 
 #include "AtomDecomposition.hpp"
 #include "CellStructureType.hpp"
+#include "HybridDecomposition.hpp"
 #include "RegularDecomposition.hpp"
 #include "grid.hpp"
 #include "lees_edwards/lees_edwards.hpp"
@@ -154,7 +155,7 @@ void CellStructure::remove_all_particles() {
 
 /* Map the data parts flags from cells to those used internally
  * by the ghost communication */
-static unsigned map_data_parts(unsigned data_parts) {
+unsigned map_data_parts(unsigned data_parts) {
   using namespace Cells;
 
   /* clang-format off */
@@ -260,3 +261,13 @@ void CellStructure::set_regular_decomposition(
   m_type = CellStructureType::CELL_STRUCTURE_REGULAR;
   local_geo.set_cell_structure_type(m_type);
 }
+
+void CellStructure::set_hybrid_decomposition(
+    boost::mpi::communicator const &comm, double cutoff_regular,
+    BoxGeometry const &box, LocalBox<double> &local_geo,
+    std::set<int> n_square_types) {
+  set_particle_decomposition(std::make_unique<HybridDecomposition>(
+      comm, cutoff_regular, box, local_geo, n_square_types));
+  m_type = CellStructureType::CELL_STRUCTURE_HYBRID;
+  local_geo.set_cell_structure_type(m_type);
+}

src/core/CellStructure.hpp

@@ -26,6 +26,7 @@
 #include "BoxGeometry.hpp"
 #include "Cell.hpp"
 #include "CellStructureType.hpp"
+#include "HybridDecomposition.hpp"
 #include "LocalBox.hpp"
 #include "Particle.hpp"
 #include "ParticleDecomposition.hpp"
@@ -74,6 +75,15 @@ enum DataPart : unsigned {
 };
 } // namespace Cells
 
+/**
+ * @brief Map the data parts flags from cells to those
+ *        used internally by the ghost communication.
+ *
+ * @param data_parts data parts flags
+ * @return ghost communication flags
+ */
+unsigned map_data_parts(unsigned data_parts);
+
 namespace Cells {
 inline ParticleRange particles(Utils::Span<Cell *> cells) {
   /* Find first non-empty cell */
@@ -522,6 +532,20 @@ struct CellStructure {
                                  double range, BoxGeometry const &box,
                                  LocalBox<double> &local_geo);
 
+  /**
+   * @brief Set the particle decomposition to HybridDecomposition.
+   *
+   * @param comm Communicator to use.
+   * @param cutoff_regular Interaction cutoff_regular.
+   * @param box Box geometry.
+   * @param local_geo Geometry of the local box.
+   * @param n_square_types Particle types to put into n_square decomposition.
+   */
+  void set_hybrid_decomposition(boost::mpi::communicator const &comm,
+                                double cutoff_regular, BoxGeometry const &box,
+                                LocalBox<double> &local_geo,
+                                std::set<int> n_square_types);
+
 public:
   template <class BondKernel> void bond_loop(BondKernel const &bond_kernel) {
     for (auto &p : local_particles()) {

src/core/CellStructureType.hpp

@@ -27,7 +27,9 @@ enum class CellStructureType : int {
   /** cell structure regular decomposition */
   CELL_STRUCTURE_REGULAR = 1,
   /** cell structure n square */
-  CELL_STRUCTURE_NSQUARE = 2
+  CELL_STRUCTURE_NSQUARE = 2,
+  /** cell structure hybrid */
+  CELL_STRUCTURE_HYBRID = 3
 };
 
 #endif // ESPRESSO_CELLSTRUCTURETYPE_HPP

src/core/HybridDecomposition.cpp

@@ -0,0 +1,164 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "HybridDecomposition.hpp"
+
+#include "CellStructure.hpp"
+#include "event.hpp"
+
+#include <utility>
+#include <utils/Vector.hpp>
+#include <utils/mpi/sendrecv.hpp>
+
+HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
+                                         double cutoff_regular,
+                                         const BoxGeometry &box_geo,
+                                         const LocalBox<double> &local_box,
+                                         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_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_regular_decomposition.get_local_cells();
+  auto local_cells_n_square = m_n_square.get_local_cells();
+  std::copy(local_cells_n_square.begin(), local_cells_n_square.end(),
+            std::back_inserter(m_local_cells));
+
+  /* Vector containing ghost cells of both child decompositions */
+  m_ghost_cells = m_regular_decomposition.get_ghost_cells();
+  auto ghost_cells_n_square = m_n_square.get_ghost_cells();
+  std::copy(ghost_cells_n_square.begin(), ghost_cells_n_square.end(),
+            std::back_inserter(m_ghost_cells));
+
+  /* Communicators that contain communications of both child decompositions */
+  m_exchange_ghosts_comm = m_regular_decomposition.exchange_ghosts_comm();
+  auto exchange_ghosts_comm_n_square = m_n_square.exchange_ghosts_comm();
+  std::copy(exchange_ghosts_comm_n_square.communications.begin(),
+            exchange_ghosts_comm_n_square.communications.end(),
+            std::back_inserter(m_exchange_ghosts_comm.communications));
+
+  m_collect_ghost_force_comm =
+      m_regular_decomposition.collect_ghost_force_comm();
+  auto collect_ghost_force_comm_n_square =
+      m_n_square.collect_ghost_force_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_regular_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);
+  }
+}
+
+void HybridDecomposition::resort(bool global,
+                                 std::vector<ParticleChange> &diff) {
+  ParticleList displaced_parts;
+
+  /* Check for n_square type particles in regular decomposition */
+  for (auto &c : m_regular_decomposition.local_cells()) {
+    for (auto it = c->particles().begin(); it != c->particles().end();) {
+      /* Particle is in the right decomposition, i.e. has no n_square type */
+      if (not is_n_square_type(it->p.type)) {
+        std::advance(it, 1);
+        continue;
+      }
+
+      /* else remove from current cell ... */
+      auto p = std::move(*it);
+      it = c->particles().erase(it);
+      diff.emplace_back(ModifiedList{c->particles()});
+      diff.emplace_back(RemovedParticle{p.identity()});
+
+      /* ... and insert into a n_square cell */
+      auto const first_local_cell = m_n_square.get_local_cells()[0];
+      first_local_cell->particles().insert(std::move(p));
+      diff.emplace_back(ModifiedList{first_local_cell->particles()});
+    }
+
+    /* Now check for regular decomposition type particles in n_square */
+    for (auto &c : m_n_square.local_cells()) {
+      for (auto it = c->particles().begin(); it != c->particles().end();) {
+        /* Particle is of n_square type */
+        if (is_n_square_type(it->p.type)) {
+          std::advance(it, 1);
+          continue;
+        }
+
+        /* else remove from current cell ... */
+        auto p = std::move(*it);
+        it = c->particles().erase(it);
+        diff.emplace_back(ModifiedList{c->particles()});
+        diff.emplace_back(RemovedParticle{p.identity()});
+
+        /* ... and insert in regular decomposition */
+        auto const target_cell = particle_to_cell(p);
+        /* if particle belongs to this node insert it into correct cell */
+        if (target_cell != nullptr) {
+          target_cell->particles().insert(std::move(p));
+          diff.emplace_back(ModifiedList{target_cell->particles()});
+        }
+        /* otherwise just put into regular decomposition */
+        else {
+          auto first_local_cell = m_regular_decomposition.get_local_cells()[0];
+          first_local_cell->particles().insert(std::move(p));
+          diff.emplace_back(ModifiedList{first_local_cell->particles()});
+        }
+      }
+    }
+  }
+
+  /* now resort into correct cells within the respective decompositions */
+  m_regular_decomposition.resort(global, diff);
+  m_n_square.resort(global, diff);
+
+  /* basically do CellStructure::ghost_count() */
+  ghost_communicator(exchange_ghosts_comm(), GHOSTTRANS_PARTNUM);
+
+  /* basically do CellStructure::ghost_update(unsigned data_parts) */
+  ghost_communicator(exchange_ghosts_comm(),
+                     map_data_parts(global_ghost_flags()));
+}
+
+Utils::Vector<std::size_t, 2>
+HybridDecomposition::parts_per_decomposition_local() const {
+  std::size_t parts_in_regular = 0;
+  std::size_t parts_in_n_square = 0;
+  for (auto &c : m_regular_decomposition.get_local_cells()) {
+    parts_in_regular += c->particles().size();
+  }
+  for (auto &c : m_n_square.get_local_cells()) {
+    parts_in_n_square += c->particles().size();
+  }
+  return {parts_in_regular, parts_in_n_square};
+}