VerletLists.h

Go to the documentation of this file.

 
#pragma once
 
#include "VerletListHelpers.h"
#include "autopas/baseFunctors/InteractionListGeneratorFunctor.h"
#include "autopas/containers/CellBasedParticleContainer.h"
#include "autopas/containers/linkedCells/LinkedCells.h"
#include "autopas/containers/linkedCells/traversals/LCC08Traversal.h"
#include "autopas/containers/verletListsCellBased/VerletListsLinkedBase.h"
#include "autopas/containers/verletListsCellBased/verletLists/traversals/VLTraversalInterface.h"
#include "autopas/options/DataLayoutOption.h"
#include "autopas/utils/WrapOpenMP.h"
 
namespace autopas {
 
template <class Particle_T>
class VerletLists : public VerletListsLinkedBase<Particle_T> {
  using ParticleType = Particle_T;
  using ParticleCellType = FullParticleCell<Particle_T>;
 
 public:
  enum BuildVerletListType {
    VerletAoS,
    VerletSoA,
  };
 
  VerletLists(const std::array<double, 3> &boxMin, const std::array<double, 3> &boxMax, const double cutoff,
              const double skin, const BuildVerletListType buildVerletListType = BuildVerletListType::VerletSoA,
              const double cellSizeFactor = 1.0)
      : VerletListsLinkedBase<Particle_T>(boxMin, boxMax, cutoff, skin, cellSizeFactor),
        _buildVerletListType(buildVerletListType) {}
 
  [[nodiscard]] ContainerOption getContainerType() const override { return ContainerOption::verletLists; }
 
  void computeInteractions(TraversalInterface *traversal) override {
    // Check if traversal is allowed for this container and give it the data it needs.
    auto *verletTraversalInterface = dynamic_cast<VLTraversalInterface<ParticleCellType> *>(traversal);
    if (verletTraversalInterface) {
      verletTraversalInterface->setCellsAndNeighborLists(this->_linkedCells.getCells(), _neighborList, _particleToIndex,
                                                         _indexToParticle);
    } else {
      utils::ExceptionHandler::exception("trying to use a traversal of wrong type in VerletLists::computeInteractions");
    }
 
    traversal->initTraversal();
    traversal->traverseParticles();
    traversal->endTraversal();
  }
 
  const VerletListHelpers<Particle_T>::NeighborListCRS &getNeighborList() const { return _neighborList; }
 
  const std::unordered_map<const Particle_T *, size_t> &getParticleIndex() const { return _particleToIndex; }
 
  void rebuildNeighborLists(TraversalInterface *traversal) override {
    this->_verletBuiltNewton3 = traversal->getUseNewton3();
    updateNeighborLists(traversal->getUseNewton3());
    // the neighbor list is now valid
    this->_neighborListIsValid.store(true, std::memory_order_relaxed);
  }
 
 private:
  size_t buildParticleIndex() {
    _particleToIndex.clear();
    _particleToIndex.reserve(this->_linkedCells.size());
    _indexToParticle.clear();
    _indexToParticle.reserve(this->_linkedCells.size());
    size_t idx = 0;
    for (auto iter = this->begin(IteratorBehavior::ownedOrHaloOrDummy); iter.isValid(); ++iter, ++idx) {
      _particleToIndex[&(*iter)] = idx;
      _indexToParticle.push_back(&*iter);
    }
    return idx;
  }
 
  virtual void updateNeighborLists(bool useNewton3) {
    const size_t N = buildParticleIndex();
    const double interactionLength = this->getCutoff() + this->getVerletSkin();
 
    DataLayoutOption dataLayout;
    if (_buildVerletListType == BuildVerletListType::VerletAoS) {
      dataLayout = DataLayoutOption::aos;
    } else if (_buildVerletListType == BuildVerletListType::VerletSoA) {
      dataLayout = DataLayoutOption::soa;
    } else {
      utils::ExceptionHandler::exception("VerletLists::updateNeighborLists(): unsupported BuildVerletListType: {}",
                                         static_cast<int>(_buildVerletListType));
    }
 
    if (autopas_get_max_threads() == 1) {
      updateNeighborListsSingleThread(N, interactionLength, dataLayout, useNewton3);
    } else {
      updateNeighborListsMultiThread(N, interactionLength, dataLayout, useNewton3);
    }
  }
 
