plane_policy.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_PLANE_POLICY_HPP
#define PECLET_VORO_PLANE_POLICY_HPP
 
#include <Kokkos_Core.hpp>
 
namespace peclet::voro {
 
struct Voronoi {
  static constexpr bool kHasWeightDof = false;
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static Real offsetFromRel(const Real r[3], Real wSelf, Real wNbr) {
    (void)wSelf;
    (void)wNbr;
    return Real(0.5) * (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static Real blockReachSq(Real rSqMax, Real wSelf, Real wMaxAll) {
    (void)wSelf;
    (void)wMaxAll;
    return Real(4) * rSqMax;
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static void planeFromNeighbour(const Real pSelf[3], const Real pNbr[3],
                                                        Real wSelf, Real wNbr, Real L, Real pdir[3],
                                                        Real& off) {
    const Real Lh = Real(0.5) * L;
    for (int d = 0; d < 3; ++d) {
      Real r = pNbr[d] - pSelf[d];
      r = r > Lh ? r - L : (r < -Lh ? r + L : r);
      pdir[d] = r;
    }
    off = offsetFromRel(pdir, wSelf, wNbr);
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static void buildNormal(const Real pSelf[3], const Real pNbr[3],
                                                 Real wSelf, Real wNbr, Real L, Real nOut[3]) {
    Real pdir[3], off;
    planeFromNeighbour(pSelf, pNbr, wSelf, wNbr, L, pdir, off);
    const Real l2 = pdir[0] * pdir[0] + pdir[1] * pdir[1] + pdir[2] * pdir[2];
    const Real a = l2 > Real(0) ? off / l2 : Real(0);
    for (int d = 0; d < 3; ++d)
      nOut[d] = a * pdir[d];
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static void chain(const Real g[3], const Real r[3], Real rho, Real c,
                                           Real fSelf[3], Real fNbr[3], Real& fwSelf, Real& fwNbr) {
    (void)r;
    (void)rho;
    (void)c;
    for (int d = 0; d < 3; ++d) {
      fSelf[d] = Real(-0.5) * g[d];
      fNbr[d] = Real(0.5) * g[d];
    }
    fwSelf = Real(0);
    fwNbr = Real(0);
  }
};
 
struct Power {
  static constexpr bool kHasWeightDof = true;
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static Real offsetFromRel(const Real r[3], Real wSelf, Real wNbr) {
    return Real(0.5) * (r[0] * r[0] + r[1] * r[1] + r[2] * r[2] + wSelf - wNbr);
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static Real blockReachSq(Real rSqMax, Real wSelf, Real wMaxAll) {
    const Real rm = Kokkos::sqrt(rSqMax);
    const Real span = wMaxAll - wSelf > Real(0) ? wMaxAll - wSelf : Real(0);
    const Real t = rm + Kokkos::sqrt(rSqMax + span);
    return t * t;
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static void planeFromNeighbour(const Real pSelf[3], const Real pNbr[3],
                                                        Real wSelf, Real wNbr, Real L, Real pdir[3],
                                                        Real& off) {
    const Real Lh = Real(0.5) * L;
    for (int d = 0; d < 3; ++d) {
      Real r = pNbr[d] - pSelf[d];
      r = r > Lh ? r - L : (r < -Lh ? r + L : r);
      pdir[d] = r;
    }
    off = offsetFromRel(pdir, wSelf, wNbr);
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static void buildNormal(const Real pSelf[3], const Real pNbr[3],
                                                 Real wSelf, Real wNbr, Real L, Real nOut[3]) {
    Real pdir[3], off;
    planeFromNeighbour(pSelf, pNbr, wSelf, wNbr, L, pdir, off);
    const Real l2 = pdir[0] * pdir[0] + pdir[1] * pdir[1] + pdir[2] * pdir[2];
    const Real a = l2 > Real(0) ? off / l2 : Real(0);
    for (int d = 0; d < 3; ++d)
      nOut[d] = a * pdir[d];
  }
 
  template <class Real>
  KOKKOS_INLINE_FUNCTION static void chain(const Real g[3], const Real r[3], Real rho, Real c,
                                           Real fSelf[3], Real fNbr[3], Real& fwSelf, Real& fwNbr) {
    const Real invRho = rho > Real(0) ? Real(1) / rho : Real(0);
    const Real alpha = Real(0.5) + Real(0.5) * c * invRho;
    const Real rg = r[0] * g[0] + r[1] * g[1] + r[2] * g[2];
    const Real kk = c * invRho * invRho * rg;  // (c/ρ²)(r·g)
    for (int d = 0; d < 3; ++d) {
      const Real jn = alpha * g[d] - kk * r[d];  // dGeom/dp_nbr component
      fNbr[d] = jn;
      fSelf[d] = -jn;
    }
    const Real w = Real(0.5) * invRho * rg;  // (g·r)/(2ρ)
    fwSelf = w;
    fwNbr = -w;
  }
};
 
template <class Policy, class Cell, class Real>
KOKKOS_INLINE_FUNCTION void chainToDofs(const Cell& c, const Real seed[3], const Real* pos,
                                        Real wSelf, const Real* weight, Real L, const Real* gx,
                                        const Real* gy, const Real* gz, Real fSelf[3], Real& fwSelf,
                                        Real* fnx, Real* fny, Real* fnz, Real* fwn) {
  fSelf[0] = fSelf[1] = fSelf[2] = Real(0);
  fwSelf = Real(0);
  const Real Lh = Real(0.5) * L;
  for (int k = 0; k < c.np; ++k) {
    const int j = c.pnbr[k];
    if (j < 0) {  // box (−1) / SDF wall (−2): no particle-DOF dependence here
      fnx[k] = fny[k] = fnz[k] = Real(0);
      fwn[k] = Real(0);
      continue;
    }
    const Real g[3] = {gx[k], gy[k], gz[k]};
    Real r[3] = {Real(0), Real(0), Real(0)}, rho = Real(0), cc = Real(0);
    if constexpr (Policy::kHasWeightDof) {
      for (int d = 0; d < 3; ++d) {
        Real rr = pos[3 * j + d] - seed[d];
        rr = rr > Lh ? rr - L : (rr < -Lh ? rr + L : rr);
        r[d] = rr;
      }
      rho = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
      cc = wSelf - weight[j];
    }
    Real fs[3], fn[3], fws = Real(0), fwn1 = Real(0);
    Policy::template chain<Real>(g, r, rho, cc, fs, fn, fws, fwn1);
    fSelf[0] += fs[0];
    fSelf[1] += fs[1];
    fSelf[2] += fs[2];
    fwSelf += fws;
    fnx[k] = fn[0];
    fny[k] = fn[1];
    fnz[k] = fn[2];
    fwn[k] = fwn1;
  }
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_PLANE_POLICY_HPP