VLCAllCellsNeighborList.h

Go to the documentation of this file.

 
#pragma once
 
#include "VLCAllCellsGeneratorFunctor.h"
#include "VLCNeighborListInterface.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/StaticBoolSelector.h"
 
namespace autopas {
 
template <class ParticleCell>
class TraversalSelector;
 
template <class Particle_T, class NeighborList>
class VLCTraversalInterface;
 
template <class Particle_T>
class VLCAllCellsNeighborList : public VLCNeighborListInterface<Particle_T> {
 public:
  using ListType = typename VerletListsCellsHelpers::AllCellsNeighborListsType<Particle_T>;
 
  using SoAPairOfParticleAndList = std::pair<size_t, std::vector<size_t, autopas::AlignedAllocator<size_t>>>;
 
  using SoAListType = typename std::vector<std::vector<SoAPairOfParticleAndList>>;
 
  [[nodiscard]] ContainerOption getContainerType() const override { return ContainerOption::verletListsCells; }
 
  void buildAoSNeighborList(LinkedCells<Particle_T> &linkedCells, bool useNewton3, double cutoff, double skin,
                            double interactionLength, const TraversalOption vlcTraversalOpt,
                            typename VerletListsCellsHelpers::VLCBuildType buildType) override {
    using namespace utils::ArrayMath::literals;
    this->_internalLinkedCells = &linkedCells;
    auto &cells = linkedCells.getCells();
    const auto cellsPerDim = this->_internalLinkedCells->getCellBlock().getCellsPerDimensionWithHalo();
 
    const auto boxSizeWithHalo = this->_internalLinkedCells->getBoxMax() - this->_internalLinkedCells->getBoxMin() +
                                 std::array<double, 3>{interactionLength, interactionLength, interactionLength} * 2.;
    const auto listLengthEstimate = VerletListsCellsHelpers::estimateListLength(
        this->_internalLinkedCells->getNumberOfParticles(IteratorBehavior::ownedOrHalo), boxSizeWithHalo,
        interactionLength, 1.3);
 
    const auto offsetsC08 = VerletListsCellsHelpers::buildC08BaseStep(utils::ArrayUtils::static_cast_copy_array<int>(
        this->_internalLinkedCells->getCellBlock().getCellsPerDimensionWithHalo()));
 
    // Helper function to estimate the number of neighbor lists for one base step
    // TODO: This is a generous and rough estimate and can probably be improved!
    const auto estimateNumLists = [&](size_t baseCellIndex) {
      // If the cell is near the end of any dimension we can reduce the estimate to the particles in the base cell,
      // because all other cells' interactions would be from halos.
      const auto cellIndex3D = utils::ThreeDimensionalMapping::oneToThreeD(baseCellIndex, cellsPerDim);
      const auto numTouchedBoundaries = [&]() {
        int acc = 0;
        for (size_t i = 0; i < 3; ++i) {
          // false == 0 ; true == 1
          acc += static_cast<int>(cellIndex3D[i] == (cellsPerDim[i] - 1));
        }
        return acc;
      }();
      if (numTouchedBoundaries > 0) {
        if (useNewton3) {
          return cells[baseCellIndex].size();
        } else {
          // In this case we have to accommodate the lists for the reverse interactions from all non-halo neighbors.
          // 1x for lists from the cell itself and once more per untouched border.
          return cells[baseCellIndex].size() * (4 - numTouchedBoundaries);
        }
      }
      size_t estimate = 0;
      std::map<int, double> offsetFactors{};
      std::map<int, double> offsetFactorsNoN3{};
      for (const auto [offsetA, offsetB, factor] : offsetsC08) {
        offsetFactors[offsetA] = std::max(offsetFactors[offsetA], factor);
        offsetFactorsNoN3[offsetB] = std::max(offsetFactors[offsetB], factor);
      }
      for (const auto &[offset, factor] : offsetFactors) {
        estimate += cells[baseCellIndex + offset].size() * factor;
      }
      if (not useNewton3) {
        for (const auto &[offset, factor] : offsetFactorsNoN3) {
          estimate += cells[baseCellIndex + offset].size() * factor;
        }
      }
      return estimate;
    };
 
    // Initialize a map of neighbor lists for each cell.
    _aosNeighborList.clear();
    const size_t numCells = cells.size();
    _aosNeighborList.resize(numCells);
    for (size_t cellIndex = 0; cellIndex < numCells; ++cellIndex) {
      const auto estimateForNumberOfLists = [&]() {
        // Usually each cell holds one list per particle, except vlc_c08, which holds all lists of a base step,
        // which involves also lists from other cells' particles.
        if (vlcTraversalOpt == TraversalOption::vlc_c08) {
          return estimateNumLists(cellIndex);
        } else {
          return cells[cellIndex].size();
        }
      }();
      auto &cell = cells[cellIndex];
      _aosNeighborList[cellIndex].reserve(estimateForNumberOfLists);
      size_t particleIndexWithinCell = 0;
      for (auto iter = cell.begin(); iter != cell.end(); ++iter, ++particleIndexWithinCell) {
        Particle_T *particle = &*iter;
        _aosNeighborList[cellIndex].emplace_back(particle, std::vector<Particle_T *>());
        _aosNeighborList[cellIndex].back().second.reserve(listLengthEstimate);
        _particleToCellMap[particle] = std::make_pair(cellIndex, particleIndexWithinCell);
      }
    }
 
    applyBuildFunctor(linkedCells, useNewton3, cutoff, skin, interactionLength, vlcTraversalOpt, buildType);
  }
 
