mac_reductions.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_MAC_REDUCTIONS_HPP
#define PECLET_FLOW_MAC_REDUCTIONS_HPP
 
#include <cstddef>
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
namespace peclet::flow {
 
using Exec = Kokkos::DefaultExecutionSpace;
using Mem = Exec::memory_space;
using DField = Kokkos::View<double*, Mem>;  // flat extended-block cell field (x-fastest)
using DConst = Kokkos::View<const double*, Mem>;
 
struct Ext3 {
  int x, y, z;
};
 
struct SumMax {
  double sum = 0.0;
  double maxabs = 0.0;
  KOKKOS_INLINE_FUNCTION SumMax& operator+=(const SumMax& o) {
    sum += o.sum;
    maxabs = (o.maxabs > maxabs) ? o.maxabs : maxabs;
    return *this;
  }
};
 
// Map an inner-cell linear index c to the extended-block flat index (x-fastest, +ghost offset).
KOKKOS_INLINE_FUNCTION std::size_t innerToExt(long c, Ext3 ext, int ghost, Ext3 inner) {
  const int ix = static_cast<int>(c % inner.x);
  const int iy = static_cast<int>((c / inner.x) % inner.y);
  const int iz = static_cast<int>(c / (static_cast<long>(inner.x) * inner.y));
  return static_cast<std::size_t>(ix + ghost) + static_cast<std::size_t>(iy + ghost) * ext.x +
         static_cast<std::size_t>(iz + ghost) * static_cast<std::size_t>(ext.x) * ext.y;
}
 
inline SumMax localSumMax(DConst f, Ext3 ext, int ghost, Ext3 inner) {
  const long n = static_cast<long>(inner.x) * inner.y * inner.z;
  SumMax r;
  if (n <= 0)
    return r;
  Kokkos::parallel_reduce(
      "peclet::flow::sum_max", Kokkos::RangePolicy<Exec>(0, n),
      KOKKOS_LAMBDA(long c, SumMax& acc) {
        const double v = f(innerToExt(c, ext, ghost, inner));
        acc.sum += v;
        const double a = Kokkos::fabs(v);
        if (a > acc.maxabs)
          acc.maxabs = a;
      },
      r);
  return r;
}
 
inline double localDot(DConst a, DConst b, Ext3 ext, int ghost, Ext3 inner) {
  const long n = static_cast<long>(inner.x) * inner.y * inner.z;
  double s = 0.0;
  if (n <= 0)
    return s;
  Kokkos::parallel_reduce(
      "peclet::flow::dot", Kokkos::RangePolicy<Exec>(0, n),
      KOKKOS_LAMBDA(long c, double& acc) {
        const std::size_t i = innerToExt(c, ext, ghost, inner);
        acc += a(i) * b(i);
      },
      s);
  return s;
}
 
inline void subtractAll(DField f, double m) {
  const std::size_t n = f.extent(0);
  Kokkos::parallel_for(
      "peclet::flow::subtract", Kokkos::RangePolicy<Exec>(0, n),
      KOKKOS_LAMBDA(std::size_t i) { f(i) -= m; });
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_MAC_REDUCTIONS_HPP