sim.hpp

Go to the documentation of this file.

 
#ifndef DEM_SIM_HPP
#define DEM_SIM_HPP
 
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <Kokkos_Core.hpp>
#include <map>
#include <memory>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
#include "broadphase_arborx.hpp"
#include "contact_preprocessing.hpp"
#include "integration.hpp"
#include "io.hpp"
#include "narrowphase.hpp"
#include "output_sdf.hpp"
#include "particles.hpp"
#include "peclet/core/common/view.hpp"  // peclet::core::toVector — single-copy device View -> host std::vector (S2a)
#include "peclet/core/geom/scene_builder.hpp"
#include "periodicity.hpp"
#include "shape_registry.hpp"  // ShapeRegistry (the base class: shapes, shells, P_)
#include "shapes_portable.hpp"
#include "sleeping.hpp"            // island sleeping / freezing (single-GPU statics)
#include "solve_driver.hpp"        // demSolveContacts + SoloSolveHooks + readInt/readFloat
#include "solve_driver_force.hpp"  // demStepForce + HertzMindlinLaw + demStepHertz
#include "solver_friction.hpp"
#include "solver_hertz.hpp"
#include "solver_multilevel.hpp"
#include "solver_position.hpp"
#include "solver_velocity.hpp"
#include "step_solve.hpp"  // demStep / computeOverlapsKokkos / growContactBuffers
 
#ifdef PECLET_DEM_MPI
#include "mpi_halo.hpp"        // ParticleHalo (gated; default module never includes it)
#include "step_solve_mpi.hpp"  // demStepMpi / demStepHertzMpi + the MPI hook policies
#endif
 
namespace peclet::dem {
 
class Simulation : public ShapeRegistry {
 public:
  explicit Simulation(int capacity) {
    registry().push_back(this);
    P_.allocate(capacity, capacity * 64, capacity * 16, /*shapes*/ 1, /*shell*/ 1, /*planes*/ 8);
    // default sphere shape (radius 1) + identity-ish defaults
    setSphereShape(1.0f);
  }
  ~Simulation() {
    auto& r = registry();
    r.erase(std::remove(r.begin(), r.end(), this), r.end());
  }
 
  // Teardown safety: the Particles SoA holds Kokkos Views, so they MUST be freed before
  // Kokkos::finalize (else "deallocated after finalize" aborts). releaseAll() (called from the
  // module's atexit, before finalize) frees every live Sim's Views, so callers need not `del sim;
  // gc.collect()` themselves.
  void releaseViews() {
    P_ = Particles{};
#ifdef PECLET_DEM_MPI
    halo_.reset();  // the halo also owns Kokkos Views (gather/forward buffers + its core
                    // sub-objects') that must be freed before Kokkos::finalize, else "deallocated
                    // after finalize" aborts. Destroying it via the unique_ptr frees them all.
#endif
  }
  static void releaseAll() {
    for (auto* s : registry())
      s->releaseViews();
  }
  static std::vector<Simulation*>& registry() {
    static std::vector<Simulation*> r;
    return r;
  }
 
  int addAnalyticWall(const std::vector<int>& nodeInts, const std::vector<float>& nodeReals,
                      int root, bool invert, float restitution, float friction) {
    namespace g = peclet::core::geom;
    if (nodeInts.size() % g::kNodeIntStride || nodeReals.size() % g::kNodeRealStride)
      throw std::runtime_error("add_analytic_wall: node arrays are not a whole number of records");
    const int n = static_cast<int>(nodeInts.size() / g::kNodeIntStride);
    if (n == 0 || static_cast<int>(nodeReals.size() / g::kNodeRealStride) != n)
      throw std::runtime_error("add_analytic_wall: node int/real counts disagree");
    if (root < 0 || root >= n)
      throw std::runtime_error("add_analytic_wall: root out of range");
    const int base = static_cast<int>(wallNodesHost_.size());
    for (int i = 0; i < n; ++i) {
      g::ShapeNode<float> nd = g::decodeNode<float>(&nodeInts[i * g::kNodeIntStride],
                                                    &nodeReals[i * g::kNodeRealStride]);
      if (nd.kind == g::kGrid)
        throw std::runtime_error("add_analytic_wall: grid leaves are not supported here");
      if (nd.kind >= g::kCsgBase) {  // rebase child indices into the shared pool
        nd.aux0 += base;
        nd.aux1 += base;
      }
      wallNodesHost_.push_back(nd);
    }
    P_.wallNodes = Kokkos::View<g::ShapeNode<float>*, CpMem>("wallNodes", wallNodesHost_.size());
    auto hn = Kokkos::create_mirror_view(P_.wallNodes);
    for (std::size_t i = 0; i < wallNodesHost_.size(); ++i)
      hn(i) = wallNodesHost_[i];
    Kokkos::deep_copy(P_.wallNodes, hn);
 
    WallSdf w{};
    w.shapeRoot = base + root;
    w.nodeCount = static_cast<int>(wallNodesHost_.size());
    // R1: the AUTHORED root transform, kept so setWallTransform composes onto it instead of onto
    // the previous frame's result (which would compound and drift).
    wallRootTf_[base + root] = wallNodesHost_[(std::size_t)(base + root)].transform;
    w.sign = invert ? -1.0f : 1.0f;
    w.restitution = restitution;
    w.friction = friction;
    const int idx = static_cast<int>(wallsHost_.size());
    wallsHost_.push_back(w);
    // Every wall's node pointer must be refreshed: the pool View was just reallocated.
    for (auto& wh : wallsHost_)
      if (wh.shapeRoot >= 0)
        wh.nodes = P_.wallNodes.data();
    P_.numWalls = static_cast<int>(wallsHost_.size());
    uploadWalls();
    ensureContactCapacity();
    return idx;
  }
 
