sleeping.hpp

Go to the documentation of this file.

 
#ifndef DEM_SLEEPING_HPP
#define DEM_SLEEPING_HPP
 
#include <Kokkos_Core.hpp>
 
#include "contact_preprocessing.hpp"  // ManifoldC, CpExec/CpMem
#include "dem_portable.hpp"
 
namespace peclet::dem {
 
inline void freezeAsleepKokkos(int numReal, Kokkos::View<const unsigned char*, CpMem> asleep,
                               Kokkos::View<const float* [3], CpMem> pos,
                               Kokkos::View<const float* [4], CpMem> quat,
                               Kokkos::View<float* [3], CpMem> posPred,
                               Kokkos::View<float* [4], CpMem> quatPred,
                               Kokkos::View<float* [3], CpMem> velPred,
                               Kokkos::View<float* [3], CpMem> angVelPred) {
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::freeze_asleep", Kokkos::RangePolicy<CpExec>(space, 0, numReal),
      KOKKOS_LAMBDA(int i) {
        if (!asleep(i))
          return;
        for (int c = 0; c < 3; ++c) {
          posPred(i, c) = pos(i, c);
          velPred(i, c) = 0.0f;
          angVelPred(i, c) = 0.0f;
        }
        for (int c = 0; c < 4; ++c)
          quatPred(i, c) = quat(i, c);
      });
  space.fence();
}
 
inline void buildInvMassEffKokkos(int numBodies, Kokkos::View<const unsigned char*, CpMem> asleep,
                                  Kokkos::View<const int*, CpMem> realIdx,
                                  Kokkos::View<const float*, CpMem> invMass,
                                  Kokkos::View<float*, CpMem> invMassEff, float sleeperFrac) {
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::inv_mass_eff", Kokkos::RangePolicy<CpExec>(space, 0, numBodies),
      KOKKOS_LAMBDA(int i) {
        invMassEff(i) = asleep(realIdx(i)) ? sleeperFrac * invMass(i) : invMass(i);
      });
  space.fence();
}
 
inline void computeManifoldSleepKokkos(Kokkos::View<const ManifoldC*, CpMem> manifolds,
                                       int numManifolds, Kokkos::View<const int*, CpMem> realIdx,
                                       Kokkos::View<const unsigned char*, CpMem> asleep,
                                       Kokkos::View<unsigned char*, CpMem> manifoldSleep) {
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::manifold_sleep", Kokkos::RangePolicy<CpExec>(space, 0, numManifolds),
      KOKKOS_LAMBDA(int idx) {
        const ManifoldC m = manifolds(idx);
        unsigned char s = 0;
        if (m.num_points > 0) {
          const bool aA = asleep(realIdx(m.bodyA)) != 0;
          const bool aB = (m.bodyB < 0) ? true : (asleep(realIdx(m.bodyB)) != 0);
          s = (aA && aB) ? 1 : 0;
        }
        manifoldSleep(idx) = s;
      });
  space.fence();
}
 
inline void computeContactSleepKokkos(Kokkos::View<const ContactC*, CpMem> contacts,
                                      int numContacts, Kokkos::View<const int*, CpMem> realIdx,
                                      Kokkos::View<const unsigned char*, CpMem> asleep,
                                      Kokkos::View<unsigned char*, CpMem> contactSleep) {
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::contact_sleep", Kokkos::RangePolicy<CpExec>(space, 0, numContacts),
      KOKKOS_LAMBDA(int idx) {
        const ContactC c = contacts(idx);
        const bool aA = asleep(realIdx(c.bodyA)) != 0;
        const bool aB = (c.bodyB < 0) ? true : (asleep(realIdx(c.bodyB)) != 0);
        contactSleep(idx) = (aA && aB) ? 1 : 0;
      });
  space.fence();
}
 
