SlicedBasedTraversal.h

Go to the documentation of this file.

 
#pragma once
 
#include "autopas/containers/TraversalInterface.h"
#include "autopas/containers/cellTraversals/CellTraversal.h"
#include "autopas/containers/verletClusterLists/VerletClusterLists.h"
#include "autopas/utils/DataLayoutConverter.h"
#include "autopas/utils/WrapOpenMP.h"
 
namespace autopas {
 
template <class ParticleCell, class Functor>
class SlicedBasedTraversal : public CellTraversal<ParticleCell>, public TraversalInterface {
 public:
  explicit SlicedBasedTraversal(const std::array<unsigned long, 3> &dims, Functor *functor,
                                const double interactionLength, const std::array<double, 3> &cellLength,
                                DataLayoutOption dataLayout, bool useNewton3, bool spaciallyForward)
      : CellTraversal<ParticleCell>(dims),
        TraversalInterface(dataLayout, useNewton3),
        _overlap{},
        _dimsSortedByLength{},
        _interactionLength(interactionLength),
        _cellLength(cellLength),
        _overlapLongestAxis(0),
        _sliceThickness{},
        _spaciallyForward(spaciallyForward),
        _dataLayoutConverter(functor, dataLayout) {
    this->init();
  }
 
  [[nodiscard]] bool isApplicable() const override {
    auto minSliceThickness = _overlapLongestAxis + 1;
    auto maxNumSlices = this->_cellsPerDimension[_dimsSortedByLength[0]] / minSliceThickness;
    return maxNumSlices > 0;
  }
 
  virtual void initSliceThickness(unsigned long minSliceThickness) {
    auto numSlices = this->_cellsPerDimension[_dimsSortedByLength[0]] / minSliceThickness;
    _sliceThickness.clear();
 
    // abort if domain is too small -> cleared _sliceThickness array indicates non applicability
    if (numSlices < 1) return;
 
    _sliceThickness.insert(_sliceThickness.begin(), numSlices, minSliceThickness);
    auto rest = this->_cellsPerDimension[_dimsSortedByLength[0]] - _sliceThickness[0] * numSlices;
    // remaining slices to distribute the remaining layers on
    auto remSlices = std::min(rest, numSlices);
    for (size_t i = 0; i < remSlices; ++i) {
      _sliceThickness[i] += rest / (remSlices - i);
      rest -= rest / (remSlices - i);
    }
    if (_spaciallyForward) {
      // decreases last _sliceThickness by _overlapLongestAxis to account for the way we handle base cells
      _sliceThickness.back() -= _overlapLongestAxis;
    }
  }
 
  void initTraversal() override {
    loadDataLayout();
    // split domain across its longest dimension
    auto minSliceThickness = _overlapLongestAxis + 1;
    initSliceThickness(minSliceThickness);
  }
 
  void endTraversal() override {
    if (this->_cells) {
      auto &cells = *(this->_cells);
      AUTOPAS_OPENMP(parallel for)
      for (size_t i = 0; i < cells.size(); ++i) {
        _dataLayoutConverter.storeDataLayout(cells[i]);
      }
    }
  }
 
 protected:
  void init();
 
  virtual void loadDataLayout() {
    if (this->_cells) {
      auto &cells = *(this->_cells);
      AUTOPAS_OPENMP(parallel for)
      for (size_t i = 0; i < cells.size(); ++i) {
        _dataLayoutConverter.loadDataLayout(cells[i]);
      }
    }
  }
 
  void reinitForVCL(const VerletClusterLists<typename ParticleCell::ParticleType> *vcl) {
    // Reinitialize the sliced traversal with up to date tower information
    const auto towerSideLength = vcl->getTowerSideLength();
    this->_cellLength = {towerSideLength[0], towerSideLength[1], vcl->getBoxMax()[2] - vcl->getBoxMin()[2]};
    const auto towersPerDim = vcl->getTowersPerDimension();
    this->_cellsPerDimension = {towersPerDim[0], towersPerDim[1], 1};
 
    // reinitialize
    init();
  }
 
  std::array<unsigned long, 3> _overlap;
 
  std::array<int, 3> _dimsSortedByLength;
 
  unsigned long _overlapLongestAxis;
 
  std::vector<unsigned long> _sliceThickness;
 
  bool _spaciallyForward;
 
  std::array<double, 3> _cellLength;
 
 private:
  double _interactionLength;
 
  utils::DataLayoutConverter<Functor> _dataLayoutConverter;
};
 
template <class ParticleCell, class Functor>
void SlicedBasedTraversal<ParticleCell, Functor>::init() {
  for (unsigned int d = 0; d < 3; d++) {
    _overlap[d] = std::ceil(_interactionLength / _cellLength[d]);
    if (not _spaciallyForward) {
      // there is potentially overlap in both directions.
      _overlap[d] *= 2;
    }
  }
 
  // find longest dimension
  auto minMaxElem = std::minmax_element(this->_cellsPerDimension.begin(), this->_cellsPerDimension.end());
  _dimsSortedByLength[0] = (int)std::distance(this->_cellsPerDimension.begin(), minMaxElem.second);
  _dimsSortedByLength[2] = (int)std::distance(this->_cellsPerDimension.begin(), minMaxElem.first);
  _dimsSortedByLength[1] = 3 - (_dimsSortedByLength[0] + _dimsSortedByLength[2]);
 
  _overlapLongestAxis = _overlap[_dimsSortedByLength[0]];
}
 
}  // namespace autopas