DirectSum.h

Go to the documentation of this file.

 
#pragma once
 
#include "autopas/cells/FullParticleCell.h"
#include "autopas/containers/CellBasedParticleContainer.h"
#include "autopas/containers/CellBorderAndFlagManager.h"
#include "autopas/containers/CompatibleTraversals.h"
#include "autopas/containers/LeavingParticleCollector.h"
#include "autopas/containers/cellTraversals/CellTraversal.h"
#include "autopas/containers/directSum/traversals/DSTraversalInterface.h"
#include "autopas/iterators/ContainerIterator.h"
#include "autopas/options/DataLayoutOption.h"
#include "autopas/particles/OwnershipState.h"
#include "autopas/utils/AutoPasMacros.h"
#include "autopas/utils/ExceptionHandler.h"
#include "autopas/utils/ParticleCellHelpers.h"
#include "autopas/utils/StringUtils.h"
#include "autopas/utils/inBox.h"
 
namespace autopas {
 
template <class Particle_T>
class DirectSum : public CellBasedParticleContainer<FullParticleCell<Particle_T>> {
 public:
  using ParticleCellType = FullParticleCell<Particle_T>;
 
  using ParticleType = Particle_T;
 
  DirectSum(const std::array<double, 3> &boxMin, const std::array<double, 3> &boxMax, double cutoff, double skin,
            const size_t aosSortingThreshold, const size_t soaSortingThreshold)
      : CellBasedParticleContainer<ParticleCellType>(boxMin, boxMax, cutoff, skin, aosSortingThreshold,
                                                     soaSortingThreshold),
        _cellBorderFlagManager() {
    using namespace autopas::utils::ArrayMath::literals;
    // 1 owned and 6 halo cells
    this->_cells.resize(7);
    this->_cells[0].setPossibleParticleOwnerships(OwnershipState::owned);
    std::for_each(++this->_cells.begin(), this->_cells.end(),
                  [&](auto &cell) { cell.setPossibleParticleOwnerships(OwnershipState::halo); });
    auto boxLength = boxMax - boxMin;
    this->_cells[0].setCellLength(boxLength);
  }
 
  [[nodiscard]] ContainerOption getContainerType() const override { return ContainerOption::directSum; }
 
  void reserve(size_t numParticles, size_t numParticlesHaloEstimate) override {
    this->getOwnedCell().reserve(numParticles);
    for (auto cellIt = ++this->_cells.begin(); cellIt != this->_cells.end(); cellIt++) {
      cellIt->reserve(numParticlesHaloEstimate);
    }
  };
 
  void addParticleImpl(const Particle_T &p) override { getOwnedCell().addParticle(p); }
 
  void addHaloParticleImpl(const Particle_T &haloParticle) override {
    const auto boxMax = this->getBoxMax();
    const auto boxMin = this->getBoxMin();
    const auto pos = haloParticle.getR();
 
    for (size_t dim = 0; dim < 3; ++dim) {
      if (pos[dim] < boxMin[dim]) {
        this->_cells[2 * dim + 1].addParticle(haloParticle);
        return;
      } else if (pos[dim] >= boxMax[dim]) {
        this->_cells[2 * dim + 2].addParticle(haloParticle);
        return;
      }
    }
  }
 
  bool updateHaloParticle(const Particle_T &haloParticle) override {
    const auto boxMax = this->getBoxMax();
    const auto boxMin = this->getBoxMin();
    const auto pos = haloParticle.getR();
    const auto skinHalf = 0.5 * this->getVerletSkin();
 
    // Look for the particle in halo cells that are within half the skin distance of its position
    for (size_t dim = 0; dim < 3; ++dim) {
      if (pos[dim] < boxMin[dim] + skinHalf) {
        if (internal::checkParticleInCellAndUpdateByIDAndPosition(this->_cells[2 * dim + 1], haloParticle, skinHalf)) {
          return true;
        }
      } else if (pos[dim] >= boxMax[dim] - skinHalf) {
        if (internal::checkParticleInCellAndUpdateByIDAndPosition(this->_cells[2 * dim + 2], haloParticle, skinHalf)) {
          return true;
        }
      }
    }
    return false;
  }
 
  void deleteHaloParticles() override {
    for (auto cellIt = ++this->_cells.begin(); cellIt != this->_cells.end(); cellIt++) {
      cellIt->clear();
    }
  }
 
  void rebuildNeighborLists(TraversalInterface *traversal) override {
    // nothing to do.
  }
 
  void computeInteractions(TraversalInterface *traversal) override {
    prepareTraversal(traversal);
 
    traversal->initTraversal();
    traversal->traverseParticles();
    traversal->endTraversal();
  }
 
