RegularDecomposition.hpp

Go to the documentation of this file.

/*
 * Copyright (C) 2010-2026 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/>.
 */
 
#pragma once
 
#include "cell_system/ParticleDecomposition.hpp"
 
#include "cell_system/Cell.hpp"
 
#include "BoxGeometry.hpp"
#include "LocalBox.hpp"
#include "Particle.hpp"
#include "ParticleList.hpp"
#include "ghosts/HaloPlan.hpp"
 
#include <utils/Vector.hpp>
 
#include <boost/mpi/communicator.hpp>
 
#include <cstddef>
#include <optional>
#include <span>
#include <utility>
#include <vector>
 
/**
 * @brief Regular decomposition cell system.
 *
 * The domain of a node is split into a 3D cell grid with dimension
 * @ref RegularDecomposition::cell_grid "cell_grid". Together with one ghost
 * cell layer on each side the overall dimension of the ghost cell grid is
 * @ref RegularDecomposition::ghost_cell_grid "ghost_cell_grid". The regular
 * decomposition enables the use of the linked cell algorithm
 * which is in turn used for setting up the Verlet list for the
 * system. You can see a 2D graphical representation of the linked
 * cell grid below.
 *
 * \image html  linked_cells.gif "Linked cells structure"
 *
 * 2D representation of a linked cell grid:
 * <tt>cell_grid = {4,4}, ghost_cell_grid = {6,6}</tt>
 *
 * Each cell has @f$ 3^D @f$ neighbor cells. Since we deal with pair forces,
 * it is sufficient to calculate only half of the interactions (Newton's law:
 * action = reaction). We have chosen the upper half e.g. all neighbor
 * cells with a higher linear index (for cell 14 they are marked in light
 * blue). Caution: This implementation needs double sided ghost
 * communication! For single sided ghost communication one would need
 * some ghost-ghost cell interaction as well, which we do not need!
 */
struct RegularDecomposition : public ParticleDecomposition {
  /** Grid dimensions per node. */
  Utils::Vector3i cell_grid = {};
  /** Cell size. */
  Utils::Vector3d cell_size = {};
  /** Offset in global grid */
  Utils::Vector3i cell_offset = {};
  /** linked cell grid with ghost frame. */
  Utils::Vector3i ghost_cell_grid = {};
  /** inverse @ref RegularDecomposition::cell_size "cell_size". */
  Utils::Vector3d inv_cell_size = {};
 
  boost::mpi::communicator m_comm;
  BoxGeometry const &m_box;
  LocalBox m_local_box;
  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;
  /**
   * Topology-agnostic direct-neighbor halo plan (see @c make_halo_plan).
   * Holds ParticleList pointers into this decomposition's cells.
   * Value-copying this object leaves these pointers dangling.
   * @todo make non-copyable or rebuild-on-copy.
   */
  GhostComm::HaloPlan m_halo_plan;
 
public:
  RegularDecomposition(boost::mpi::communicator comm, double range,
                       BoxGeometry const &box_geo, LocalBox const &local_geo,
                       std::optional<std::pair<int, int>> fully_connected);
 
  GhostComm::HaloPlan const *halo_plan() const override { return &m_halo_plan; }
 
  std::span<Cell *const> local_cells() const override { return m_local_cells; }
  std::span<Cell *const> ghost_cells() const override { return m_ghost_cells; }
 
  /* Getter needed for HybridDecomposition */
  auto const &get_local_cells() const { return m_local_cells; }
  auto const &get_ghost_cells() const { return m_ghost_cells; }
 
  Cell *particle_to_cell(Particle const &p) override {
    return position_to_cell(p.pos());
  }
 
  Cell const *particle_to_cell(Particle const &p) const override {
    return position_to_cell(p.pos());
  }
 
  void resort(bool global, std::vector<ParticleChange> &diff) override;
  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};
  }
 
  BoxGeometry const &box() const override { return m_box; }
 
private:
  /** Fill @c m_local_cells list and @c m_ghost_cells list for use with regular
   *  decomposition.
   */
  void mark_cells();
 
  int calc_processor_min_num_cells() const;
 
  int position_to_cell_index(Utils::Vector3d const &pos) const;
 
  /**
   * @brief Get pointer to the cell which corresponds to the position if the
   * position is in the node's spatial domain, otherwise a nullptr.
   */
  Cell *position_to_cell(Utils::Vector3d const &pos) {
    auto const index = position_to_cell_index(pos);
    return (index < 0) ? nullptr : &(cells.at(static_cast<std::size_t>(index)));
  }
  Cell const *position_to_cell(Utils::Vector3d const &pos) const {
    auto const index = position_to_cell_index(pos);
    return (index < 0) ? nullptr : &(cells.at(static_cast<std::size_t>(index)));
  }
 
  /**
   * @brief Move particles into the cell system if it belongs to this node.
   *
   * Moves all particles from src into the local cell
   * system if they do belong here. Otherwise the
   * particles are moved into rest.
   *
   * @param src Particles to move.
   * @param rest Output list for left-over particles.
   * @param modified_cells Local cells that were touched.
   */
  void move_if_local(ParticleList &src, ParticleList &rest,
                     std::vector<ParticleChange> &modified_cells);
 
  /**
   * @brief Split particle list by direction.
   *
   * Moves all particles from @p src into @p left
   * or @p right depending on whether they belong
   * to the left or right side of the local node
   * in direction @p dir.
   *
   * @param src Particles to sort.
   * @param left Particles that should go to the left
   * @param right Particles that should go to the right
   * @param dir Direction to consider.
   */
  void move_left_or_right(ParticleList &src, ParticleList &left,
                          ParticleList &right, int dir) const;
 
  /**
   * @brief One round of particle exchange with the next neighbors.
   *
   * @param[in] pl Particle on the move
   * @param[out] modified_cells Cells that got touched.
   */
  void exchange_neighbors(ParticleList &pl,
                          std::vector<ParticleChange> &modified_cells);
 
  /**
   *  @brief Calculate cell grid dimensions, cell sizes and number of cells.
   *
   *  Calculates the cell grid, based on the local box size and the range.
   *  If the number of cells is larger than @c max_num_cells,
   *  it increases @c max_range until the number of cells is
   *  smaller or equal to @c max_num_cells. It sets:
   *  @c cell_grid,
   *  @c ghost_cell_grid,
   *  @c cell_size, and
   *  @c inv_cell_size.
   *
   *  @param range interaction range. All pairs closer
   *               than this distance are found.
   */
  void create_cell_grid(double range);
 
  /** Init cell interactions for cell system regular decomposition.
   *  Initializes the interacting neighbor cell list of a cell.
   *  This list of interacting neighbor cells is used by the Verlet
   *  algorithm.
   */
  void init_cell_interactions();
 
  /** @brief Build a direct-neighbor @ref GhostComm::HaloPlan.
   *
   *  Every one of the up-to-26 stencil neighbors is addressed as an explicit
   *  peer (no axis-by-axis relay). The result is cached in @ref m_halo_plan.
   */
  GhostComm::HaloPlan make_halo_plan();
 
  /** Maximal number of cells per node. In order to avoid memory
   *  problems due to the cell grid, one has to specify the maximal
   *  number of cells. If the number of cells is larger
   *  than @c max_num_cells, the cell grid is reduced.
   *  @c max_num_cells has to be larger than 27, e.g. one inner cell.
   */
  static constexpr int max_num_cells = 32768;
};