block_container.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_BLOCK_CONTAINER_HPP
#define PECLET_FLOW_VOF_BLOCK_CONTAINER_HPP
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdio>
#include <Kokkos_Core.hpp>
#include <map>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
 
#include "mac_stencils.hpp"  // peclet::flow::SExec, SField, I3, L3
#include "peclet/core/decomp/block_decomposer.hpp"
#include "vof/advect_wy.hpp"
#include "vof/block_batch.hpp"
#include "vof/colour_field.hpp"
#include "vof/curvature_field.hpp"
#include "vof/surface_tension.hpp"
#include "policy.hpp"
 
namespace peclet::flow::vof {
 
// ---------------------------------------------------------------------------------------------
// index boxes
// ---------------------------------------------------------------------------------------------
 
struct VofBox {
  int lo[3]{0, 0, 0};
  int hi[3]{0, 0, 0};
 
  int n(int d) const { return hi[d] - lo[d]; }
  long cells() const { return static_cast<long>(n(0)) * n(1) * n(2); }
  bool empty() const { return n(0) <= 0 || n(1) <= 0 || n(2) <= 0; }
  I3 size() const { return I3{n(0), n(1), n(2)}; }
  I3 origin() const { return I3{lo[0], lo[1], lo[2]}; }
 
  bool contains(const VofBox& o) const {
    for (int d = 0; d < 3; ++d)
      if (o.lo[d] < lo[d] || o.hi[d] > hi[d])
        return false;
    return true;
  }
  bool operator==(const VofBox& o) const {
    for (int d = 0; d < 3; ++d)
      if (lo[d] != o.lo[d] || hi[d] != o.hi[d])
        return false;
    return true;
  }
  bool operator!=(const VofBox& o) const { return !(*this == o); }
 
  static VofBox grown(const VofBox& b, int m) {
    VofBox r = b;
    for (int d = 0; d < 3; ++d) {
      r.lo[d] -= m;
      r.hi[d] += m;
    }
    return r;
  }
  static VofBox intersect(const VofBox& a, const VofBox& b) {
    VofBox r;
    for (int d = 0; d < 3; ++d) {
      r.lo[d] = a.lo[d] > b.lo[d] ? a.lo[d] : b.lo[d];
      r.hi[d] = a.hi[d] < b.hi[d] ? a.hi[d] : b.hi[d];
    }
    return r;
  }
};
 
struct VofRun {
  int g0 = 0;    
  int len = 0;   
  int loc0 = 0;  
};
 
inline int vofAxisRuns(int a, int b, int gs, bool periodic, VofRun out[5], int outside = 0) {
  int n = 0;
  if (b <= a)
    return 0;
  if (!periodic) {
    const int s = a < -outside ? -outside : a;
    const int e = b > gs + outside ? gs + outside : b;
    if (e > s)
      out[n++] = VofRun{s, e - s, s - a};
    return n;
  }
  if (b - a > 3 * gs)  // guarded by vofClampBox(); a block box is at most the whole axis
    throw std::runtime_error("peclet::flow::vof: block range exceeds three periodic images");
  int k = a;
  while (k < b && n < 5) {
    const int w = ((k % gs) + gs) % gs;
    int len = gs - w;
    if (len > b - k)
      len = b - k;
    out[n++] = VofRun{w, len, k - a};
    k += len;
  }
  return n;
}
 
inline VofBox vofClampBox(VofBox b, I3 gs, const std::array<bool, 3>& per) {
  const int g[3] = {gs.x, gs.y, gs.z};
  for (int d = 0; d < 3; ++d) {
    if (per[d]) {
      if (b.n(d) > g[d]) {
        b.lo[d] = 0;
        b.hi[d] = g[d];
      }
    } else {
      if (b.lo[d] < 0)
        b.lo[d] = 0;
      if (b.hi[d] > g[d])
        b.hi[d] = g[d];
      if (b.hi[d] <= b.lo[d]) {
        b.lo[d] = 0;
        b.hi[d] = g[d] < 1 ? 1 : g[d];
      }
    }
  }
  return b;
}
 
// ---------------------------------------------------------------------------------------------
// master assignment (rung W1 item a)
// ---------------------------------------------------------------------------------------------
 
enum class VofMasterAssign {
  RoundRobin = 0,  
  Lpt = 1,         
  WeightedOrb = 2  
};
 
inline void vofAssignLpt(const std::vector<long>& w, int size, std::vector<int>& master) {
  const std::size_t n = w.size();
  master.assign(n, 0);
  if (n == 0 || size <= 1)
    return;
  std::vector<std::size_t> ord(n);
  for (std::size_t i = 0; i < n; ++i)
    ord[i] = i;
  // stable, and the comparator is a strict weak ordering on the weight alone, so equal weights
  // keep block-id order -> with equal blocks LPT reproduces the round robin exactly.
  std::stable_sort(ord.begin(), ord.end(),
                   [&](std::size_t a, std::size_t b) { return w[a] > w[b]; });
  std::vector<long> load(static_cast<std::size_t>(size), 0);
  for (std::size_t k = 0; k < n; ++k) {
    int best = 0;
    for (int r = 1; r < size; ++r)
      if (load[static_cast<std::size_t>(r)] < load[static_cast<std::size_t>(best)])
        best = r;
    master[ord[k]] = best;
    load[static_cast<std::size_t>(best)] += w[ord[k]];
  }
}
 
inline void vofAssignOrb(const std::vector<long>& w, int size, std::vector<int>& master) {
  const std::size_t n = w.size();
  master.assign(n, 0);
  if (n == 0 || size <= 1)
    return;
  std::vector<peclet::core::Real> wr(n);
  for (std::size_t i = 0; i < n; ++i)
    wr[i] = static_cast<peclet::core::Real>(w[i]);
  peclet::core::decomp::BlockDecomposer<1> dec(
      static_cast<std::size_t>(size),
      peclet::core::IVec<1>{static_cast<peclet::core::Index>(n)}, wr);
  for (int r = 0; r < size; ++r) {
    const auto b = dec.block(static_cast<std::size_t>(r));
    for (auto i = b.origin[0]; i < b.origin[0] + b.size[0]; ++i)
      master[static_cast<std::size_t>(i)] = r;
  }
}
 
// ---------------------------------------------------------------------------------------------
// one block
// ---------------------------------------------------------------------------------------------
 
struct VofBlockStats {
  long id = 0;
  int master = 0;
  int lo[3]{0, 0, 0};
  int hi[3]{0, 0, 0};
  double volume = 0.0;  
  double centroid[3]{0, 0, 0};
  double velocity[3]{0, 0, 0};  
  double moment[6]{0, 0, 0, 0, 0, 0};
  bool recentred = false;  
  long cells = 0;          
  double discarded = 0.0;
  double area = 0.0;
  long debrisCells = 0;         
  double debrisVolume = 0.0;    
  double debrisReturned = 0.0;  
  double debrisLost = 0.0;      
  long debrisUnresolved = 0;    
  long fullAxis = 0;
  double residueReturned = 0.0;
};
 
class VofBlock {
 public:
  long id = 0;
  int master = 0;
  VofBox box;  
  VofBox seed_{};
 
  bool mine() const { return mine_; }
  WyAdvector& advector() { return adv_; }
  const WyAdvector& advector() const { return adv_; }
  VofBox extended(int ghost) const { return VofBox::grown(box, ghost); }
  const VofBlockStats& stats() const { return st_; }
  VofCurvature& curvature() { return curv_; }
  const VofCurvature& curvature() const { return curv_; }
  SField csfForce(int c) const { return f_[c]; }
  static constexpr int kAuxLen = 10;
  void serializeAux(double out[kAuxLen]) const {
    for (int d = 0; d < 3; ++d)
      out[d] = prevCentroid_[d];
    out[3] = hasPrev_ ? 1.0 : 0.0;
    out[4] = st_.debrisReturned;
    out[5] = st_.debrisLost;
    out[6] = static_cast<double>(st_.debrisUnresolved);  // a cell count: exact in a double
    out[7] = st_.discarded;  // §12.2: the re-centring ledger migrates too
    out[8] = static_cast<double>(st_.fullAxis);
    out[9] = st_.residueReturned;
  }
  void deserializeAux(const double in[kAuxLen]) {
    for (int d = 0; d < 3; ++d)
      prevCentroid_[d] = in[d];
    hasPrev_ = (in[3] != 0.0);
    st_.debrisReturned = in[4];
    st_.debrisLost = in[5];
    st_.debrisUnresolved = static_cast<long>(in[6]);
    st_.discarded = in[7];
    st_.fullAxis = static_cast<long>(in[8]);
    st_.residueReturned = in[9];
  }
 
  friend class VofBlockSet;
 
 private:
  bool mine_ = false;
  bool allocated_ = false;
  WyAdvector adv_;
  VofCurvature curv_;
  SField f_[3];
  VofBlockStats st_;
  double prevCentroid_[3]{0, 0, 0};
  bool hasPrev_ = false;
};
 
// ---------------------------------------------------------------------------------------------
// the exchange interface (implemented in vof/block_exchange.hpp)
// ---------------------------------------------------------------------------------------------
 
class VofBlockSet;
 
struct VofBlockExchangeBase {
  virtual ~VofBlockExchangeBase() = default;
  virtual void gatherFaceVel(std::vector<VofBlock>& blocks, int ghost) = 0;
  virtual void scatterColourMax(std::vector<VofBlock>& blocks, SField cLocal) = 0;
  virtual void syncTable(std::vector<VofBlock>& blocks) = 0;
  virtual void gatherColour(std::vector<VofBlock>& blocks, SField cLocal) = 0;
  virtual void scatterForceSum(std::vector<VofBlock>& blocks, SField fx, SField fy, SField fz) = 0;
  virtual void scatterColourSum(std::vector<VofBlock>&, SField) {
    throw std::runtime_error("peclet::flow::vof: this block exchange has no colour SUM scatter");
  }
  virtual void migrateColour(std::vector<VofBlock>& blocks, const std::vector<int>& oldMaster) = 0;
  virtual long gatherBytes() const { return 0; }
  virtual long scatterBytes() const { return 0; }
};
 
// ---------------------------------------------------------------------------------------------
// the block set
// ---------------------------------------------------------------------------------------------
 
class VofBlockSet {
 public:
  void init(I3 gs, std::array<bool, 3> per, int rank, int size, double h = 1.0) {
    h_ = h;
    gs_ = gs;
    per_ = per;
    rank_ = rank;
    size_ = size < 1 ? 1 : size;
    blocks_.clear();
    statsPend_.clear();
    step_ = 0;
  }
 