  [[nodiscard]] std::vector<Particle_T> updateContainer(bool keepNeighborListsValid) override {
    if (keepNeighborListsValid) {
      return LeavingParticleCollector::collectParticlesAndMarkNonOwnedAsDummy(*this);
    }
    // first we delete halo particles, as we don't want them here.
    deleteHaloParticles();
    getOwnedCell().deleteDummyParticles();
 
    std::vector<Particle_T> invalidParticles{};
    auto &particleVec = getOwnedCell()._particles;
    for (auto iter = particleVec.begin(); iter != particleVec.end();) {
      if (utils::notInBox(iter->getR(), this->getBoxMin(), this->getBoxMax())) {
        invalidParticles.push_back(*iter);
        // swap-delete
        *iter = particleVec.back();
        particleVec.pop_back();
      } else {
        ++iter;
      }
    }
    return invalidParticles;
  }
 
  [[nodiscard]] TraversalSelectorInfo getTraversalSelectorInfo() const override {
    using namespace autopas::utils::ArrayMath::literals;
 
    // direct sum consists of seven cells (owned + two halo cells in each dimension)
    return TraversalSelectorInfo(
        // DS container can be viewed as a 3x3x3 grid with some halo cells being combined.
        {3, 3, 3},
        // intentionally use cutoff here, as ds_sequential should be using the cutoff.
        this->getCutoff(), this->getBoxMax() - this->getBoxMin(), 0);
  }
 
  [[nodiscard]] ContainerIterator<Particle_T, true, false> begin(
      IteratorBehavior behavior = IteratorBehavior::ownedOrHalo,
      utils::optRef<typename ContainerIterator<Particle_T, true, false>::ParticleVecType> additionalVectors =
          std::nullopt) override {
    return ContainerIterator<Particle_T, true, false>(*this, behavior, additionalVectors);
  }
 
  [[nodiscard]] ContainerIterator<Particle_T, false, false> begin(
      IteratorBehavior behavior = IteratorBehavior::ownedOrHalo,
      utils::optRef<typename ContainerIterator<Particle_T, false, false>::ParticleVecType> additionalVectors =
          std::nullopt) const override {
    return ContainerIterator<Particle_T, false, false>(*this, behavior, additionalVectors);
  }
 
  template <typename Lambda>
  void forEach(Lambda forEachLambda, IteratorBehavior behavior) {
    if (behavior & IteratorBehavior::owned) {
      getOwnedCell().forEach(forEachLambda);
    }
    if (behavior & IteratorBehavior::halo) {
      for (auto cellIt = ++this->_cells.begin(); cellIt != this->_cells.end(); cellIt++) {
        cellIt->forEach(forEachLambda);
      }
    }
    // sanity check
    if (not(behavior & IteratorBehavior::ownedOrHalo)) {
      utils::ExceptionHandler::exception("Encountered invalid iterator behavior!");
    }
  }
 
  template <typename Lambda, typename A>
  void reduce(Lambda reduceLambda, A &result, IteratorBehavior behavior) {
    if (behavior & IteratorBehavior::owned) {
      getOwnedCell().reduce(reduceLambda, result);
    }
    if (behavior & IteratorBehavior::halo) {
      for (auto cellIt = ++this->_cells.begin(); cellIt != this->_cells.end(); cellIt++) {
        cellIt->reduce(reduceLambda, result);
      }
    }
    // sanity check
    if (not(behavior & IteratorBehavior::ownedOrHalo)) {
      utils::ExceptionHandler::exception("Encountered invalid iterator behavior!");
    }
  }
 
  [[nodiscard]] ContainerIterator<Particle_T, true, true> getRegionIterator(
      const std::array<double, 3> &lowerCorner, const std::array<double, 3> &higherCorner, IteratorBehavior behavior,
      utils::optRef<typename ContainerIterator<Particle_T, true, true>::ParticleVecType> additionalVectors =
          std::nullopt) override {
    return ContainerIterator<Particle_T, true, true>(*this, behavior, additionalVectors, lowerCorner, higherCorner);
  }
 
  [[nodiscard]] ContainerIterator<Particle_T, false, true> getRegionIterator(
      const std::array<double, 3> &lowerCorner, const std::array<double, 3> &higherCorner, IteratorBehavior behavior,
      utils::optRef<typename ContainerIterator<Particle_T, false, true>::ParticleVecType> additionalVectors =
          std::nullopt) const override {
    return ContainerIterator<Particle_T, false, true>(*this, behavior, additionalVectors, lowerCorner, higherCorner);
  }
 