  size_t getNumberOfPartners(const Particle_T *particle) const override {
    const auto &[cellIndex, particleIndexInCell] = _particleToCellMap.at(const_cast<Particle_T *>(particle));
    return _aosNeighborList.at(cellIndex).at(particleIndexInCell).second.size();
  }
 
  typename VerletListsCellsHelpers::AllCellsNeighborListsType<Particle_T> &getAoSNeighborList() {
    return _aosNeighborList;
  }
 
  auto &getSoANeighborList() { return _soaNeighborList; }
 
  void generateSoAFromAoS(LinkedCells<Particle_T> &linkedCells) override {
    _soaNeighborList.clear();
 
    // particle pointer to global index of particle
    std::unordered_map<Particle_T *, size_t> particlePtrToIndex;
    particlePtrToIndex.reserve(linkedCells.size());
    size_t i = 0;
    for (auto iter = linkedCells.begin(IteratorBehavior::ownedOrHaloOrDummy); iter.isValid(); ++iter, ++i) {
      particlePtrToIndex[&(*iter)] = i;
    }
 
    _soaNeighborList.resize(linkedCells.getCells().size());
 
    // iterate over cells and for each create the soa lists from the aos lists
    for (size_t firstCellIndex = 0; firstCellIndex < _aosNeighborList.size(); ++firstCellIndex) {
      const auto &aosLists = _aosNeighborList[firstCellIndex];
      auto &soaLists = _soaNeighborList[firstCellIndex];
      soaLists.reserve(aosLists.size());
 
      // iterate over pairs of particle and neighbor list
      for (const auto &[particlePtr, neighbors] : aosLists) {
        // global index of current particle
        size_t currentParticleGlobalIndex = particlePtrToIndex.at(particlePtr);
 
        // create SoA neighbor list for current particle
        std::vector<size_t, autopas::AlignedAllocator<size_t>> currentSoANeighborList{};
        currentSoANeighborList.reserve(neighbors.size());
 
        // fill the SoA neighbor list with the indices of the particles from the corresponding AoS neighbor list
        for (const auto &neighborOfCurrentParticle : neighbors) {
          currentSoANeighborList.emplace_back(particlePtrToIndex.at(neighborOfCurrentParticle));
        }
 
        // add the newly constructed pair of particle index and SoA neighbor list to cell
        soaLists.emplace_back(currentParticleGlobalIndex, currentSoANeighborList);
      }
    }
  }
 