  void setMetric(const VofMetric& g) {
    metric_ = g;
    curvProto.metric = g;
    for (auto& b : blocks_) {
      b.adv_.metric = g;
      b.curv_.metric = g;
    }
  }
  const VofMetric& metric() const { return metric_; }
 
  void setExchange(std::shared_ptr<VofBlockExchangeBase> x) { exch_ = std::move(x); }
  VofBlockExchangeBase* exchange() const { return exch_.get(); }
 
  int margin() const { return margin_; }
  int ghost() const { return ghost_; }
  int recentrePad = 2;
  double cflLimit = 0.25;
  bool batched = true;
  VofMasterAssign assignMode = VofMasterAssign::RoundRobin;
  long reassignEvery = 0;
  bool usePool = true;
  int poolCapacity = 8;  
  bool csfEnabled = false;
  double sigma = 0.0;  
  VofCurvature curvProto = blockCurvProto();
  static VofCurvature blockCurvProto() {
    VofCurvature c;
    c.kappaMax = -1.0;
    return c;
  }
  double bubbleEps = 1e-12;
  double wispEps = 0.0;
  void setWispEps(double eps) {
    wispEps = eps;
    curvProto.pureEps = eps;  // the cascade's purity test follows the advector (review finding 1)
    for (auto& b : blocks_)
      if (b.allocated_) {
        b.adv_.wispEps = eps;
        b.curv_.pureEps = eps;
      }
  }
  double areaEps = 1e-8;
  bool allowShrink = true;
  static constexpr double kDebrisMaxC = 0.5;
  bool debrisCensus = true;
  bool debrisRemove = false;
  bool timingOn = false;
  double debrisSeconds = 0.0;
  long debrisCalls = 0;
  static double tick_() {
    Kokkos::fence();
    return WyAdvector::wallSeconds();
  }
  static constexpr double kDebrisMinAttached = 1.0;
  double cellSize() const { return h_; }
 
  std::size_t count() const { return blocks_.size(); }
  std::vector<VofBlock>& blocks() {
    flushStats();
    return blocks_;
  }
  const std::vector<VofBlock>& blocks() const { return blocks_; }
  long step() const { return step_; }
  void setStep(long s) { step_ = s; }
 
  int masterOf(long id) const { return static_cast<int>(id % size_); }
 
  void seedSphere(double cx, double cy, double cz, double r, int subLevels = 4, int extraPad = 0) {
    flushStats();  // C2: deferred statistics are indexed by block; apply them before the table grows
    VofBox bb;
    const double c[3] = {cx, cy, cz};
    for (int d = 0; d < 3; ++d) {
      bb.lo[d] = static_cast<int>(std::floor((c[d] - r) / h_)) - 1;
      bb.hi[d] = static_cast<int>(std::ceil((c[d] + r) / h_)) + 1;
    }
    VofBlock b;
    b.id = static_cast<long>(blocks_.size());
    b.master = masterOf(b.id);
    b.box = vofClampBox(VofBox::grown(bb, margin_ + extraPad), gs_, per_);
    b.mine_ = (b.master == rank_);
    blocks_.push_back(std::move(b));
    const std::size_t idx = blocks_.size() - 1;
    // Re-install every earlier block's hook: `push_back` may have reallocated the table, and the
    // hooks capture `this` + an index (so they survive it), but the advectors themselves moved.
    for (std::size_t k = 0; k + 1 < blocks_.size(); ++k)
      if (blocks_[k].mine_ && blocks_[k].allocated_)
        installHook(k);
    VofBlock& nb = blocks_[idx];
    if (nb.mine_) {
      allocate(idx);
      const I3 e = nb.adv_.extent();
      const int g = ghost_;
      const I3 o = nb.box.origin();
      const double hh = h_;
      SField cf = nb.adv_.colour();
      Kokkos::parallel_for(
          "vof::block::seed_sphere",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const double gx = (x - g + o.x) * hh, gy = (y - g + o.y) * hh, gz = (z - g + o.z) * hh;
            cf(L3(x, y, z, e)) = sphereCellFraction(cx, cy, cz, r, gx, gy, gz, hh, subLevels);
          });
      Kokkos::fence();
      fillBlockGhosts(blocks_[idx]);
      measure(blocks_[idx], 0.0);
    }
  }
 
  void advect(double dt, SField cLocal) {
    if (!exch_)
      throw std::runtime_error("peclet::flow::vof::VofBlockSet: no exchange installed");
    exch_->gatherFaceVel(blocks_, ghost_);
    std::vector<VofBlock*> deb;
    const bool batch = batchEligible();
    if (batch)
      advectBatched(dt);  // read #1 also applies the previous step's deferred statistics (C2)
    else
      flushStats();
    for (auto& b : blocks_) {
      if (!b.mine_)
        continue;
      if (!batch)
        b.adv_.advect(dt, step_);
      b.st_.recentred = false;
      // Before the re-centring, so the bubble box a re-centre measures is the one the census saw.
      // §13: also whenever the advectors run a wisp guard -- the sub-wispEps residue must be
      // returned on the kinematic path too (the leak is not a CSF matter).
      if ((csfEnabled && debrisCensus) || wispEps > 0.0)
        deb.push_back(&b);
    }
    // After every block has advected: each block's pass touches only its own colour, so running
    // them together (their sums in one launch) is the per-block sequence reordered, bit for bit.
    if (!deb.empty())
      debrisPassBatch(deb, /*deferLedger=*/batch);
    if (batch) {
      // §5.9 stage 5: every bubble box in one launch, then read #2 -- the boxes, the debris sums
      // and the WO-R ledger in ONE device->host copy; the recentre decisions exactly as before.
      recentreBatched(!deb.empty());
    } else {
      for (std::size_t k = 0; k < blocks_.size(); ++k)
        if (blocks_[k].mine_)
          recentre(k);
    }
    exch_->syncTable(blocks_);
    for (auto& b : blocks_)
      b.mine_ = (b.master == rank_);
    // Rung W1: re-balance the master assignment on the CURRENT boxes, before the scatter, so the
    // step's union is produced by the new owners and there is no half-migrated state anywhere.
    // The boxes have just been replicated, so `plannedMasters()` is the same on every rank.
    lastReassigned_ = 0;
    if (reassignEvery > 0 && ((step_ + 1) % reassignEvery) == 0)
      lastReassigned_ = assignMasters();
    exch_->scatterColourMax(blocks_, cLocal);
    if (batch) {
      measureBatched(dt);  // C2 (§5.10): one launch, no host read; applied at the next read #1
    } else {
      for (auto& b : blocks_)
        if (b.mine_)
          measure(b, dt);
    }
    ++step_;
  }
 
  long lastReassigned() const { return lastReassigned_; }
 
  // ---- rung W2: per-block curvature + the CSF face force ---------------------------------------
 
  void enableCsf(double sigmaValue) {
    csfEnabled = true;
    sigma = sigmaValue;
    for (std::size_t i = 0; i < blocks_.size(); ++i)
      if (blocks_[i].mine_ && blocks_[i].allocated_)
        allocateCsf(i);
  }
 
  void computeCsf(SField fx, SField fy, SField fz) {
    if (!exch_)
      throw std::runtime_error("peclet::flow::vof::VofBlockSet: no exchange installed");
    curvStats_ = VofCurvature::Stats{};
    if (csfBatchEligible()) {  // C1 §5.9 stage 6: every pass for all blocks in one launch
      computeCsfBatched();
      exch_->scatterForceSum(blocks_, fx, fy, fz);
      return;
    }
    // Tier 3 of every local block in ONE batched launch (VofCurvature::computeBegin): bit-identical
    // to the per-block cascade, which ran the 5^3 fit's full latency once per block, serially.
    bool batch = true;
    for (auto& b : blocks_)
      if (b.mine_ && !b.curv_.useWorklist)
        batch = false;
    if (batch) {
      std::vector<VofCurvFallbackJob> jobs;
      for (auto& b : blocks_)
        if (b.mine_) {
          b.curv_.computeBegin(b.adv_.colour());
          jobs.push_back(b.curv_.fallbackJob(b.adv_.colour()));
        }
      VofCurvature::fallbackBatch(jobs);
    }
    for (auto& b : blocks_) {
      if (!b.mine_)
        continue;
      const VofCurvature::Stats st =
          batch ? b.curv_.computeEnd() : b.curv_.compute(b.adv_.colour());
      curvStats_.interfacial += st.interfacial;
      curvStats_.hf += st.hf;
      curvStats_.hfMixed += st.hfMixed;
      curvStats_.hfFit += st.hfFit;
      curvStats_.pv += st.pv;
      curvStats_.pvReduced += st.pvReduced;
      curvStats_.noEstimate += st.noEstimate;
      curvStats_.clipped += st.clipped;
      buildCsfForce(b);
    }
    exch_->scatterForceSum(blocks_, fx, fy, fz);
  }
  VofCurvature::Stats csfCurvatureStats() const { return curvStats_; }
 
  void setKappaMax(double km) {
    curvProto.kappaMax = km;
    for (auto& b : blocks_)
      b.curv_.kappaMax = km;
  }
 
  void allocateCsf(std::size_t idx) {
    VofBlock& b = blocks_[idx];
    const I3 n = b.adv_.inner();
    b.curv_.init(n.x, n.y, n.z, ghost_);
    // EVERY tunable of the prototype, the metric included (review finding 1: this list once
    // dropped `pureEps`, so a block WY at wispEps 1e-8 was judged at the 1e-10 purity floor and the
    // HF tier decayed into the PV fallback silently, 63.8 % -> 52.3 %; and it never carried the
    // metric, so a cascade re-allocated by a re-centring ran an anisotropic block at unit metric).
    b.curv_.copyTunablesFrom(curvProto);
    const long len = static_cast<long>(b.adv_.extent().x) * b.adv_.extent().y * b.adv_.extent().z;
    for (int c = 0; c < 3; ++c)
      b.f_[c] = SField("vof::block::csf", len);
  }
 
  void buildCsfForce(VofBlock& b) {
    const I3 e = b.adv_.extent(), n = b.adv_.inner();
    const int g = ghost_;
    // Phase 3 (V3.1): the same per-axis pressure-gradient weight `1/h_a'^2` the global
    // `Solver::addCsfRhs` uses. `h_` stays the (index) block spacing and is not that weight.
    const double sig = sigma;
    const double wgt[3] = {metric_.h[0] * metric_.h[0], metric_.h[1] * metric_.h[1],
                           metric_.h[2] * metric_.h[2]};
    SField cv = b.adv_.colour(), kp = b.curv_.kappa(), kb = b.curv_.branch();
    const long sy = e.x, sz = static_cast<long>(e.x) * e.y;
    for (int c = 0; c < 3; ++c) {
      SField ff = b.f_[c];
      const double wc = wgt[c];
      const long strd = (c == 0) ? 1 : (c == 1 ? sy : sz);
      Kokkos::parallel_for(
          "vof::block::csf_force",
          MDRange3<SExec>(SExec(), {g, g, g},
                                                        {g + n.x, g + n.y, g + n.z}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const long i = L3(x, y, z, e);
            const double dC = cv(i) - cv(i - strd);
            double f = 0.0;
            if (dC != 0.0) {
              double kf = 0.0;
              vof::csfFaceCurvature(kp(i - strd), kb(i - strd), kp(i), kb(i), kf);
              f = vof::csfFaceForce(sig, kf, dC, wc);
            }
            ff(i) = f;
          });
    }
  }
 
