CellFunctor.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
 
#include "autopas/baseFunctors/PairwiseFunctor.h"
#include "autopas/cells/SortedCellView.h"
#include "autopas/options/DataLayoutOption.h"
#include "autopas/utils/SortedSoAView.h"
#include "autopas/utils/SortingThresholdInfo2B.h"
#include "autopas/utils/SortingThresholdInfoInterface.h"
#include "autopas/utils/SortingThresholdInfoSingle.h"
#include "autopas/utils/WrapOpenMP.h"
 
namespace autopas::internal {
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional = true>
class CellFunctor {
 public:
  explicit CellFunctor(ParticleFunctor_T &f, const double sortingCutoff, DataLayoutOption dataLayout, bool useNewton3)
      : _functor(f), _sortingCutoff(sortingCutoff), _dataLayout(dataLayout), _useNewton3(useNewton3) {
    if constexpr (ParticleFunctor_T::supportsSoASorting) {
      if (dataLayout == DataLayoutOption::soa) {
        _soaThreadData.resize(autopas::autopas_get_max_threads());
      }
    }
  }
 
  void processCell(ParticleCell_T &cell);
 
  void processCellPair(ParticleCell_T &cell1, ParticleCell_T &cell2,
                       const std::array<double, 3> &sortingDirection = {0., 0., 0.});
 
  [[nodiscard]] DataLayoutOption::Value getDataLayout() const { return _dataLayout; }
 
  [[nodiscard]] bool getNewton3() const { return _useNewton3; }
 
  [[nodiscard]] bool getBidirectional() const { return bidirectional; }
 
  void setAoSSortingThresholds(const SortingThresholdInfoInterface &aosSortingThreshold);
 
  void setSoASortingThresholds(const SortingThresholdInfoInterface &soaSortingThreshold);
 
  [[nodiscard]] SoASortingData computeSortingData(const std::vector<std::pair<double, size_t>> &projIdxI,
                                                  const std::vector<std::pair<double, size_t>> &projIdxJ,
                                                  std::vector<size_t> &maxIndexCache,
                                                  std::vector<size_t> &minIndexCache) const;
 
 private:
  [[nodiscard]] bool shouldUseSorting(size_t particleCount, const std::array<double, 3> &sortingDirection,
                                      bool useSoA) const {
    if (sortingDirection[0] != 0.0 or sortingDirection[1] != 0.0 or sortingDirection[2] != 0.0) {
      const auto &thresholds = useSoA ? _soaSortingThresholds : _aosSortingThresholds;
      return particleCount >= thresholds.getThresholdByConfig(_useNewton3, sortingDirection, bidirectional);
    }
    return false;
  }
 
  [[nodiscard]] static SortingThresholdInfo2B resolveThreshold(const SortingThresholdInfoInterface &info) {
    if (const auto *info2B = dynamic_cast<const SortingThresholdInfo2B *>(&info)) {
      return *info2B;
    }
    if (const auto *infoSingle = dynamic_cast<const SortingThresholdInfoSingle *>(&info)) {
      return SortingThresholdInfo2B(infoSingle->threshold);
    }
    throw utils::ExceptionHandler::AutoPasException(
        "CellFunctor was called with wrong SortingThresholdInfo Type. (Supported are single and 2B)");
  }
 
  void processCellAoSImpl(ParticleCell_T &cell);
 
  void processCellPairAoSImpl(ParticleCell_T &cell1, ParticleCell_T &cell2,
                              const std::array<double, 3> &sortingDirection);
 
  void processCellPairSoAImpl(ParticleCell_T &cell1, ParticleCell_T &cell2,
                              const std::array<double, 3> &sortingDirection);
 
  ParticleFunctor_T &_functor;
 
  const double _sortingCutoff;
 
  SortingThresholdInfo2B _aosSortingThresholds{8};
  SortingThresholdInfo2B _soaSortingThresholds{100};
 