  template <typename Lambda>
  void forEachInRegion(Lambda forEachLambda, const std::array<double, 3> &lowerCorner,
                       const std::array<double, 3> &higherCorner, IteratorBehavior behavior) {
    if (behavior & IteratorBehavior::owned) {
      getOwnedCell().forEach(forEachLambda, lowerCorner, higherCorner, behavior);
    }
    if (behavior & IteratorBehavior::halo) {
      for (auto cellIt = ++this->_cells.begin(); cellIt != this->_cells.end(); cellIt++) {
        cellIt->forEach(forEachLambda, lowerCorner, higherCorner, behavior);
      }
    }
    // sanity check
    if (not(behavior & IteratorBehavior::ownedOrHalo)) {
      utils::ExceptionHandler::exception("Encountered invalid iterator behavior!");
    }
  }
 
  template <typename Lambda, typename A>
  void reduceInRegion(Lambda reduceLambda, A &result, const std::array<double, 3> &lowerCorner,
                      const std::array<double, 3> &higherCorner, IteratorBehavior behavior) {
    if (behavior & IteratorBehavior::owned) {
      getOwnedCell().reduce(reduceLambda, result, lowerCorner, higherCorner, behavior);
    }
    if (behavior & IteratorBehavior::halo) {
      for (auto cellIt = ++this->_cells.begin(); cellIt != this->_cells.end(); cellIt++) {
        cellIt->reduce(reduceLambda, result, lowerCorner, higherCorner, behavior);
      }
    }
    // sanity check
    if (not(behavior & IteratorBehavior::ownedOrHalo)) {
      utils::ExceptionHandler::exception("Encountered invalid iterator behavior!");
    }
  }
 
  std::tuple<const Particle_T *, size_t, size_t> getParticle(size_t cellIndex, size_t particleIndex,
                                                             IteratorBehavior iteratorBehavior,
                                                             const std::array<double, 3> &boxMin,
                                                             const std::array<double, 3> &boxMax) const override {
    return getParticleImpl<true>(cellIndex, particleIndex, iteratorBehavior, boxMin, boxMax);
  }
  std::tuple<const Particle_T *, size_t, size_t> getParticle(size_t cellIndex, size_t particleIndex,
                                                             IteratorBehavior iteratorBehavior) const override {
    // this is not a region iter hence we stretch the bounding box to the numeric max
    constexpr std::array<double, 3> boxMin{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
                                           std::numeric_limits<double>::lowest()};
 
    constexpr std::array<double, 3> boxMax{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(),
                                           std::numeric_limits<double>::max()};
    return getParticleImpl<false>(cellIndex, particleIndex, iteratorBehavior, boxMin, boxMax);
  }
 
  bool deleteParticle(Particle_T &particle) override {
    // swap-delete helper function
    auto swapDelFromCell = [&](auto &particleCell) -> bool {
      auto &particleVec = particleCell._particles;
      const bool isRearParticle = &particle == &particleVec.back();
      particle = particleVec.back();
      particleVec.pop_back();
      return isRearParticle;
    };
 
    // deduce into which vector the reference points
    if (particle.isOwned()) {
      return swapDelFromCell(getOwnedCell());
    } else if (particle.isHalo()) {
      const auto boxMin = this->getBoxMin();
      const auto boxMax = this->getBoxMax();
      const auto skinHalf = 0.5 * this->getVerletSkin();
      const auto pos = particle.getR();
 
      // Look for the particle in halo cells that are within half the skinHalf distance of its position
      for (size_t dim = 0; dim < 3; ++dim) {
        if (pos[dim] < boxMin[dim] + skinHalf) {
          if (swapDelFromCell(this->_cells[2 * dim + 1])) {
            return true;
          }
        } else if (pos[dim] >= boxMax[dim] - skinHalf) {
          if (swapDelFromCell(this->_cells[2 * dim + 2])) {
            return true;
          }
        }
      }
    }
    return false;
  }
 
  bool deleteParticle(size_t cellIndex, size_t particleIndex) override {
    auto &particleVec = this->_cells[cellIndex]._particles;
    auto &particle = particleVec[particleIndex];
    // swap-delete
    particle = particleVec.back();
    particleVec.pop_back();
    return particleIndex < particleVec.size();
  }
 
 private:
  class DirectSumCellBorderAndFlagManager : public internal::CellBorderAndFlagManager {
    using index_t = std::size_t;
 
   public:
    [[nodiscard]] bool cellCanContainHaloParticles(index_t index1d) const override {
      return index1d >= 1 and index1d <= 6;
    }
 
    [[nodiscard]] bool cellCanContainOwnedParticles(index_t index1d) const override { return index1d == 0; }
 
  } _cellBorderFlagManager;
 
