verlet_skin.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_VERLET_SKIN_HPP
#define PECLET_VORO_VERLET_SKIN_HPP
 
#include <Kokkos_Core.hpp>
#include <string>
 
#include "peclet/core/common/view.hpp"
 
namespace peclet::voro {
 
enum SkinTrigger : int {
  kSkinNone = 0,
  kSkinMover = 1,
  // kSkinBoundary = 2,  // DEFERRED (Risk 1d): |SDF(p)| < security radius contact-status change.
};
 
template <class Real>
struct VerletSkin {
  using Mem = peclet::core::MemSpace;
  using Exec = peclet::core::ExecSpace;
  Kokkos::View<Real*, Mem> xRef;  
  Real skin = 0;                  
  int N = 0;
 
  void alloc(int n, Real skinWidth) {
    N = n;
    skin = skinWidth;
    xRef = Kokkos::View<Real*, Mem>(
        Kokkos::view_alloc(std::string("verlet.xRef"), Kokkos::WithoutInitializing), (size_t)n * 3);
  }
 
  void reset(const Kokkos::View<Real*, Mem>& pos) { Kokkos::deep_copy(xRef, pos); }
};
 
template <class Real>
int flagSkinMovers(const Kokkos::View<Real*, peclet::core::MemSpace>& pos,
                   const Kokkos::View<Real*, peclet::core::MemSpace>& xRef, Real skin,
                   const Real L[3], const Kokkos::View<int*, peclet::core::MemSpace>& outFlags) {
  using Exec = peclet::core::ExecSpace;
  const int N = (int)outFlags.extent(0);
  const Real half2 = Real(0.25) * skin * skin;  // (skin/2)^2
  const Real Lx = L[0], Ly = L[1], Lz = L[2];
  const Real Lxh = Real(0.5) * Lx, Lyh = Real(0.5) * Ly, Lzh = Real(0.5) * Lz;
  int count = 0;
  Kokkos::parallel_reduce(
      "verlet.flag", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, int& acc) {
        Real dx = pos(3 * i + 0) - xRef(3 * i + 0);
        Real dy = pos(3 * i + 1) - xRef(3 * i + 1);
        Real dz = pos(3 * i + 2) - xRef(3 * i + 2);
        dx = dx > Lxh ? dx - Lx : (dx < -Lxh ? dx + Lx : dx);
        dy = dy > Lyh ? dy - Ly : (dy < -Lyh ? dy + Ly : dy);
        dz = dz > Lzh ? dz - Lz : (dz < -Lzh ? dz + Lz : dz);
        const Real d2 = dx * dx + dy * dy + dz * dz;
        const bool mover = d2 > half2;
        outFlags(i) = mover ? kSkinMover : kSkinNone;
        if (mover)
          ++acc;
      },
      count);
  return count;
}
 
template <class Real>
Real maxDisplacement(const Kokkos::View<Real*, peclet::core::MemSpace>& pos,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& xRef, const Real L[3]) {
  using Exec = peclet::core::ExecSpace;
  const int N = (int)(xRef.extent(0) / 3);
  const Real Lx = L[0], Ly = L[1], Lz = L[2];
  const Real Lxh = Real(0.5) * Lx, Lyh = Real(0.5) * Ly, Lzh = Real(0.5) * Lz;
  Real m2 = 0;
  Kokkos::parallel_reduce(
      "verlet.maxDisp", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, Real& lm) {
        Real dx = pos(3 * i + 0) - xRef(3 * i + 0);
        Real dy = pos(3 * i + 1) - xRef(3 * i + 1);
        Real dz = pos(3 * i + 2) - xRef(3 * i + 2);
        dx = dx > Lxh ? dx - Lx : (dx < -Lxh ? dx + Lx : dx);
        dy = dy > Lyh ? dy - Ly : (dy < -Lyh ? dy + Ly : dy);
        dz = dz > Lzh ? dz - Lz : (dz < -Lzh ? dz + Lz : dz);
        const Real d2 = dx * dx + dy * dy + dz * dz;
        if (d2 > lm)
          lm = d2;
      },
      Kokkos::Max<Real>(m2));
  return Kokkos::sqrt(m2);
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_VERLET_SKIN_HPP