SlicedC02BasedTraversal.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
 
#include "SlicedBasedTraversal.h"
#include "autopas/utils/DataLayoutConverter.h"
#include "autopas/utils/ThreeDimensionalMapping.h"
#include "autopas/utils/WrapOpenMP.h"
 
namespace autopas {
 
template <class ParticleCell, class Functor>
class SlicedC02BasedTraversal : public SlicedBasedTraversal<ParticleCell, Functor> {
 public:
  explicit SlicedC02BasedTraversal(const std::array<unsigned long, 3> &dims, Functor *functor,
                                   const double interactionLength, const std::array<double, 3> &cellLength,
                                   DataLayoutOption dataLayout, bool useNewton3, bool spaciallyForward)
      : SlicedBasedTraversal<ParticleCell, Functor>(dims, functor, interactionLength, cellLength, dataLayout,
                                                    useNewton3, spaciallyForward) {}
 
  template <typename LoopBody>
  inline void cSlicedTraversal(LoopBody &&loopBody);
 
  [[nodiscard]] bool isApplicable() const override {
    return this->_cellsPerDimension[this->_dimsSortedByLength[0]] >= this->_overlapLongestAxis;
  }
 
  void initTraversal() override {
    this->loadDataLayout();
    // split domain across its longest dimension
    auto minSliceThickness = this->_overlapLongestAxis;
    this->initSliceThickness(minSliceThickness);
  }
};
 
template <class ParticleCell, class Functor>
template <typename LoopBody>
void SlicedC02BasedTraversal<ParticleCell, Functor>::cSlicedTraversal(LoopBody &&loopBody) {
  using std::array;
 
  auto numSlices = this->_sliceThickness.size();
  // check if applicable
 
  std::array<size_t, 2> overLapps23{this->_overlap[this->_dimsSortedByLength[1]],
                                    this->_overlap[this->_dimsSortedByLength[2]]};
 
  if (not this->_spaciallyForward) {
    overLapps23 = {0ul, 0ul};
  }
 
  for (size_t offset = 0; offset < 2; offset++) {
    // although every thread gets exactly one iteration (=slice) this is faster than a normal parallel region
    AUTOPAS_OPENMP(parallel for schedule(dynamic, 1))
    for (size_t slice = offset; slice < numSlices; slice += 2) {
      array<unsigned long, 3> myStartArray{0, 0, 0};
      for (size_t i = 0; i < slice; ++i) {
        myStartArray[this->_dimsSortedByLength[0]] += this->_sliceThickness[i];
      }
 
      const auto lastLayer = myStartArray[this->_dimsSortedByLength[0]] + this->_sliceThickness[slice];
      for (unsigned long dimSlice = myStartArray[this->_dimsSortedByLength[0]]; dimSlice < lastLayer; ++dimSlice) {
        for (unsigned long dimMedium = 0;
             dimMedium < this->_cellsPerDimension[this->_dimsSortedByLength[1]] - overLapps23[0]; ++dimMedium) {
          for (unsigned long dimShort = 0;
               dimShort < this->_cellsPerDimension[this->_dimsSortedByLength[2]] - overLapps23[1]; ++dimShort) {
            array<unsigned long, 3> idArray = {};
            idArray[this->_dimsSortedByLength[0]] = dimSlice;
            idArray[this->_dimsSortedByLength[1]] = dimMedium;
            idArray[this->_dimsSortedByLength[2]] = dimShort;
            loopBody(idArray[0], idArray[1], idArray[2]);
          }
        }
      }
    }
  }
}
 
}  // namespace autopas