  template <bool regionIter>
  std::tuple<const Particle_T *, size_t, size_t> getParticleImpl(size_t cellIndex, size_t particleIndex,
                                                                 IteratorBehavior iteratorBehavior,
                                                                 const std::array<double, 3> &boxMin,
                                                                 const std::array<double, 3> &boxMax) const {
    // first and last relevant cell index
    const auto [startCellIndex, endCellIndex] = [&]() -> std::tuple<size_t, size_t> {
      // shortcuts to limit the iterator to only part of the domain.
      if (not(iteratorBehavior & IteratorBehavior::halo)) {
        // only owned region
        return {0, 0};
      }
      if (not(iteratorBehavior & IteratorBehavior::owned)) {
        // only halo region
        return {1, 6};
      }
      if constexpr (regionIter) {
        // if the region lies fully within the container only look at the owned cell
        if (utils::ArrayMath::less(this->getBoxMin(), boxMin) and utils::ArrayMath::less(boxMax, this->getBoxMax())) {
          return {0, 0};
        }
      }
      // all cells
      return {0, 6};
    }();
 
    // if we are at the start of an iteration determine this thread's cell index
    if (cellIndex == 0 and particleIndex == 0) {
      cellIndex =
          startCellIndex + ((iteratorBehavior & IteratorBehavior::forceSequential) ? 0 : autopas_get_thread_num());
    }
    // abort if the index is out of bounds
    if (cellIndex >= this->_cells.size()) {
      return {nullptr, 0, 0};
    }
    // check the data behind the indices
    if (particleIndex >= this->_cells[cellIndex].size() or
        not containerIteratorUtils::particleFulfillsIteratorRequirements<regionIter>(
            this->_cells[cellIndex][particleIndex], iteratorBehavior, boxMin, boxMax)) {
      // either advance them to something interesting or invalidate them.
      std::tie(cellIndex, particleIndex) =
          advanceIteratorIndices<regionIter>(cellIndex, particleIndex, iteratorBehavior, boxMin, boxMax);
    }
 
    // shortcut if the given index doesn't exist
    if (cellIndex >= this->_cells.size()) {
      return {nullptr, 0, 0};
    }
    const Particle_T *retPtr = &this->_cells[cellIndex][particleIndex];
 
    return {retPtr, cellIndex, particleIndex};
  }
 
  template <bool regionIter>
  std::tuple<size_t, size_t> advanceIteratorIndices(size_t cellIndex, size_t particleIndex,
                                                    IteratorBehavior iteratorBehavior,
                                                    const std::array<double, 3> &boxMin,
                                                    const std::array<double, 3> &boxMax) const {
    // Find the indices for the next particle
    const size_t stride = (iteratorBehavior & IteratorBehavior::forceSequential) ? 1 : autopas_get_num_threads();
 
    do {
      // advance to the next particle
      ++particleIndex;
      // If this breaches the end of a cell, find the next non-empty cell and reset particleIndex.
      while (particleIndex >= this->_cells[cellIndex].size()) {
        cellIndex += stride;
        particleIndex = 0;
        // If there are no more reasonable cells return invalid indices.
        if (cellIndex > ((not(iteratorBehavior & IteratorBehavior::halo)) ? 0 : (this->_cells.size() - 1))) {
          return {std::numeric_limits<decltype(cellIndex)>::max(), std::numeric_limits<decltype(particleIndex)>::max()};
        }
      }
    } while (not containerIteratorUtils::particleFulfillsIteratorRequirements<regionIter>(
        this->_cells[cellIndex][particleIndex], iteratorBehavior, boxMin, boxMax));
 
    // the indices returned at this point should always be valid
    return {cellIndex, particleIndex};
  }
 
  template <typename Traversal>
  void prepareTraversal(Traversal &traversal) {
    auto *dsTraversal = dynamic_cast<DSTraversalInterface *>(traversal);
    auto *cellTraversal = dynamic_cast<CellTraversal<ParticleCellType> *>(traversal);
    if (dsTraversal && cellTraversal) {
      cellTraversal->setAoSSortingThreshold(this->_aosSortingThreshold);
      cellTraversal->setSoASortingThreshold(this->_soaSortingThreshold);
      cellTraversal->setCellsToTraverse(this->_cells);
    } else {
      utils::ExceptionHandler::exception(
          "The selected traversal is not compatible with the DirectSum container. TraversalID: {}",
          traversal->getTraversalType());
    }
  }
 
  ParticleCellType &getOwnedCell() { return this->_cells[0]; };
 
  ParticleCellType &getHaloCell() { return this->_cells[1]; };
};
 
}  // namespace autopas