VerletClusterLists.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <iterator>
 
#include "ClusterTowerBlock2D.h"
#include "autopas/cells/FullParticleCell.h"
#include "autopas/containers/CompatibleTraversals.h"
#include "autopas/containers/LeavingParticleCollector.h"
#include "autopas/containers/NeighborListsBuffer.h"
#include "autopas/containers/ParticleContainerInterface.h"
#include "autopas/containers/ParticleDeletedObserver.h"
#include "autopas/containers/cellTraversals/BalancedTraversal.h"
#include "autopas/containers/verletClusterLists/Cluster.h"
#include "autopas/containers/verletClusterLists/ClusterTower.h"
#include "autopas/containers/verletClusterLists/VerletClusterListsRebuilder.h"
#include "autopas/containers/verletClusterLists/traversals/VCLTraversalInterface.h"
#include "autopas/iterators/ContainerIterator.h"
#include "autopas/options/LoadEstimatorOption.h"
#include "autopas/particles/OwnershipState.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/Timer.h"
#include "autopas/utils/WrapOpenMP.h"
#include "autopas/utils/inBox.h"
#include "autopas/utils/markParticleAsDeleted.h"
 
namespace autopas {
 
template <class Particle_T>
class VerletClusterLists : public ParticleContainerInterface<Particle_T>, public internal::ParticleDeletedObserver {
 public:
  using ParticleType = Particle_T;
  using ParticleCellType = internal::ClusterTower<Particle_T>;
 
  struct ClusterRange {
    size_t startTowerIndex{};
    size_t startIndexInTower{};
    size_t numClusters{};
  };
 
  VerletClusterLists(const std::array<double, 3> &boxMin, const std::array<double, 3> &boxMax, double cutoff,
                     double skin, size_t clusterSize, LoadEstimatorOption loadEstimator = LoadEstimatorOption::none)
      : ParticleContainerInterface<Particle_T>(skin),
        _towerBlock{boxMin, boxMax, cutoff + skin},
        _clusterSize{clusterSize},
        _particlesToAdd(autopas_get_max_threads()),
        _cutoff{cutoff},
        _loadEstimator(loadEstimator) {
    // always have at least one tower.
    _towerBlock.addTower(_clusterSize);
  }
 
  [[nodiscard]] ContainerOption getContainerType() const override { return ContainerOption::verletClusterLists; }
 
  BalancedTraversal::EstimatorFunction getLoadEstimatorFunction() {
    // (Explicit) static cast required for Apple Clang (last tested version: 17.0.0)
    switch (static_cast<LoadEstimatorOption::Value>(this->_loadEstimator)) {
      case LoadEstimatorOption::neighborListLength: {
        return [&](const std::array<unsigned long, 3> &cellsPerDimension,
                   const std::array<unsigned long, 3> &lowerCorner, const std::array<unsigned long, 3> &upperCorner) {
          // the neighborListLength function defined for verletListsCells in not compatible with this container.
          unsigned long sum = 0;
          for (unsigned long x = lowerCorner[0]; x <= upperCorner[0]; x++) {
            for (unsigned long y = lowerCorner[1]; y <= upperCorner[1]; y++) {
              unsigned long cellLoad = 0;
              auto &tower = _towerBlock.getTowerByIndex2D(x, y);
              for (auto &cluster : tower.getClusters()) {
                if (cluster.getNeighbors()) {
                  cellLoad += cluster.getNeighbors()->size();
                }
              }
              sum += cellLoad;
            }
          }
          return sum;
        };
      }
      case LoadEstimatorOption::none:
        [[fallthrough]];
      default: {
        return
            [&](const std::array<unsigned long, 3> &cellsPerDimension, const std::array<unsigned long, 3> &lowerCorner,
                const std::array<unsigned long, 3> &upperCorner) { return 1; };
      }
    }
  }
 
  void computeInteractions(TraversalInterface *traversal) override {
    if (_isValid == ValidityState::cellsAndListsValid) {
      utils::ExceptionHandler::exception(
          "VerletClusterLists::computeInteractions(): Trying to do a pairwise iteration, even though verlet lists are "
          "not valid.");
    }
    auto *traversalInterface = dynamic_cast<VCLTraversalInterface<Particle_T> *>(traversal);
    if (traversalInterface) {
      traversalInterface->setClusterLists(*this);
      traversalInterface->setTowers(_towerBlock.getTowersRef());
    } else {
      utils::ExceptionHandler::exception(
          "Trying to use a traversal of wrong type in VerletClusterLists::computeInteractions. TraversalID: {}",
          traversal->getTraversalType());
    }
    if (auto *balancedTraversal = dynamic_cast<BalancedTraversal *>(traversal)) {
      balancedTraversal->setLoadEstimator(getLoadEstimatorFunction());
    }
 
    traversal->initTraversal();
    traversal->traverseParticles();
    traversal->endTraversal();
  }
 
  void setAoSSortingThresholds(std::shared_ptr<const SortingThresholdInfoInterface> aosSortingThresholds) override {}
 
  void setSoASortingThresholds(std::shared_ptr<const SortingThresholdInfoInterface> soaSortingThresholds) override {}
 
  void reserve(size_t numParticles, size_t numParticlesHaloEstimate) override {
    const auto particlesPerTower = (numParticles + numParticlesHaloEstimate) / _towerBlock.size();
    for (auto &tower : _towerBlock) {
      tower.reserve(particlesPerTower);
    }
  }
 
  void addParticleImpl(const Particle_T &p) override {
    _isValid.store(ValidityState::invalid, std::memory_order_relaxed);
    _particlesToAdd[autopas_get_thread_num()].push_back(p);
  }
 