  void setDomain(float lx, float ly, float lz, bool px, bool py, bool pz) {
    P_.domain = Domain{F3{0, 0, 0}, F3{lx, ly, lz}, F3{lx, ly, lz}, px, py, pz};
    P_.skin = 0.1f * P_.globalScale;
  }
  // CUDA Simulation::set_domain(min, max): arbitrary origin; keeps the current periodicity flags.
  void setDomainMinMax(F3 mn, F3 mx) {
    P_.domain = Domain{mn,
                       mx,
                       F3{mx.x - mn.x, mx.y - mn.y, mx.z - mn.z},
                       P_.domain.periodic_x,
                       P_.domain.periodic_y,
                       P_.domain.periodic_z};
    P_.skin = 0.1f * P_.globalScale;
  }
  void enablePeriodicity(bool x, bool y, bool z) {
    P_.domain.periodic_x = x;
    P_.domain.periodic_y = y;
    P_.domain.periodic_z = z;
  }
  // The suite-canonical domain quartet (suite/docs/NAMING.md 1.1): `origin` is the lower corner,
  // `extent` the SIZE, `periodic` the per-axis flags — the same four names flow, voro and the AMR
  // octree use. `setDomainCanonical` is what the keyword form of `set_domain` binds to.
  void setDomainCanonical(F3 extent, F3 origin, bool px, bool py, bool pz) {
    P_.domain = Domain{origin, F3{origin.x + extent.x, origin.y + extent.y, origin.z + extent.z},
                       extent, px,
                       py,     pz};
    P_.skin = 0.1f * P_.globalScale;
  }
  std::tuple<float, float, float> domainOrigin() const {
    return {P_.domain.min.x, P_.domain.min.y, P_.domain.min.z};
  }
  std::tuple<float, float, float> domainExtent() const {
    return {P_.domain.size.x, P_.domain.size.y, P_.domain.size.z};
  }
  std::tuple<bool, bool, bool> domainPeriodic() const {
    return {P_.domain.periodic_x, P_.domain.periodic_y, P_.domain.periodic_z};
  }
  void setGravity(float gx, float gy, float gz) { P_.gravity = F3{gx, gy, gz}; }
  std::tuple<float, float, float> gravity() const {
    return {P_.gravity.x, P_.gravity.y, P_.gravity.z};
  }
  void setThermostat(float temperature, float tau, float kB) {  // Berendsen; tau=0 disables
    P_.thermostatTemp = temperature;
    P_.thermostatTau = tau;
    P_.thermostatKB = kB;
  }
  void setSolverIterations(int pos, int vel) {
    P_.positionIterations = pos;
    P_.velocityIterations = vel;
  }
  // Select the single-GPU collision solves: true (default) = colored Gauss–Seidel for both the
  // restitution and the overlap solve, false = count-averaged Jacobi (legacy). For A/B validation.
  void setVelocityUseGS(bool useGS) { P_.velocityUseGS = useGS; }
  bool velocityUseGS() const { return P_.velocityUseGS; }
  void setStabilizationMode(const std::string& mode) {
    if (mode == "off")
      P_.stabilizationMode = 0;
    else if (mode == "onesided")
      P_.stabilizationMode = 1;
    else if (mode == "multilevel")
      P_.stabilizationMode = 2;
    else if (mode == "escalate")
      P_.stabilizationMode = 3;
    else if (mode == "ordered")
      P_.stabilizationMode = 4;
    else
      throw std::invalid_argument(
          "set_stabilization: expected 'off', 'onesided', 'multilevel', 'escalate' or "
          "'ordered'");
  }
  std::string stabilizationMode() const {
    static const char* const names[] = {"off", "onesided", "multilevel", "escalate", "ordered"};
    return names[P_.stabilizationMode];
  }
  void setRestitutionModel(const std::string& model) {
    if (model == "newton")
      P_.restitutionModel = 0;
    else if (model == "poisson")
      P_.restitutionModel = 1;
    else
      throw std::invalid_argument("set_restitution_model: expected 'newton' or 'poisson'");
  }
  std::string restitutionModel() const { return P_.restitutionModel == 1 ? "poisson" : "newton"; }
  void setSleeping(bool enabled, float threshold_scale = 2.0f, int consecutive = 64,
                   float wake_scale = 40.0f, bool wake_on_lost_contact = false,
                   float immovable_frac = 0.01f) {
    P_.sleepingEnabled = enabled;
    if (threshold_scale > 0.0f)
      P_.sleepScale = threshold_scale;
    if (consecutive > 0)
      P_.sleepK = consecutive;
    if (wake_scale > 0.0f)
      P_.wakeScale = wake_scale;  // hysteresis: wake only well above the residual settling jitter
    P_.sleepWakeLostContact = wake_on_lost_contact;
    if (immovable_frac >= 0.0f)
      P_.sleepImmovableFrac = immovable_frac;
  }
  bool sleeping() const { return P_.sleepingEnabled; }
  void setVerletSkin(float skin_frac) {
    P_.verletSkinFrac = skin_frac;
    P_.impNumPairs = -1;  // invalidate the cache
  }
  float verletSkin() const { return P_.verletSkinFrac; }
  void setCudaGraphs(bool enabled) { P_.cudaGraphs = enabled; }
  bool cudaGraphs() const { return P_.cudaGraphs; }
  void setFusedSweeps(const std::string& mode) {
    if (mode == "auto")
      P_.fusedSweeps = -1;
    else if (mode == "off")
      P_.fusedSweeps = 0;
    else if (mode == "on")
      P_.fusedSweeps = 1;
    else
      throw std::invalid_argument("set_fused_sweeps: expected 'auto', 'on' or 'off'");
  }
  std::string fusedSweeps() const {
    return P_.fusedSweeps < 0 ? "auto" : (P_.fusedSweeps == 1 ? "on" : "off");
  }
  void setIncrementalColoring(bool enabled) { P_.incrementalColoring = enabled; }
  bool incrementalColoring() const { return P_.incrementalColoring; }
  void setShellDetection(const std::string& mode) {
    if (mode == "one_way")
      P_.shellDetection = 0;
    else if (mode == "two_way")
      P_.shellDetection = 1;
    else
      throw std::invalid_argument("set_shell_detection: expected 'one_way' or 'two_way'");
    ensureContactCapacity();  // the setup-time bound doubles for two-way
  }
  std::string shellDetection() const { return P_.shellDetection == 1 ? "two_way" : "one_way"; }
  int numAsleep() {
    int n = 0;
    auto a = P_.asleep;
    Kokkos::parallel_reduce(
        "peclet::dem::count_asleep", Kokkos::RangePolicy<CpExec>(0, P_.numReal),
        KOKKOS_LAMBDA(int i, int& acc) { acc += a(i) ? 1 : 0; }, n);
    return n;
  }
  void setHertzMaterial(int mat, float youngs, float poisson) {
    if (mat < 0 || mat >= 8)
      throw std::invalid_argument("material id out of range");
    auto he = Kokkos::create_mirror_view(P_.hertzE);
    auto hn = Kokkos::create_mirror_view(P_.hertzNu);
    Kokkos::deep_copy(he, P_.hertzE);
    Kokkos::deep_copy(hn, P_.hertzNu);
    he(mat) = youngs;
    hn(mat) = poisson;
    Kokkos::deep_copy(P_.hertzE, he);
    Kokkos::deep_copy(P_.hertzNu, hn);
  }
  void stepHertz(int substeps, float skin_frac) {
    requireDt("step_hertz");
    demStepHertz(P_, P_.dt, substeps, skin_frac);
  }
  void setGlobalScale(float s) {
    P_.globalScale = s;
    P_.skin = 0.1f * s;
  }
  float globalScale() const { return P_.globalScale; }
  void setDt(float dt) {
    if (!(dt > 0.0f))
      throw std::invalid_argument("set_dt: dt must be > 0 (relax() is the dynamics-free substep)");
    P_.dt = dt;
    dtSet_ = true;
  }
  float dt() const { return P_.dt; }
  int capacity() const { return P_.capacity; }
  // (restitution_normal, restitution_tangent, friction) to match CUDA set_material_params; the
  // Kokkos pipeline currently carries normal restitution + dynamic friction (tangential restitution
  // unused).
  // NOTE (found via peclet-examples/stirred-column, 2026-08-30): the DEFAULT body-body material is
  // frictionless. add_analytic_wall / add_sdf_wall set the particle-WALL material only, so a bed
  // more than a few layers deep with the default body-body friction behaves like a liquid -- it
  // transmits full hydrostatic pressure to the container and the position solve squeezes grains
  // through the boundary. The failure is silent and looks exactly like a solver-convergence bug
  // (raising the position iterations 4,4 -> 24,12 and halving dt changed the measured leakage not
  // at all). Set a non-zero friction here for any deep bed.
  void setMaterialParams(float restitution_normal, float restitution_tangent, float friction) {
    P_.restitutionNormal = restitution_normal;
    P_.frictionDynamic = friction;
    P_.restitutionTangent = restitution_tangent;  // Walton beta (impact law, cone-clamped)
  }
  void setMaterialIds(const std::vector<int>& ids) {
    auto h = Kokkos::create_mirror_view(P_.materialId);
    Kokkos::deep_copy(h, P_.materialId);
    for (int i = 0; i < P_.numReal && i < (int)ids.size(); ++i)
      h(i) = static_cast<unsigned char>(ids[i]);
    Kokkos::deep_copy(P_.materialId, h);
  }
  void setPairMaterial(int a, int b, float restitution, float friction) {
    if (a < 0 || b < 0 || a >= kMaxMaterials || b >= kMaxMaterials)
      throw std::invalid_argument("material id out of range");
    if (P_.pairMaterials.extent(0) == 0) {
      P_.pairMaterials =
          Kokkos::View<float*, CpMem>("pairMaterials", kMaxMaterials * kMaxMaterials * 2);
      auto h0 = Kokkos::create_mirror_view(P_.pairMaterials);
      for (int i = 0; i < kMaxMaterials * kMaxMaterials; ++i) {
        h0(2 * i) = P_.restitutionNormal;
        h0(2 * i + 1) = P_.frictionDynamic;
      }
      Kokkos::deep_copy(P_.pairMaterials, h0);
    }
    auto h = Kokkos::create_mirror_view(P_.pairMaterials);
    Kokkos::deep_copy(h, P_.pairMaterials);
    for (auto xy : {std::pair<int, int>{a, b}, std::pair<int, int>{b, a}}) {
      h((xy.first * kMaxMaterials + xy.second) * 2) = restitution;
      h((xy.first * kMaxMaterials + xy.second) * 2 + 1) = friction;
    }
    Kokkos::deep_copy(P_.pairMaterials, h);
  }
  void setWallMaterialId(int wid, int mat) {
    if (wid < 0 || wid >= P_.numWalls)
      throw std::invalid_argument("wall index out of range");
    auto h = Kokkos::create_mirror_view(P_.walls);
    Kokkos::deep_copy(h, P_.walls);
    h(wid).materialId = mat;
    Kokkos::deep_copy(P_.walls, h);
  }
  void addPlane(float px, float py, float pz, float nx, float ny, float nz) {
    auto h = Kokkos::create_mirror_view(P_.planes);
    Kokkos::deep_copy(h, P_.planes);
    if (P_.numPlanes < static_cast<int>(P_.planes.extent(0)))
      h(P_.numPlanes++) = PlaneP{F3{px, py, pz}, F3{nx, ny, nz}};
    Kokkos::deep_copy(P_.planes, h);
  }
 