  void scatterSum(SField sLocal) {
    if (!exch_)
      throw std::runtime_error("peclet::flow::vof::VofBlockSet: no exchange installed");
    exch_->scatterColourSum(blocks_, sLocal);
  }
 
  void scatter(SField cLocal) {
    if (!exch_)
      throw std::runtime_error("peclet::flow::vof::VofBlockSet: no exchange installed");
    exch_->syncTable(blocks_);
    exch_->scatterColourMax(blocks_, cLocal);
  }
 
  // ---- rung W1 item (a): master assignment and re-assignment ---------------------------------
 
  std::vector<int> plannedMasters() const {
    std::vector<long> w(blocks_.size());
    for (std::size_t i = 0; i < blocks_.size(); ++i)
      w[i] = blocks_[i].box.cells();
    std::vector<int> m;
    switch (assignMode) {
      case VofMasterAssign::Lpt:
        vofAssignLpt(w, size_, m);
        break;
      case VofMasterAssign::WeightedOrb:
        vofAssignOrb(w, size_, m);
        break;
      case VofMasterAssign::RoundRobin:
      default:
        m.assign(blocks_.size(), 0);
        for (std::size_t i = 0; i < blocks_.size(); ++i)
          m[i] = static_cast<int>(blocks_[i].id % size_);
        break;
    }
    return m;
  }
 
  long assignMasters() {
    flushStats();  // C2: a migrating block carries its previous centroid (serializeAux)
    if (blocks_.empty() || size_ <= 1)
      return 0;
    const std::vector<int> want = plannedMasters();
    std::vector<int> old(blocks_.size());
    long moved = 0;
    for (std::size_t i = 0; i < blocks_.size(); ++i) {
      old[i] = blocks_[i].master;
      if (old[i] != want[i])
        ++moved;
    }
    if (moved == 0)
      return 0;
    for (std::size_t i = 0; i < blocks_.size(); ++i) {
      blocks_[i].master = want[i];
      blocks_[i].st_.master = want[i];
    }
    // gaining ranks allocate FIRST (zeroed), so the migration has somewhere to land
    for (std::size_t i = 0; i < blocks_.size(); ++i)
      if (want[i] == rank_ && old[i] != rank_) {
        blocks_[i].mine_ = true;
        allocate(i);
      }
    if (exch_)
      exch_->migrateColour(blocks_, old);
    for (std::size_t i = 0; i < blocks_.size(); ++i) {
      if (old[i] == rank_ && want[i] != rank_) {
        release(i);
        blocks_[i].mine_ = false;
      } else if (want[i] == rank_ && old[i] != rank_) {
        // A GAINED block only. A block this rank keeps is left completely alone -- `measure` would
        // reset its previous centroid and blank one step of its reported velocity.
        blocks_[i].mine_ = true;
        fillBlockGhosts(blocks_[i]);
        double aux[VofBlock::kAuxLen];
        blocks_[i].serializeAux(aux);  // the migrated previous centroid survives the re-measure
        measure(blocks_[i], 0.0);
        blocks_[i].deserializeAux(aux);
      }
    }
    return moved;
  }
 
  // ---- general seeding: a block from an arbitrary global colour field -------------------------
 
  void seedBox(const VofBox& bb, int extraPad = 0) {
    flushStats();  // C2, as in seedSphere
    VofBlock b;
    b.id = static_cast<long>(blocks_.size());
    b.master = static_cast<int>(b.id % size_);
    b.box = vofClampBox(VofBox::grown(bb, margin_ + extraPad), gs_, per_);
    b.seed_ = bb;
    b.mine_ = (b.master == rank_);
    blocks_.push_back(std::move(b));
    const std::size_t idx = blocks_.size() - 1;
    for (std::size_t k = 0; k + 1 < blocks_.size(); ++k)
      if (blocks_[k].mine_ && blocks_[k].allocated_)
        installHook(k);
    if (blocks_[idx].mine_)
      allocate(idx);
  }
 
  void finishSeeding(SField cLocal) {
    if (!exch_)
      throw std::runtime_error("peclet::flow::vof::VofBlockSet: no exchange installed");
    flushStats();  // C2: its per-block measure() must not be overwritten by an older packet
    exch_->gatherColour(blocks_, cLocal);
    for (auto& b : blocks_)
      if (b.mine_) {
        clipToSeed(b);
        fillBlockGhosts(b);
        measure(b, 0.0);
      }
  }
 
  // ---- rung W3: checkpoint / restart of the block state ---------------------------------------
  //
  // A run longer than one job's wall clock has to write its markers out and read them back. The
  // union colour field is NOT enough to do that: `seedBox`/`finishSeeding` gather each marker out
  // of the UNION and clip to the seed extent, so two markers that touch each adopt a slice of the
  // other (WO-W12 open item 5, measured -2.7 % / +7.1 %). A block's ONLY state is its own inner
  // colour (everything else is either replicated -- the table -- or recomputed every step), so a
  // checkpoint is exactly {box, colour} per block and the restart below is EXACT whatever the
  // markers are doing.
 
  std::vector<double> blockColourHost(std::size_t idx) {
    VofBlock& b = blocks_.at(idx);
    if (!b.mine_ || !b.allocated_)
      return {};
    const I3 n = b.adv_.inner();
    SField c = b.adv_.colour();
    auto hc = Kokkos::create_mirror_view(c);
    Kokkos::deep_copy(hc, c);
    std::vector<double> out(static_cast<std::size_t>(n.x) * static_cast<std::size_t>(n.y) *
                            static_cast<std::size_t>(n.z));
    std::size_t q = 0;
    for (int z = 0; z < n.z; ++z)
      for (int y = 0; y < n.y; ++y)
        for (int x = 0; x < n.x; ++x)
          out[q++] = hc(b.adv_.index(x, y, z));
    return out;
  }
 
  std::vector<double> blockKappaHost(std::size_t idx) {
    VofBlock& b = blocks_.at(idx);
    if (!b.mine_ || !b.allocated_ || !csfEnabled)
      return {};
    const I3 n = b.adv_.inner();
    SField k = b.curv_.kappa();
    auto hk = Kokkos::create_mirror_view(k);
    Kokkos::deep_copy(hk, k);
    std::vector<double> out(static_cast<std::size_t>(n.x) * static_cast<std::size_t>(n.y) *
                            static_cast<std::size_t>(n.z));
    std::size_t q = 0;
    for (int z = 0; z < n.z; ++z)
      for (int y = 0; y < n.y; ++y)
        for (int x = 0; x < n.x; ++x)
          out[q++] = hk(b.adv_.index(x, y, z));
    return out;
  }
 
  void seedBoxWithColour(const VofBox& bb, const std::vector<double>& colour) {
    flushStats();  // C2, as in seedSphere
    VofBlock b;
    b.id = static_cast<long>(blocks_.size());
    b.master = masterOf(b.id);
    b.box = vofClampBox(bb, gs_, per_);
    b.mine_ = (b.master == rank_);
    blocks_.push_back(std::move(b));
    const std::size_t idx = blocks_.size() - 1;
    for (std::size_t k = 0; k + 1 < blocks_.size(); ++k)
      if (blocks_[k].mine_ && blocks_[k].allocated_)
        installHook(k);
    if (!blocks_[idx].mine_)
      return;
    allocate(idx);
    VofBlock& nb = blocks_[idx];
    const I3 n = nb.adv_.inner();
    const std::size_t want = static_cast<std::size_t>(n.x) * static_cast<std::size_t>(n.y) *
                             static_cast<std::size_t>(n.z);
    if (colour.size() != want)
      throw std::runtime_error(
          "peclet::flow::vof::VofBlockSet::seedBoxWithColour: colour size does not match the box");
    SField c = nb.adv_.colour();
    auto hc = Kokkos::create_mirror_view(c);
    Kokkos::deep_copy(hc, c);  // ghosts as allocated (zero); the fill below re-does them anyway
    std::size_t q = 0;
    for (int z = 0; z < n.z; ++z)
      for (int y = 0; y < n.y; ++y)
        for (int x = 0; x < n.x; ++x)
          hc(nb.adv_.index(x, y, z)) = colour[q++];
    Kokkos::deep_copy(c, hc);
    Kokkos::fence();
    fillBlockGhosts(nb);
    measure(nb, 0.0);
  }
 
  void clipToSeed(VofBlock& b) {
    const I3 e = b.adv_.extent(), n = b.adv_.inner(), o = b.box.origin();
    const int g = ghost_;
    const int lx = b.seed_.lo[0], hx = b.seed_.hi[0], ly = b.seed_.lo[1], hy = b.seed_.hi[1],
              lz = b.seed_.lo[2], hz = b.seed_.hi[2];
    if (hx <= lx || hy <= ly || hz <= lz)
      return;  // no seed box recorded (a sphere seed): nothing to clip
    SField c = b.adv_.colour();
    Kokkos::parallel_for(
        "vof::block::clip_seed",
        MDRange3<SExec>(SExec(), {0, 0, 0}, {n.x, n.y, n.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const int gx = x + o.x, gy = y + o.y, gz = z + o.z;
          if (gx < lx || gx >= hx || gy < ly || gy >= hy || gz < lz || gz >= hz)
            c(L3(x + g, y + g, z + g, e)) = 0.0;
        });
    Kokkos::fence();
  }
 
  void fillBlockGhosts(VofBlock& b) {
    const I3 e = b.adv_.extent(), n = b.adv_.inner(), o = b.box.origin();
    const int g = ghost_;
    SField f = b.adv_.colour();
    // 1. every ghost cell -> 0 (the block's far field: C = 0, the continuous phase, by the margin).
    Kokkos::parallel_for(
        "vof::block::ghost_zero",
        MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          if (x >= g && x < g + n.x && y >= g && y < g + n.y && z >= g && z < g + n.z)
            return;
          f(L3(x, y, z, e)) = 0.0;
        });
    // 2. an axis the block spans entirely, and which is globally periodic, wraps within the block.
    const bool sx = per_[0] && b.box.n(0) == gs_.x;
    const bool sy = per_[1] && b.box.n(1) == gs_.y;
    const bool sz = per_[2] && b.box.n(2) == gs_.z;
    if (sx || sy || sz)
      vof::periodicFill(f, e, g, sx, sy, sz);
    // 3. the part of the extended box outside a NON-periodic domain takes the globally clamped
    //    (zero-gradient) value — the same rule the structured colour field uses.
    if (!per_[0] || !per_[1] || !per_[2])
      vof::clampFill(f, e, g, o, gs_, per_[0], per_[1], per_[2]);
  }
 
