VerletListHelpers.h

Go to the documentation of this file.

 
#pragma once
 
#include <atomic>
#include <unordered_map>
#include <vector>
 
#include "autopas/baseFunctors/PairwiseFunctor.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/SoA.h"
namespace autopas {
 
template <class Particle_T>
class VerletListHelpers {
 public:
  struct NeighborListCRS {
    std::vector<size_t> offsets;                            
    std::vector<size_t, AlignedAllocator<size_t>> indices;  
 
    [[nodiscard]] size_t size() const { return offsets.empty() ? 0u : offsets.size() - 1u; }
    [[nodiscard]] size_t count(const size_t i) const { return offsets[i + 1] - offsets[i]; }
    [[nodiscard]] const size_t *begin(const size_t i) const { return indices.data() + offsets[i]; }
    [[nodiscard]] size_t *begin(const size_t i) { return indices.data() + offsets[i]; }
    [[nodiscard]] std::span<const size_t> getNeighbors(const size_t i) const { return {begin(i), count(i)}; }
  };
 
  class CRSNeighborListPolicy {
   public:
    CRSNeighborListPolicy(std::vector<std::vector<size_t>> &neighborLists,
                          std::unordered_map<const Particle_T *, size_t> &particleToIndex)
        : _neighborLists(neighborLists), _particleToIndex(particleToIndex) {}
 
    void add(Particle_T *i, Particle_T *j) { _neighborLists[_particleToIndex[i]].push_back(_particleToIndex[j]); }
 
   private:
    std::vector<std::vector<size_t>> &_neighborLists;
    std::unordered_map<const Particle_T *, size_t> &_particleToIndex;
  };
 
  class VerletListCounterFunctor : public PairwiseFunctor<Particle_T, VerletListCounterFunctor> {
   public:
    using SoAArraysType = typename Particle_T::SoAArraysType;
 
    struct alignas(64) PaddedAtomic {
      std::atomic<size_t> value{0};
    };
 
    VerletListCounterFunctor(std::vector<PaddedAtomic> &counts,
                             const std::unordered_map<const Particle_T *, size_t> &particleToIndex,
                             double interactionLength)
        : PairwiseFunctor<Particle_T, VerletListCounterFunctor>(interactionLength),
          _counts(counts),
          _particleToIndex(particleToIndex),
          _interactionLengthSquared(interactionLength * interactionLength) {}
 
    std::string getName() override { return "VerletListCounterFunctor"; }
    bool isRelevantForTuning() override { return false; }
    bool allowsNewton3() override { return true; }
    bool allowsNonNewton3() override { return true; }
 
    void AoSFunctor(Particle_T &i, Particle_T &j, bool /*newton3*/) override {
      using namespace autopas::utils::ArrayMath::literals;
      if (i.isDummy() or j.isDummy()) return;
      const auto displacement = i.getR() - j.getR();
      if (utils::ArrayMath::dot(displacement, displacement) < _interactionLengthSquared) {
        _counts[_particleToIndex.at(&i)].value.fetch_add(1, std::memory_order_relaxed);
        // newton3 ignored: AoSFunctor(j, i) is also called for newton3=false.
      }
    }
 
    void SoAFunctorSingle(SoAView<SoAArraysType> soa, const bool newton3) override {
      if (soa.size() == 0) return;
      auto **const __restrict ptrptr = soa.template begin<Particle_T::AttributeNames::ptr>();
      const double *const __restrict xptr = soa.template begin<Particle_T::AttributeNames::posX>();
      const double *const __restrict yptr = soa.template begin<Particle_T::AttributeNames::posY>();
      const double *const __restrict zptr = soa.template begin<Particle_T::AttributeNames::posZ>();
      const size_t n = soa.size();
      for (size_t i = 0; i < n; ++i) {
        const size_t iIdx = _particleToIndex.at(ptrptr[i]);
        size_t localCount = 0;
        for (size_t j = i + 1; j < n; ++j) {
          const double dx = xptr[i] - xptr[j];
          const double dy = yptr[i] - yptr[j];
          const double dz = zptr[i] - zptr[j];
          if (dx * dx + dy * dy + dz * dz < _interactionLengthSquared) {
            ++localCount;
            if (not newton3) {
              _counts[_particleToIndex.at(ptrptr[j])].value.fetch_add(1, std::memory_order_relaxed);
            }
          }
        }
        _counts[iIdx].value.fetch_add(localCount, std::memory_order_relaxed);
      }
    }
 