  void addHaloParticleImpl(const Particle_T &haloParticle) override {
    _isValid.store(ValidityState::invalid, std::memory_order_relaxed);
    _particlesToAdd[autopas_get_thread_num()].push_back(haloParticle);
  }
 
  bool updateHaloParticle(const Particle_T &haloParticle) override {
    using namespace autopas::utils::ArrayMath::literals;
 
    const auto &haloPos = haloParticle.getR();
    // this might be called from a parallel region so force this iterator to be sequential
    for (auto it = getRegionIterator(haloPos - (this->getVerletSkin() / 2.), haloPos + (this->getVerletSkin() / 2.),
                                     IteratorBehavior::halo | IteratorBehavior::forceSequential, std::nullopt);
         it.isValid(); ++it) {
      if (haloParticle.getID() == it->getID()) {
        // don't simply copy haloParticle over iter. This would trigger a dataRace with other regionIterators that
        // overlap with this region.
        it->setR(haloPos);
        it->setV(haloParticle.getV());
        it->setF(haloParticle.getF());
        return true;
      }
    }
    return false;
  }
 
  void deleteHaloParticles() override {
    // Step 1: Remove particles from _particlesToAdd
    for (auto &particleVec : _particlesToAdd) {
      for (size_t j = 0; j < particleVec.size();) {
        if (particleVec[j].isHalo()) {
          particleVec[j] = particleVec[particleVec.size() - 1];
          particleVec.pop_back();
        } else {
          ++j;
        }
      }
    }
    // Step 2: Remove particles from _towers
    bool deletedSomething = false;
    AUTOPAS_OPENMP(parallel for reduction(|| : deletedSomething))
    // Thanks to clang 13 this has to be a index based loop instead of a range based
    for (size_t i = 0; i < _towerBlock.size(); ++i) {
      auto &tower = _towerBlock[i];
      const auto towerSize = tower.size();
      auto numTailDummies = tower.getNumTailDummyParticles();
      // iterate over all non-tail dummies. Avoid underflows.
      for (size_t j = 0; numTailDummies < towerSize and j < towerSize - numTailDummies;) {
        if (tower[j].isHalo()) {
          // swap-"delete"
          tower[j] = tower[towerSize - 1 - numTailDummies];
          // Since we can't pop the moved particle here mark it for deletion.
          internal::markParticleAsDeleted(tower[towerSize - 1 - numTailDummies]);
          ++numTailDummies;
          deletedSomething = true;
        } else {
          ++j;
        }
      }
      // if anything was marked for deletion actually delete it now.
      if (deletedSomething) {
        tower.deleteDummyParticles();
      }
    }
    if (deletedSomething) {
      _isValid.store(ValidityState::invalid, std::memory_order_relaxed);
    }
  }
 
  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);
  }
 
  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 {
    using namespace autopas::utils::ArrayMath::literals;
 
    // in this context cell == tower
    // catching the edge case that towers are not yet built -> Iterator jumps to additional vectors
    if (_towerBlock.empty()) {
      return {nullptr, 0, 0};
    }
 
    std::array<double, 3> boxMinWithSafetyMargin = boxMin;
    std::array<double, 3> boxMaxWithSafetyMargin = boxMax;
    if constexpr (regionIter) {
      // We extend the search box for cells here since particles might have moved
      boxMinWithSafetyMargin -= 0.5 * this->getVerletSkin();
      boxMaxWithSafetyMargin += 0.5 * this->getVerletSkin();
    }
 
    // first and last relevant cell index
    const auto [startCellIndex, endCellIndex] = [&]() -> std::tuple<size_t, size_t> {
      if constexpr (regionIter) {
        // We extend the search box for cells here since particles might have moved
        return {_towerBlock.getTowerIndex1DAtPosition(boxMinWithSafetyMargin),
                _towerBlock.getTowerIndex1DAtPosition(boxMaxWithSafetyMargin)};
      } else {
        if (not(iteratorBehavior & IteratorBehavior::halo)) {
          // only potentially owned region
          return {_towerBlock.getFirstOwnedTowerIndex(), _towerBlock.getLastOwnedTowerIndex()};
        } else {
          // whole range of cells
          return {0, _towerBlock.size() - 1};
        }
      }
    }();
 
    // if we are at the start of an iteration ...
    if (cellIndex == 0 and particleIndex == 0) {
      cellIndex =
          startCellIndex + ((iteratorBehavior & IteratorBehavior::forceSequential) ? 0 : autopas_get_thread_num());
    }
    // abort if the start index is already out of bounds
    if (cellIndex >= _towerBlock.size()) {
      return {nullptr, 0, 0};
    }
    // check the data behind the indices
    if (particleIndex >= _towerBlock[cellIndex].getNumActualParticles() or
        not containerIteratorUtils::particleFulfillsIteratorRequirements<regionIter>(
            _towerBlock[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,
                                             boxMinWithSafetyMargin, boxMaxWithSafetyMargin, endCellIndex);
    }
    // shortcut if the given index doesn't exist
    if (cellIndex > endCellIndex) {
      return {nullptr, 0, 0};
    }
    const Particle_T *retPtr = &_towerBlock[cellIndex][particleIndex];
 
    return {retPtr, cellIndex, particleIndex};
  }
 
  bool deleteParticle(Particle_T &particle) override {
    // This function doesn't actually delete anything as it would mess up the clusters' references.
    internal::markParticleAsDeleted(particle);
    return false;
  }
 
  bool deleteParticle(size_t cellIndex, size_t particleIndex) override {
    // This function doesn't actually delete anything as it would mess up the clusters' references.
    internal::markParticleAsDeleted(_towerBlock[cellIndex][particleIndex]);
    return false;
  }
 
  [[nodiscard]] std::vector<Particle_T> updateContainer(bool keepNeighborListsValid) override {
    if (keepNeighborListsValid) {
      return LeavingParticleCollector::collectParticlesAndMarkNonOwnedAsDummy(*this);
    }
    // First delete all halo particles.
    this->deleteHaloParticles();
    // Delete dummy particles.
    AUTOPAS_OPENMP(parallel for)
    for (size_t i = 0ul; i < _towerBlock.size(); ++i) {
      _towerBlock[i].deleteDummyParticles();
    }
    // Delete dummy particles also from the _particlesToAdd vector
    AUTOPAS_OPENMP(parallel for)
    for (size_t i = 0ul; i < _particlesToAdd.size(); ++i) {
      _particlesToAdd[i].erase(std::remove_if(_particlesToAdd[i].begin(), _particlesToAdd[i].end(),
                                              [](const auto &p) { return p.isDummy(); }),
                               _particlesToAdd[i].end());
    }
 
    // next find invalid particles
    std::vector<Particle_T> invalidParticles;
 
    // custom openmp reduction to concatenate all local vectors to one at the end of a parallel region
    AUTOPAS_OPENMP(declare reduction(
        vecMergeParticle : std::vector<Particle_T> : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())))
    AUTOPAS_OPENMP(parallel reduction(vecMergeParticle : invalidParticles)) {
      for (auto iter = this->begin(IteratorBehavior::owned); iter.isValid(); ++iter) {
        if (not utils::inBox(iter->getR(), this->getBoxMin(), this->getBoxMax())) {
          invalidParticles.push_back(*iter);
          internal::deleteParticle(iter);
        }
      }
    }
    _isValid.store(ValidityState::invalid, std::memory_order_relaxed);
    return invalidParticles;
  }
 
  [[nodiscard]] TraversalSelectorInfo getTraversalSelectorInfo() const override {
    using namespace autopas::utils::ArrayMath::literals;
    // Here, the towers might not yet be built, hence do not use members like _towerBlock.getTowersPerDim
    const auto boxSizeWithHalo = this->getHaloBoxMax() - this->getHaloBoxMin();
    const auto [towerSideLength, towersPerDim] = _towerBlock.estimateOptimalGridSideLength(
        this->getNumberOfParticles(IteratorBehavior::ownedOrHalo), _clusterSize);
    const std::array<double, 3> towerSize = {towerSideLength[0], towerSideLength[1],
                                             this->getHaloBoxMax()[2] - this->getHaloBoxMin()[2]};
    const std::array<unsigned long, 3> towerDimensions = {towersPerDim[0], towersPerDim[1], 1};
    return TraversalSelectorInfo(towerDimensions, this->getInteractionLength(), towerSize, _clusterSize);
  }
 
  [[nodiscard]] ContainerIterator<Particle_T, true, false> begin(
      IteratorBehavior behavior = IteratorBehavior::ownedOrHalo,
      utils::optRef<typename ContainerIterator<Particle_T, true, false>::ParticleVecType> additionalVectors =
          std::nullopt) override {
    // Note: particlesToAddEmpty() can only be called if the container status is not invalid. If the status is set to
    // invalid, we do writing operations on _particlesToAdd and can not read from it without race conditions.
    if (_isValid != ValidityState::invalid) {
      // we call particlesToAddEmpty() as a sanity check to ensure there are actually no particles in _particlesToAdd if
      // the status is not invalid
      if (not particlesToAddEmpty(autopas_get_thread_num())) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::begin(): Error: particle container is valid, but _particlesToAdd isn't empty!");
      }
      // If the particles are sorted into the towers, we can simply use the iteration over towers + additionalVectors
      // from LogicHandler.
      return ContainerIterator<Particle_T, true, false>(*this, behavior, additionalVectors);
    } else {
      // if the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      // store all pointers in a temporary which is passed to the ParticleIterator constructor.
      typename ContainerIterator<Particle_T, true, false>::ParticleVecType additionalVectorsToPass;
      appendBuffersHelper(additionalVectors, additionalVectorsToPass);
      return ContainerIterator<Particle_T, true, false>(*this, behavior, std::ref(additionalVectorsToPass));
    }
  }
 
  [[nodiscard]] ContainerIterator<Particle_T, false, false> begin(
      IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo,
      utils::optRef<typename ContainerIterator<Particle_T, false, false>::ParticleVecType> additionalVectors =
          std::nullopt) const override {
    // Note: particlesToAddEmpty() can only be called if the container status is not invalid. If the status is set to
    // invalid, we do writing operations on _particlesToAdd and can not read from from it without race conditions.
    if (_isValid != ValidityState::invalid) {
      // we call particlesToAddEmpty() as a sanity check to ensire there are actually no particles in _particlesToAdd if
      // the status is not invalid
      if (not particlesToAddEmpty(autopas_get_thread_num())) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::begin() const: Error: particle container is valid, but _particlesToAdd isn't empty!");
      }
      // If the particles are sorted into the towers, we can simply use the iteration over towers + additionalVectors
      // from LogicHandler.
      return ContainerIterator<Particle_T, false, false>(*this, behavior, additionalVectors);
    } else {
      // if the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      // store all pointers in a temporary which is passed to the ParticleIterator constructor.
      typename ContainerIterator<Particle_T, false, false>::ParticleVecType additionalVectorsToPass;
      appendBuffersHelper(additionalVectors, additionalVectorsToPass);
      return ContainerIterator<Particle_T, false, false>(*this, behavior, std::ref(additionalVectorsToPass));
    }
  }
 
  template <typename Lambda>
  void forEach(Lambda forEachLambda, IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo) {
    for (auto &tower : _towerBlock) {
      tower.forEach(forEachLambda, behavior);
    }
    for (auto &vector : this->_particlesToAdd) {
      for (auto &particle : vector) {
        if (behavior.contains(particle)) {
          forEachLambda(particle);
        }
      }
    }
  }
 
  template <typename Lambda>
  void forEach(Lambda forEachLambda, IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo) const {
    if (_isValid != ValidityState::invalid) {
      if (not particlesToAddEmpty()) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::forEach() const: Error: particle container is valid, but _particlesToAdd isn't "
            "empty!");
      }
    }
 
    // If the particles are sorted into the towers, we can simply use the iteration over towers.
    for (auto &tower : _towerBlock) {
      tower.forEach(forEachLambda, behavior);
    }
 
    // if the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
    if (_isValid == ValidityState::invalid) {
      for (auto &particlesToAddPerThread : _particlesToAdd) {
        for (auto &particle : particlesToAddPerThread) {
          if (behavior.contains(particle)) {
            forEachLambda(particle);
          }
        }
      }
    }
  }
 
  template <typename Lambda, typename A>
  void reduce(Lambda reduceLambda, A &result, IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo) {
    for (auto &tower : _towerBlock) {
      tower.reduce(reduceLambda, result, behavior);
    }
    for (auto &vector : this->_particlesToAdd) {
      for (auto &p : vector) {
        if (behavior.contains(p)) {
          reduceLambda(p, result);
        }
      }
    }
  }
 
  template <typename Lambda, typename A>
  void reduce(Lambda reduceLambda, A &result,
              IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo) const {
    if (_isValid != ValidityState::invalid) {
      if (not particlesToAddEmpty()) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::reduce() const: Error: particle container is valid, but _particlesToAdd isn't empty!");
      }
    }
 
    // If the particles are sorted into the towers, we can simply use the iteration over towers.
    for (auto tower : _towerBlock) {
      tower.reduce(reduceLambda, result, behavior);
    }
 
    if (_isValid == ValidityState::invalid) {
      // if the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      for (auto &particlesToAddPerThread : _particlesToAdd) {
        for (auto &particle : particlesToAddPerThread) {
          if (behavior.contains(particle)) {
            reduceLambda(particle, result);
          }
        }
      }
    }
  }
 
  [[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 {
    // Note: particlesToAddEmpty() can only be called if the container status is not invalid. If the status is set to
    // invalid, we do writing operations on _particlesToAdd and can not read from from it without race conditions.
    if (_isValid != ValidityState::invalid) {
      // we call particlesToAddEmpty() as a sanity check to ensure there are actually no particles in _particlesToAdd
      // if the status is not invalid
      if (not particlesToAddEmpty(autopas_get_thread_num())) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::reduce() const: Error: particle container is valid, but _particlesToAdd isn't empty!");
      }
      // If the particles are sorted into the towers, we can simply use the iteration over towers + additionalVectors
      // from LogicHandler.
      return ContainerIterator<Particle_T, true, true>(*this, behavior, additionalVectors, lowerCorner, higherCorner);
    } else {
      // if the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      // store all pointers in a temporary which is passed to the ParticleIterator constructor.
      typename ContainerIterator<Particle_T, true, true>::ParticleVecType additionalVectorsToPass;
      appendBuffersHelper(additionalVectors, additionalVectorsToPass);
      return ContainerIterator<Particle_T, true, true>(*this, behavior, std::ref(additionalVectorsToPass), 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 {
    // Note: particlesToAddEmpty() can only be called if the container status is not invalid. If the status is set to
    // invalid, we do writing operations on _particlesToAdd and can not read from from it without race conditions.
    if (_isValid != ValidityState::invalid) {
      // we call particlesToAddEmpty() as a sanity check to ensire there are actually no particles in _particlesToAdd if
      // the status is not invalid
      if (not particlesToAddEmpty(autopas_get_thread_num())) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::getRegionIterator() const: Error: particle container is valid, but _particlesToAdd "
            "isn't empty!");
      }
      // If the particles are sorted into the towers, we can simply use the iteration over towers + additionalVectors
      // from LogicHandler.
      return ContainerIterator<Particle_T, false, true>(*this, behavior, additionalVectors, lowerCorner, higherCorner);
    } else {
      // if the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      // store all pointers in a temporary which is passed to the ParticleIterator constructor.
      typename ContainerIterator<Particle_T, false, true>::ParticleVecType additionalVectorsToPass;
      appendBuffersHelper(additionalVectors, additionalVectorsToPass);
      return ContainerIterator<Particle_T, false, true>(*this, behavior, std::ref(additionalVectorsToPass), lowerCorner,
                                                        higherCorner);
    }
  }
 
  template <typename Lambda>
  void forEachInRegion(Lambda forEachLambda, const std::array<double, 3> &lowerCorner,
                       const std::array<double, 3> &higherCorner,
                       IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo) {
    for (size_t i = 0; i < _towerBlock.size(); ++i) {
      if (_towerBlock.ignoreCellForIteration(i, behavior)) {
        continue;
      }
      auto &tower = _towerBlock[i];
      const auto [towerLowCorner, towerHighCorner] = _towerBlock.getTowerBoundingBox(i);
      // particles can move over cell borders. Calculate the volume this cell's particles can be.
      const auto towerLowCornerSkin = utils::ArrayMath::subScalar(towerLowCorner, this->getVerletSkin() * 0.5);
      const auto towerHighCornerSkin = utils::ArrayMath::addScalar(towerHighCorner, this->getVerletSkin() * 0.5);
      if (utils::boxesOverlap(towerLowCornerSkin, towerHighCornerSkin, lowerCorner, higherCorner)) {
        tower.forEach(forEachLambda, lowerCorner, higherCorner, behavior);
      }
    }
    for (auto &vector : _particlesToAdd) {
      for (auto &particle : vector) {
        if (behavior.contains(particle)) {
          if (utils::inBox(particle.getR(), lowerCorner, higherCorner)) {
            forEachLambda(particle);
          }
        }
      }
    }
  }
 
  template <typename Lambda>
  void forEachInRegion(Lambda forEachLambda, const std::array<double, 3> &lowerCorner,
                       const std::array<double, 3> &higherCorner,
                       IteratorBehavior behavior = autopas::IteratorBehavior::ownedOrHalo) const {
    if (_isValid != ValidityState::invalid) {
      if (not particlesToAddEmpty()) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::forEachInRegion() const: Error: particle container is valid, but _particlesToAdd "
            "isn't empty!");
      }
    }
 
    // If the particles are sorted into the towers, we can simply use the iteration over towers.
    for (size_t i = 0; i < _towerBlock.size(); ++i) {
      if (_towerBlock.ignoreCellForIteration(i, behavior)) {
        continue;
      }
      auto &tower = _towerBlock[i];
      const auto [towerLowCorner, towerHighCorner] = _towerBlock.getTowerBoundingBox(i);
      // particles can move over cell borders. Calculate the volume this cell's particles can be.
      const auto towerLowCornerSkin = utils::ArrayMath::subScalar(towerLowCorner, this->getVerletSkin() * 0.5);
      const auto towerHighCornerSkin = utils::ArrayMath::addScalar(towerHighCorner, this->getVerletSkin() * 0.5);
      if (utils::boxesOverlap(towerLowCornerSkin, towerHighCornerSkin, lowerCorner, higherCorner)) {
        tower.forEach(forEachLambda, lowerCorner, higherCorner, behavior);
      }
    }
 
    if (_isValid == ValidityState::invalid) {
      // If the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      for (auto &particlesToAddPerThread : _particlesToAdd) {
        for (auto &particle : particlesToAddPerThread) {
          if (behavior.contains(particle)) {
            if (utils::inBox(particle.getR(), lowerCorner, higherCorner)) {
              forEachLambda(particle);
            }
          }
        }
      }
    }
  }
 
  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 = autopas::IteratorBehavior::ownedOrHalo) {
    for (size_t i = 0; i < _towerBlock.size(); ++i) {
      if (_towerBlock.ignoreCellForIteration(i, behavior)) {
        continue;
      }
      auto &tower = _towerBlock[i];
      const auto [towerLowCorner, towerHighCorner] = _towerBlock.getTowerBoundingBox(i);
      // particles can move over cell borders. Calculate the volume this cell's particles can be.
      const auto towerLowCornerSkin = utils::ArrayMath::subScalar(towerLowCorner, this->getVerletSkin() * 0.5);
      const auto towerHighCornerSkin = utils::ArrayMath::addScalar(towerHighCorner, this->getVerletSkin() * 0.5);
      if (utils::boxesOverlap(towerLowCornerSkin, towerHighCornerSkin, lowerCorner, higherCorner)) {
        tower.reduce(reduceLambda, result, lowerCorner, higherCorner, behavior);
      }
    }
    for (auto &vector : _particlesToAdd) {
      for (auto &particle : vector) {
        if (behavior.contains(particle)) {
          if (utils::inBox(particle.getR(), lowerCorner, higherCorner)) {
            reduceLambda(particle, result);
          }
        }
      }
    }
  }
 
  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 = autopas::IteratorBehavior::ownedOrHalo) const {
    if (_isValid != ValidityState::invalid) {
      if (not particlesToAddEmpty()) {
        autopas::utils::ExceptionHandler::exception(
            "VerletClusterLists::reduceInRegion() const: Error: particle container is valid, but _particlesToAdd isn't "
            "empty!");
      }
    }
    // If the particles are sorted into the towers, we can simply use the iteration over towers.
    for (size_t i = 0; i < _towerBlock.size(); ++i) {
      if (_towerBlock.ignoreCellForIteration(i, behavior)) {
        continue;
      }
      auto &tower = _towerBlock[i];
      const auto [towerLowCorner, towerHighCorner] = _towerBlock.getTowerBoundingBox(i);
      // particles can move over cell borders. Calculate the volume this cell's particles can be.
      const auto towerLowCornerSkin = utils::ArrayMath::subScalar(towerLowCorner, this->getVerletSkin() * 0.5);
      const auto towerHighCornerSkin = utils::ArrayMath::addScalar(towerHighCorner, this->getVerletSkin() * 0.5);
      if (utils::boxesOverlap(towerLowCornerSkin, towerHighCornerSkin, lowerCorner, higherCorner)) {
        tower.reduce(reduceLambda, result, lowerCorner, higherCorner, behavior);
      }
    }
 
    if (_isValid == ValidityState::invalid) {
      // If the particles are not sorted into the towers, we have to also iterate over _particlesToAdd.
      for (auto &particlesToAddPerThread : _particlesToAdd) {
        for (auto &particle : particlesToAddPerThread) {
          if (behavior.contains(particle)) {
            if (utils::inBox(particle.getR(), lowerCorner, higherCorner)) {
              reduceLambda(particle, result);
            }
          }
        }
      }
    }
  }
 
  void rebuildNeighborLists(TraversalInterface *traversal) override {
    // The builder might have a different newton3 choice than the traversal. This typically only happens in unit tests
    // when rebuildTowersAndClusters() was not called explicitly.
    if (_isValid == ValidityState::invalid or traversal->getUseNewton3() != _builder->getNewton3()) {
      // clear the lists buffer because clusters will be recreated
      _neighborLists.clear();
      rebuildTowersAndClusters(traversal->getUseNewton3());
    }
    _builder->rebuildNeighborListsAndFillClusters();
 
    auto *clusterTraversalInterface = dynamic_cast<VCLTraversalInterface<Particle_T> *>(traversal);
    if (clusterTraversalInterface) {
      if (clusterTraversalInterface->needsStaticClusterThreadPartition()) {
        calculateClusterThreadPartition();
      }
    } else {
      utils::ExceptionHandler::exception(
          "Trying to use a traversal of wrong type in VerletClusterLists::rebuildNeighborLists. TraversalID: {}",
          traversal->getTraversalType());
    }
  }
 
  template <bool inParallel, class LoopBody>
  void traverseClusters(LoopBody &&loopBody) {
    if (inParallel) {
      traverseClustersParallel<LoopBody>(std::forward<LoopBody>(loopBody));
    } else {
      traverseClustersSequential<LoopBody>(std::forward<LoopBody>(loopBody));
    }
  }
 
  [[nodiscard]] size_t size() const override {
    size_t sum = std::accumulate(_towerBlock.begin(), _towerBlock.end(), 0,
                                 [](size_t acc, const auto &tower) { return acc + tower.size(); });
    sum = std::accumulate(_particlesToAdd.begin(), _particlesToAdd.end(), sum,
                          [](size_t acc, const auto &buffer) { return acc + buffer.size(); });
    return sum;
  }
 
  [[nodiscard]] size_t getNumberOfParticles(IteratorBehavior behavior) const override {
    // sum up all particles in towers that fulfill behavior
    size_t sum = std::accumulate(_towerBlock.begin(), _towerBlock.end(), 0, [&behavior](size_t acc, const auto &tower) {
      return acc + tower.getNumberOfParticles(behavior);
    });
 
    // Since we can not directly insert particles into towers without a rebuild of the whole data structure,
    // _particlesToAdd is used to store all these particles temporarily until the next rebuild inserts them into the
    // towers data structure. However, these particles already belong to the respective tower, so we have to count them
    // as well.
    sum = std::accumulate(
        _particlesToAdd.begin(), _particlesToAdd.end(), sum, [&behavior](size_t acc, const auto &buffer) {
          return acc +
                 (std::count_if(buffer.begin(), buffer.end(), [&behavior](auto p) { return behavior.contains(p); }));
        });
 
    return sum;
  }
 
  const auto &getClusterThreadPartition() const { return _clusterThreadPartition; }
 
  auto getNumClusters() const { return _numClusters; }
 
  auto getTowerSideLength() const { return _towerBlock.getTowerSideLength(); }
 
  auto getTowersPerDimension() const { return _towerBlock.getTowersPerDim(); }
 
  auto getClusterSize() const { return _clusterSize; }
 
  auto getNumTowersPerInteractionLength() const { return _towerBlock.getNumTowersPerInteractionLength(); }
 
  template <class Functor>
  void loadParticlesIntoSoAs(Functor &functor) {
    const auto numTowers = _towerBlock.size();
    AUTOPAS_OPENMP(parallel for schedule(dynamic))
    for (size_t index = 0; index < numTowers; index++) {
      _towerBlock[index].loadSoA(functor);
    }
  }
 
  template <class Functor>
  void extractParticlesFromSoAs(Functor &functor) {
    const auto numTowers = _towerBlock.size();
    AUTOPAS_OPENMP(parallel for schedule(dynamic))
    for (size_t index = 0; index < numTowers; index++) {
      _towerBlock[index].extractSoA(functor);
    }
  }
 
  internal::ClusterTower<Particle_T> &getTowerByIndex(size_t x, size_t y) {
    return _towerBlock.getTowerByIndex2D(x, y);
  }
 
  internal::ClusterTowerBlock2D<Particle_T> &getTowerBlock() { return _towerBlock; }
 
  [[nodiscard]] const std::array<double, 3> &getBoxMax() const override { return _towerBlock.getBoxMax(); }
 
  [[nodiscard]] const std::array<double, 3> &getHaloBoxMax() const { return _towerBlock.getHaloBoxMax(); }
 
  [[nodiscard]] const std::array<double, 3> &getBoxMin() const override { return _towerBlock.getBoxMin(); }
 
  [[nodiscard]] const std::array<double, 3> &getHaloBoxMin() const { return _towerBlock.getHaloBoxMin(); }
 
  [[nodiscard]] double getCutoff() const override { return _cutoff; }
 
  void setCutoff(double cutoff) override { _cutoff = cutoff; }
 
  [[nodiscard]] double getVerletSkin() const override { return this->_skin; }
 
  void setSkin(double skin) { this->_skin = skin; }
 
  [[nodiscard]] double getInteractionLength() const override { return _cutoff + this->getVerletSkin(); }
 
  void deleteAllParticles() override {
    _isValid.store(ValidityState::invalid, std::memory_order_relaxed);
    std::for_each(_particlesToAdd.begin(), _particlesToAdd.end(), [](auto &buffer) { buffer.clear(); });
    std::for_each(_towerBlock.begin(), _towerBlock.end(), [](auto &tower) { tower.clear(); });
  }
 
  const typename internal::VerletClusterListsRebuilder<Particle_T>::NeighborListsBuffer_T &getNeighborLists() const {
    return _neighborLists;
  }
 
  void rebuildTowersAndClusters(bool newton3) {
    using namespace utils::ArrayMath::literals;
    // collect all particles to add from across the thread buffers
    typename decltype(_particlesToAdd)::value_type particlesToAdd;
    const size_t numParticlesToAdd =
        std::accumulate(_particlesToAdd.begin(), _particlesToAdd.end(), 0,
                        [](size_t acc, const auto &buffer) { return acc + buffer.size(); });
    particlesToAdd.reserve(numParticlesToAdd);
    std::for_each(_particlesToAdd.begin(), _particlesToAdd.end(), [&](auto &particlesBuffer) {
      particlesToAdd.insert(particlesToAdd.end(), particlesBuffer.begin(), particlesBuffer.end());
      particlesBuffer.clear();
    });
 
    const double interactionLength = _cutoff + this->_skin;
    _builder = std::make_unique<internal::VerletClusterListsRebuilder<Particle_T>>(
        _towerBlock, particlesToAdd, _neighborLists, _clusterSize, interactionLength * interactionLength, newton3);
 
    _numClusters = _builder->rebuildTowersAndClusters();
 
    _isValid.store(ValidityState::cellsValidListsInvalid, std::memory_order_relaxed);
    for (auto &tower : _towerBlock) {
      tower.setParticleDeletionObserver(this);
    }
  }
 
  template <class LoopBody>
  void traverseClustersSequential(LoopBody &&loopBody) {
    for (size_t x = 0; x < _towerBlock.getTowersPerDim()[0]; x++) {
      for (size_t y = 0; y < _towerBlock.getTowersPerDim()[1]; y++) {
        auto &tower = _towerBlock.getTowerByIndex2D(x, y);
        for (auto clusterIter = tower.getFirstOwnedCluster(); clusterIter < tower.getFirstTailHaloCluster();
             ++clusterIter) {
          loopBody(*clusterIter);
        }
      }
    }
  }
 
  template <class LoopBody>
  void traverseClustersParallel(LoopBody &&loopBody) {
    const auto towersPerDimX = _towerBlock.getTowersPerDim()[0];
    const auto towersPerDimY = _towerBlock.getTowersPerDim()[1];
    AUTOPAS_OPENMP(parallel for schedule(dynamic) collapse(2))
    for (size_t x = 0; x < towersPerDimX; x++) {
      for (size_t y = 0; y < towersPerDimY; y++) {
        auto &tower = _towerBlock.getTowerByIndex2D(x, y);
 
        for (auto clusterIter = tower.getFirstOwnedCluster(); clusterIter < tower.getFirstTailHaloCluster();
             ++clusterIter) {
          loopBody(*clusterIter);
        }
      }
    }
  }
 
 protected:
  void calculateClusterThreadPartition() {
    size_t numClusterPairs = 0;
    this->template traverseClusters<false>(
        [&numClusterPairs](auto &cluster) { numClusterPairs += cluster.getNeighbors()->size(); });
 
    constexpr int minNumClusterPairsPerThread = 1000;
    auto numThreads =
        std::clamp(static_cast<int>(numClusterPairs / minNumClusterPairsPerThread), 1, autopas_get_max_threads());
 
    size_t numClusterPairsPerThread =
        std::max(static_cast<unsigned long>(std::ceil(static_cast<double>(numClusterPairs) / numThreads)), 1ul);
    if (numClusterPairsPerThread * numThreads < numClusterPairs) {
      autopas::utils::ExceptionHandler::exception(
          "VerletClusterLists::calculateClusterThreadPartition(): numClusterPairsPerThread ({}) * numThreads ({})={} "
          "should always "
          "be at least the amount of Cluster Pairs ({})!",
          numClusterPairsPerThread, numThreads, numClusterPairsPerThread * numThreads, numClusterPairs);
    }
    fillClusterRanges(numClusterPairsPerThread, numThreads);
  }
 
  void fillClusterRanges(size_t numClusterPairsPerThread, int numThreads) {
    if (numClusterPairsPerThread < 1) {
      autopas::utils::ExceptionHandler::exception(
          "VerletClusterLists::fillClusterRanges(): numClusterPairsPerThread({}) is less than one, this is not "
          "supported "
          "and will lead to errors!",
          numClusterPairsPerThread);
    }
    _clusterThreadPartition.resize(numThreads);
 
    size_t currentThread = 0;
    size_t currentNumClustersToAdd = 0;
    size_t numClusterPairsTotal = 0;
    bool threadIsInitialized = false;
    // Iterate over the clusters of all towers
    for (size_t currentTowerIndex = 0; currentTowerIndex < _towerBlock.size(); currentTowerIndex++) {
      auto &currentTower = _towerBlock[currentTowerIndex];
      const auto firstOwnedClusterIndex =
          std::distance(currentTower.getClusters().begin(), currentTower.getFirstOwnedCluster());
      const auto firstTailHaloClusterIndex =
          std::distance(currentTower.getClusters().begin(), currentTower.getFirstTailHaloCluster());
      for (size_t currentClusterInTower = firstOwnedClusterIndex; currentClusterInTower < firstTailHaloClusterIndex;
           ++currentClusterInTower) {
        auto &currentCluster = currentTower.getCluster(currentClusterInTower);
 
        // If on a new thread, start with the clusters for this thread here.
        if (not threadIsInitialized) {
          _clusterThreadPartition[currentThread] = {currentTowerIndex, currentClusterInTower, 0};
          threadIsInitialized = true;
        }
 
        currentNumClustersToAdd++;
        numClusterPairsTotal += currentCluster.getNeighbors()->size();
 
        // If the thread is finished, write number of clusters and start new thread.
        if (numClusterPairsTotal >= numClusterPairsPerThread * (currentThread + 1)) {
          // Add the number of clusters for the finished thread.
          _clusterThreadPartition[currentThread].numClusters += currentNumClustersToAdd;
          currentNumClustersToAdd = 0;
          // Go to next thread!
          currentThread++;
          // if we are already at the end of all threads, go back to last thread!
          // this is a safety precaution and should not really matter.
          if (currentThread >= numThreads) {
            --currentThread;
            threadIsInitialized = true;
          } else {
            threadIsInitialized = false;
          }
        }
      }
    }
    if (not threadIsInitialized) {
      _clusterThreadPartition[currentThread] = {0, 0, 0};
    }
    // Make sure the last cluster range contains the rest of the clusters, even if there is not the perfect number left.
    if (currentNumClustersToAdd != 0) {
      _clusterThreadPartition[currentThread].numClusters += currentNumClustersToAdd;
    }
    // Theoretically, some threads may still remain. This ensures that their numClusters are set to 0.
    while (++currentThread < numThreads) {
      _clusterThreadPartition[currentThread] = {0, 0, 0};
    }
  }
 
  void notifyParticleDeleted() override {
    // this is potentially called from a threaded environment, so we have to make this atomic here!
    _isValid.store(ValidityState::invalid, std::memory_order_relaxed);
  }
 
 private:
  template <bool regionIter>
  [[nodiscard]] 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 std::array<double, 3> &boxMinWithSafetyMargin,
      const std::array<double, 3> &boxMaxWithSafetyMargin, size_t endCellIndex) const {
    // Finding the indices for the next particle
    const size_t stride = (iteratorBehavior & IteratorBehavior::forceSequential) ? 1 : autopas_get_num_threads();
 
    // helper function to determine if the cell can even contain particles of interest to the iterator
    auto towerIsRelevant = [&]() -> bool {
      // special case: Towers are not yet built, then the tower acting as buffer is always relevant.
      if (_towerBlock.size() == 1) {
        return true;
      }
      bool isRelevant = true;
      if constexpr (regionIter) {
        // is the cell in the region?
        const auto [towerLowCorner, towerHighCorner] = _towerBlock.getTowerBoundingBox(cellIndex);
        isRelevant =
            utils::boxesOverlap(towerLowCorner, towerHighCorner, boxMinWithSafetyMargin, boxMaxWithSafetyMargin);
      }
      return isRelevant;
    };
 
    do {
      // advance to the next particle
      ++particleIndex;
      // If this breaches the end of a cell, find the next non-empty cell and reset particleIndex.
 
      // If cell has wrong type, or there are no more particles in this cell jump to the next
      while (not towerIsRelevant() or particleIndex >= _towerBlock[cellIndex].getNumActualParticles()) {
        cellIndex += stride;
        particleIndex = 0;
 
        // If we notice that there is nothing else to look at set invalid values, so we get a nullptr next time and
        // break.
        if (cellIndex > endCellIndex) {
          return {std::numeric_limits<size_t>::max(), particleIndex};
        }
      }
    } while (not containerIteratorUtils::particleFulfillsIteratorRequirements<regionIter>(
        _towerBlock[cellIndex][particleIndex], iteratorBehavior, boxMin, boxMax));
 
    // the indices returned at this point should always be valid
    return {cellIndex, particleIndex};
  }
 
  internal::ClusterTowerBlock2D<Particle_T> _towerBlock;
 
  autopas::LoadEstimatorOption _loadEstimator;
 
  size_t _clusterSize;
 
  size_t _numClusters{0};
 
  mutable std::vector<std::vector<Particle_T>> _particlesToAdd;
 
  [[nodiscard]] bool particlesToAddEmpty(int bufferID = -1) const {
    if (bufferID == -1) {
      for (auto &threadBuffer : _particlesToAdd) {
        if (not threadBuffer.empty()) {
          return false;
        }
      }
      return true;
    } else {
      return _particlesToAdd[bufferID].empty();
    }
  }
 
  template <class VecVec>
  void appendBuffersHelper(utils::optRef<VecVec> additionalVectors, VecVec &outVec) const {
    if (additionalVectors.has_value()) {
      outVec.reserve(_particlesToAdd.size() + additionalVectors->get().size());
      outVec.insert(outVec.end(), additionalVectors->get().begin(), additionalVectors->get().end());
    } else {
      outVec.reserve(_particlesToAdd.size());
    }
    for (auto &vec : _particlesToAdd) {
      outVec.push_back(&vec);
    }
  }
 
  std::vector<ClusterRange> _clusterThreadPartition;
 
  double _cutoff{};
 
  enum class ValidityState : unsigned char {
    invalid = 0,                 // nothing is valid.
    cellsValidListsInvalid = 1,  // only the cell structure is valid, but the lists are not.
    cellsAndListsValid = 2       // the cells and lists are valid
  };
 
  std::atomic<ValidityState> _isValid{ValidityState::invalid};
 
  std::unique_ptr<internal::VerletClusterListsRebuilder<Particle_T>> _builder;
 
  typename internal::VerletClusterListsRebuilder<Particle_T>::NeighborListsBuffer_T _neighborLists{};
};
 
}  // namespace autopas