  void updateNeighborListsSingleThread(size_t N, double interactionLength, DataLayoutOption dataLayout,
                                       bool useNewton3) {
    std::vector<std::vector<size_t>> tempLists(N);
 
    typename VerletListHelpers<Particle_T>::CRSNeighborListPolicy policy(tempLists, _particleToIndex);
 
    InteractionListGeneratorFunctor<Particle_T, typename VerletListHelpers<Particle_T>::CRSNeighborListPolicy> f(
        policy, interactionLength, useNewton3);
    auto traversal = LCC08Traversal<
        ParticleCellType,
        InteractionListGeneratorFunctor<Particle_T, typename VerletListHelpers<Particle_T>::CRSNeighborListPolicy>>(
        this->_linkedCells.getCellBlock().getCellsPerDimensionWithHalo(), f, this->getInteractionLength(),
        this->_linkedCells.getCellBlock().getCellLength(), dataLayout, useNewton3);
    this->_linkedCells.computeInteractions(&traversal);
 
    // Prefix-sum + copy:
    _neighborList.offsets.resize(N + 1);
    _neighborList.offsets[0] = 0;
    for (size_t i = 0; i < N; ++i) {
      _neighborList.offsets[i + 1] = _neighborList.offsets[i] + tempLists[i].size();
    }
    const size_t totalNeighbors = _neighborList.offsets[N];
    _neighborList.indices.resize(totalNeighbors);
    for (size_t i = 0; i < N; ++i) {
      std::copy(tempLists[i].begin(), tempLists[i].end(),
                _neighborList.indices.begin() + static_cast<std::ptrdiff_t>(_neighborList.offsets[i]));
    }
 
    AutoPasLog(DEBUG, "VerletLists::updateNeighborLists (1T): {} particles, {} neighbors, avg {:.2f}", N,
               totalNeighbors, N > 0 ? static_cast<double>(totalNeighbors) / static_cast<double>(N) : 0.0);
  }
 
  void updateNeighborListsMultiThread(size_t N, double interactionLength, DataLayoutOption dataLayout,
                                      bool useNewton3) {
    using PaddedAtomic = VerletListHelpers<Particle_T>::VerletListCounterFunctor::PaddedAtomic;
 
    // Pass 1: Count neighbors per particle
    std::vector<PaddedAtomic> counts(N);
    {
      typename VerletListHelpers<Particle_T>::VerletListCounterFunctor counter(counts, _particleToIndex,
                                                                               interactionLength);
      auto traversal =
          LCC08Traversal<ParticleCellType, typename VerletListHelpers<Particle_T>::VerletListCounterFunctor>(
              this->_linkedCells.getCellBlock().getCellsPerDimensionWithHalo(), counter, this->getInteractionLength(),
              this->_linkedCells.getCellBlock().getCellLength(), dataLayout, useNewton3);
      this->_linkedCells.computeInteractions(&traversal);
    }
 
    // Prefix sum: counts -> CRS offsets, allocate flat indices
    _neighborList.offsets.resize(N + 1);
    _neighborList.offsets[0] = 0;
    for (size_t i = 0; i < N; ++i) {
      _neighborList.offsets[i + 1] = _neighborList.offsets[i] + counts[i].value.load(std::memory_order_relaxed);
    }
    const size_t totalNeighbors = _neighborList.offsets[N];
    _neighborList.indices.resize(totalNeighbors);
 
    AutoPasLog(DEBUG, "VerletLists::updateNeighborLists (MT): {} particles, {} neighbors, avg {:.2f}", N,
               totalNeighbors, N > 0 ? static_cast<double>(totalNeighbors) / static_cast<double>(N) : 0.0);
 
    // Pass 2: Fill CRS indices directly
    // Reuse the PaddedAtomic array as fill cursors, seeding each with offsets[i].
    for (size_t i = 0; i < N; ++i) {
      counts[i].value.store(_neighborList.offsets[i], std::memory_order_relaxed);
    }
    {
      typename VerletListHelpers<Particle_T>::VerletListFillerFunctor filler(_neighborList, counts, _particleToIndex,
                                                                             interactionLength);
      auto traversal =
          LCC08Traversal<ParticleCellType, typename VerletListHelpers<Particle_T>::VerletListFillerFunctor>(
              this->_linkedCells.getCellBlock().getCellsPerDimensionWithHalo(), filler, this->getInteractionLength(),
              this->_linkedCells.getCellBlock().getCellLength(), dataLayout, useNewton3);
      this->_linkedCells.computeInteractions(&traversal);
    }
  }
 
  std::unordered_map<const Particle_T *, size_t> _particleToIndex;
 
  std::vector<ParticleType *> _indexToParticle;
 
  VerletListHelpers<Particle_T>::NeighborListCRS _neighborList;
 
  BuildVerletListType _buildVerletListType;
};
 
}  // namespace autopas