    void SoAFunctorPair(SoAView<SoAArraysType> soa1, SoAView<SoAArraysType> soa2, bool /*newton3*/) override {
      if (soa1.size() == 0 || soa2.size() == 0) return;
      auto **const __restrict ptr1ptr = soa1.template begin<Particle_T::AttributeNames::ptr>();
      const double *const __restrict x1ptr = soa1.template begin<Particle_T::AttributeNames::posX>();
      const double *const __restrict y1ptr = soa1.template begin<Particle_T::AttributeNames::posY>();
      const double *const __restrict z1ptr = soa1.template begin<Particle_T::AttributeNames::posZ>();
      auto **const __restrict ptr2ptr = soa2.template begin<Particle_T::AttributeNames::ptr>();
      const double *const __restrict x2ptr = soa2.template begin<Particle_T::AttributeNames::posX>();
      const double *const __restrict y2ptr = soa2.template begin<Particle_T::AttributeNames::posY>();
      const double *const __restrict z2ptr = soa2.template begin<Particle_T::AttributeNames::posZ>();
      // newton3 ignored: SoAFunctorPair(soa2,soa1) is also called for newton3=false.
      const size_t n1 = soa1.size(), n2 = soa2.size();
      for (size_t i = 0; i < n1; ++i) {
        const size_t iIdx = _particleToIndex.at(ptr1ptr[i]);
        size_t localCount = 0;
        for (size_t j = 0; j < n2; ++j) {
          const double dx = x1ptr[i] - x2ptr[j];
          const double dy = y1ptr[i] - y2ptr[j];
          const double dz = z1ptr[i] - z2ptr[j];
          if (dx * dx + dy * dy + dz * dz < _interactionLengthSquared) {
            ++localCount;
          }
        }
        _counts[iIdx].value.fetch_add(localCount, std::memory_order_relaxed);
      }
    }
 
    constexpr static std::array<typename Particle_T::AttributeNames, 4> getNeededAttr() {
      return {Particle_T::AttributeNames::ptr, Particle_T::AttributeNames::posX, Particle_T::AttributeNames::posY,
              Particle_T::AttributeNames::posZ};
    }
 
    constexpr static std::array<typename Particle_T::AttributeNames, 0> getComputedAttr() { return {}; }
 
   private:
    std::vector<PaddedAtomic> &_counts;
    const std::unordered_map<const Particle_T *, size_t> &_particleToIndex;
    double _interactionLengthSquared;
  };
 
  class VerletListFillerFunctor : public PairwiseFunctor<Particle_T, VerletListFillerFunctor> {
   public:
    using SoAArraysType = typename Particle_T::SoAArraysType;
    using PaddedAtomic = typename VerletListCounterFunctor::PaddedAtomic;
 
    VerletListFillerFunctor(NeighborListCRS &neighborList, std::vector<PaddedAtomic> &fillPos,
                            const std::unordered_map<const Particle_T *, size_t> &particleToIndex,
                            double interactionLength)
        : PairwiseFunctor<Particle_T, VerletListFillerFunctor>(interactionLength),
          _neighborList(neighborList),
          _fillPos(fillPos),
          _particleToIndex(particleToIndex),
          _interactionLengthSquared(interactionLength * interactionLength) {}
 
    std::string getName() override { return "VerletListFillerFunctor"; }
    bool isRelevantForTuning() override { return false; }
    bool allowsNewton3() override { return true; }
    bool allowsNonNewton3() override { return true; }
 
    void AoSFunctor(Particle_T &i, Particle_T &j, bool /*newton3*/) override {
      using namespace autopas::utils::ArrayMath::literals;
      if (i.isDummy() or j.isDummy()) return;
      const auto displacement = i.getR() - j.getR();
      if (utils::ArrayMath::dot(displacement, displacement) < _interactionLengthSquared) {
        const size_t iIdx = _particleToIndex.at(&i);
        const size_t jIdx = _particleToIndex.at(&j);
        _neighborList.indices[_fillPos[iIdx].value.fetch_add(1, std::memory_order_relaxed)] = jIdx;
        // newton3 ignored: AoSFunctor(j, i) is also called for newton3=false.
      }
    }
 