  std::vector<VofBlockStats> statsAll() {
    flushStats();  // C2: the last batched measure's values reach the host here or at read #1
    // `id`, `master` and the box come from the REPLICATED table, so they are reported on every
    // rank; only the measured entries (volume, centroid, velocity, moments, area) are the
    // master's.  Filling them here rather than in `measure()` is what makes that promise true on
    // a rank that masters nothing.
    std::vector<VofBlockStats> v;
    v.reserve(blocks_.size());
    for (const auto& b : blocks_) {
      VofBlockStats q = b.st_;
      q.id = b.id;
      q.master = b.master;
      for (int d = 0; d < 3; ++d) {
        q.lo[d] = b.box.lo[d];
        q.hi[d] = b.box.hi[d];
      }
      q.cells = b.box.cells();
      v.push_back(q);
    }
    return v;
  }
 
  void masterCensus(std::vector<long>& masters, std::vector<long>& cells) const {
    masters.assign(size_, 0);
    cells.assign(size_, 0);
    for (const auto& b : blocks_) {
      masters[b.master] += 1;
      cells[b.master] += VofBox::grown(b.box, 0).cells();
    }
  }
  double cellImbalance() const {
    std::vector<long> m, c;
    masterCensus(m, c);
    long tot = 0, mx = 0;
    for (long v : c) {
      tot += v;
      mx = v > mx ? v : mx;
    }
    if (tot == 0)
      return 1.0;
    return static_cast<double>(mx) * size_ / static_cast<double>(tot);
  }
 
  I3 globalSize() const { return gs_; }
  std::array<bool, 3> periodic() const { return per_; }
  int rank() const { return rank_; }
  int size() const { return size_; }
 
 private:
  void installHook(std::size_t idx) {
    VofBlockSet* self = this;
    blocks_[idx].adv_.exchange = [self, idx](SField) { self->fillBlockGhosts(self->blocks_[idx]); };
  }
 
  void allocate(std::size_t idx) {
    VofBlock& b = blocks_[idx];
    b.adv_ = acquireAdvector(b.box.n(0), b.box.n(1), b.box.n(2));
    installHook(idx);
    b.allocated_ = true;
    b.st_.id = b.id;
    b.st_.master = b.master;
    if (csfEnabled)
      allocateCsf(idx);
  }
 
  void release(std::size_t idx) {
    VofBlock& b = blocks_[idx];
    retireAdvector(std::move(b.adv_));
    b.adv_ = WyAdvector();
    b.curv_ = VofCurvature();
    for (int c = 0; c < 3; ++c)
      b.f_[c] = SField();
    b.allocated_ = false;
    b.hasPrev_ = false;
  }
 
  WyAdvector acquireAdvector(int nx, int ny, int nz) {
    const std::array<int, 3> key{nx, ny, nz};
    if (usePool) {
      auto it = pool_.find(key);
      if (it != pool_.end() && !it->second.empty()) {
        WyAdvector a = std::move(it->second.back());
        it->second.pop_back();
        ++poolHits_;
        Kokkos::deep_copy(a.colour(), 0.0);
        for (int d = 0; d < 3; ++d)
          Kokkos::deep_copy(a.faceVel(d), 0.0);
        Kokkos::fence();
        a.cflLimit = cflLimit;
        a.wispEps = wispEps;
        a.metric = metric_;  // a pooled advector carries whatever metric it last had
        a.globalMax = nullptr;
        a.exchange = nullptr;  // re-installed by installHook
        return a;
      }
    }
    ++poolMisses_;
    WyAdvector a;
    a.init(nx, ny, nz, h_, ghost_);
    a.cflLimit = cflLimit;
    a.wispEps = wispEps;
    a.metric = metric_;  // a re-centred block's fresh advector keeps the anisotropic metric
    a.globalMax = nullptr;  // the block IS the whole domain of its own advector
    return a;
  }
 
  void retireAdvector(WyAdvector&& a) {
    if (!usePool || a.colour().extent(0) == 0)
      return;
    const I3 n = a.inner();
    const std::array<int, 3> key{n.x, n.y, n.z};
    auto& v = pool_[key];
    if (static_cast<int>(v.size()) >= poolCapacity)
      return;
    a.exchange = nullptr;  // the hook captures a block index; never carry it into the pool
    v.push_back(std::move(a));
  }
 
 public:
  long poolHits() const { return poolHits_; }
  long poolMisses() const { return poolMisses_; }
  void clearPool() { pool_.clear(); }
 