  // Add a static, world-space SDF wall/container the grains collide against (a drum barrel, hopper,
  // vibrating tray). `grid` is a flat [nx*ny*nz] signed-distance field, x-fastest (idx = x + y*nx +
  // z*nx*ny), sampled at world nodes origin + (x,y,z)*spacing — POSITIVE in the void where the
  // grains live, NEGATIVE inside the solid wall (so a grain surface point reads the penetration
  // depth and the outward gradient is the push-out normal). restitution/friction are the binary
  // particle–wall material (independent of the body-body material). The wall is motionless but
  // carries a rigid-body surface-velocity field (set via setWallVelocity) so a grain touching it
  // feels the wall's motion. Returns the wall's index (for setWallVelocity). Add walls before
  // stepping.
  int addSdfWall(const std::vector<float>& grid, int nx, int ny, int nz, F3 origin, F3 spacing,
                 float restitution, float friction) {
    if (nx < 2 || ny < 2 || nz < 2)
      throw std::runtime_error("addSdfWall: grid dims must be >= 2 on each axis");
    if (static_cast<long>(nx) * ny * nz != static_cast<long>(grid.size()))
      throw std::runtime_error("addSdfWall: grid.size() must equal nx*ny*nz");
 
    WallSdf w{};
    w.grid.nx = nx;
    w.grid.ny = ny;
    w.grid.nz = nz;
    w.grid.offset = static_cast<int>(wallGridHost_.size());
    w.grid.origin = toCoreVec(origin);
    w.grid.invSpacing = peclet::core::Vec3<float>{spacing.x > 0 ? 1.0f / spacing.x : 0.0f,
                                                  spacing.y > 0 ? 1.0f / spacing.y : 0.0f,
                                                  spacing.z > 0 ? 1.0f / spacing.z : 0.0f};
    // A CONTAINER: beyond the stored box is wall-side, so the off-grid residual is SUBTRACTED.
    w.grid.extension = peclet::core::geom::GridExtension::kContainer;
    w.restitution = restitution;
    w.friction = friction;
    const int idx = static_cast<int>(wallsHost_.size());
    wallsHost_.push_back(w);
    wallGridHost_.insert(wallGridHost_.end(), grid.begin(), grid.end());
    // Upload the concatenated grid samples (only changes when a wall is added, not per step).
    P_.wallGrid = Kokkos::View<float*, CpMem>("wallGrid", wallGridHost_.size());
    auto hg = Kokkos::create_mirror_view(P_.wallGrid);
    for (size_t i = 0; i < wallGridHost_.size(); ++i)
      hg(i) = wallGridHost_[i];
    Kokkos::deep_copy(P_.wallGrid, hg);
    uploadWalls();
    P_.numWalls = static_cast<int>(wallsHost_.size());
    if (friction > P_.wallFrictionMax)
      P_.wallFrictionMax = friction;
    ensureContactCapacity();
    return idx;
  }
 
  // Set a wall's rigid-body surface-velocity field v(x) = linVel + angVel × (x − center). A grain
  // in contact feels this velocity even though the geometry never moves: set angVel for a rotating
  // drum (about `center` on the axis), or drive linVel sinusoidally each step for a vibrating wall.
  // Cheap (a few host scalars); safe to call every step.
  void setWallVelocity(int wallIndex, F3 linVel, F3 angVel, F3 center) {
    if (wallIndex < 0 || wallIndex >= static_cast<int>(wallsHost_.size()))
      throw std::runtime_error("setWallVelocity: wall index out of range");
    wallsHost_[wallIndex].linVel = linVel;
    wallsHost_[wallIndex].angVel = angVel;
    wallsHost_[wallIndex].center = center;
    uploadWalls();
  }
 
  // R1: RIGID-BODY PLACEMENT of an analytic wall -- rotate/translate the GEOMETRY, not just its
  // surface-velocity field. setWallVelocity alone is enough for an axisymmetric wall (a drum
  // barrel looks the same at every angle), but a stirrer blade has to actually move. The world
  // transform W is composed ONTO THE AUTHORED root transform (kept from add_analytic_wall), so
  // repeated calls place the wall absolutely and never compound; the KB-sized node pool is
  // re-uploaded, which is the whole cost.
  //
  // W is the placement of the authored frame in the world: a probe point p is evaluated as
  // eval_authored(toCanonical(W, p)). Drive it together with setWallVelocity (v = linVel +
  // angVel x (p - center)) so the surface velocity a grain feels matches the geometry it touches;
  // nothing here infers one from the other.
  void setWallTransform(int wallIndex, F3 translation, float qx, float qy, float qz, float qw) {
    namespace g = peclet::core::geom;
    if (wallIndex < 0 || wallIndex >= static_cast<int>(wallsHost_.size()))
      throw std::runtime_error("setWallTransform: wall index out of range");
    const int root = wallsHost_[(std::size_t)wallIndex].shapeRoot;
    if (root < 0)
      throw std::runtime_error(
          "setWallTransform: this wall is a sampled GRID SDF, which has no shape tree to place -- "
          "only analytic walls (add_analytic_wall) can be moved");
    const auto it = wallRootTf_.find(root);
    if (it == wallRootTf_.end())
      throw std::runtime_error("setWallTransform: authored root transform missing (internal)");
    g::Transform<float> W;
    W.translation = peclet::core::Vec3<float>{translation.x, translation.y, translation.z};
    W.rotation = peclet::core::Quat<float>{qx, qy, qz, qw};
    W.scale = 1.0f;
    wallNodesHost_[(std::size_t)root].transform =
        g::SceneBuilder<float>::composeTransform(W, it->second);
    auto hn = Kokkos::create_mirror_view(P_.wallNodes);
    for (std::size_t i = 0; i < wallNodesHost_.size(); ++i)
      hn(i) = wallNodesHost_[i];
    Kokkos::deep_copy(P_.wallNodes, hn);
  }
 
