InteractionListGeneratorFunctor.h

Go to the documentation of this file.

 
#pragma once
 
#include <atomic>
 
#include "autopas/baseFunctors/PairwiseFunctor.h"
#include "autopas/particles/OwnershipState.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/SoA.h"
namespace autopas {
 
template <class Particle_T, bool isInternal = false>
class InteractionListGeneratorFunctor
    : public PairwiseFunctor<Particle_T, InteractionListGeneratorFunctor<Particle_T, isInternal>> {
 public:
  using SoAArraysType = Particle_T::SoAArraysType;
 
  using NeighborListAoSType = std::unordered_map<Particle_T *, std::vector<Particle_T *>>;
 
  InteractionListGeneratorFunctor(NeighborListAoSType &neighborListsAoS, double interactionLength,
                                  bool gatherNewton3Lists)
      : PairwiseFunctor<Particle_T, InteractionListGeneratorFunctor>(interactionLength),
        _neighborListsAoS(neighborListsAoS),
        _interactionLengthSquared(interactionLength * interactionLength),
        _gatherNewton3Lists(gatherNewton3Lists) {}
 
  std::string getName() override { return "InteractionListGeneratorFunctor"; }
 
  template <class ParticleIterator_T>
  void initializeNeighborList(ParticleIterator_T particlesBegin) {
    _neighborListsAoS.clear();
    for (auto iter = particlesBegin; iter.isValid(); ++iter) {
      _neighborListsAoS[&(*iter)];
    }
  }
 
  bool isRelevantForTuning() override { return not isInternal; }
 
  bool allowsNewton3() override { return true; }
 
  bool allowsNonNewton3() override { return not _gatherNewton3Lists; }
 
  void AoSFunctor(Particle_T &i, Particle_T &j, bool newton3) override {
    using namespace autopas::utils::ArrayMath::literals;
 
    if (_gatherNewton3Lists and not newton3) [[unlikely]] {
      utils::ExceptionHandler::exception(
          "InteractionListGeneratorFunctor should not be used with newton3=false and gatherNewton3Lists=true.");
    }
 
    if (i.isDummy() or j.isDummy()) {
      return;
    }
    auto dist = i.getR() - j.getR();
 
    double distsquare = utils::ArrayMath::dot(dist, dist);
    if (distsquare < _interactionLengthSquared) {
      // Assuming this functor is used like any other functor, this is thread safe: _neighborListsAoS is an
      // unordered_map, meaning we can push_back to particle i's list with the same thread safety as for writing to its
      // force buffer in e.g. a LJ functor.
 
      // This is only thread-safe if all keys are inserted prior to applying this functor, as a rehash might lead to
      // dangling references, but at() protects against this: if the key exists, we can push_back safely; if a key
      // doesn't, an exception is thrown anyway.
 
      // - If newton3=false & gatherNewton3Lists=false, we only need to add the i->j interaction to i's list as the j->i
      // interaction is handled in another call.
      // - If newton3=true & gatherNewton3Lists=false, we need to add the i->j interaction to j's list and the j->i
      // interaction to j's list.
      // - If newton3=true & gatherNewton3Lists=true, we only need to add the i->j interaction to i's list as we don't
      // want the j->i interaction in the lists. (Could also be the other way around)
 
      _neighborListsAoS.at(&i).push_back(&j);
      if (newton3 and not _gatherNewton3Lists) {
        _neighborListsAoS.at(&j).push_back(&i);
      }
    }
  }
 
  void SoAFunctorSingle(SoAView<SoAArraysType> soa, 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 auto *const __restrict ownedStatePtr = soa.template begin<Particle_T::AttributeNames::ownershipState>();
 
    size_t numPart = soa.size();
    for (size_t i = 0; i < numPart; ++i) {
      if (ownedStatePtr[i] == OwnershipState::dummy) continue;
      auto &iList = _neighborListsAoS.at(ptrptr[i]);
      for (size_t j = i + 1; j < numPart; ++j) {
        if (ownedStatePtr[j] == OwnershipState::dummy) continue;
        const double drx = xptr[i] - xptr[j];
        const double dry = yptr[i] - yptr[j];
        const double drz = zptr[i] - zptr[j];
 
        const double drx2 = drx * drx;
        const double dry2 = dry * dry;
        const double drz2 = drz * drz;
 
        const double dr2 = drx2 + dry2 + drz2;
 
        if (dr2 < _interactionLengthSquared) {
          // This SoAFunctorSingle implementation always used newton3 in practice, regardless of the option.
          // This is thread safe (see comment in AoS functor)
          iList.push_back(ptrptr[j]);
          if (not _gatherNewton3Lists) {
            _neighborListsAoS.at(ptrptr[j]).push_back(ptrptr[i]);
          }
        }
      }
    }
  }
 
  void SoAFunctorPair(SoAView<SoAArraysType> soa1, SoAView<SoAArraysType> soa2, bool newton3) override {
    if (_gatherNewton3Lists and not newton3) [[unlikely]] {
      utils::ExceptionHandler::exception(
          "InteractionListGeneratorFunctor should not be used with newton3=false and gatherNewton3Lists=true.");
    }
 
    if (soa1.size() == 0 or 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>();
    const auto *const __restrict ownedState1Ptr = soa1.template begin<Particle_T::AttributeNames::ownershipState>();
 
    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>();
    const auto *const __restrict ownedState2Ptr = soa2.template begin<Particle_T::AttributeNames::ownershipState>();
 
    const size_t numPart1 = soa1.size();
    for (size_t i = 0; i < numPart1; ++i) {
      auto &iList = _neighborListsAoS.at(ptr1ptr[i]);
      if (ownedState1Ptr[i] == OwnershipState::dummy) continue;
      const size_t numPart2 = soa2.size();
 
      for (size_t j = 0; j < numPart2; ++j) {
        if (ownedState2Ptr[j] == OwnershipState::dummy) continue;
        const double drx = x1ptr[i] - x2ptr[j];
        const double dry = y1ptr[i] - y2ptr[j];
        const double drz = z1ptr[i] - z2ptr[j];
 
        const double drx2 = drx * drx;
        const double dry2 = dry * dry;
        const double drz2 = drz * drz;
 
        const double dr2 = drx2 + dry2 + drz2;
 
        if (dr2 < _interactionLengthSquared) {
          // This is thread safe (see comment in AoS functor)
          iList.push_back(ptr2ptr[j]);
          if (newton3 and not _gatherNewton3Lists) {
            _neighborListsAoS.at(ptr2ptr[j]).push_back(ptr1ptr[i]);
          }
        }
      }
    }
  }
 
  void SoAFunctorVerlet(SoAView<SoAArraysType> soa, const size_t indexFirst,
                        const std::vector<size_t, AlignedAllocator<size_t>> &verletList, bool newton3) override {
    if (_gatherNewton3Lists and not newton3) [[unlikely]] {
      utils::ExceptionHandler::exception(
          "InteractionListGeneratorFunctor should not be used with newton3=false and gatherNewton3Lists=true.");
    }
 
    if (soa.size() == 0 or verletList.empty()) 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 auto *const __restrict ownedStatePtr = soa.template begin<Particle_T::AttributeNames::ownershipState>();
 
    if (ownedStatePtr[indexFirst] == OwnershipState::dummy) return;
 
    auto &firstList = _neighborListsAoS.at(ptrptr[indexFirst]);
    const double xFirst = xptr[indexFirst];
    const double yFirst = yptr[indexFirst];
    const double zFirst = zptr[indexFirst];
 
    for (const size_t j : verletList) {
      if (ownedStatePtr[j] == OwnershipState::dummy) continue;
      const double drx = xFirst - xptr[j];
      const double dry = yFirst - yptr[j];
      const double drz = zFirst - zptr[j];
 
      const double drx2 = drx * drx;
      const double dry2 = dry * dry;
      const double drz2 = drz * drz;
 
      const double dr2 = drx2 + dry2 + drz2;
 
      if (dr2 < _interactionLengthSquared) {
        // This is thread safe (see comment in AoS functor)
        firstList.push_back(ptrptr[j]);
        if (newton3 and not _gatherNewton3Lists) {
          _neighborListsAoS.at(ptrptr[j]).push_back(ptrptr[indexFirst]);
        }
      }
    }
  }
 
  constexpr static std::array<typename Particle_T::AttributeNames, 5> getNeededAttr() {
    return std::array<typename Particle_T::AttributeNames, 5>{
        Particle_T::AttributeNames::ptr, Particle_T::AttributeNames::posX, Particle_T::AttributeNames::posY,
        Particle_T::AttributeNames::posZ, Particle_T::AttributeNames::ownershipState};
  }
 
  constexpr static std::array<typename Particle_T::AttributeNames, 5> getNeededAttr(std::false_type) {
    return getNeededAttr();
  }
 
  constexpr static std::array<typename Particle_T::AttributeNames, 0> getComputedAttr() {
    return std::array<typename Particle_T::AttributeNames, 0>{/*Nothing*/};
  }
 
 private:
  NeighborListAoSType &_neighborListsAoS;
  double _interactionLengthSquared;
 
  bool _gatherNewton3Lists{false};
};
 
}  // namespace autopas