  void setUpTraversal(TraversalInterface *traversal) override {
    auto vTraversal = dynamic_cast<VLCTraversalInterface<Particle_T, VLCAllCellsNeighborList<Particle_T>> *>(traversal);
 
    if (vTraversal) {
      vTraversal->setVerletList(*this);
    } else {
      autopas::utils::ExceptionHandler::exception(
          "Trying to use a traversal of wrong type in VerletListCells.h. TraversalID: {}",
          traversal->getTraversalType());
    }
  }
 
 private:
  void applyBuildFunctor(LinkedCells<Particle_T> &linkedCells, bool useNewton3, double cutoff, double skin,
                         double interactionLength, const TraversalOption &vlcTraversalOpt,
                         typename VerletListsCellsHelpers::VLCBuildType buildType) override {
    // Generate the build traversal with the traversal selector and apply the build functor with it.
    TraversalSelector<FullParticleCell<Particle_T>> traversalSelector;
    // Argument "cluster size" does not matter here.
    TraversalSelectorInfo traversalSelectorInfo(linkedCells.getCellBlock().getCellsPerDimensionWithHalo(),
                                                interactionLength, linkedCells.getCellBlock().getCellLength(), 0);
 
    const auto dataLayout =
        buildType == VerletListsCellsHelpers::VLCBuildType::aosBuild ? DataLayoutOption::aos : DataLayoutOption::soa;
 
    // LC traversal that will be used to build the List
    const TraversalOption lcBuildTraversalOpt =
        vlcTraversalOpt == TraversalOption::vlc_c08 ? TraversalOption::lc_c08 : TraversalOption::lc_c18;
 
    using namespace utils::ArrayMath::literals;
    const auto boxSizeWithHalo = linkedCells.getBoxMax() - linkedCells.getBoxMin() +
                                 std::array<double, 3>{interactionLength, interactionLength, interactionLength} * 2.;
    const auto listLengthEstimate = VerletListsCellsHelpers::estimateListLength(
        linkedCells.getNumberOfParticles(IteratorBehavior::ownedOrHalo), boxSizeWithHalo, interactionLength, 1.3);
 
    if (vlcTraversalOpt == TraversalOption::vlc_c08) {
      VLCAllCellsGeneratorFunctor<Particle_T, TraversalOption::vlc_c08> f(
          _aosNeighborList, _particleToCellMap, cutoff + skin, listLengthEstimate,
          linkedCells.getCellBlock().getCellsPerDimensionWithHalo());
      f.setCells(&this->_internalLinkedCells->getCells());
      // Build the AoS list using the AoS or SoA functor depending on buildType
      auto buildTraversal = traversalSelector.template generateTraversal<std::remove_reference_t<decltype(f)>>(
          lcBuildTraversalOpt, f, traversalSelectorInfo, dataLayout, useNewton3);
      linkedCells.computeInteractions(buildTraversal.get());
    } else {
      VLCAllCellsGeneratorFunctor<Particle_T, TraversalOption::vlc_c18> f(
          _aosNeighborList, _particleToCellMap, cutoff + skin, listLengthEstimate,
          linkedCells.getCellBlock().getCellsPerDimensionWithHalo());
      // Build the AoS list using the AoS or SoA functor depending on buildType
      auto buildTraversal = traversalSelector.template generateTraversal<std::remove_reference_t<decltype(f)>>(
          lcBuildTraversalOpt, f, traversalSelectorInfo, dataLayout, useNewton3);
      linkedCells.computeInteractions(buildTraversal.get());
    }
  }
 
  typename VerletListsCellsHelpers::AllCellsNeighborListsType<Particle_T> _aosNeighborList{};
 
  std::unordered_map<Particle_T *, std::pair<size_t, size_t>> _particleToCellMap{};
 
  SoAListType _soaNeighborList{};
};
}  // namespace autopas