  // Diagnostic/authoring probe: the wall's own SDF at world points (flat [m*3]), sign included --
  // POSITIVE in the void where the grains live. This is exactly what the narrow phase reads
  // (sampleWallSdf), so it is the honest way to check a placement or draw a stirrer.
  std::vector<float> wallSdfAt(int wallIndex, const std::vector<float>& pts) const {
    if (wallIndex < 0 || wallIndex >= static_cast<int>(wallsHost_.size()))
      throw std::runtime_error("wall_sdf_at: wall index out of range");
    const int m = static_cast<int>(pts.size() / 3);
    std::vector<float> out((std::size_t)m, 0.0f);
    const WallSdf w = wallsHost_[(std::size_t)wallIndex];
    Kokkos::View<float*, CpMem> d("wallProbe", std::max(1, m));
    Kokkos::View<float*, CpMem> q("wallProbePts", std::max(1, 3 * m));
    auto hq = Kokkos::create_mirror_view(q);
    for (int i = 0; i < 3 * m; ++i)
      hq(i) = pts[(std::size_t)i];
    Kokkos::deep_copy(q, hq);
    auto grid = P_.wallGrid;
    Kokkos::parallel_for(
        "peclet::dem::wall_sdf_probe", Kokkos::RangePolicy<CpExec>(0, m), KOKKOS_LAMBDA(int i) {
          d(i) = sampleWallSdf(F3{q(3 * i), q(3 * i + 1), q(3 * i + 2)}, w, grid);
        });
    Kokkos::fence();
    auto hd = Kokkos::create_mirror_view(d);
    Kokkos::deep_copy(hd, d);
    for (int i = 0; i < m; ++i)
      out[(std::size_t)i] = hd(i);
    return out;
  }
 