  const DataLayoutOption::Value _dataLayout;
 
  const bool _useNewton3;
 
  struct alignas(64) SoAThreadData {
    SoA<typename ParticleCell_T::ParticleType::SoAArraysType> sortedSoa1;
    SoA<typename ParticleCell_T::ParticleType::SoAArraysType> sortedSoa2;
    std::vector<std::pair<double, size_t>> projIdx1;
    std::vector<std::pair<double, size_t>> projIdx2;
    std::vector<size_t> maxIndex;
    std::vector<size_t> minIndex;
  };
 
  std::vector<SoAThreadData> _soaThreadData{};
};
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::setAoSSortingThresholds(
    const SortingThresholdInfoInterface &aosSortingThreshold) {
  _aosSortingThresholds = resolveThreshold(aosSortingThreshold);
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::setSoASortingThresholds(
    const SortingThresholdInfoInterface &soaSortingThreshold) {
  _soaSortingThresholds = resolveThreshold(soaSortingThreshold);
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::processCell(ParticleCell_T &cell) {
  const bool isAoS = _dataLayout == DataLayoutOption::aos;
  const bool isSoA = _dataLayout == DataLayoutOption::soa;
 
  // Return early if the cell is empty.
  if ((isSoA and cell._particleSoABuffer.size() == 0) or (isAoS and cell.isEmpty())) {
    return;
  }
  // Avoid force calculations if the cell contains only halo particles or if the cell is empty (=dummy)
  if (not cell.canHaveOwnedParticles()) {
    return;
  }
 
  if (isAoS) {
    processCellAoSImpl(cell);
  } else if (isSoA) {
    _functor.SoAFunctorSingle(cell._particleSoABuffer, _useNewton3);
  }
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::processCellPair(
    ParticleCell_T &cell1, ParticleCell_T &cell2, const std::array<double, 3> &sortingDirection) {
  const bool isAoS = _dataLayout == DataLayoutOption::aos;
  const bool isSoA = _dataLayout == DataLayoutOption::soa;
 
  // Return early if a cell is empty.
  if ((isSoA and (cell1._particleSoABuffer.size() == 0 or cell2._particleSoABuffer.size() == 0)) or
      (isAoS and (cell1.isEmpty() or cell2.isEmpty()))) {
    return;
  }
 
  if (not cell1.canHaveOwnedParticles()) {
    // Nothing to do if cell1 has no owned particles and we don't write to cell2 particles.
    if constexpr (not bidirectional) {
      if (not _useNewton3) {
        return;
      }
    }
    // Nothing to do if both cells cannot have owned particles.
    if (not cell2.canHaveOwnedParticles()) {
      return;
    }
  }
 
  if (isAoS) {
    processCellPairAoSImpl(cell1, cell2, sortingDirection);
  } else if (isSoA) {
    processCellPairSoAImpl(cell1, cell2, sortingDirection);
  }
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::processCellAoSImpl(ParticleCell_T &cell) {
  // helper function
  const auto interactParticles = [this](auto &p1, auto &p2) {
    this->_functor.AoSFunctor(p1, p2, this->_useNewton3);
 
    if (not this->_useNewton3) {
      this->_functor.AoSFunctor(p2, p1, false);
    }
  };
  // @todo: Which parameter from the Thresholds is used for inner cell sorting is currently arbitrary, think of a good
  // way to handle this
  if (cell.size() >= _aosSortingThresholds.getThresholdByConfig(_useNewton3, {1, 0, 0}, bidirectional)) {
    SortedCellView<ParticleCell_T> cellSorted(cell, utils::ArrayMath::normalize(cell.getCellLength()));
 
    for (auto cellIter1 = cellSorted._particles.begin(); cellIter1 != cellSorted._particles.end(); ++cellIter1) {
      auto &[p1Projection, p1Ptr] = *cellIter1;
      // start inner loop ahead of the outer loop
      for (auto cellIter2 = std::next(cellIter1); cellIter2 != cellSorted._particles.end(); ++cellIter2) {
        auto &[p2Projection, p2Ptr] = *cellIter2;
        if (std::abs(p2Projection - p1Projection) > _sortingCutoff) {
          break;
        }
        interactParticles(*p1Ptr, *p2Ptr);
      }
    }
  } else {
    for (auto p1Ptr = cell.begin(); p1Ptr != cell.end(); ++p1Ptr) {
      auto p2Ptr = p1Ptr;
      ++p2Ptr;
      for (; p2Ptr != cell.end(); ++p2Ptr) {
        interactParticles(*p1Ptr, *p2Ptr);
      }
    }
  }
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::processCellPairAoSImpl(
    ParticleCell_T &cell1, ParticleCell_T &cell2, const std::array<double, 3> &sortingDirection) {
  const auto interactParticles = [this](auto &p1, auto &p2) {
    this->_functor.AoSFunctor(p1, p2, this->_useNewton3);
    if constexpr (bidirectional) {
      if (not this->_useNewton3) {
        this->_functor.AoSFunctor(p2, p1, false);
      }
    }
  };
 
  if (shouldUseSorting(cell1.size() + cell2.size(), sortingDirection, /*useSoA*/ false)) {
    // Use sorted cell views
    SortedCellView<ParticleCell_T> cell1Sorted(cell1, sortingDirection);
    SortedCellView<ParticleCell_T> cell2Sorted(cell2, sortingDirection);
 
    for (auto &[p1Projection, p1Ptr] : cell1Sorted._particles) {
      for (auto &[p2Projection, p2Ptr] : cell2Sorted._particles) {
        if (std::abs(p2Projection - p1Projection) > _sortingCutoff) {
          break;
        }
 
        interactParticles(*p1Ptr, *p2Ptr);
      }
    }
  } else {
    // Without sorting
    for (auto &p1 : cell1) {
      for (auto &p2 : cell2) {
        interactParticles(p1, p2);
      }
    }
  }
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
SoASortingData CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::computeSortingData(
    const std::vector<std::pair<double, size_t>> &projIdxI, const std::vector<std::pair<double, size_t>> &projIdxJ,
    std::vector<size_t> &maxIndexCache, std::vector<size_t> &minIndexCache) const {
  const size_t nI = projIdxI.size();
  const size_t nJ = projIdxJ.size();
 
  // No j-particles to interact with: skip all i-particles. Guards the projIdxJ[0] access below, which would
  // otherwise be an out-of-bounds read if projIdxJ were empty.
  if (nJ == 0) {
    maxIndexCache.assign(nI, 0);
    minIndexCache.assign(nI, 0);
    return {nI, nI, maxIndexCache, minIndexCache};
  }
 
  // Compute startI: the first i-particle that can interact with any j-particle.
  // Any i with `projI[i] < projJ[0] - cutoff` is strictly farther than cutoff from every j along the
  // sorting axis, so it cannot contribute an interaction and is skipped. Particles exactly at the cutoff
  // distance are kept, consistent with the inclusive cutoff test used below for minIndex/maxIndex.
  const double startThreshold = projIdxJ[0].first - _sortingCutoff;
  auto startIter = std::lower_bound(projIdxI.begin(), projIdxI.end(), startThreshold,
                                    [](const auto &elem, double val) { return elem.first < val; });
  const size_t startI = static_cast<size_t>(startIter - projIdxI.begin());
 
  // Compute endI: the exclusive upper bound of i-particles that can interact with any j-particle. Symmetric
  // to startI: any i with `projI[i] > projJ.back() + cutoff` is strictly farther than cutoff from every j
  // along the sorting axis, so it and all following i (projIdxI is sorted) can be skipped unconditionally.
  const double endThreshold = projIdxJ[nJ - 1].first + _sortingCutoff;
  auto endIter = std::upper_bound(projIdxI.begin(), projIdxI.end(), endThreshold,
                                  [](double val, const auto &elem) { return val < elem.first; });
  const size_t endI = static_cast<size_t>(endIter - projIdxI.begin());
 
  // Compute maxIndexCache and minIndexCache in a single O(nI + nJ) sweep.
  // Both bounds are monotonically non-decreasing with i because projIdxI is sorted, so each pointer only
  // advances, this includes some particles redundant for calculation but avoids a O(n^2) scan.
  // For each i, the j-particles in [minIndex, maxIndex) are the candidates satisfying the 1-D cutoff check.
  maxIndexCache.assign(nI, 0);
  minIndexCache.assign(nI, 0);
  size_t jUpper = 0, jLower = 0;
  for (size_t i = startI; i < endI; ++i) {
    // jUpper: first j where projJ > projI[i] + cutoff (exclusive upper bound).
    while (jUpper < nJ and projIdxJ[jUpper].first <= projIdxI[i].first + _sortingCutoff) {
      ++jUpper;
    }
    maxIndexCache[i] = jUpper;
    // jLower: first j where projJ >= projI[i] - cutoff (lower bound).
    while (jLower < nJ and projIdxJ[jLower].first < projIdxI[i].first - _sortingCutoff) {
      ++jLower;
    }
    minIndexCache[i] = jLower;
  }
 
  return {startI, endI, maxIndexCache, minIndexCache};
}
 
template <class ParticleCell_T, class ParticleFunctor_T, bool bidirectional>
void CellFunctor<ParticleCell_T, ParticleFunctor_T, bidirectional>::processCellPairSoAImpl(
    ParticleCell_T &cell1, ParticleCell_T &cell2, const std::array<double, 3> &sortingDirection) {
  if constexpr (ParticleFunctor_T::supportsSoASorting) {
    if (shouldUseSorting(cell1._particleSoABuffer.size() + cell2._particleSoABuffer.size(), sortingDirection,
                         /*useSoA*/ true)) {
      using Particle_T = ParticleCell_T::ParticleType;
      auto &threadData = _soaThreadData[autopas::autopas_get_thread_num()];
 
      SortedSoAView<Particle_T, ParticleFunctor_T> view1(cell1._particleSoABuffer, sortingDirection,
                                                         threadData.sortedSoa1, threadData.projIdx1);
      SortedSoAView<Particle_T, ParticleFunctor_T> view2(cell2._particleSoABuffer, sortingDirection,
                                                         threadData.sortedSoa2, threadData.projIdx2);
 
      _functor.SoAFunctorPairSorted(
          view1.getView(), view2.getView(),
          computeSortingData(threadData.projIdx1, threadData.projIdx2, threadData.maxIndex, threadData.minIndex),
          _useNewton3);
 
      if constexpr (bidirectional) {
        if (not _useNewton3) {
          // Roles swapped relative to the call above: startI/maxIndex are now trivial (0 / nJ), endI/minIndex do
          // the pruning instead.
          _functor.SoAFunctorPairSorted(
              view2.getView(), view1.getView(),
              computeSortingData(threadData.projIdx2, threadData.projIdx1, threadData.maxIndex, threadData.minIndex),
              false);
        }
      }
 
      view1.scatterBack();
      // if we are not bidirectional and don't use newton3 we don't need to scatter back the second view
      if constexpr (bidirectional) {
        view2.scatterBack();
      } else if (_useNewton3) {
        view2.scatterBack();
      }
      return;
    }
  }
  _functor.SoAFunctorPair(cell1._particleSoABuffer, cell2._particleSoABuffer, _useNewton3);
  if constexpr (bidirectional) {
    if (not _useNewton3) {
      _functor.SoAFunctorPair(cell2._particleSoABuffer, cell1._particleSoABuffer, false);
    }
  }
}
}  // namespace autopas::internal