periodicity.hpp

Go to the documentation of this file.

 
#ifndef DEM_PERIODICITY_HPP
#define DEM_PERIODICITY_HPP
 
#include <Kokkos_Core.hpp>
 
#include "dem_portable.hpp"
#include "integration.hpp"  // Domain, V3/V4/Vf/Vi, detail::st3/st4
 
namespace peclet::dem {
 
inline int calculateGhostCapacity(int nReal, Domain dom, float skin) {
  const float sx = dom.periodic_x ? skin : 0.0f;
  const float sy = dom.periodic_y ? skin : 0.0f;
  const float sz = dom.periodic_z ? skin : 0.0f;
  const double innerX = (dom.size.x - 2.0f * sx > 0.0f) ? (dom.size.x - 2.0f * sx) : 0.0;
  const double innerY = (dom.size.y - 2.0f * sy > 0.0f) ? (dom.size.y - 2.0f * sy) : 0.0;
  const double innerZ = (dom.size.z - 2.0f * sz > 0.0f) ? (dom.size.z - 2.0f * sz) : 0.0;
  const double volTotal = static_cast<double>(dom.size.x) * dom.size.y * dom.size.z;
  if (volTotal <= 0.0)
    return nReal * 8;  // degenerate box (domain not set yet) -> generous fallback
  const double ghostFraction = (volTotal - innerX * innerY * innerZ) / volTotal;
  const double estGhosts = (ghostFraction > 0.5) ? nReal * 32.0 : nReal * ghostFraction * 4.0;
  return nReal + static_cast<int>(estGhosts) + 4096;  // + fixed buffer (matches CUDA)
}
 
inline void generateGhostsKokkos(int numReal, int capacity, Domain dom, float skin, V3 pos,
                                 Vf invMass, V3 posPred, V3 vel, V3 velPred, V4 quat, V4 quatPred,
                                 V3 angVel, V3 angVelPred, Vf scale, Vi shapeId, Vi realIndices,
                                 Kokkos::View<int, CpMem> topGhost, Vi gid = Vi{},
                                 Kokkos::View<unsigned char*, CpMem> materialId = {}) {
  using detail::st3;
  using detail::st4;
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::generate_ghosts", Kokkos::RangePolicy<CpExec>(space, 0, numReal),
      KOKKOS_LAMBDA(int i) {
        const F3 p = ldF3(pos, i);
        const int sx = dom.periodic_x
                           ? ((p.x < dom.min.x + skin) ? 1 : ((p.x > dom.max.x - skin) ? -1 : 0))
                           : 0;
        const int sy = dom.periodic_y
                           ? ((p.y < dom.min.y + skin) ? 1 : ((p.y > dom.max.y - skin) ? -1 : 0))
                           : 0;
        const int sz = dom.periodic_z
                           ? ((p.z < dom.min.z + skin) ? 1 : ((p.z > dom.max.z - skin) ? -1 : 0))
                           : 0;
        const int ax[2] = {0, sx}, ay[2] = {0, sy}, az[2] = {0, sz};
        const int nx = (sx == 0) ? 1 : 2, ny = (sy == 0) ? 1 : 2, nz = (sz == 0) ? 1 : 2;
 
        const F3 pPred = ldF3(posPred, i);
        for (int a = 0; a < nx; ++a)
          for (int b = 0; b < ny; ++b)
            for (int c = 0; c < nz; ++c) {
              const int ix = ax[a], iy = ay[b], iz = az[c];
              if (ix == 0 && iy == 0 && iz == 0)
                continue;
              const int slot = Kokkos::atomic_fetch_add(&topGhost(), 1);
              if (slot >= capacity) {
                Kokkos::atomic_add(&topGhost(), -1);
                continue;
              }
              const F3 shift{ix * dom.size.x, iy * dom.size.y, iz * dom.size.z};
              st3(pos, slot, add3(p, shift));
              st3(posPred, slot, add3(pPred, shift));
              invMass(slot) = invMass(i);
              st3(vel, slot, ldF3(vel, i));
              st3(velPred, slot, ldF3(velPred, i));
              st4(quat, slot, ldF4(quat, i));
              st4(quatPred, slot, ldF4(quatPred, i));
              st3(angVel, slot, ldF3(angVel, i));
              st3(angVelPred, slot, ldF3(angVelPred, i));
              scale(slot) = scale(i);
              shapeId(slot) = shapeId(i);
              realIndices(slot) = i;
              if (gid.extent(0) > 0)
                gid(slot) = gid(i);
              // Ghost slots must carry the owner's material id: the narrowphase reads matId by RAW
              // slot index, so a stale/zero ghost entry silently applied material 0's pair row to
              // every cross-periodic-boundary contact.
              if (materialId.extent(0) > 0)
                materialId(slot) = materialId(i);
            }
      });
  space.fence();
}
 
}  // namespace peclet::dem
 
#endif  // DEM_PERIODICITY_HPP