  // positions: flat [n*3]; (re)sets the real-particle count and default state.
  void setPositions(const std::vector<float>& xyz) {
    const int n = static_cast<int>(xyz.size() / 3);
    if (n > P_.capacity)
      throw std::invalid_argument(
          "set_positions: " + std::to_string(n) +
          " particles exceed Simulation(capacity=" + std::to_string(P_.capacity) + ")");
    P_.numReal = n;
    P_.numParticles = n;
    P_.hertzNumPairs = -1;  // particle indices changed: invalidate the hertz pair cache
    P_.hertzPrevCount = 0;
    P_.prevPairCount = 0;  // stale persistent-pair ledger must not warm-start the new set
#ifdef PECLET_DEM_MPI
    mpiGidsGlobal_ = false;  // new particle set -> re-base the global ids at the next stepMpi
#endif
    auto pos = Kokkos::create_mirror_view(P_.pos);
    auto q = Kokkos::create_mirror_view(P_.quat);
    auto im = Kokkos::create_mirror_view(P_.invMass);
    auto sc = Kokkos::create_mirror_view(P_.scale);
    auto ii = Kokkos::create_mirror_view(P_.invInertia);
    auto sid = Kokkos::create_mirror_view(P_.shapeId);
    auto vel = Kokkos::create_mirror_view(P_.vel);
    auto av = Kokkos::create_mirror_view(P_.angVel);
    auto gi = Kokkos::create_mirror_view(P_.gid);
    for (int i = 0; i < n; ++i) {
      gi(i) = i;  // identity global id; the MPI enable re-bases it to a global Exscan offset
      pos(i, 0) = xyz[3 * i];
      pos(i, 1) = xyz[3 * i + 1];
      pos(i, 2) = xyz[3 * i + 2];
      q(i, 0) = 0;
      q(i, 1) = 0;
      q(i, 2) = 0;
      q(i, 3) = 1;
      im(i) = 1.0f;
      sc(i) = 1.0f;
      sid(i) = 0;
      ii(i, 0) = defaultInvI_.x;
      ii(i, 1) = defaultInvI_.y;
      ii(i, 2) = defaultInvI_.z;
      vel(i, 0) = vel(i, 1) = vel(i, 2) = 0;
      av(i, 0) = av(i, 1) = av(i, 2) = 0;
    }
    Kokkos::deep_copy(P_.pos, pos);
    Kokkos::deep_copy(P_.quat, q);
    Kokkos::deep_copy(P_.invMass, im);
    Kokkos::deep_copy(P_.scale, sc);
    Kokkos::deep_copy(P_.invInertia, ii);
    Kokkos::deep_copy(P_.shapeId, sid);
    Kokkos::deep_copy(P_.vel, vel);
    Kokkos::deep_copy(P_.angVel, av);
    Kokkos::deep_copy(P_.gid, gi);
    Kokkos::deep_copy(P_.targetScale, P_.scale);  // unscaled growth target = the set scale
  }
  void setScalesUniform(float s) {
    auto sc = Kokkos::create_mirror_view(P_.scale);
    for (int i = 0; i < P_.numReal; ++i)
      sc(i) = s;
    Kokkos::deep_copy(P_.scale, sc);
    Kokkos::deep_copy(P_.targetScale, P_.scale);
  }
  // per-particle scales (growth target): flat [n]. scale starts at target (growth factor applies in
  // step).
  void setScales(const std::vector<float>& s) {
    auto tsc = Kokkos::create_mirror_view(P_.targetScale);
    for (int i = 0; i < P_.numReal && i < (int)s.size(); ++i)
      tsc(i) = s[i];
    Kokkos::deep_copy(P_.targetScale, tsc);
    Kokkos::deep_copy(P_.scale, P_.targetScale);
  }
  void setVelocities(const std::vector<float>& v) {
    auto vel = Kokkos::create_mirror_view(P_.vel);
    for (int i = 0; i < P_.numReal && 3 * i + 2 < (int)v.size(); ++i) {
      vel(i, 0) = v[3 * i];
      vel(i, 1) = v[3 * i + 1];
      vel(i, 2) = v[3 * i + 2];
    }
    Kokkos::deep_copy(P_.vel, vel);
  }
  // Per-particle external FORCE (fluid drag etc.), an (N,3) flat array. Applied in the next
  // step()'s velocity predict as dv = F*invMass*dt. Persists across steps until re-set or cleared.
  void setExternalForces(const std::vector<float>& f) {
    auto ef = Kokkos::create_mirror_view(P_.extForce);
    for (int i = 0; i < P_.numReal && 3 * i + 2 < (int)f.size(); ++i) {
      ef(i, 0) = f[3 * i];
      ef(i, 1) = f[3 * i + 1];
      ef(i, 2) = f[3 * i + 2];
    }
    Kokkos::deep_copy(P_.extForce, ef);
  }
  void clearExternalForces() { Kokkos::deep_copy(P_.extForce, 0.0f); }
  // Per-particle external TORQUE in the WORLD frame, an (N,3) flat array. Applied in the next
  // step()'s angular predict as the body-frame Euler update dw = invI*(tau_body - w x I w)*dt.
  // Persists across steps until re-set or cleared. Only bodies with a finite inertia (invInertia
  // > 0, i.e. a registered non-degenerate shape) respond -- a torque on a point mass is silently
  // inert, exactly as the gyroscopic term is.
  void setExternalTorques(const std::vector<float>& t) {
    auto et = Kokkos::create_mirror_view(P_.extTorque);
    for (int i = 0; i < P_.numReal && 3 * i + 2 < (int)t.size(); ++i) {
      et(i, 0) = t[3 * i];
      et(i, 1) = t[3 * i + 1];
      et(i, 2) = t[3 * i + 2];
    }
    Kokkos::deep_copy(P_.extTorque, et);
    P_.extTorqueActive = true;
  }
  void clearExternalTorques() {
    Kokkos::deep_copy(P_.extTorque, 0.0f);
    P_.extTorqueActive = false;
  }
  const V3& externalTorquesView() const { return P_.extTorque; }
  const V3& externalForcesView() const { return P_.extForce; }
  const Vf& invMassView() const { return P_.invMass; }
  // rigid-body rotation state (the pipeline integrates the gyroscopic Euler term + quaternion
  // already)
  void setQuaternions(const std::vector<float>& q) {
    auto h = Kokkos::create_mirror_view(P_.quat);
    for (int i = 0; i < P_.numReal && 4 * i + 3 < (int)q.size(); ++i) {
      h(i, 0) = q[4 * i];
      h(i, 1) = q[4 * i + 1];
      h(i, 2) = q[4 * i + 2];
      h(i, 3) = q[4 * i + 3];
    }
    Kokkos::deep_copy(P_.quat, h);
  }
  void setAngularVelocities(const std::vector<float>& w) {
    auto h = Kokkos::create_mirror_view(P_.angVel);
    for (int i = 0; i < P_.numReal && 3 * i + 2 < (int)w.size(); ++i) {
      h(i, 0) = w[3 * i];
      h(i, 1) = w[3 * i + 1];
      h(i, 2) = w[3 * i + 2];
    }
    Kokkos::deep_copy(P_.angVel, h);
  }
  void setInvInertia(const std::vector<float>& ii) {
    auto h = Kokkos::create_mirror_view(P_.invInertia);
    for (int i = 0; i < P_.numReal && 3 * i + 2 < (int)ii.size(); ++i) {
      h(i, 0) = ii[3 * i];
      h(i, 1) = ii[3 * i + 1];
      h(i, 2) = ii[3 * i + 2];
    }
    Kokkos::deep_copy(P_.invInertia, h);
  }
  void setInvMass(const std::vector<float>& im) {
    auto h = Kokkos::create_mirror_view(P_.invMass);
    for (int i = 0; i < P_.numReal && i < (int)im.size(); ++i)
      h(i) = im[i];
    Kokkos::deep_copy(P_.invMass, h);
  }
  std::vector<float> getAngularVelocities() const {
    return peclet::core::toVector(
        Kokkos::subview(P_.angVel, Kokkos::make_pair(0, P_.numReal), Kokkos::ALL));
  }
  std::vector<float> getInvInertia() const {
    return peclet::core::toVector(
        Kokkos::subview(P_.invInertia, Kokkos::make_pair(0, P_.numReal), Kokkos::ALL));
  }
  // growth: factor *= exp(rate*dt) per step (capped at 1); new_factor<0 keeps/initialises (0.01 if
  // inactive).
  void setGrowthParams(float rate, float new_factor) {
    if (P_.growthFactor == -1.0f)
      P_.growthFactor = (new_factor > 0.0f) ? new_factor : 0.01f;
    else if (new_factor > 0.0f)
      P_.growthFactor = new_factor;
    P_.growthRate = rate;
    if (P_.growthFactor > 0.0f)
      updateGrowthScalesKokkos(P_.numReal, P_.scale, P_.targetScale, P_.growthFactor);
  }
  float growthFactor() const { return P_.growthFactor; }
  float getGrowthRate() const { return P_.growthRate; }
  // per-particle mass = 1/invMass (0 for fixed/infinite-mass particles), CUDA
  // Simulation::get_masses.
  std::vector<float> getMasses() const {
    auto im = Kokkos::create_mirror_view(P_.invMass);
    Kokkos::deep_copy(im, P_.invMass);
    std::vector<float> out(P_.numReal);
    for (int i = 0; i < P_.numReal; ++i)
      out[i] = (im(i) > 0.0f) ? (1.0f / im(i)) : 0.0f;
    return out;
  }
 
  std::vector<float> getPositions() const {
    return peclet::core::toVector(
        Kokkos::subview(P_.pos, Kokkos::make_pair(0, P_.numReal), Kokkos::ALL));
  }
  std::vector<float> getVelocities() const {
    return peclet::core::toVector(
        Kokkos::subview(P_.vel, Kokkos::make_pair(0, P_.numReal), Kokkos::ALL));
  }
  std::vector<float> getQuaternions() const {
    return peclet::core::toVector(
        Kokkos::subview(P_.quat, Kokkos::make_pair(0, P_.numReal), Kokkos::ALL));
  }
  std::vector<float> getScales() const {
    return peclet::core::toVector(Kokkos::subview(P_.scale, Kokkos::make_pair(0, P_.numReal)));
  }
 
  void step(int n = 1) {
    requireDt("step");
    P_.splitStats.mlHubAggregated = 0;  // split_stats: max over this call's substeps
    for (int i = 0; i < n; ++i)
      demStep(P_);
  }
  void relax(int n = 1) {
    const float dt = P_.dt;
    P_.dt = 0.0f;
    P_.splitStats.mlHubAggregated = 0;
    for (int i = 0; i < n; ++i)
      demStep(P_);
    P_.dt = dt;
  }
 
  // Max pair interpenetration on the current committed state (CUDA Simulation::compute_overlaps).
  float computeOverlaps() { return computeOverlapsKokkos(P_); }
 
  // LAMMPS "dump custom" of the current committed state (CUDA Simulation::export_lammps). Radius =
  // scale*globalScale*baseRadius; bounds computed from the particle AABBs.
  void exportLammps(const std::string& filename, int step) const {
    const std::vector<float> pos = getPositions(), vel = getVelocities(), quat = getQuaternions();
    auto sc = Kokkos::create_mirror_view(P_.scale);
    Kokkos::deep_copy(sc, P_.scale);
    std::vector<float> radii(P_.numReal);
    for (int i = 0; i < P_.numReal; ++i)
      radii[i] = sc(i) * P_.globalScale * baseRadius_;
    const bool pbc = P_.domain.periodic_x || P_.domain.periodic_y || P_.domain.periodic_z;
    peclet::dem::writeLammpsDump(filename, step, pos, vel, quat, radii, nullptr, nullptr, pbc);
  }
 