    void SoAFunctorSingle(SoAView<SoAArraysType> soa, const bool newton3) override {
      if (soa.size() == 0) return;
      auto **const __restrict ptrptr = soa.template begin<Particle_T::AttributeNames::ptr>();
      const double *const __restrict xptr = soa.template begin<Particle_T::AttributeNames::posX>();
      const double *const __restrict yptr = soa.template begin<Particle_T::AttributeNames::posY>();
      const double *const __restrict zptr = soa.template begin<Particle_T::AttributeNames::posZ>();
      const size_t n = soa.size();
      for (size_t i = 0; i < n; ++i) {
        const size_t iIdx = _particleToIndex.at(ptrptr[i]);
        for (size_t j = i + 1; j < n; ++j) {
          const double dx = xptr[i] - xptr[j];
          const double dy = yptr[i] - yptr[j];
          const double dz = zptr[i] - zptr[j];
          if (dx * dx + dy * dy + dz * dz < _interactionLengthSquared) {
            const size_t jIdx = _particleToIndex.at(ptrptr[j]);
            _neighborList.indices[_fillPos[iIdx].value.fetch_add(1, std::memory_order_relaxed)] = jIdx;
            if (not newton3) {
              _neighborList.indices[_fillPos[jIdx].value.fetch_add(1, std::memory_order_relaxed)] = iIdx;
            }
          }
        }
      }
    }
 
    void SoAFunctorPair(SoAView<SoAArraysType> soa1, SoAView<SoAArraysType> soa2, bool /*newton3*/) override {
      if (soa1.size() == 0 || soa2.size() == 0) return;
      auto **const __restrict ptr1ptr = soa1.template begin<Particle_T::AttributeNames::ptr>();
      const double *const __restrict x1ptr = soa1.template begin<Particle_T::AttributeNames::posX>();
      const double *const __restrict y1ptr = soa1.template begin<Particle_T::AttributeNames::posY>();
      const double *const __restrict z1ptr = soa1.template begin<Particle_T::AttributeNames::posZ>();
      auto **const __restrict ptr2ptr = soa2.template begin<Particle_T::AttributeNames::ptr>();
      const double *const __restrict x2ptr = soa2.template begin<Particle_T::AttributeNames::posX>();
      const double *const __restrict y2ptr = soa2.template begin<Particle_T::AttributeNames::posY>();
      const double *const __restrict z2ptr = soa2.template begin<Particle_T::AttributeNames::posZ>();
      // newton3 ignored: SoAFunctorPair(soa2,soa1) is also called for newton3=false.
      const size_t n1 = soa1.size(), n2 = soa2.size();
      for (size_t i = 0; i < n1; ++i) {
        const size_t iIdx = _particleToIndex.at(ptr1ptr[i]);
        for (size_t j = 0; j < n2; ++j) {
          const double dx = x1ptr[i] - x2ptr[j];
          const double dy = y1ptr[i] - y2ptr[j];
          const double dz = z1ptr[i] - z2ptr[j];
          if (dx * dx + dy * dy + dz * dz < _interactionLengthSquared) {
            _neighborList.indices[_fillPos[iIdx].value.fetch_add(1, std::memory_order_relaxed)] =
                _particleToIndex.at(ptr2ptr[j]);
          }
        }
      }
    }
 
    constexpr static std::array<typename Particle_T::AttributeNames, 4> getNeededAttr() {
      return {Particle_T::AttributeNames::ptr, Particle_T::AttributeNames::posX, Particle_T::AttributeNames::posY,
              Particle_T::AttributeNames::posZ};
    }
 
    constexpr static std::array<typename Particle_T::AttributeNames, 0> getComputedAttr() { return {}; }
 
   private:
    NeighborListCRS &_neighborList;
    std::vector<PaddedAtomic> &_fillPos;
    const std::unordered_map<const Particle_T *, size_t> &_particleToIndex;
    double _interactionLengthSquared;
  };
};  // class VerletListHelpers
}  // namespace autopas