inline void wakeDisturbedKokkos(Kokkos::View<const ManifoldC*, CpMem> manifolds, int numManifolds,
                                Kokkos::View<const int*, CpMem> realIdx,
                                Kokkos::View<const float* [3], CpMem> velPred, float wakeSpeed,
                                Kokkos::View<unsigned char*, CpMem> asleep,
                                Kokkos::View<unsigned char*, CpMem> sleepCounter,
                                Kokkos::View<unsigned char*, CpMem> movingWall,
                                Kokkos::View<int*, CpMem> curCount, int numReal) {
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::wake_reset", Kokkos::RangePolicy<CpExec>(space, 0, numReal),
      KOKKOS_LAMBDA(int i) {
        movingWall(i) = 0;
        curCount(i) = 0;
      });
  const float wakeSpeed2 = wakeSpeed * wakeSpeed;
  Kokkos::parallel_for(
      "peclet::dem::wake_disturbed", Kokkos::RangePolicy<CpExec>(space, 0, numManifolds),
      KOKKOS_LAMBDA(int idx) {
        const ManifoldC m = manifolds(idx);
        if (m.num_points <= 0)
          return;
        const int ra = realIdx(m.bodyA);
        Kokkos::atomic_add(&curCount(ra), 1);
        int rb = -1;
        if (m.bodyB >= 0) {
          rb = realIdx(m.bodyB);
          if (rb == ra)
            return;  // periodic self-pair
          Kokkos::atomic_add(&curCount(rb), 1);
        }
        // Moving wall (boundary manifold with nonzero wall velocity): never sleeps.
        if (m.bodyB < 0) {
          const F4 wv = m.wallVel_sum;
          if (wv.x != 0.0f || wv.y != 0.0f || wv.z != 0.0f) {
            movingWall(ra) = 1;
            if (asleep(ra)) {
              asleep(ra) = 0;
              sleepCounter(ra) = 0;
            }
          }
          return;
        }
        const bool sA = asleep(ra) != 0, sB = asleep(rb) != 0;
        if (sA == sB)
          return;  // both awake or both asleep: no cross-wake here
        // Exactly one asleep: wake it if the awake side moves faster than wakeSpeed.
        const int awake = sA ? rb : ra, sleeper = sA ? ra : rb;
        const float vx = velPred(awake, 0), vy = velPred(awake, 1), vz = velPred(awake, 2);
        if (vx * vx + vy * vy + vz * vz > wakeSpeed2) {
          asleep(sleeper) = 0;
          sleepCounter(sleeper) = 0;
        }
      });
  space.fence();
}
 
inline void wakeContactChangeKokkos(int numReal, Kokkos::View<const int*, CpMem> curCount,
                                    Kokkos::View<int*, CpMem> prevCount,
                                    Kokkos::View<const unsigned char*, CpMem> movingWall,
                                    Kokkos::View<unsigned char*, CpMem> asleep,
                                    Kokkos::View<unsigned char*, CpMem> sleepCounter,
                                    bool wakeOnChange) {
  CpExec space;
  Kokkos::parallel_for(
      "peclet::dem::wake_contact_change", Kokkos::RangePolicy<CpExec>(space, 0, numReal),
      KOKKOS_LAMBDA(int i) {
        // Rule (b): a change in the live contact count. In a near-threshold settled bed contacts
        // FLICKER in and out every substep, so waking on any change stops the bed ever freezing;
        // waking only on a LOST contact (support removed), and only when enabled, is the robust
        // form. A gained contact is extra support and never a reason to wake.
        if (wakeOnChange && asleep(i) && curCount(i) < prevCount(i)) {
          asleep(i) = 0;
          sleepCounter(i) = 0;
        }
        if (movingWall(i))
          sleepCounter(i) = 0;
        prevCount(i) = curCount(i);
      });
  space.fence();
}
 
inline void updateSleepKokkos(int numReal, Kokkos::View<const float* [3], CpMem> vel,
                              Kokkos::View<const float* [3], CpMem> angVel,
                              Kokkos::View<const float*, CpMem> rad,
                              Kokkos::View<const unsigned char*, CpMem> grounded,
                              Kokkos::View<const unsigned char*, CpMem> movingWall,
                              Kokkos::View<unsigned char*, CpMem> asleep,
                              Kokkos::View<unsigned char*, CpMem> sleepCounter, float sleepSpeed,
                              int K, Kokkos::View<float* [3], CpMem> velOut,
                              Kokkos::View<float* [3], CpMem> angVelOut) {
  CpExec space;
  const float s2 = sleepSpeed * sleepSpeed;
  Kokkos::parallel_for(
      "peclet::dem::update_sleep", Kokkos::RangePolicy<CpExec>(space, 0, numReal),
      KOKKOS_LAMBDA(int i) {
        if (asleep(i))
          return;  // stays asleep until an explicit wake
        const float vx = vel(i, 0), vy = vel(i, 1), vz = vel(i, 2);
        const float wx = angVel(i, 0), wy = angVel(i, 1), wz = angVel(i, 2);
        const float r = rad(i);
        const bool lowLin = (vx * vx + vy * vy + vz * vz) < s2;
        const bool lowAng = (wx * wx + wy * wy + wz * wz) * r * r < s2;
        if (lowLin && lowAng && grounded(i) > 0 && !movingWall(i)) {
          int c = static_cast<int>(sleepCounter(i)) + 1;
          if (c >= K) {
            asleep(i) = 1;
            for (int k = 0; k < 3; ++k) {
              velOut(i, k) = 0.0f;
              angVelOut(i, k) = 0.0f;
            }
            c = K;  // saturate the counter
          }
          sleepCounter(i) = static_cast<unsigned char>(c > 255 ? 255 : c);
        } else {
          sleepCounter(i) = 0;
        }
      });
  space.fence();
}
 
}  // namespace peclet::dem
 
#endif  // DEM_SLEEPING_HPP