  // SDF field over the domain -> ImageData VTI (CUDA Simulation::export_sdf).
  void exportSdf(const std::string& filename, int rx, int ry, int rz) {
    const std::vector<float> grid = getSdfGrid(rx, ry, rz);
    const float mn[3] = {P_.domain.min.x, P_.domain.min.y, P_.domain.min.z};
    const float mx[3] = {P_.domain.max.x, P_.domain.max.y, P_.domain.max.z};
    peclet::dem::writeSdfVti(filename, grid, rx, ry, rz, mn, mx);
  }
 
#ifdef PECLET_DEM_MPI
  // Block decomposition over the GLOBAL domain (once); the per-block solver stays non-periodic, the
  // halo supplies the periodic wrap. gsize is the ORB cell grid. Mirror of Simulation::mpi_init.
  void initMpi(std::tuple<double, double, double> origin, std::tuple<double, double, double> size,
               std::tuple<long, long, long> gsize, std::tuple<bool, bool, bool> periodic,
               MPI_Comm comm) {
    halo_->initMpi({std::get<0>(origin), std::get<1>(origin), std::get<2>(origin)},
                   {std::get<0>(size), std::get<1>(size), std::get<2>(size)},
                   {std::get<0>(gsize), std::get<1>(gsize), std::get<2>(gsize)},
                   {std::get<0>(periodic), std::get<1>(periodic), std::get<2>(periodic)}, comm);
  }
  // Enable the distributed step. rcut is a LOWER bound on the ghost-band width: the XPBD step
  // uses max(rcut, 2.1 R_max) with R_max the GLOBAL maximum radius, re-evaluated every step
  // (growth included -- the same band as the single-rank periodic ghosts), so rcut <= 0 means
  // exactly the contact reach; sync_every: owner/ghost reconciliation interval; conservation is
  // exact at any value, larger = more lag, fewer messages. rebalance_every: re-decompose by
  // particle count + migrate ownership every
  // N distributed steps to keep the per-rank load even as a packing densifies (0 = never; the
  // partition is then fixed at the initial decomposition, as before). A pure redistribution — the
  // physics result is unchanged.
  void enableMpiStep(double rcut, int sync_every = 1, bool forward_rotation = true,
                     int rebalance_every = 0, double verlet_skin = 0.0) {
    mpiRcut_ = rcut;
    mpiSyncEvery_ = sync_every < 1 ? 1 : sync_every;
    mpiForwardRotation_ = forward_rotation;
    mpiRebalanceEvery_ = rebalance_every < 0 ? 0 : rebalance_every;
    // Verlet-skin ghost reuse (D2): rebuild the halo topology only when a particle has moved >
    // skin, instead of every substep. 0 (default) keeps the exact per-substep rebuild.
    halo_->setVerletSkin(static_cast<float>(verlet_skin));
  }
  // Halo rebuild stats (D2): topology rebuilds vs total gather() calls (for benchmarking).
  long mpiRebuilds() const { return halo_->numRebuilds(); }
  long mpiGathers() const { return halo_->numGathers(); }
  // Migrate ownership now so each rank holds a near-equal particle count. Safe to call at a step
  // boundary; returns this rank's new owned count. Exposed for manual / adaptive balancing.
  int rebalance() { return halo_->rebalance(P_); }
  // Co-rebalance: migrate ownership onto the weighted ORB of per-cell weights `w` (the SAME
  // partition the coupled flow solver redistributes onto from the same weight field), its split
  // planes on multiples of `align` cells (1 = unaligned, the pre-`align` behaviour bit for bit;
  // a coupled run passes what flow's rebalance_by_weights returned). Returns new owned count.
  // Validated here because core's aligned weighted init checks its contract by assert only.
  int migrateToWeights(const std::vector<peclet::core::Real>& w, int align = 1) {
    if (!halo_->inited())
      throw std::runtime_error("migrate_to_weights: call init_mpi first");
    const auto& G = halo_->decomposer().globalSize();
    int size = 1;
    MPI_Comm_size(halo_->comm(), &size);
    if (align < 1 || (align & (align - 1)) != 0)
      throw std::invalid_argument("migrate_to_weights: align must be a power of two >= 1, got " +
                                  std::to_string(align));
    double coarseCells = 1.0;
    std::size_t cells = 1;
    for (int k = 0; k < 3; ++k) {
      if (G[k] % align != 0)
        throw std::invalid_argument("migrate_to_weights: align " + std::to_string(align) +
                                    " does not divide the ORB grid (" + std::to_string(G[0]) +
                                    ", " + std::to_string(G[1]) + ", " + std::to_string(G[2]) +
                                    ")");
      coarseCells *= static_cast<double>(G[k] / align);
      cells *= static_cast<std::size_t>(G[k]);
    }
    if (coarseCells < static_cast<double>(size))
      throw std::invalid_argument("migrate_to_weights: align " + std::to_string(align) +
                                  " leaves fewer aligned boxes than ranks");
    if (w.size() != cells)
      throw std::invalid_argument("migrate_to_weights: weights has " + std::to_string(w.size()) +
                                  " entries, the ORB grid has " + std::to_string(cells) + " cells");
    return halo_->migrateToWeights(P_, w, align);
  }
  // The current block decomposition (C++ only: tests compare it with the partition flow builds).
  const peclet::core::decomp::BlockDecomposer<3>& decomposer() const { return halo_->decomposer(); }
  // Globally-unique particle ids (persistent-pair and Mindlin-history keys are gid-based):
  // re-base each rank's identity ids by an exclusive scan of the owned counts, once per particle
  // set. Migration and rebalance carry gids, so the ids stay stable afterwards; set_positions
  // resets the flag. Any pre-MPI ledger was keyed with the identity ids, which the re-base makes
  // stale — cold-start both engines' histories (one soft restart, negligible).
  void ensureGlobalGids() {
    if (mpiGidsGlobal_)
      return;
    long base = 0, mine = P_.numReal;
    MPI_Exscan(&mine, &base, 1, MPI_LONG, MPI_SUM, halo_->comm());
    if (halo_->rank() == 0)
      base = 0;  // MPI_Exscan leaves rank 0's recvbuf undefined
    fillGidBaseKokkos(P_.gid, P_.numReal, static_cast<int>(base));
    P_.prevPairCount = 0;
    P_.hertzPrevCount = 0;
    P_.hertzNumPairs = -1;
    mpiGidsGlobal_ = true;
  }
  void stepMpi(int nsteps) {
    requireDt("step_mpi");
    const double rcut = mpiRcut_;  // demStepMpi widens it to the global contact reach
    ensureGlobalGids();
    P_.splitStats.mlHubAggregated = 0;  // split_stats: max over this call's substeps
    if (P_.orphanClampCount.data() != nullptr)
      Kokkos::deep_copy(P_.orphanClampCount, 0);  // split_stats.orphanClamps: this call's sum
    for (int s = 0; s < nsteps; ++s) {
      if (mpiRebalanceEvery_ > 0 && mpiStepCount_ % mpiRebalanceEvery_ == 0)
        halo_->rebalance(P_);
      demStepMpi(P_, *halo_, rcut, mpiSyncEvery_, mpiForwardRotation_);
      ++mpiStepCount_;
    }
  }
  void stepHertzMpi(int substeps, float skin_frac) {
    requireDt("step_hertz_mpi");
    ensureGlobalGids();
    if (mpiRebalanceEvery_ > 0 && mpiHertzCalls_ % mpiRebalanceEvery_ == 0)
      halo_->rebalance(P_);
    ++mpiHertzCalls_;
    demStepHertzMpi(P_, *halo_, P_.dt, substeps, skin_frac);
  }
  int rank() const { return halo_->rank(); }
  int numGhost() const { return halo_->numGhost(); }
  void debugCaptureContacts(bool on) { P_.debugCapture = on; }
  void debugIterationCounters(bool on) { P_.iterCounters = on; }
  void debugNoAdaptiveStop(bool on) { P_.noAdaptiveStop = on; }
  const DebugContactCapture& debugCapturedContacts() const { return P_.debugCaptured; }
#endif  // PECLET_DEM_MPI
 
