LCC01Traversal.h

Go to the documentation of this file.

 
#pragma once
 
#include "LCTraversalInterface.h"
#include "autopas/baseFunctors/CellFunctor.h"
#include "autopas/baseFunctors/CellFunctor3B.h"
#include "autopas/containers/cellTraversals/C01BasedTraversal.h"
#include "autopas/options/DataLayoutOption.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/ArrayUtils.h"
#include "autopas/utils/WrapOpenMP.h"
#include "autopas/utils/checkFunctorType.h"
 
namespace autopas {
 
template <class ParticleCell, class Functor, bool combineSoA = false>
class LCC01Traversal : public C01BasedTraversal<ParticleCell, Functor, (combineSoA ? 2 : 3)>,
                       public LCTraversalInterface {
 public:
  explicit LCC01Traversal(const std::array<unsigned long, 3> &dims, Functor &functor, const double interactionLength,
                          const std::array<double, 3> &cellLength, DataLayoutOption dataLayout, bool useNewton3)
      : C01BasedTraversal<ParticleCell, Functor, (combineSoA ? 2 : 3)>(dims, functor, interactionLength, cellLength,
                                                                       dataLayout, useNewton3),
        _cellFunctor(functor, interactionLength /*should use cutoff here, if not used to build verlet-lists*/,
                     dataLayout, useNewton3),
        _functor(functor),
        _cacheOffset(DEFAULT_CACHE_LINE_SIZE / sizeof(unsigned int)) {
    computeOffsets();
  }
 
  void computeOffsets();
 
  void traverseParticles() override;
 
  [[nodiscard]] bool isApplicableToDomain() const override { return true; }
 
  [[nodiscard]] TraversalOption getTraversalType() const override {
    return (combineSoA) ? TraversalOption::lc_c01_combined_SoA : TraversalOption::lc_c01;
  }
 
  void setAoSSortingThresholds(const SortingThresholdInfoInterface &aosSortingThreshold) override {
    _cellFunctor.setAoSSortingThresholds(aosSortingThreshold);
  }
 
  void setSoASortingThresholds(const SortingThresholdInfoInterface &soaSortingThreshold) override {
    _cellFunctor.setSoASortingThresholds(soaSortingThreshold);
  }
 
 private:
  // CellOffsets needs to store interaction pairs or triplets depending on the Functor type.
  using CellOffsetsType = std::conditional_t<decltype(utils::isPairwiseFunctor<Functor>())::value,
                                             std::vector<std::vector<std::pair<long, std::array<double, 3>>>>,
                                             std::vector<std::tuple<long, long, std::array<double, 3>>>>;
 
  // CellFunctor type for either Pairwise or Triwise Functors.
  using CellFunctorType = std::conditional_t<decltype(utils::isPairwiseFunctor<Functor>())::value,
                                             internal::CellFunctor<ParticleCell, Functor, /*bidirectional*/ false>,
                                             internal::CellFunctor3B<ParticleCell, Functor, /*bidirectional*/ false>>;
 
  inline void processBaseCell(std::vector<ParticleCell> &cells, unsigned long x, unsigned long y, unsigned long z);
 
  inline void processBaseCellPairwise(std::vector<ParticleCell> &cells, unsigned long x, unsigned long y,
                                      unsigned long z);
 
  inline void processBaseCellTriwise(std::vector<ParticleCell> &cells, unsigned long x, unsigned long y,
                                     unsigned long z);
 
  void computePairwiseOffsets();
 
  void computeTriwiseOffsets();
 
  template <std::size_t... I>
  constexpr void appendNeeded(ParticleCell &cell, ParticleCell &appendCell, std::index_sequence<I...>) {
    cell._particleSoABuffer.template append<std::get<I>(Functor::getNeededAttr(std::false_type()))...>(
        appendCell._particleSoABuffer);
  }
 
  void resizeBuffers();
 
  CellOffsetsType _cellOffsets;
 
  CellFunctorType _cellFunctor;
 
  Functor &_functor;
 
  std::vector<std::vector<ParticleCell>> _combinationSlices;
 
  std::vector<unsigned int> _currentSlices;
 
  const unsigned int _cacheOffset;
};
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::computeOffsets() {
  if constexpr (utils::isPairwiseFunctor<Functor>()) {
    computePairwiseOffsets();
  } else if (utils::isTriwiseFunctor<Functor>()) {
    computeTriwiseOffsets();
  } else {
    utils::ExceptionHandler::exception("LCC01Traversal::computeOffsets(): Functor is not valid.");
  }
}
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::computePairwiseOffsets() {
  _cellOffsets.resize(2 * this->_overlap[0] + 1);
 
  const auto interactionLengthSquare{this->_interactionLength * this->_interactionLength};
 
  for (long x = -this->_overlap[0]; x <= 0l; ++x) {
    for (long y = -this->_overlap[1]; y <= static_cast<long>(this->_overlap[1]); ++y) {
      for (long z = -this->_overlap[2]; z <= static_cast<long>(this->_overlap[2]); ++z) {
        const std::array<double, 3> pos = {
            std::max(0l, (std::abs(x) - 1l)) * this->_cellLength[0],
            std::max(0l, (std::abs(y) - 1l)) * this->_cellLength[1],
            std::max(0l, (std::abs(z) - 1l)) * this->_cellLength[2],
        };
        const double distSquare = utils::ArrayMath::dot(pos, pos);
        if (distSquare <= interactionLengthSquare) {
          const long currentOffset = utils::ThreeDimensionalMapping::threeToOneD(
              x, y, z, utils::ArrayUtils::static_cast_copy_array<long>(this->_cellsPerDimension));
          const bool containCurrentOffset =
              std::any_of(_cellOffsets[x + this->_overlap[0]].cbegin(), _cellOffsets[x + this->_overlap[0]].cend(),
                          [currentOffset](const auto &e) { return e.first == currentOffset; });
          if (containCurrentOffset) {
            continue;
          }
          for (long ix = x; ix <= std::abs(x); ++ix) {
            const long offset = utils::ThreeDimensionalMapping::threeToOneD(
                ix, y, z, utils::ArrayUtils::static_cast_copy_array<long>(this->_cellsPerDimension));
            const size_t index = ix + this->_overlap[0];
 
            // Calculate the sorting direction from the base cell (x, y, z) and the other cell by use of the offset (ix,
            // y, z).
            std::array<double, 3> sortingDir = {static_cast<double>(ix) * this->_cellLength[0],
                                                static_cast<double>(y) * this->_cellLength[1],
                                                static_cast<double>(z) * this->_cellLength[2]};
 
            // the offset to the current cell itself is zero.
            if (ix == 0 and y == 0 and z == 0) {
              sortingDir = {1., 1., 1.};
            }
            sortingDir = utils::ArrayMath::normalize(sortingDir);
 
            if (y == 0l and z == 0l) {
              // make sure center of slice is always at the beginning
              _cellOffsets[index].insert(_cellOffsets[index].cbegin(), std::make_pair(offset, sortingDir));
            } else {
              _cellOffsets[index].emplace_back(offset, sortingDir);
            }
          }
        }
      }
    }
  }
}
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::computeTriwiseOffsets() {
  using namespace utils::ArrayMath::literals;
  // Reserve approximately. Overestimates more for larger overlap.
  const int cubeSize = this->_overlap[0] * this->_overlap[1] * this->_overlap[2];
  _cellOffsets.reserve(cubeSize * cubeSize / 4);
 
  // Helper function to get minimal distance between two cells
  auto cellDistance = [&](long x1, long y1, long z1, long x2, long y2, long z2) {
    return std::array<double, 3>{std::max(0l, (std::abs(x1 - x2) - 1l)) * this->_cellLength[0],
                                 std::max(0l, (std::abs(y1 - y2) - 1l)) * this->_cellLength[1],
                                 std::max(0l, (std::abs(z1 - z2) - 1l)) * this->_cellLength[2]};
  };
 
  const auto interactionLengthSquare{this->_interactionLength * this->_interactionLength};
  _cellOffsets.emplace_back(0, 0, std::array<double, 3>{1., 1., 1.});
 
  // offsets for the first cell
  for (long x1 = -this->_overlap[0]; x1 <= static_cast<long>(this->_overlap[0]); ++x1) {
    for (long y1 = -this->_overlap[1]; y1 <= static_cast<long>(this->_overlap[1]); ++y1) {
      for (long z1 = -this->_overlap[2]; z1 <= static_cast<long>(this->_overlap[2]); ++z1) {
        // check distance between base cell and cell 1
        const auto dist01 = cellDistance(0l, 0l, 0l, x1, y1, z1);
 
        const double distSquare = utils::ArrayMath::dot(dist01, dist01);
        if (distSquare > interactionLengthSquare) continue;
 
        // offsets for cell 2
        for (long x2 = -this->_overlap[0]; x2 <= static_cast<long>(this->_overlap[0]); ++x2) {
          for (long y2 = -this->_overlap[1]; y2 <= static_cast<long>(this->_overlap[1]); ++y2) {
            for (long z2 = -this->_overlap[2]; z2 <= static_cast<long>(this->_overlap[2]); ++z2) {
              // check distance between cell 1 and cell 2
              const auto dist12 = cellDistance(x1, y1, z1, x2, y2, z2);
 
              const double dist12Squared = utils::ArrayMath::dot(dist12, dist12);
              if (dist12Squared > interactionLengthSquare) continue;
 
              // check distance between base cell and cell 2
              const auto dist02 = cellDistance(0l, 0l, 0l, x2, y2, z2);
 
              const double dist02Squared = utils::ArrayMath::dot(dist02, dist02);
              if (dist02Squared > interactionLengthSquare) continue;
 
              const long offset1 = utils::ThreeDimensionalMapping::threeToOneD(
                  x1, y1, z1, utils::ArrayUtils::static_cast_copy_array<long>(this->_cellsPerDimension));
 
              const long offset2 = utils::ThreeDimensionalMapping::threeToOneD(
                  x2, y2, z2, utils::ArrayUtils::static_cast_copy_array<long>(this->_cellsPerDimension));
 
              // Only add unique combinations. E.g.: (5, 8) == (8, 5)
              if (offset2 <= offset1) continue;
 
              // sorting direction from base cell to the first different cell
              std::array<double, 3> sortDirection{};
              if (offset1 == 0) {
                sortDirection = {x2 * this->_cellLength[0], y2 * this->_cellLength[1], z2 * this->_cellLength[2]};
              } else {
                sortDirection = {x1 * this->_cellLength[0], y1 * this->_cellLength[1], z1 * this->_cellLength[2]};
              }
              _cellOffsets.emplace_back(offset1, offset2, utils::ArrayMath::normalize(sortDirection));
            }
          }
        }
      }
    }
  }
}
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::processBaseCell(std::vector<ParticleCell> &cells,
                                                                               unsigned long x, unsigned long y,
                                                                               unsigned long z) {
  if constexpr (utils::isPairwiseFunctor<Functor>()) {
    processBaseCellPairwise(cells, x, y, z);
  } else if constexpr (utils::isTriwiseFunctor<Functor>()) {
    processBaseCellTriwise(cells, x, y, z);
  } else {
    utils::ExceptionHandler::exception(
        "LCC01Traversal::processBaseCell(): Functor {} is not of type PairwiseFunctor or TriwiseFunctor.",
        _functor.getName());
  }
}
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::processBaseCellPairwise(std::vector<ParticleCell> &cells,
                                                                                       unsigned long x, unsigned long y,
                                                                                       unsigned long z) {
  unsigned long baseIndex = utils::ThreeDimensionalMapping::threeToOneD(x, y, z, this->_cellsPerDimension);
  ParticleCell &baseCell = cells[baseIndex];
  const size_t cOffSize = _cellOffsets.size();
 
  if constexpr (combineSoA) {
    // Iteration along x
 
    const auto threadID = static_cast<size_t>(autopas_get_thread_num());
    auto &currentSlice = _currentSlices[threadID * _cacheOffset];
    auto &combinationSlice = _combinationSlices[threadID];
 
    // First cell needs to initialize whole buffer
    if (x == this->_overlap[0]) {
      currentSlice = 0;
      for (unsigned int offsetSlice = 0; offsetSlice < cOffSize; offsetSlice++) {
        combinationSlice[offsetSlice]._particleSoABuffer.clear();
        for (const auto &offset : _cellOffsets[offsetSlice]) {
          const unsigned long otherIndex = baseIndex + offset.first;
          ParticleCell &otherCell = cells[otherIndex];
          appendNeeded(combinationSlice[offsetSlice], otherCell,
                       std::make_index_sequence<Functor::getNeededAttr(std::false_type()).size()>{});
        }
      }
    } else {
      // reduce size
      size_t i = 0;
      const size_t midSlice = (currentSlice + this->_overlap[0] + 1) % cOffSize;
      for (size_t slice = (currentSlice + 1) % cOffSize; slice != midSlice; ++slice %= cOffSize, ++i) {
        size_t newSize = 0;
        for (const auto &offset : _cellOffsets[i]) {
          const unsigned long otherIndex = baseIndex + offset.first;
          ParticleCell &otherCell = cells[otherIndex];
          newSize += otherCell.size();
        }
        combinationSlice[slice]._particleSoABuffer.resizeArrays(newSize);
      }
      // append buffers
      for (size_t slice = midSlice; slice != currentSlice; ++slice %= cOffSize, ++i) {
        for (auto offsetIndex = _cellOffsets[(i + 1) % cOffSize].size(); offsetIndex < _cellOffsets[i].size();
             ++offsetIndex) {
          const unsigned long otherIndex = baseIndex + _cellOffsets[i][offsetIndex].first;
          ParticleCell &otherCell = cells[otherIndex];
          appendNeeded(combinationSlice[slice], otherCell,
                       std::make_index_sequence<Functor::getNeededAttr(std::false_type()).size()>{});
        }
      }
 
      combinationSlice[currentSlice]._particleSoABuffer.clear();
 
      for (const auto &offset : _cellOffsets.back()) {
        const unsigned long otherIndex = baseIndex + offset.first;
        ParticleCell &otherCell = cells[otherIndex];
        appendNeeded(combinationSlice[currentSlice], otherCell,
                     std::make_index_sequence<Functor::getNeededAttr(std::false_type()).size()>{});
      }
 
      ++currentSlice %= cOffSize;
    }
 
    // calculate all interactions
    for (unsigned int slice = 0; slice < cOffSize; slice++) {
      if (slice == (currentSlice + this->_overlap[0]) % cOffSize) {
        // slice contains base cell -> skip particles of base cell. This is not supported by CellFunctor, so call
        // pairwise functor directly.
        auto startIndex = baseCell.size();
        auto endIndex = combinationSlice[slice]._particleSoABuffer.size();
        _functor.SoAFunctorPair(baseCell._particleSoABuffer,
                                {&(combinationSlice[slice]._particleSoABuffer), startIndex, endIndex}, false);
        // compute base cell
        this->_cellFunctor.processCell(baseCell);
      } else {
        this->_cellFunctor.processCellPair(baseCell, combinationSlice[slice]);
      }
    }
  } else {
    for (const auto &slice : _cellOffsets) {
      for (auto const &[offset, r] : slice) {
        const unsigned long otherIndex = baseIndex + offset;
        ParticleCell &otherCell = cells[otherIndex];
 
        if (baseIndex == otherIndex) {
          this->_cellFunctor.processCell(baseCell);
        } else {
          this->_cellFunctor.processCellPair(baseCell, otherCell, r);
        }
      }
    }
  }
}
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::processBaseCellTriwise(std::vector<ParticleCell> &cells,
                                                                                      unsigned long x, unsigned long y,
                                                                                      unsigned long z) {
  unsigned long baseIndex = utils::ThreeDimensionalMapping::threeToOneD(x, y, z, this->_cellsPerDimension);
  ParticleCell &baseCell = cells[baseIndex];
 
  for (auto const &[offset1, offset2, r] : _cellOffsets) {
    const unsigned long otherIndex1 = baseIndex + offset1;
    const unsigned long otherIndex2 = baseIndex + offset2;
    ParticleCell &otherCell1 = cells[otherIndex1];
    ParticleCell &otherCell2 = cells[otherIndex2];
 
    if (baseIndex == otherIndex1 and baseIndex == otherIndex2) {
      this->_cellFunctor.processCell(baseCell);
    } else if (baseIndex == otherIndex1 and baseIndex != otherIndex2) {
      this->_cellFunctor.processCellPair(baseCell, otherCell2);
    } else if (baseIndex != otherIndex1 and baseIndex == otherIndex2) {
      this->_cellFunctor.processCellPair(baseCell, otherCell1);
    } else if (baseIndex != otherIndex1 and otherIndex1 == otherIndex2) {
      this->_cellFunctor.processCellPair(baseCell, otherCell1);
    } else {
      this->_cellFunctor.processCellTriple(baseCell, otherCell1, otherCell2, r);
    }
  }
}
 
template <class ParticleCell, class PairwiseFunctor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, PairwiseFunctor, combineSoA>::resizeBuffers() {
  const auto numThreads = static_cast<size_t>(autopas_get_max_threads());
  if (_combinationSlices.size() != numThreads) {
    _combinationSlices.resize(numThreads);
    const auto cellOffsetsSize = _cellOffsets.size();
    std::for_each(_combinationSlices.begin(), _combinationSlices.end(),
                  [cellOffsetsSize](auto &e) { e.resize(cellOffsetsSize); });
    _currentSlices.resize(numThreads * _cacheOffset);
  }
}
 
template <class ParticleCell, class Functor, bool combineSoA>
inline void LCC01Traversal<ParticleCell, Functor, combineSoA>::traverseParticles() {
  auto &cells = *(this->_cells);
  if (not this->isApplicableToDomain()) {
    if constexpr (combineSoA) {
      utils::ExceptionHandler::exception(
          "The C01 traversal with combined SoA buffers cannot work with data layout AoS and enabled newton3 (unless "
          "only one thread is used)!");
    } else {
      utils::ExceptionHandler::exception(
          "The C01 traversal cannot work with enabled newton3 (unless only one thread is used)!");
    }
  }
  if constexpr (combineSoA) {
    resizeBuffers();
  }
  this->c01Traversal([&](unsigned long x, unsigned long y, unsigned long z) { this->processBaseCell(cells, x, y, z); });
}
 
}  // namespace autopas