  // NOTE: everything below is public ONLY because nvcc refuses an extended `__host__ __device__`
  // lambda inside a private or protected member function ("The enclosing parent function ... cannot
  // have private or protected access within its class"). These are implementation details of the
  // step; treat them as such.
 public:
  bool bubbleBox(const VofBlock& b, VofBox& out) const {
    const I3 e = b.adv_.extent(), n = b.adv_.inner();
    const int g = ghost_;
    const double eps = bubbleEps > wispEps ? bubbleEps : wispEps;  // §12.2
    SField c = b.adv_.colour();
    int lo[3] = {n.x, n.y, n.z}, hi[3] = {-1, -1, -1};
    using MD = MDRange3<SExec>;
    MD pol(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z});
    Kokkos::parallel_reduce(
        "vof::block::bbox_lo", pol,
        KOKKOS_LAMBDA(int x, int y, int z, int& lx, int& ly, int& lz) {
          if (!(Kokkos::fabs(c(L3(x, y, z, e))) > eps))
            return;
          const int px = x - g, py = y - g, pz = z - g;
          lx = px < lx ? px : lx;
          ly = py < ly ? py : ly;
          lz = pz < lz ? pz : lz;
        },
        Kokkos::Min<int>(lo[0]), Kokkos::Min<int>(lo[1]), Kokkos::Min<int>(lo[2]));
    Kokkos::parallel_reduce(
        "vof::block::bbox_hi", pol,
        KOKKOS_LAMBDA(int x, int y, int z, int& hx, int& hy, int& hz) {
          if (!(Kokkos::fabs(c(L3(x, y, z, e))) > eps))
            return;
          const int px = x - g, py = y - g, pz = z - g;
          hx = px > hx ? px : hx;
          hy = py > hy ? py : hy;
          hz = pz > hz ? pz : hz;
        },
        Kokkos::Max<int>(hi[0]), Kokkos::Max<int>(hi[1]), Kokkos::Max<int>(hi[2]));
    Kokkos::fence();
    if (hi[0] < lo[0])
      return false;  // the block is empty
    for (int d = 0; d < 3; ++d) {
      out.lo[d] = lo[d] + b.box.lo[d];
      out.hi[d] = hi[d] + 1 + b.box.lo[d];
    }
    return true;
  }
 
  void recentre(std::size_t idx) {
    VofBox bb;
    if (!bubbleBox(blocks_[idx], bb))
      return;
    recentreTo(idx, bb);
  }
  void recentreTo(std::size_t idx, const VofBox& bb) {
    VofBlock& b = blocks_[idx];
    const VofBox need = vofClampBox(VofBox::grown(bb, margin_), gs_, per_);
    bool wasteful = false;
    if (allowShrink)
      for (int d = 0; d < 3; ++d)
        if (b.box.n(d) > need.n(d) + 2 * (margin_ + recentrePad))
          wasteful = true;
    if (b.box.contains(need) && !wasteful)
      return;
    const VofBox want = VofBox::grown(bb, margin_ + recentrePad);
    const VofBox nb = vofClampBox(want, gs_, per_);
    if (nb == b.box)
      return;
    // §12.3: an axis on which the box spans the whole periodic domain is copied WITH the wrap
    // (source index = destination global index mod L) -- lossless, a marker as long as the domain
    // is legal. Every other axis keeps the plain copy by unwrapped global index.
    const int L[3] = {gs_.x, gs_.y, gs_.z};
    bool wrap[3];
    bool anyWrap = false;
    for (int d = 0; d < 3; ++d) {
      wrap[d] = per_[d] && nb.n(d) == L[d];
      anyWrap = anyWrap || (wrap[d] && want.n(d) >= L[d]);
    }
    // What the move DROPS: the colour of old-box cells that fall outside the new box. Measured
    // directly rather than as `sum(old) - sum(new)`, which at |sum| ~ 1e2 is 1e-13 of summation
    // rounding and says nothing about the wisps (measured -5.7e-14 that way on the G1 scene).
    const double dropped = outOfBoxSum(b, nb, wrap);
    // exact copy by global index (both boxes are in the same unwrapped global frame)
    WyAdvector fresh = acquireAdvector(nb.n(0), nb.n(1), nb.n(2));
    const I3 se = b.adv_.extent(), de = fresh.extent();
    const I3 dn = fresh.inner();
    const int g = ghost_;
    const int sx = b.box.lo[0], sy = b.box.lo[1], sz = b.box.lo[2];
    const int snx = b.box.n(0), sny = b.box.n(1), snz = b.box.n(2);
    const int dx = nb.lo[0], dy = nb.lo[1], dz = nb.lo[2];
    const bool wx = wrap[0], wy = wrap[1], wz = wrap[2];
    const int Lx = L[0], Ly = L[1], Lz = L[2];
    SField src = b.adv_.colour(), dst = fresh.colour();
    Kokkos::parallel_for(
        "vof::block::recentre_copy",
        MDRange3<SExec>(SExec(), {0, 0, 0}, {dn.x, dn.y, dn.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const int gx = x + dx, gy = y + dy, gz = z + dz;
          int lx = gx - sx, ly = gy - sy, lz = gz - sz;
          if (wx)
            lx = ((lx % Lx) + Lx) % Lx;
          if (wy)
            ly = ((ly % Ly) + Ly) % Ly;
          if (wz)
            lz = ((lz % Lz) + Lz) % Lz;
          double v = 0.0;
          if (lx >= 0 && lx < snx && ly >= 0 && ly < sny && lz >= 0 && lz < snz)
            v = src(L3(lx + g, ly + g, lz + g, se));
          dst(L3(x + g, y + g, z + g, de)) = v;
        });
    Kokkos::fence();
    b.box = nb;
    retireAdvector(std::move(b.adv_));
    b.adv_ = std::move(fresh);
    b.adv_.cflLimit = cflLimit;
    b.adv_.globalMax = nullptr;
    installHook(idx);
    if (csfEnabled)
      allocateCsf(idx);
    fillBlockGhosts(blocks_[idx]);
    blocks_[idx].st_.recentred = true;
    blocks_[idx].st_.discarded += dropped;
    if (anyWrap)
      ++blocks_[idx].st_.fullAxis;
  }
 
  void debrisPass(VofBlock& b) {
    std::vector<VofBlock*> one{&b};
    debrisPassBatch(one);
  }
  void debrisPassBatch(std::vector<VofBlock*>& bs, bool deferLedger = false) {
    const double t0 = timingOn ? tick_() : 0.0;
    if (batchEligible()) {
      debrisBatched(bs);
      if (!deferLedger) {
        Kokkos::deep_copy(pkHost_, pk_);
        applyDebrisLedger(bs, csfEnabled && debrisCensus);
      }
    } else {
      debrisPassBatchImpl(bs);
    }
    if (timingOn) {
      Kokkos::fence();
      debrisSeconds += tick_() - t0;
    }
    debrisCalls += static_cast<long>(bs.size());
  }
  struct DebrisJob {
    double* c;
    const long *listD, *listR, *listA;
    long nDc, nDr, nR, nA;
  };
  static constexpr int kDebrisBatch = 16;
  struct DebrisTable {
    DebrisJob job[kDebrisBatch];
    int nj;
  };
  void debrisPassBatchImpl(std::vector<VofBlock*>& bs) {
    const bool doDebris = csfEnabled && debrisCensus, doResidue = wispEps > 0.0;
    const int g = ghost_;
    // The interfacial predicate of D and A is `max(interfaceEps, wispEps)` (review finding 7): R
    // is `0 < |C| <= wispEps`, so the three lists are DISJOINT whatever the two thresholds are --
    // on a kinematic block run curvProto.interfaceEps is 0 while wispEps is 1e-8, and an A/R or
    // D/R overlap would drop a W share or count one cell in both ledgers.
    const double weps = wispEps, cfull = kDebrisMaxC;
    const double ieps = curvProto.interfaceEps > weps ? curvProto.interfaceEps : weps;
    // per-block list offsets into the shared scratch
    std::vector<long> off(bs.size() + 1, 0);
    for (std::size_t k = 0; k < bs.size(); ++k) {
      const I3 n = bs[k]->adv_.inner();
      off[k + 1] = off[k] + static_cast<long>(n.x) * n.y * n.z;
    }
    auto ensure = [&](LField& f, const char* name) {
      if (static_cast<long>(f.extent(0)) < off.back())
        f = LField(Kokkos::view_alloc(std::string(name), Kokkos::WithoutInitializing), off.back());
    };
    if (doDebris)
      ensure(listD_, "vof::block::debris_list");
    if (doResidue)
      ensure(listR_, "vof::block::residue_list");
    struct Pending {
      VofBlock* b;
      long k, nD, nR, nA, nDr;
      bool removeD;
    };
    std::vector<Pending> pend;
    std::vector<DebrisJob> jobs;
    for (std::size_t k = 0; k < bs.size(); ++k) {
      VofBlock& b = *bs[k];
      b.st_.debrisCells = 0;
      b.st_.debrisVolume = 0.0;
      const I3 e = b.adv_.extent(), n = b.adv_.inner();
      const long region = static_cast<long>(n.x) * n.y * n.z, o = off[k];
      const long sy = e.x, sz = static_cast<long>(e.x) * e.y;
      const int nx = n.x, ny = n.y;
      SField c = b.adv_.colour();
      long nD = 0, nR = 0;
      if (doDebris) {
        LField listD = listD_;
        Kokkos::parallel_scan(
            "vof::block::debris_mark", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
            KOKKOS_LAMBDA(const long r, long& upd, const bool final) {
              const int ix = static_cast<int>(r % nx);
              const int iy = static_cast<int>((r / nx) % ny);
              const int iz = static_cast<int>(r / (static_cast<long>(nx) * ny));
              const long i = L3(g + ix, g + iy, g + iz, e);
              if (!vofDebrisCell(c, i, sy, sz, ieps, cfull))
                return;
              if (final)
                listD(o + upd) = i;
              ++upd;
            },
            nD);
      }
      if (doResidue) {
        LField listR = listR_;
        Kokkos::parallel_scan(
            "vof::block::residue_mark", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
            KOKKOS_LAMBDA(const long r, long& upd, const bool final) {
              const int ix = static_cast<int>(r % nx);
              const int iy = static_cast<int>((r / nx) % ny);
              const int iz = static_cast<int>(r / (static_cast<long>(nx) * ny));
              const double ci = c(L3(g + ix, g + iy, g + iz, e));
              if (!vofResidueCell(ci, weps))
                return;
              if (final)
                listR(o + upd) = L3(g + ix, g + iy, g + iz, e);
              ++upd;
            },
            nR);
      }
      const bool removeD = doDebris && debrisRemove && nD > 0;
      if (!removeD && nR == 0) {
        if (doDebris && nD > 0) {  // census only: report what WOULD be removed
          LField listD = listD_;
          double dVd = 0.0;
          Kokkos::parallel_reduce(
              "vof::block::debris_census", Kokkos::RangePolicy<SExec>(SExec(), 0, 1),
              KOKKOS_LAMBDA(const long, double& acc) {
                double a = 0.0;
                for (long t = 0; t < nD; ++t)
                  a += c(listD(o + t));
                acc += a;
              },
              dVd);
          b.st_.debrisCells = nD;
          b.st_.debrisVolume = dVd;
        }
        continue;  // nothing to remove: nothing else runs
      }
      ensure(listA_, "vof::block::attached_list");
      LField listA = listA_;
      long nA = 0;
      Kokkos::parallel_scan(
          "vof::block::debris_attached", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
          KOKKOS_LAMBDA(const long r, long& upd, const bool final) {
            const int ix = static_cast<int>(r % nx);
            const int iy = static_cast<int>((r / nx) % ny);
            const int iz = static_cast<int>(r / (static_cast<long>(nx) * ny));
            const long i = L3(g + ix, g + iy, g + iz, e);
            if (!vofAttachedCell(c, i, sy, sz, ieps, cfull))
              return;
            if (final)
              listA(o + upd) = i;
            ++upd;
          },
          nA);
      Pending pb{&b, static_cast<long>(k), nD, nR, nA, removeD ? nD : 0, removeD};
      pend.push_back(pb);
      DebrisJob j;
      j.c = c.data();
      j.listD = doDebris ? listD_.data() + o : nullptr;
      j.listR = doResidue ? listR_.data() + o : nullptr;
      j.listA = listA_.data() + o;
      j.nDc = doDebris ? nD : 0;
      j.nDr = pb.nDr;
      j.nR = nR;
      j.nA = nA;
      jobs.push_back(j);
    }
    if (pend.empty())
      return;
    // [0] sum_D, [1] sum_R, [2] W, [3] VA, [4] removed volume (then `lost`), [5] acted?
    // The SUMS run on one thread per block in list order (a fixed summation order: the state stays
    // decomposition-independent bitwise); the per-cell writes are independent once dV and W are
    // known -- D, R and A are disjoint -- so they run in parallel and every cell receives exactly
    // the double the serial loop gave it (review finding 4). `lost` (a cap at 1, expected never to
    // fire) is summed in list order from a per-cell scratch only when some cell was capped.
    const long np = static_cast<long>(pend.size());
    if (static_cast<long>(debrisOut_.extent(0)) < 6 * np) {
      debrisOut_ = SField("vof::block::debris_out", 6 * np);
      debrisOutHost_ = Kokkos::create_mirror_view(debrisOut_);
    }
    SField out = debrisOut_;
    const double vaMin = kDebrisMinAttached;
    for (long j0 = 0; j0 < np; j0 += kDebrisBatch) {
      DebrisTable T;
      T.nj = static_cast<int>(std::min<long>(kDebrisBatch, np - j0));
      for (int q = 0; q < T.nj; ++q)
        T.job[q] = jobs[j0 + q];
      const long base = 6 * j0;
      Kokkos::parallel_for(
          "vof::block::debris_sums", Kokkos::RangePolicy<SExec>(SExec(), 0, T.nj),
          KOKKOS_LAMBDA(const long q) {
            const DebrisJob& J = T.job[q];
            const VofRawField c{J.c};
            const long nDc = J.nDc, nDr = J.nDr, nR = J.nR, nA = J.nA;
            double sD = 0.0, sR = 0.0, W = 0.0, VA = 0.0;
            for (long t = 0; t < nDc; ++t)
              sD += c(J.listD[t]);
            for (long t = 0; t < nR; ++t)
              sR += c(J.listR[t]);
            for (long t = 0; t < nA; ++t) {
              const double ca = c(J.listA[t]);
              W += ca * (1.0 - ca);
              VA += ca;
            }
            const long ob = base + 6 * q;
            out(ob + 0) = sD;
            out(ob + 1) = sR;
            out(ob + 2) = W;
            out(ob + 3) = VA;
            out(ob + 4) = 0.0;
            out(ob + 5) = (VA < vaMin || !(W > 0.0)) ? 0.0 : 1.0;
            // the removed volume, in the fixed order D (if removing) then R -- the serial order
            double dV = 0.0;
            for (long t = 0; t < nDr; ++t)
              dV += c(J.listD[t]);
            for (long t = 0; t < nR; ++t)
              dV += c(J.listR[t]);
            out(ob + 4) = dV;  // carried to the write kernel; `lost` replaces it below
          });
    }
    Kokkos::deep_copy(debrisOutHost_, out);  // ONE device->host copy for every block
    auto hoAll = debrisOutHost_;
    for (long p = 0; p < np; ++p) {
      const Pending& pb = pend[p];
      VofBlock& b = *pb.b;
      const long o = off[pb.k];
      const I3 n = b.adv_.inner();
      const long region = static_cast<long>(n.x) * n.y * n.z;
      double ho[6];
      for (int q = 0; q < 6; ++q)
        ho[q] = hoAll(6 * p + q);
      const long nD = pb.nD, nR = pb.nR, nA = pb.nA, nDr = pb.nDr;
      const bool removeD = pb.removeD;
      if (ho[5] != 0.0) {
        const double dV = ho[4], W = ho[2];
        if (static_cast<long>(debrisEx_.extent(0)) < region)
          debrisEx_ = SField(Kokkos::view_alloc("vof::block::debris_excess",
                                                Kokkos::WithoutInitializing),
                             region);
        SField ex = debrisEx_;
        SField c = b.adv_.colour();
        LField listD = listD_, listR = listR_, listA = listA_;
        const long nZ = nDr + nR, nTot = nZ + nA;
        long capped = 0;
        Kokkos::parallel_reduce(
            "vof::block::debris_write", Kokkos::RangePolicy<SExec>(SExec(), 0, nTot),
            KOKKOS_LAMBDA(const long t, long& cap) {
              if (t < nDr) {
                c(listD(o + t)) = 0.0;
                return;
              }
              if (t < nZ) {
                c(listR(o + t - nDr)) = 0.0;
                return;
              }
              const long ta = t - nZ;
              const long i = listA(o + ta);
              const double ca = c(i);
              const double d = dV * (ca * (1.0 - ca)) / W;
              double cn = ca + d;
              double e = 0.0;
              if (cn > 1.0) {
                e = cn - 1.0;
                cn = 1.0;
                ++cap;
              }
              ex(ta) = e;
              c(i) = cn;
            },
            capped);
        double lost = 0.0;
        if (capped > 0) {
          Kokkos::parallel_reduce(
              "vof::block::debris_lost", Kokkos::RangePolicy<SExec>(SExec(), 0, 1),
              KOKKOS_LAMBDA(const long, double& acc) {
                double a = 0.0;  // list order, as the serial loop accumulated it
                for (long t = 0; t < nA; ++t)
                  a += ex(t);
                acc += a;
              },
              lost);
        }
        ho[4] = lost;
      }
      if (doDebris && !debrisRemove && nD > 0) {  // census of what WOULD be removed
        b.st_.debrisCells = nD;
        b.st_.debrisVolume = ho[0];
      }
      if (ho[5] == 0.0) {  // step 3: no attached interface worth a cell -- the clip bounds it
        if (removeD)
          b.st_.debrisUnresolved += nD;
        continue;
      }
      fillBlockGhosts(b);  // the inner colour moved; the next reader expects valid ghosts
      if (removeD) {
        b.st_.debrisCells = nD;
        b.st_.debrisVolume = ho[0];
        b.st_.debrisReturned += ho[0];
      }
      b.st_.residueReturned += ho[1];
      b.st_.debrisLost += ho[4];
    }
  }
 
  double outOfBoxSum(const VofBlock& b, const VofBox& nb) const {
    const bool nowrap[3] = {false, false, false};
    return outOfBoxSum(b, nb, nowrap);
  }
  double outOfBoxSum(const VofBlock& b, const VofBox& nb, const bool wrap[3]) const {
    const I3 e = b.adv_.extent(), n = b.adv_.inner();
    const int g = ghost_;
    const int sx = b.box.lo[0], sy = b.box.lo[1], sz = b.box.lo[2];
    const int lx = nb.lo[0], ly = nb.lo[1], lz = nb.lo[2];
    const int hx = nb.hi[0], hy = nb.hi[1], hz = nb.hi[2];
    const bool wx = wrap[0], wy = wrap[1], wz = wrap[2];
    SField c = b.adv_.colour();
    double v = 0.0;
    Kokkos::parallel_reduce(
        "vof::block::out_of_box",
        MDRange3<SExec>(SExec(), {g, g, g},
                                                      {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z, double& a) {
          const int gx = x - g + sx, gy = y - g + sy, gz = z - g + sz;
          const bool inX = wx || (gx >= lx && gx < hx), inY = wy || (gy >= ly && gy < hy),
                     inZ = wz || (gz >= lz && gz < hz);
          if (inX && inY && inZ)
            return;
          a += c(L3(x, y, z, e));
        },
        v);
    Kokkos::fence();
    return v;
  }
 
  double innerSum(const VofBlock& b) const {
    const I3 e = b.adv_.extent(), n = b.adv_.inner();
    const int g = ghost_;
    SField c = b.adv_.colour();
    double v = 0.0;
    Kokkos::parallel_reduce(
        "vof::block::inner_sum",
        MDRange3<SExec>(SExec(), {g, g, g},
                                                      {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z, double& a) { a += c(L3(x, y, z, e)); }, v);
    Kokkos::fence();
    return v;
  }
 
  double interfaceArea(const VofBlock& b) const {
    const I3 e = b.adv_.extent(), n = b.adv_.inner();
    const int g = ghost_;
    const double hh = h_, aeps = areaEps;
    SField c = b.adv_.colour();
    const long sy = e.x, sz = static_cast<long>(e.x) * e.y;
    double a = 0.0;
    Kokkos::parallel_reduce(
        "vof::block::area",
        MDRange3<SExec>(SExec(), {g, g, g},
                                                      {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
          const long i = L3(x, y, z, e);
          const double ci = c(i);
          if (!(ci > aeps) || !(ci < 1.0 - aeps))
            return;
          double st[27];
          for (int dz = -1; dz <= 1; ++dz)
            for (int dy = -1; dy <= 1; ++dy)
              for (int dx = -1; dx <= 1; ++dx)
                st[vof::plicSt(dx + 1, dy + 1, dz + 1)] = c(i + dx + dy * sy + dz * sz);
          double m[3];
          vof::mycNormal(st, m);
          const double al = vof::plicAlpha(m[0], m[1], m[2], ci);
          double v[8][3], ctr[3], ar = 0.0;
          const int nv = vof::plicPolygon(m[0], m[1], m[2], al, v);
          vof::polygonAreaCentroid(v, nv, ctr, ar);
          acc += ar;
        },
        a);
    Kokkos::fence();
    (void)hh;
    return a;  // CELL units squared, matching `volume`'s cell volumes
  }
 
  void measure(VofBlock& b, double dt) {
    const I3 e = b.adv_.extent(), n = b.adv_.inner(), o = b.box.origin();
    const int g = ghost_;
    const double hh = h_;
    SField c = b.adv_.colour();
    using MD = MDRange3<SExec>;
    MD pol(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z});
    double v = 0.0, mx = 0.0, my = 0.0, mz = 0.0;
    Kokkos::parallel_reduce(
        "vof::block::moments1", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& av, double& ax, double& ay, double& az) {
          const double q = c(L3(x, y, z, e));
          av += q;
          ax += q * (x - g + o.x + 0.5) * hh;
          ay += q * (y - g + o.y + 0.5) * hh;
          az += q * (z - g + o.z + 0.5) * hh;
        },
        v, mx, my, mz);
    Kokkos::fence();
    b.st_.id = b.id;
    b.st_.master = b.master;
    for (int d = 0; d < 3; ++d) {
      b.st_.lo[d] = b.box.lo[d];
      b.st_.hi[d] = b.box.hi[d];
    }
    b.st_.cells = b.box.cells();
    b.st_.volume = v;
    if (v <= 0.0) {
      b.st_.area = 0.0;
      for (int d = 0; d < 3; ++d)
        b.st_.centroid[d] = b.st_.velocity[d] = 0.0;
      for (int d = 0; d < 6; ++d)
        b.st_.moment[d] = 0.0;
      return;
    }
    const double cx = mx / v, cy = my / v, cz = mz / v;
    double sxx = 0, syy = 0, szz = 0, sxy = 0, sxz = 0, syz = 0;
    Kokkos::parallel_reduce(
        "vof::block::moments2", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& axx, double& ayy, double& azz, double& axy,
                      double& axz, double& ayz) {
          const double q = c(L3(x, y, z, e));
          const double px = (x - g + o.x + 0.5) * hh - cx, py = (y - g + o.y + 0.5) * hh - cy,
                       pz = (z - g + o.z + 0.5) * hh - cz;
          axx += q * px * px;
          ayy += q * py * py;
          azz += q * pz * pz;
          axy += q * px * py;
          axz += q * px * pz;
          ayz += q * py * pz;
        },
        sxx, syy, szz, sxy, sxz, syz);
    Kokkos::fence();
    const double nc[3] = {cx, cy, cz};
    for (int d = 0; d < 3; ++d) {
      b.st_.velocity[d] = (b.hasPrev_ && dt > 0.0) ? (nc[d] - b.prevCentroid_[d]) / dt : 0.0;
      b.prevCentroid_[d] = nc[d];
      b.st_.centroid[d] = nc[d];
    }
    b.hasPrev_ = true;
    b.st_.area = interfaceArea(b);
    const double iv = 1.0 / v;
    b.st_.moment[0] = sxx * iv;
    b.st_.moment[1] = syy * iv;
    b.st_.moment[2] = szz * iv;
    b.st_.moment[3] = sxy * iv;
    b.st_.moment[4] = sxz * iv;
    b.st_.moment[5] = syz * iv;
  }
 
  // ---- C1: the batched container stages (`vof/block_batch.hpp`) -------------------------------
 
  bool csfBatchEligible() const {
    if (!batchEligible())
      return false;
    for (const auto& b : blocks_) {
      if (!b.mine_)
        continue;
      if (!b.curv_.ready() || !b.curv_.useWorklist || b.curv_.timingOn || !b.f_[0].extent(0))
        return false;
    }
    return true;
  }
 
  void computeCsfBatched() {
    std::vector<VofBlock*> mb;
    for (auto& b : blocks_)
      if (b.mine_)
        mb.push_back(&b);
    const std::size_t nb = mb.size();
    if (nb == 0)
      return;
    std::vector<VofCurvJob> jobs(nb);
    bool anyClip = false;
    for (std::size_t k = 0; k < nb; ++k) {
      VofBlock& b = *mb[k];
      const VofCurvature& cv = b.curv_;
      const VofCurvFallbackJob fj = cv.fallbackJob(b.adv_.colour());
      VofCurvJob J;
      J.c = fj.c;
      J.mx = fj.mx;
      J.my = fj.my;
      J.mz = fj.mz;
      J.al = fj.al;
      J.kap = fj.kap;
      J.br = fj.br;
      J.slot = cv.slotMap().data();
      for (int c = 0; c < 3; ++c)
        J.f[c] = b.f_[c].data();
      J.e = cv.extent();
      J.n = cv.inner();
      J.g = cv.ghost();
      J.ieps = cv.interfaceEps;
      J.mtol = cv.monoTol;
      J.ptW = cv.ptWeightWidth * cv.metric.maxH();  // heightPass's `ptW`
      J.peps = cv.pureEps;
      J.dW = fj.dW;  // fallbackPass's `weightWidth * metric.maxH()`
      J.cmin = fj.cmin;
      J.km = cv.kappaClipValue();
      J.forceFb = cv.debugForceFallback;
      J.oneDir = cv.debugSingleDirection;
      J.useFit = cv.useMixedHeightFit;
      J.gm = cv.metric;
      anyClip = anyClip || (J.km > 0.0);
      jobs[k] = J;
    }
    if (static_cast<std::size_t>(cStartG_.extent(0)) < nb) {
      cStartG_ = LField("vof::block::batch_curv_sG", nb);
      cEndG_ = LField("vof::block::batch_curv_eG", nb);
      cStartI_ = LField("vof::block::batch_curv_sI", nb);
      cEndI_ = LField("vof::block::batch_curv_eI", nb);
      cCnt_ = LField("vof::block::batch_curv_counts", 8 * nb);
      cCntHost_ = Kokkos::create_mirror_view(cCnt_);
    }
    Kokkos::deep_copy(SExec(), cCnt_, 0L);
    const double wgt[3] = {metric_.h[0] * metric_.h[0], metric_.h[1] * metric_.h[1],
                           metric_.h[2] * metric_.h[2]};  // buildCsfForce's `1/h_a'^2` weights
    for (std::size_t j0 = 0; j0 < nb; j0 += kVofBlockBatch) {
      VofCurvTable T;
      T.nj = static_cast<int>(std::min<std::size_t>(kVofBlockBatch, nb - j0));
      T.base = static_cast<int>(j0);
      for (int k = 0; k < T.nj; ++k)
        T.job[k] = jobs[j0 + k];
      const long rg = vofSetOffsets(T, vofCurvGrown);
      if (static_cast<long>(cListG_.extent(0)) < rg)
        cListG_ = LField(
            Kokkos::view_alloc(std::string("vof::block::batch_curv_listG"), Kokkos::WithoutInitializing),
            rg);
      // design G §5.2: the chunk's moment cache, sized and reallocated exactly like cListG_ (the
      // counts stay on the device); its positions are the slots, hence the int guard
      if (rg >= (1L << 31))
        throw std::runtime_error(
            "peclet::flow::vof::VofBlockSet: grown list capacity >= 2^31 (the slot map is int)");
      if (static_cast<long>(cCacheG_.extent(0)) < rg)
        cCacheG_ = Kokkos::View<PvMoments*, SMem>(
            Kokkos::view_alloc(std::string("vof::block::batch_curv_cacheG"), Kokkos::WithoutInitializing),
            rg);
      T.cache = cCacheG_.data();
      vofCurvCompact(T, true, cListG_, cStartG_, cEndG_);
      vofCurvPlanesZero(T);
      vofCurvPlanesList(T, cListG_, cStartG_, cEndG_);
      const long ri = vofSetOffsets(T, vofCurvInner);
      if (static_cast<long>(cListI_.extent(0)) < ri)
        cListI_ = LField(
            Kokkos::view_alloc(std::string("vof::block::batch_curv_listI"), Kokkos::WithoutInitializing),
            ri);
      vofCurvCompact(T, false, cListI_, cStartI_, cEndI_);
      vofCurvHfReset(T);
      vofCurvListPasses(T, anyClip, cListI_, cStartI_, cEndI_, cCnt_);  // tiers 1-3, clip, census
      vofCsfForceBatch(T, sigma, wgt[0], wgt[1], wgt[2]);
    }
    Kokkos::deep_copy(cCntHost_, cCnt_);  // the one read of the curvature stage
    for (std::size_t k = 0; k < nb; ++k) {
      const long* ct = &cCntHost_(8 * k);
      curvStats_.interfacial += ct[0];
      curvStats_.hf += ct[1];
      curvStats_.hfMixed += ct[2];
      curvStats_.hfFit += ct[3];
      curvStats_.pv += ct[4];
      curvStats_.pvReduced += ct[5];
      curvStats_.noEstimate += ct[6];
      curvStats_.clipped += ct[7];
    }
  }
 
  bool batchEligible() const {
    if (!batched)
      return false;
    for (const auto& b : blocks_) {
      if (!b.mine_)
        continue;
      const WyAdvector& a = b.adv_;
      if (!b.allocated_ || a.hasGeometry() || a.debugRecomputeDilation || !a.useWorklist ||
          a.timingOn || a.globalMax || !a.exchange || a.ghost() != ghost_)
        return false;
    }
    return true;
  }
 
  VofBlockJob makeJob(VofBlock& b, double dt) const {
    WyAdvector& a = b.adv_;
    VofBlockJob J;
    J.c = a.colour().data();
    J.mx = a.planeM(0).data();
    J.my = a.planeM(1).data();
    J.mz = a.planeM(2).data();
    J.al = a.planeAlpha().data();
    J.fl = a.faceFlux().data();
    for (int d = 0; d < 3; ++d)
      J.u[d] = a.faceVel(d).data();
    J.cc = a.dilationFlag().data();
    J.outside = a.hasOutsideMask() ? a.outsideMask().data() : nullptr;
    J.e = a.extent();
    J.n = a.inner();
    J.o = b.box.origin();
    J.g = a.ghost();
    J.dth = dt / a.h();  // WyAdvector::advect's `dth`
    J.weps = a.wispEps;
    J.interfaceCfl = a.interfaceLocalCfl;
    const int L[3] = {gs_.x, gs_.y, gs_.z};
    for (int d = 0; d < 3; ++d)
      J.span[d] = per_[d] && b.box.n(d) == L[d];
    for (int f = 0; f < 6; ++f)
      J.bcOwn[f] = a.bcFaceOwned(f) ? 1u : 0u;
    return J;
  }
 
  static VofBlockTable batchTable(const std::vector<VofBlockJob>& jobs, std::size_t j0) {
    VofBlockTable T;
    T.nj = static_cast<int>(std::min<std::size_t>(kVofBlockBatch, jobs.size() - j0));
    T.base = static_cast<int>(j0);
    for (int k = 0; k < T.nj; ++k)
      T.job[k] = jobs[j0 + k];
    return T;
  }
 
  void batchGhostFill(VofBlockTable T) {
    vofSetOffsets(T, vofJobLen);
    vofBatchGhostZero(T);
    for (int a = 0; a < 3; ++a) {
      VofBlockTable P;
      P.nj = 0;
      P.base = T.base;
      for (int k = 0; k < T.nj; ++k)
        if (T.job[k].span[a])
          P.job[P.nj++] = T.job[k];
      if (P.nj == 0)
        continue;
      const int b = (a + 1) % 3, c = (a + 2) % 3;
      vofSetOffsets(P, [b, c](const VofBlockJob& J) {
        const int dims[3] = {J.e.x, J.e.y, J.e.z};
        return static_cast<long>(dims[b]) * dims[c];
      });
      vofBatchPeriodic(P, a);
    }
    if (!per_[0] || !per_[1] || !per_[2])
      vofBatchClamp(T, gs_, per_[0], per_[1], per_[2]);
  }
 
  void ensureRead1(std::size_t len) {
    if (static_cast<std::size_t>(bRead1_.extent(0)) >= len)
      return;
    if (bRead1_.extent(0) == 0)
      bRead1_ = SField("vof::block::batch_read1", len);
    else
      Kokkos::resize(bRead1_, len);
    bRead1Host_ = Kokkos::create_mirror_view(bRead1_);
  }
 
  void measureBatched(double dt) {
    flushStats();
    std::vector<std::size_t> idx;
    for (std::size_t k = 0; k < blocks_.size(); ++k)
      if (blocks_[k].mine_)
        idx.push_back(k);
    const std::size_t nb = idx.size();
    if (nb == 0)
      return;
    std::vector<VofBlockJob> jobs(nb);
    for (std::size_t k = 0; k < nb; ++k) {
      VofBlock& b = blocks_[idx[k]];
      jobs[k] = makeJob(b, 0.0);
      b.st_.id = b.id;
      b.st_.master = b.master;
      for (int d = 0; d < 3; ++d) {
        b.st_.lo[d] = b.box.lo[d];
        b.st_.hi[d] = b.box.hi[d];
      }
      b.st_.cells = b.box.cells();
    }
    ensureRead1(kVofStats * nb);
    for (std::size_t j0 = 0; j0 < nb; j0 += kVofBlockBatch) {
      const VofBlockTable T = batchTable(jobs, j0);
      VofStatsPrev P;
      for (int k = 0; k < T.nj; ++k) {
        const VofBlock& b = blocks_[idx[j0 + k]];
        for (int d = 0; d < 3; ++d)
          P.c[k][d] = b.prevCentroid_[d];
        P.valid[k] = b.hasPrev_ ? 1u : 0u;
      }
      vofBatchStats(T, P, h_, areaEps, dt, bRead1_, 0);
    }
    statsPend_ = std::move(idx);
  }
 
  void applyStats(const double* hs) {
    for (std::size_t k = 0; k < statsPend_.size(); ++k) {
      VofBlock& b = blocks_[statsPend_[k]];
      const double* q = hs + kVofStats * k;
      b.st_.volume = q[0];
      if (q[0] <= 0.0) {
        b.st_.area = 0.0;
        for (int d = 0; d < 3; ++d)
          b.st_.centroid[d] = b.st_.velocity[d] = 0.0;
        for (int d = 0; d < 6; ++d)
          b.st_.moment[d] = 0.0;
        continue;
      }
      for (int d = 0; d < 3; ++d) {
        b.st_.velocity[d] = q[4 + d];
        b.prevCentroid_[d] = q[1 + d];
        b.st_.centroid[d] = q[1 + d];
      }
      b.hasPrev_ = true;
      b.st_.area = q[7];
      for (int d = 0; d < 6; ++d)
        b.st_.moment[d] = q[8 + d];
    }
    statsPend_.clear();
  }
 
  void flushStats() {
    if (statsPend_.empty())
      return;
    const std::size_t len = kVofStats * statsPend_.size();
    Kokkos::deep_copy(Kokkos::subview(bRead1Host_, std::make_pair(std::size_t(0), len)),
                      Kokkos::subview(bRead1_, std::make_pair(std::size_t(0), len)));
    applyStats(bRead1Host_.data());
  }
 
  void advectBatched(double dt) {
    std::vector<VofBlock*> mb;
    for (auto& b : blocks_)
      if (b.mine_)
        mb.push_back(&b);
    if (mb.empty()) {
      flushStats();
      return;
    }
    const std::size_t nb = mb.size();
    std::vector<VofBlockJob> jobs(nb);
    bool anyMask = false;
    for (std::size_t k = 0; k < nb; ++k) {
      jobs[k] = makeJob(*mb[k], dt);
      anyMask = anyMask || (jobs[k].outside != nullptr);
    }
    if (static_cast<std::size_t>(bStart_.extent(0)) < nb) {
      bStart_ = LField("vof::block::batch_start", nb);
      bEnd_ = LField("vof::block::batch_end", nb);
    }
    // read #1's packet: the deferred statistics of the last batched measure (C2) first, then this
    // step's Courant numbers -- one device->host copy for both
    const std::size_t off = kVofStats * statsPend_.size();
    ensureRead1(off + nb);
    auto cfl = Kokkos::subview(bRead1_, std::make_pair(off, off + nb));
    ensurePacket(nb);
    // (1) Courant numbers: one launch, one read, the throw before any block moves (R4)
    for (std::size_t j0 = 0; j0 < nb; j0 += kVofBlockBatch)
      vofBatchCfl(batchTable(jobs, j0), cfl);
    Kokkos::deep_copy(Kokkos::subview(bRead1Host_, std::make_pair(std::size_t(0), off + nb)),
                      Kokkos::subview(bRead1_, std::make_pair(std::size_t(0), off + nb)));
    applyStats(bRead1Host_.data());
    for (std::size_t k = 0; k < nb; ++k)
      mb[k]->adv_.enforceCfl(bRead1Host_(off + k), dt);
    Kokkos::deep_copy(SExec(), pk_, 0.0);  // the WO-R ledger slots accumulate over the sweeps
    const int* perm = kWySweepPerm[static_cast<int>(step_ % 6)];
    for (std::size_t j0 = 0; j0 < nb; j0 += kVofBlockBatch) {
      VofBlockTable T = batchTable(jobs, j0);
      // (2) THE dilation flag, frozen once from C^n
      vofSetOffsets(T, vofJobLen);
      vofBatchFreeze(T);
      // (3) three directional sweeps
      for (int s = 0; s < 3; ++s) {
        const int d = perm[s];
        const long region = vofSetOffsets(T, vofJobGrown1);
        if (static_cast<long>(bList_.extent(0)) < region)
          bList_ = LField(Kokkos::view_alloc(std::string("vof::block::batch_list"),
                                             Kokkos::WithoutInitializing),
                          region);
        vofBatchWorklist(T, bList_, bStart_, bEnd_);
        vofBatchPlic(T, bList_, bStart_, bEnd_);
        vofSetOffsets(T, [d](const VofBlockJob& J) { return vofJobFaces(J, d); });
        vofBatchFlux(T, d);
        if (anyMask)
          vofBatchBcLedger(T, d, pk_, kPk, kPkBcv);
        vofSetOffsets(T, vofJobInner);
        vofBatchUpdate(T, d);
        batchGhostFill(T);
      }
    }
    for (std::size_t k = 0; k < nb; ++k)
      mb[k]->adv_.noteBatchedStep();
  }
 
  void debrisBatched(std::vector<VofBlock*>& bs) {
    const std::size_t nb = bs.size();
    if (nb == 0)
      return;
    VofDebrisFlags F;
    F.doDebris = csfEnabled && debrisCensus;
    F.doResidue = wispEps > 0.0;
    F.remove = debrisRemove;
    F.vaMin = kDebrisMinAttached;
    // the predicates of debrisPassBatchImpl (review finding 7 there: the three lists are disjoint)
    const double weps = wispEps, cfull = kDebrisMaxC;
    const double ieps = curvProto.interfaceEps > weps ? curvProto.interfaceEps : weps;
    std::vector<VofBlockJob> jobs(nb);
    for (std::size_t k = 0; k < nb; ++k) {
      bs[k]->st_.debrisCells = 0;
      bs[k]->st_.debrisVolume = 0.0;
      jobs[k] = makeJob(*bs[k], 0.0);
    }
    if (static_cast<std::size_t>(dCap_.extent(0)) < nb) {
      for (int q = 0; q < 6; ++q)
        dRange_[q] = LField("vof::block::batch_debris_range", nb);
      dCap_ = Kokkos::View<int*, SMem>("vof::block::batch_debris_cap", nb);
    }
    ensurePacket(nb);
    Kokkos::deep_copy(SExec(), dCap_, 0);
    const LField sD = dRange_[0], eD = dRange_[1], sR = dRange_[2], eR = dRange_[3],
                 sA = dRange_[4], eA = dRange_[5];
    for (std::size_t j0 = 0; j0 < nb; j0 += kVofBlockBatch) {
      VofBlockTable T = batchTable(jobs, j0);
      const long region = vofSetOffsets(T, vofJobInner);
      auto ensure = [&](LField& f, const char* name) {
        if (static_cast<long>(f.extent(0)) < region)
          f = LField(Kokkos::view_alloc(std::string(name), Kokkos::WithoutInitializing), region);
      };
      if (F.doDebris) {
        ensure(listD_, "vof::block::debris_list");
        vofBatchMark(T, 0, ieps, cfull, weps, listD_, sD, eD);
      }
      if (F.doResidue) {
        ensure(listR_, "vof::block::residue_list");
        vofBatchMark(T, 1, ieps, cfull, weps, listR_, sR, eR);
      }
      ensure(listA_, "vof::block::attached_list");
      vofBatchMark(T, 2, ieps, cfull, weps, listA_, sA, eA);
      vofBatchDebrisSums(T, F, listD_, listR_, listA_, sD, eD, sR, eR, sA, eA, pk_, kPk, 0);
      if (static_cast<long>(debrisEx_.extent(0)) < region)
        debrisEx_ = SField(
            Kokkos::view_alloc("vof::block::debris_excess", Kokkos::WithoutInitializing), region);
      vofBatchDebrisWrite(T, F, listD_, listR_, listA_, sD, sR, sA, eA, pk_, kPk, 0, debrisEx_,
                          dCap_);
      vofBatchDebrisLost(T, sA, eA, pk_, kPk, 0, debrisEx_, dCap_);
      batchGhostFill(T);  // the inner colour moved; the next reader expects valid ghosts
    }
  }
 
  void applyDebrisLedger(std::vector<VofBlock*>& bs, bool doDebris) {
    for (std::size_t k = 0; k < bs.size(); ++k) {
      VofBlockStats& st = bs[k]->st_;
      const double* ho = &pkHost_(kPk * k);
      const long nD = static_cast<long>(ho[kDbND]);
      const bool removeD = doDebris && debrisRemove && nD > 0;
      const bool acted = ho[kDbActed] != 0.0;
      if (ho[kDbPending] == 0.0) {
        if (doDebris && nD > 0) {  // census only: report what WOULD be removed
          st.debrisCells = nD;
          st.debrisVolume = ho[kDbSumD];
        }
        continue;  // nothing to remove: nothing else ran
      }
      if (doDebris && !debrisRemove && nD > 0) {  // census of what WOULD be removed
        st.debrisCells = nD;
        st.debrisVolume = ho[kDbSumD];
      }
      if (!acted) {  // step 3: no attached interface worth a cell -- the clip bounds it
        if (removeD)
          st.debrisUnresolved += nD;
        continue;
      }
      if (removeD) {
        st.debrisCells = nD;
        st.debrisVolume = ho[kDbSumD];
        st.debrisReturned += ho[kDbSumD];
      }
      st.residueReturned += ho[kDbSumR];
      st.debrisLost += ho[kDbDV];  // `lost`, written over the dV slot by the act
    }
  }
 
  static constexpr int kPkBox = kDbSlots, kPkBcv = kDbSlots + 6, kPk = kDbSlots + 12;
  void ensurePacket(std::size_t nb) {
    if (static_cast<std::size_t>(pk_.extent(0)) < kPk * nb) {
      pk_ = SField("vof::block::batch_packet", kPk * nb);
      pkHost_ = Kokkos::create_mirror_view(pk_);
    }
  }
 
  void recentreBatched(bool ranDebris) {
    std::vector<std::size_t> idx;
    std::vector<VofBlock*> mb;
    for (std::size_t k = 0; k < blocks_.size(); ++k)
      if (blocks_[k].mine_) {
        idx.push_back(k);
        mb.push_back(&blocks_[k]);
      }
    const std::size_t nb = idx.size();
    if (nb == 0)
      return;
    std::vector<VofBlockJob> jobs(nb);
    for (std::size_t k = 0; k < nb; ++k)
      jobs[k] = makeJob(*mb[k], 0.0);
    ensurePacket(nb);
    const double eps = bubbleEps > wispEps ? bubbleEps : wispEps;  // bubbleBox's threshold (§12.2)
    for (std::size_t j0 = 0; j0 < nb; j0 += kVofBlockBatch)
      vofBatchBox(batchTable(jobs, j0), eps, pk_, kPk, kPkBox);
    Kokkos::deep_copy(pkHost_, pk_);  // read #2
    if (ranDebris)
      applyDebrisLedger(mb, csfEnabled && debrisCensus);
    for (std::size_t k = 0; k < nb; ++k) {  // the WO-R ledger (no block installs a mask today)
      if (jobs[k].outside == nullptr)
        continue;
      for (int f = 0; f < 6; ++f)
        if (jobs[k].bcOwn[f])
          mb[k]->adv_.addBcFaceVolume(f, pkHost_(kPk * k + kPkBcv + f));
    }
    for (std::size_t k = 0; k < nb; ++k) {
      const double* hb = &pkHost_(kPk * k + kPkBox);
      const int lo[3] = {static_cast<int>(hb[0]), static_cast<int>(hb[1]), static_cast<int>(hb[2])};
      const int hi[3] = {static_cast<int>(hb[3]), static_cast<int>(hb[4]), static_cast<int>(hb[5])};
      if (hi[0] < lo[0])
        continue;  // the block is empty (bubbleBox returned false)
      VofBlock& b = *mb[k];
      VofBox bb;
      for (int d = 0; d < 3; ++d) {
        bb.lo[d] = lo[d] + b.box.lo[d];
        bb.hi[d] = hi[d] + 1 + b.box.lo[d];
      }
      recentreTo(idx[k], bb);
    }
  }
 
 private:
  VofMetric metric_;  
  double h_ = 1.0;
  I3 gs_{0, 0, 0};
  std::array<bool, 3> per_{true, true, true};
  int rank_ = 0, size_ = 1;
  int margin_ =
      3;           
  int ghost_ = 3;  
  long step_ = 0;
  std::vector<VofBlock> blocks_;
  std::shared_ptr<VofBlockExchangeBase> exch_;
  std::map<std::array<int, 3>, std::vector<WyAdvector>> pool_;
  long poolHits_ = 0, poolMisses_ = 0;
  long lastReassigned_ = 0;
  VofCurvature::Stats curvStats_{};
  LField listD_, listA_, listR_;  
  SField debrisEx_;               
  SField debrisOut_;  
  SField::host_mirror_type debrisOutHost_;
  // C1 batched stages: per-job device scalars (indexed by the job's place among the master blocks)
  SField bRead1_;  
  SField::host_mirror_type bRead1Host_;
  std::vector<std::size_t> statsPend_;  
  LField bList_, bStart_, bEnd_;  
  LField dRange_[6];              
  Kokkos::View<int*, SMem> dCap_; 
  LField cListG_, cListI_, cStartG_, cEndG_, cStartI_, cEndI_;  
  Kokkos::View<PvMoments*, SMem> cCacheG_;  
  LField cCnt_;                   
  LField::host_mirror_type cCntHost_;
  SField pk_;                     
  SField::host_mirror_type pkHost_;
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_BLOCK_CONTAINER_HPP