  // SDF grid (get_sdf_grid): Eikonal reconstruction over the domain, flat x-fastest, negative
  // inside solid.
  std::vector<float> getSdfGrid(int rx, int ry, int rz) {
    return peclet::dem::generateSdfKokkos(rx, ry, rz, P_.domain.min, P_.domain.max, P_.numReal,
                                          P_.pos, P_.quat, P_.scale, P_.shapeId, P_.shapes,
                                          P_.domain.periodic_x, P_.domain.periodic_y,
                                          P_.domain.periodic_z, P_.sdfGrid, P_.globalScale);
  }
 
  int numParticles() const { return P_.numReal; }
  // Live device Views of the owned particle state, for the zero-copy device-array export (H2): the
  // binding wraps a [0,numReal) subview as a DLPack/__cuda_array_interface__ array (or a NumPy view
  // on a host backend) referencing this memory — no device->host copy.
  const V3& positionsView() const { return P_.pos; }
  const V3& velocitiesView() const { return P_.vel; }
  int numContacts() { return readInt(P_.contactCount); }
  int numManifolds() { return readInt(P_.manifoldCount); }
  // TEST-ONLY colouring self-check: over the LAST solved substep's colourings, count how many
  // (manifold, contact) pairs violate the graph-colouring invariant "no two same-colour items share
  // a body" — must be exactly (0, 0) for a valid colouring (the incremental warm-start path must
  // never import a conflict). Returns {velocity-colour conflicts, position-colour conflicts}.
  std::pair<int, int> debugColoringConflicts() {
    using peclet::dem::CpExec;
    using peclet::dem::CpMem;
    CpExec space;
    const int nm = readInt(P_.manifoldCount);
    const int nc = readInt(P_.contactCount);
    // Colouring vertices include the hub copy slots (docs/contact_solve_framework.md §4.4), which
    // sit above the SoA's live span; under MPI the contacts also reference ghost slots.
    const int nb =
        std::max({P_.numParticles, P_.numReal, P_.velCopies.slotBase, P_.posCopies.slotBase}) +
        std::max(P_.velCopies.nCopies, P_.posCopies.nCopies) + 1;
    Kokkos::View<std::uint64_t*, CpMem> seen("dbg_color_seen", std::max(nb, P_.capacity + 1));
    int velConf = 0, posConf = 0;
    if (nm > 0) {
      Kokkos::deep_copy(space, seen, std::uint64_t(0));
      auto manifolds = P_.manifolds;
      auto realIdx = P_.realIndices;
      auto mColor = P_.manifoldColor;
      const peclet::dem::SlotOverride ov =
          P_.velCopies.nHubs > 0
              ? peclet::dem::SlotOverride{P_.velCopies.slotA, P_.velCopies.slotB, {}, 1.0f}
              : peclet::dem::SlotOverride{};
      Kokkos::parallel_reduce(
          "peclet::dem::dbg_vel_color", Kokkos::RangePolicy<CpExec>(space, 0, nm),
          KOKKOS_LAMBDA(int idx, int& acc) {
            const int c = mColor(idx);
            if (c < 0)
              return;
            const auto m = manifolds(idx);
            const std::uint64_t bit = std::uint64_t(1) << c;
            if ((Kokkos::atomic_fetch_or(&seen(ov.slotA(idx, realIdx(m.bodyA))), bit) >> c) & 1)
              acc += 1;
            if (m.bodyB >= 0 &&
                ((Kokkos::atomic_fetch_or(&seen(ov.slotB(idx, realIdx(m.bodyB))), bit) >> c) & 1))
              acc += 1;
          },
          velConf);
    }
    // Position graph: its edges are the position units (docs/contact_solve_framework.md §4.3),
    // each the leader contact's (bodyA, bodyB) -- a unit's contacts share one colour by design.
    const int nu = P_.numPosUnits;
    if (nc > 0 && nu > 0) {
      Kokkos::deep_copy(space, seen, std::uint64_t(0));
      auto contacts = P_.contacts;
      auto uColor = P_.unitColor;
      const peclet::dem::PosUnits units{P_.unitStart, P_.unitContacts};
      const peclet::dem::SlotOverride ov =
          P_.posCopies.nHubs > 0
              ? peclet::dem::SlotOverride{P_.posCopies.slotA, P_.posCopies.slotB, {}, 1.0f}
              : peclet::dem::SlotOverride{};
      Kokkos::parallel_reduce(
          "peclet::dem::dbg_pos_color", Kokkos::RangePolicy<CpExec>(space, 0, nu),
          KOKKOS_LAMBDA(int idx, int& acc) {
            const int c = uColor(idx);
            if (c < 0)
              return;
            const auto ct = contacts(units.leader(idx));
            const std::uint64_t bit = std::uint64_t(1) << c;
            if ((Kokkos::atomic_fetch_or(&seen(ov.slotA(idx, ct.bodyA)), bit) >> c) & 1)
              acc += 1;
            if (ct.bodyB >= 0 &&
                ((Kokkos::atomic_fetch_or(&seen(ov.slotB(idx, ct.bodyB)), bit) >> c) & 1))
              acc += 1;
          },
          posConf);
    }
    space.fence();
    return {velConf, posConf};
  }
  // TEST-ONLY: (velocity, position) edges the LAST substep's colourings left uncoloured (-1, or
  // -3 = no free colour): manifolds and position units. After hub splitting this is 0 or the step
  // threw (the colouring invariant, docs/contact_solve_framework.md §4.2, §12 S3). Single-rank:
  // under MPI the non-owned tail of the manifold list keeps stale colours.
  std::pair<int, int> debugColoringLeftovers() {
    using peclet::dem::CpExec;
    CpExec space;
    const int nm = readInt(P_.manifoldCount);
    const int nu = P_.numPosUnits;
    int velLeft = 0, posLeft = 0;
    auto mColor = P_.manifoldColor;
    auto uColor = P_.unitColor;
    if (nm > 0 && P_.velocityUseGS)
      Kokkos::parallel_reduce(
          "peclet::dem::dbg_vel_left", Kokkos::RangePolicy<CpExec>(space, 0, nm),
          KOKKOS_LAMBDA(int i, int& acc) {
            acc += (mColor(i) == -1 || mColor(i) == peclet::dem::kColorUncolourable) ? 1 : 0;
          },
          velLeft);
    if (nu > 0)
      Kokkos::parallel_reduce(
          "peclet::dem::dbg_pos_left", Kokkos::RangePolicy<CpExec>(space, 0, nu),
          KOKKOS_LAMBDA(int u, int& acc) {
            acc += (uColor(u) == -1 || uColor(u) == peclet::dem::kColorUncolourable) ? 1 : 0;
          },
          posLeft);
    space.fence();
    return {velLeft, posLeft};
  }
  // TEST-ONLY (docs/contact_solve_framework.md §WO-4 item 3): same-colour pairs at an AGGREGATE
  // in the last multilevel hierarchy, summed over levels -- a valid coarse colouring is
  // group-disjoint per (level, colour). 0 when no hierarchy was built.
  int debugMultilevelColoringConflicts() {
    using peclet::dem::CpExec;
    using peclet::dem::CpMem;
    const auto& L = P_.mlLast;
    if (L.numLevels <= 0 || L.numManifolds <= 0)
      return 0;
    CpExec space;
    const int nm = L.numManifolds, nb = L.numBodies;
    Kokkos::View<int*, CpMem> grp("dbg_ml_grp", nb);
    auto realIdx = P_.realIndices;
    auto manifolds = P_.manifolds;
    auto cp = P_.mlColorPacked;
    auto parent = P_.mlParent;
    Kokkos::parallel_for(
        "peclet::dem::dbg_ml_grp0", Kokkos::RangePolicy<CpExec>(space, 0, nb),
        KOKKOS_LAMBDA(int i) { grp(i) = i; });
    int total = 0;
    for (int lvl = 1; lvl <= L.numLevels; ++lvl) {
      const int pOff = L.parentOff[lvl - 1], ng = L.numGroups[lvl - 1];
      const int shift = 6 * (lvl - 1);
      Kokkos::parallel_for(
          "peclet::dem::dbg_ml_compose", Kokkos::RangePolicy<CpExec>(space, 0, nb),
          KOKKOS_LAMBDA(int i) { grp(i) = parent(pOff + grp(i)); });
      Kokkos::View<std::uint64_t*, CpMem> seen("dbg_ml_seen", ng + 1);
      int conf = 0;
      Kokkos::parallel_reduce(
          "peclet::dem::dbg_ml_color", Kokkos::RangePolicy<CpExec>(space, 0, nm),
          KOKKOS_LAMBDA(int idx, int& acc) {
            const int c = static_cast<int>((cp(idx) >> shift) & 63);
            if (c == peclet::dem::kMlSlotSkip)
              return;
            const auto m = manifolds(idx);
            const int gA = grp(realIdx(m.bodyA));
            const int gB = m.bodyB >= 0 ? grp(realIdx(m.bodyB)) : -1;
            const std::uint64_t bit = std::uint64_t(1) << c;
            if ((Kokkos::atomic_fetch_or(&seen(gA), bit) >> c) & 1)
              acc += 1;
            if (gB >= 0 && gB != gA && ((Kokkos::atomic_fetch_or(&seen(gB), bit) >> c) & 1))
              acc += 1;
          },
          conf);
      total += conf;
    }
    space.fence();
    return total;
  }
  // diagnostics.split_stats() (docs/contact_solve_framework.md §WO-4 item 4): the last substep's
  // body copies.
  peclet::dem::SplitStats debugSplitStats() const {
    peclet::dem::SplitStats st = P_.splitStats;
    if (P_.orphanClampCount.data() != nullptr)  // device counter (§13.3): read on demand only
      Kokkos::deep_copy(st.orphanClamps, P_.orphanClampCount);
    return st;
  }
  float maxOverlap() {
    float h;
    Kokkos::deep_copy(h, P_.maxOverlap);
    return h;
  }
  std::tuple<double, float, int> restBankStats() {
    return restBankStatsKokkos(P_.prevRestBank, P_.prevPairCount);
  }
  std::tuple<double, float, int> restOrphanStats() {
    return restBankStatsKokkos(P_.bodyOrphan, P_.numReal);
  }
 
  // ParaView PolyData (points + Radius + Velocity), faithful to CUDA Simulation::write_vtp:
  // Radius = scale * globalScale * baseRadius.
  void writeVtp(const std::string& filename) const {
    auto pos = Kokkos::create_mirror_view(P_.pos);
    Kokkos::deep_copy(pos, P_.pos);
    auto sc = Kokkos::create_mirror_view(P_.scale);
    Kokkos::deep_copy(sc, P_.scale);
    auto vel = Kokkos::create_mirror_view(P_.vel);
    Kokkos::deep_copy(vel, P_.vel);
    const int n = P_.numReal;
 
    std::ofstream out(filename);
    if (!out)
      throw std::runtime_error("Could not open file for writing: " + filename);
    out << "<?xml version=\"1.0\"?>\n";
    out << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    out << "  <PolyData>\n";
    out << "    <Piece NumberOfPoints=\"" << n << "\" NumberOfVerts=\"0\" "
        << "NumberOfLines=\"0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">\n";
    out << "      <Points>\n";
    out << "        <DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" "
           "format=\"ascii\">\n";
    for (int i = 0; i < n; ++i)
      out << pos(i, 0) << " " << pos(i, 1) << " " << pos(i, 2) << " ";
    out << "\n        </DataArray>\n";
    out << "      </Points>\n";
    out << "      <PointData Scalars=\"Radius\">\n";
    out << "        <DataArray type=\"Float32\" Name=\"Radius\" NumberOfComponents=\"1\" "
           "format=\"ascii\">\n";
    for (int i = 0; i < n; ++i)
      out << sc(i) * P_.globalScale * baseRadius_ << " ";
    out << "\n        </DataArray>\n";
    out << "        <DataArray type=\"Float32\" Name=\"Velocity\" NumberOfComponents=\"3\" "
           "format=\"ascii\">\n";
    for (int i = 0; i < n; ++i)
      out << vel(i, 0) << " " << vel(i, 1) << " " << vel(i, 2) << " ";
    out << "\n        </DataArray>\n";
    out << "      </PointData>\n";
    out << "    </Piece>\n";
    out << "  </PolyData>\n";
    out << "</VTKFile>\n";
    out.close();
    std::printf("Exported VTP: %s\n", filename.c_str());
  }
 
 private:
  bool dtSet_ = false;  // set_dt was called: every stepper checks it (no silent default)
  void requireDt(const char* who) const {
    if (!dtSet_)
      throw std::runtime_error(std::string(who) +
                               ": call set_dt(dt) first (there is no default time step)");
  }
  // (Re)upload just the small WallSdf array (velocity fields change every step for a vibrating
  // wall; the grid samples are uploaded once in addSdfWall).
  void uploadWalls() {
    const int n = std::max<int>(1, static_cast<int>(wallsHost_.size()));
    if (static_cast<int>(P_.walls.extent(0)) < n)
      P_.walls = Kokkos::View<WallSdf*, CpMem>("walls", n);
    auto h = Kokkos::create_mirror_view(P_.walls);
    for (size_t i = 0; i < wallsHost_.size(); ++i)
      h(i) = wallsHost_[i];
    Kokkos::deep_copy(P_.walls, h);
  }
 
  std::vector<WallSdf> wallsHost_;
  std::vector<float> wallGridHost_;
  std::vector<peclet::core::geom::ShapeNode<float>> wallNodesHost_;
  // R1: authored root transform per analytic-wall root node (see setWallTransform).
  std::map<int, peclet::core::geom::Transform<float>> wallRootTf_;
#ifdef PECLET_DEM_MPI
  std::unique_ptr<ParticleHalo> halo_ = std::make_unique<ParticleHalo>();
  double mpiRcut_ = 0.0;
  int mpiSyncEvery_ = 1;
  bool mpiForwardRotation_ = true;
  int mpiRebalanceEvery_ = 0;
  long mpiStepCount_ = 0;
  long mpiHertzCalls_ =
      0;  // step_hertz_mpi call count (rebalance_every cadence for the force path)
  bool mpiGidsGlobal_ = false;  // gids re-based to a global Exscan offset (once per particle set)
#endif
};
 
}  // namespace peclet::dem
 
#endif  // DEM_SIM_HPP