curvature_field.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_CURVATURE_FIELD_HPP
#define PECLET_FLOW_VOF_CURVATURE_FIELD_HPP
 
#include <algorithm>
#include <chrono>
#include <Kokkos_Core.hpp>
#include <stdexcept>
#include <vector>
 
#include "mac_stencils.hpp"  // peclet::flow::SExec, SField, I3, L3
#include "policy.hpp"
#include "vof/advect_wy.hpp"  // wyIsMixed, wyReconstructCell
#include "vof/curvature.hpp"
#include "vof/surface_tension.hpp"  // csfKappaDefined
 
namespace peclet::flow::vof {
 
KOKKOS_INLINE_FUNCTION bool vofIsInterface(double c, double eps) {
  return c > eps && c < 1.0 - eps;
}
 
template <class SF>
KOKKOS_INLINE_FUNCTION void curvHeightCell(long i, SF c, SF mx, SF my, SF mz, SF al, SF kap, SF br,
                                           long s0, long s1, long s2, double mtol, double ptW,
                                           double ieps, bool forceFb, bool oneDir, bool useFit,
                                           VofMetric g, double peps) {
  const long st[3] = {s0, s1, s2};
  kap(i) = 0.0;
  if (!vofIsInterface(c(i), ieps)) {
    br(i) = static_cast<double>(kCurvNone);
    return;
  }
  if (forceFb) {
    br(i) = -1.0;
    return;
  }
  // Order the three column directions by the cell's own |n_d|. The normal is used as a
  // CATEGORICAL selector only — it is never differenced (see curvature.hpp).
  const double am[3] = {Kokkos::fabs(mx(i)), Kokkos::fabs(my(i)), Kokkos::fabs(mz(i))};
  int ord[3] = {0, 1, 2};
  for (int p = 1; p < 3; ++p)  // insertion sort, descending
    for (int q = p; q > 0 && am[ord[q]] > am[ord[q - 1]]; --q) {
      const int t = ord[q];
      ord[q] = ord[q - 1];
      ord[q - 1] = t;
    }
  const int ntry = oneDir ? 1 : 3;
  for (int t = 0; t < ntry; ++t) {
    const int d = ord[t];
    const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
    const long sd = st[d], s1 = st[d1], s2 = st[d2];
    double hh[9];
    int orient0 = 0;
    bool ok = true;
    for (int q = 0; q < 3 && ok; ++q)
      for (int p = 0; p < 3 && ok; ++p) {
        const long base = i + (p - 1) * s1 + (q - 1) * s2;
        double col[kHfColumn];
        for (int k = 0; k < kHfColumn; ++k)
          col[k] = c(base + (k - kHfColumn / 2) * sd);
        double h;
        int orient;
        if (!hfColumnHeight(col, kHfColumn, h, orient, mtol, peps)) {
          ok = false;
          break;
        }
        if (orient0 == 0)
          orient0 = orient;
        else if (orient != orient0) {
          ok = false;  // the nine columns do not describe one single-valued surface
          break;
        }
        hh[p + 3 * q] = h;
      }
    if (!ok)
      continue;
    // Phase 3 (V2.3): the heights are in cells along `d` and the patch spans cells along
    // `d1`/`d2`, so the metric turns them into the physical graph. Unit metric == today.
    kap(i) = hfPatchKappa(hh, g, d);
    br(i) = static_cast<double>(t == 0 ? kCurvHf : kCurvHfMixed);
    return;
  }
 
  // Tier 2b — the MIXED height function. No direction gives nine consistent columns; pool
  // the interface positions of whichever of the 27 columns DO close and fit a paraboloid
  // through them in the target's own frame. Every cell it reads is inside the same
  // 3x3x7-per-direction footprint tier 1 already used, so it costs no extra halo.
  if (useFit) {
    const double m0 = mx(i), m1 = my(i), m2 = mz(i);
    const double mi[3] = {m0, m1, m2};
    double nn[3] = {0.0, 0.0, 0.0};
    // PHYSICAL frame (Phase 3, V2.4); identity at the unit metric.
    if (vofPhysNormalInv(mi, g, nn) > 0.0) {
      double t1[3], t2[3];
      curvFrame(nn, t1, t2);
      double vtx[8][3], ctr[3], area;
      const int nvt = plicPolygon(m0, m1, m2, al(i), vtx);
      polygonAreaCentroid(vtx, nvt, ctr, area);
      const double org[3] = {ctr[0] - 0.5, ctr[1] - 0.5, ctr[2] - 0.5};
      PtFit pf;
      ptFitInit(pf);
      for (int d = 0; d < 3; ++d) {
        const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
        const long sd = st[d], s1 = st[d1], s2 = st[d2];
        for (int q = -1; q <= 1; ++q)
          for (int p = -1; p <= 1; ++p) {
            const long base = i + p * s1 + q * s2;
            double col[kHfColumn];
            for (int k = 0; k < kHfColumn; ++k)
              col[k] = c(base + (k - kHfColumn / 2) * sd);
            double hv;
            int orient;
            if (!hfColumnHeight(col, kHfColumn, hv, orient, mtol, peps))
              continue;
            double X[3];
            X[d1] = static_cast<double>(p);
            X[d2] = static_cast<double>(q);
            X[d] = orient * hv;  // the interface POSITION along d (h is the signed height)
            ptFitAdd(pf, X, org, t1, t2, nn, ptW, g);
          }
      }
      double a[6];
      bool red = false;
      if (pf.npt >= 6 && ptFitSolve(pf, a, red) && !red) {
        kap(i) = paraboloidKappa(a);
        br(i) = static_cast<double>(kCurvHfFit);
        return;
      }
    }
  }
  br(i) = -1.0;  // to the fallback
}
 
template <class SF>
KOKKOS_INLINE_FUNCTION void curvFallbackStore(long i, const PvFit& fit, SF kap, SF br) {
  double a[6];
  bool red = false;
  if (!pvFitSolve(fit, a, red)) {
    kap(i) = 0.0;
    br(i) = static_cast<double>(kCurvNoEstimate);
    return;
  }
  kap(i) = paraboloidKappa(a);
  br(i) = static_cast<double>(red ? kCurvPvReduced : kCurvPv);
}
 
// ---- tier 3 at cost (design G, `doc/vof_curvature_cost_design.md` §5.2-§5.4) -----------------
 
#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
#define PECLET_VOF_DEVICE_NOINLINE __noinline__
#else
#define PECLET_VOF_DEVICE_NOINLINE
#endif
 
KOKKOS_INLINE_FUNCTION PECLET_VOF_DEVICE_NOINLINE void curvMomentsBuild(double mx, double my,
                                                                        double mz, double alpha,
                                                                        const VofMetric& g,
                                                                        PvMoments& P) {
  pvMomentsBuild(mx, my, mz, alpha, g, P);
}
//
// Each interfacial cell's PLIC polygon data is built ONCE per curvature pass, in the planes pass,
// into a per-cell cache indexed by a slot map, and every target reads it (the polygon used to be
// rebuilt by each of the ~16 targets whose 5^3 stencil holds the cell: 63 % of tier 3, and then
// re-integrated in each target's frame: another 30 %). The cache holds the polygon's 3-D area
// moments (core `PvMoments`), and each target TRANSFORMS them into its frame (`pvTermMoments`):
// the same integrals as the Green's-theorem terms, a RECORDED change at round-off (WO-3). The
// tier-3 body is one template over an ENTRY ACCESSOR, so the compacted path (`VofPvCached`: the
// slot map + the cache) and the dense compaction oracle (`VofPvOnTheFly`: the entry built from the
// planes on the spot) share one body; both hand out bit-identical entries because both run
// `pvMomentsBuild` on the same plane and metric.
 
struct VofPvCached {
  const int* slot;
  const PvMoments* cache;
  KOKKOS_INLINE_FUNCTION bool present(long j) const { return slot[j] >= 0; }
  KOKKOS_INLINE_FUNCTION const PvMoments& entry(long j) const { return cache[slot[j]]; }
};
 
template <class SF>
struct VofPvOnTheFly {
  SF c, mx, my, mz, al;
  double ieps;
  VofMetric g;
  KOKKOS_INLINE_FUNCTION bool present(long j) const { return vofIsInterface(c(j), ieps); }
  KOKKOS_INLINE_FUNCTION PvMoments entry(long j) const {
    PvMoments P;
    curvMomentsBuild(mx(j), my(j), mz(j), al(j), g, P);
    return P;
  }
};
 
template <class SF, class Acc>
KOKKOS_INLINE_FUNCTION void curvFallbackCell(long i, const Acc& acc, SF kap, SF br, long sy,
                                             long sz, int gr, double dW, double cmin, VofMetric g) {
  if (br(i) >= 0.0)
    return;
 
  double nn[3], t1[3], t2[3], org[3];
  if (!pvFrameMoments(acc.entry(i), g, nn, t1, t2, org)) {
    kap(i) = 0.0;
    br(i) = static_cast<double>(kCurvNoEstimate);
    return;
  }
 
  PvFit fit;
  pvFitInit(fit);
  for (int oz = -gr; oz <= gr; ++oz)
    for (int oy = -gr; oy <= gr; ++oy)
      for (int ox = -gr; ox <= gr; ++ox) {
        const long j = i + (long)ox + (long)oy * sy + (long)oz * sz;
        if (!acc.present(j))
          continue;
        const double off[3] = {static_cast<double>(ox), static_cast<double>(oy),
                               static_cast<double>(oz)};
        PvTerm t;
        if (pvTermMoments(t, acc.entry(j), off, org, t1, t2, nn, dW, cmin, g))
          pvFitAccumLower(fit, t);
      }
  curvFallbackStore(i, fit, kap, br);
}
 
using VofPvTerms = Kokkos::View<PvTerm*, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged>;
using VofPvInts = Kokkos::View<int*, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged>;
using VofPvDbls = Kokkos::View<double*, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged>;
struct VofPvTeamScratch {
  VofPvInts kidx;
  VofPvTerms term;
  VofPvDbls ent;
  static std::size_t shmem(int nk) {
    return VofPvInts::shmem_size(nk) + VofPvTerms::shmem_size(nk) +
           VofPvDbls::shmem_size(kPvEntries);
  }
  template <class Member>
  KOKKOS_INLINE_FUNCTION VofPvTeamScratch(const Member& tm, int nk)
      : kidx(tm.team_scratch(0), nk)
      , term(tm.team_scratch(0), nk)
      , ent(tm.team_scratch(0), kPvEntries) {}
};
 
template <class Member, class SF, class Acc>
KOKKOS_INLINE_FUNCTION PECLET_VOF_DEVICE_NOINLINE void curvFallbackTeam(
    const Member& tm, const VofPvTeamScratch& sc, long i, const Acc& acc, SF kap, SF br, long sy,
    long sz, int gr, double dW, double cmin, VofMetric g) {
  if (br(i) >= 0.0)  // uniform across the team: nothing writes br(i) before the barriers below
    return;
  double nn[3], t1[3], t2[3], org[3];
  const bool framed = pvFrameMoments(acc.entry(i), g, nn, t1, t2, org);
  const int side = 2 * gr + 1, nk = side * side * side;
  int n = 0;
  if (framed) {
    Kokkos::parallel_scan(
        Kokkos::TeamThreadRange(tm, nk),
        [&](const int k, int& upd, const bool final) {
          const int ox = k % side - gr, oy = (k / side) % side - gr, oz = k / (side * side) - gr;
          const long j = i + (long)ox + (long)oy * sy + (long)oz * sz;
          if (acc.present(j)) {
            if (final)
              sc.kidx(upd) = k;
            ++upd;
          }
        },
        n);
    tm.team_barrier();
    Kokkos::parallel_for(Kokkos::TeamThreadRange(tm, n), [&](const int p) {
      const int k = sc.kidx(p);
      const int ox = k % side - gr, oy = (k / side) % side - gr, oz = k / (side * side) - gr;
      const long j = i + (long)ox + (long)oy * sy + (long)oz * sz;
      const double off[3] = {static_cast<double>(ox), static_cast<double>(oy),
                             static_cast<double>(oz)};
      PvTerm t;
      pvTermMoments(t, acc.entry(j), off, org, t1, t2, nn, dW, cmin, g);
      sc.term(p) = t;
    });
    tm.team_barrier();
    Kokkos::parallel_for(Kokkos::TeamThreadRange(tm, kPvEntries), [&](const int e) {
      double a = 0.0;
      for (int p = 0; p < n; ++p)
        if (sc.term(p).ok)
          pvFitAccumEntry(a, sc.term(p), e);
      sc.ent(e) = a;
    });
  }
  tm.team_barrier();  // the entries are in scratch, and every lane has read br(i)
  Kokkos::single(Kokkos::PerTeam(tm), [&]() {
    if (!framed) {
      kap(i) = 0.0;
      br(i) = static_cast<double>(kCurvNoEstimate);
      return;
    }
    int npoly = 0;
    for (int p = 0; p < n; ++p)
      npoly += sc.term(p).ok ? 1 : 0;
    double ent[kPvEntries];
    for (int e = 0; e < kPvEntries; ++e)
      ent[e] = sc.ent(e);
    PvFit fit;
    pvFitFromEntries(fit, ent, npoly);
    curvFallbackStore(i, fit, kap, br);
  });
}
 
#if defined(KOKKOS_ENABLE_HIP)
inline constexpr int kVofWarp = 64;
#else
inline constexpr int kVofWarp = 32;
#endif
 
struct VofRawField {
  double* p;
  KOKKOS_INLINE_FUNCTION double& operator()(long i) const { return p[i]; }
};
 
struct VofCurvFallbackJob {
  double *c, *mx, *my, *mz, *al, *kap, *br;
  const int* slot;         
  const PvMoments* cache;  
  const long* list;
  long n, sy, sz;
  int gr;
  double dW, cmin, ieps;
  VofMetric gm;
};
 
template <int N>
struct VofCurvFallbackTable {
  VofCurvFallbackJob job[N];
  long off[N + 1];
  int nj;
};
 
class VofCurvature {
 public:
  using IField = Kokkos::View<int*, SMem>;
  using MomCache = Kokkos::View<PvMoments*, SMem>;
  struct Stats {
    long interfacial = 0;  
    long hf = 0;           
    long hfMixed = 0;      
    long hfFit = 0;        
    long pv = 0;           
    long pvReduced = 0;    
    long noEstimate = 0;   
    long clipped = 0;
  };
 
  void init(int nx, int ny, int nz, int ghost) {
    if (ghost < kHfColumn / 2)
      throw std::invalid_argument(
          "peclet::flow::vof::VofCurvature: ghost width must be >= 3 (the height-function column "
          "reaches 3 cells and the 5^3 PV stencil's MYC normals reach 3 cells)");
    if (nx < 1 || ny < 1 || nz < 1)
      throw std::invalid_argument("peclet::flow::vof::VofCurvature: empty block");
    n_ = I3{nx, ny, nz};
    g_ = ghost;
    e_ = I3{nx + 2 * ghost, ny + 2 * ghost, nz + 2 * ghost};
    len_ = static_cast<long>(e_.x) * e_.y * e_.z;
    mx_ = SField("vof::curv::mx", len_);
    my_ = SField("vof::curv::my", len_);
    mz_ = SField("vof::curv::mz", len_);
    alpha_ = SField("vof::curv::alpha", len_);
    kappa_ = SField("vof::curv::kappa", len_);
    branch_ = SField("vof::curv::branch", len_);
    // WO-V9: the compaction lists. The grown one covers every cell `reconstructPlanes` writes; the
    // inner one every cell the cascade runs on. Sized at the full region (worst case: an entirely
    // interfacial block) so a scan can never overflow them.
    const long grown = (long)(n_.x + 2 * kPvHalf) * (n_.y + 2 * kPvHalf) * (n_.z + 2 * kPvHalf);
    listG_ = LField("vof::curv::listG", grown);
    listI_ = LField("vof::curv::listI", (long)n_.x * n_.y * n_.z);
    // Design G §5.2: the slot map over the extended block (read only inside the grown region,
    // where the planes pass writes it) and the per-cell moment cache, sized to listG_'s capacity
    // on the first worklist planes pass (`ensureCache_`; a block-container cascade whose passes are
    // batched never runs one and never allocates it).
    slot_ = IField("vof::curv::slot", len_);
    cache_ = MomCache();
  }
 
  bool ready() const { return kappa_.extent(0) != 0; }
 
  // ---- WO-V9: compaction, and the per-pass timers --------------------------------------------
  //
  // The cascade is the most divergent kernel in the VoF pipeline: on a resolved droplet the
  // interfacial cells are ~0.5 % of the block, so a dense `parallel_for` puts at most one or two
  // active lanes in a warp and the other thirty run the guard and idle for the whole height
  // function.  Compacting the interfacial cells into a contiguous list first makes those warps
  // full.  It is a pure re-ordering — every cell reads the same neighbours and writes the same
  // value — so the two paths are BIT-IDENTICAL and the ctest gates that.
  //
  // The timers follow `WyAdvector`'s rule exactly: fence at the boundaries only when armed.
  struct Timing {
    double planes = 0.0;    
    double height = 0.0;    
    double fallback = 0.0;  
    double census = 0.0;    
    double compact = 0.0;   
    double clip = 0.0;      
    long calls = 0;
  };
  bool useWorklist = true;
  bool timingOn = false;
  Timing tm;
  void resetTiming() { tm = Timing(); }
  const Timing& timing() const { return tm; }
  double tick_() const {
    if (!timingOn)
      return 0.0;
    Kokkos::fence();
    return std::chrono::duration<double>(std::chrono::steady_clock::now().time_since_epoch())
        .count();
  }
  void addT_(double& acc, double t0) {
    if (!timingOn)
      return;
    Kokkos::fence();
    acc +=
        std::chrono::duration<double>(std::chrono::steady_clock::now().time_since_epoch()).count() -
        t0;
  }
  I3 inner() const { return n_; }
  I3 extent() const { return e_; }
  int ghost() const { return g_; }
  SField kappa() const { return kappa_; }
  SField branch() const { return branch_; }
  SField planeM(int d) const { return d == 0 ? mx_ : (d == 1 ? my_ : mz_); }
  SField planeAlpha() const { return alpha_; }
  IField slotMap() const { return slot_; }
 
  // ---- tunables (all measured knobs, defaults are the literature values) ---------------------
  VofMetric metric;
  double weightWidth = kPvWeightWidth;
  double pureEps = 0.0;
  double monoTol = 1e-6;
  double ptWeightWidth = kPtWeightWidth;
  double cosMin = 0.2;
  double interfaceEps = 0.0;
  bool debugForceFallback = false;
  bool debugSingleDirection = false;
  bool useMixedHeightFit = false;
  double kappaMax = 0.0;
  void copyTunablesFrom(const VofCurvature& src) {
    metric = src.metric;
    weightWidth = src.weightWidth;
    pureEps = src.pureEps;
    monoTol = src.monoTol;
    ptWeightWidth = src.ptWeightWidth;
    cosMin = src.cosMin;
    interfaceEps = src.interfaceEps;
    debugForceFallback = src.debugForceFallback;
    debugSingleDirection = src.debugSingleDirection;
    useMixedHeightFit = src.useMixedHeightFit;
    kappaMax = src.kappaMax;
    useWorklist = src.useWorklist;
  }
  double kappaClipValue() const { return kappaMax < 0.0 ? 1.0 / metric.minH() : kappaMax; }
 
  Stats compute(SField c) {
    computeBegin(c);
    const double t3 = tick_();
    fallbackPass(c);
    addT_(tm.fallback, t3);
    return computeEnd();
  }
 
  void computeBegin(SField c) {
    if (!ready())
      throw std::runtime_error("peclet::flow::vof::VofCurvature::compute: init() not called");
    const double t0 = tick_();
    if (useWorklist)
      compact(c);
    addT_(tm.compact, t0);
    const double t1 = tick_();
    reconstructPlanes(c);
    addT_(tm.planes, t1);
    const double t2 = tick_();
    heightPass(c);
    addT_(tm.height, t2);
  }
  Stats computeEnd() {
    const double t5 = tick_();
    const long nclip = clipPass();
    addT_(tm.clip, t5);
    const double t4 = tick_();
    Stats s = census();
    addT_(tm.census, t4);
    s.clipped = nclip;
    ++tm.calls;
    return s;
  }
  VofCurvFallbackJob fallbackJob(SField c) const {
    VofCurvFallbackJob j;
    j.c = c.data();
    j.mx = mx_.data();
    j.my = my_.data();
    j.mz = mz_.data();
    j.al = alpha_.data();
    j.kap = kappa_.data();
    j.br = branch_.data();
    j.slot = slot_.data();
    j.cache = cache_.data();
    j.list = listI_.data();
    j.n = nI_;
    j.sy = e_.x;
    j.sz = static_cast<long>(e_.x) * e_.y;
    j.gr = kPvHalf;
    j.dW = weightWidth * metric.maxH();
    j.cmin = cosMin;
    j.ieps = interfaceEps;
    j.gm = metric;
    return j;
  }
  static constexpr int kFallbackBatch = 16;
  static void fallbackBatch(const std::vector<VofCurvFallbackJob>& jobs) {
    using Table = VofCurvFallbackTable<kFallbackBatch>;
    for (std::size_t j0 = 0; j0 < jobs.size(); j0 += kFallbackBatch) {
      Table T;
      T.nj = static_cast<int>(std::min<std::size_t>(kFallbackBatch, jobs.size() - j0));
      T.off[0] = 0;
      for (int k = 0; k < T.nj; ++k) {
        T.job[k] = jobs[j0 + k];
        T.off[k + 1] = T.off[k] + T.job[k].n;
      }
      for (int k = T.nj; k < kFallbackBatch; ++k)
        T.off[k + 1] = T.off[T.nj];
      if (T.off[T.nj] == 0)
        continue;
      if constexpr (!Kokkos::SpaceAccessibility<Kokkos::HostSpace,
                                                SField::memory_space>::accessible) {
        // C3 (§5.11): one team (a warp) per target on a device; the host keeps one thread per
        // target below.
        using Policy = Kokkos::TeamPolicy<SExec>;
        int gr = 0;
        for (int k = 0; k < T.nj; ++k)
          gr = std::max(gr, T.job[k].gr);
        const int nk = (2 * gr + 1) * (2 * gr + 1) * (2 * gr + 1);
        Policy pol(SExec(), T.off[T.nj], kVofWarp);
        pol.set_scratch_size(0, Kokkos::PerTeam(VofPvTeamScratch::shmem(nk)));
        Kokkos::parallel_for(
            "vof::curv::pv_batch_team", pol, KOKKOS_LAMBDA(const typename Policy::member_type& tm) {
              const long t = tm.league_rank();
              int b = 0;
              while (t >= T.off[b + 1])
                ++b;
              const VofCurvFallbackJob& J = T.job[b];
              const int side = 2 * J.gr + 1;
              const VofPvTeamScratch sc(tm, side * side * side);
              curvFallbackTeam(tm, sc, J.list[t - T.off[b]], VofPvCached{J.slot, J.cache},
                               VofRawField{J.kap}, VofRawField{J.br}, J.sy, J.sz, J.gr, J.dW,
                               J.cmin, J.gm);
            });
        continue;
      }
      Kokkos::parallel_for(
          "vof::curv::pv_batch", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
          KOKKOS_LAMBDA(long t) {
            int b = 0;
            while (t >= T.off[b + 1])
              ++b;
            const VofCurvFallbackJob& J = T.job[b];
            curvFallbackCell(J.list[t - T.off[b]], VofPvCached{J.slot, J.cache}, VofRawField{J.kap},
                             VofRawField{J.br}, J.sy, J.sz, J.gr, J.dW, J.cmin, J.gm);
          });
    }
  }
 
  void compact(SField c) {
    const I3 e = e_, n = n_;
    const int g = g_, gr = kPvHalf;
    const double ieps = interfaceEps;
    {
      const int rx = n.x + 2 * gr, ry = n.y + 2 * gr, rz = n.z + 2 * gr;
      const long region = (long)rx * ry * rz;
      LField list = listG_;
      long cnt = 0;
      Kokkos::parallel_scan(
          "vof::curv::compactG", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
          KOKKOS_LAMBDA(const long r, long& upd, const bool final) {
            const int ix = (int)(r % rx);
            const int iy = (int)((r / rx) % ry);
            const int iz = (int)(r / ((long)rx * ry));
            const long i = L3(g - gr + ix, g - gr + iy, g - gr + iz, e);
            if (vofIsInterface(c(i), ieps)) {
              if (final)
                list(upd) = i;
              ++upd;
            }
          },
          cnt);
      Kokkos::fence();
      nG_ = cnt;
    }
    {
      const long region = (long)n.x * n.y * n.z;
      const int nx = n.x, ny = n.y;
      LField list = listI_;
      long cnt = 0;
      Kokkos::parallel_scan(
          "vof::curv::compactI", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
          KOKKOS_LAMBDA(const long r, long& upd, const bool final) {
            const int ix = (int)(r % nx);
            const int iy = (int)((r / nx) % ny);
            const int iz = (int)(r / ((long)nx * ny));
            const long i = L3(g + ix, g + iy, g + iz, e);
            if (vofIsInterface(c(i), ieps)) {
              if (final)
                list(upd) = i;
              ++upd;
            }
          },
          cnt);
      Kokkos::fence();
      nI_ = cnt;
    }
  }
 
  // The four passes below are PUBLIC only because nvcc forbids 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"). They are implementation detail; call
  // `compute()`. `WyAdvector` is public for the same reason.
  //
  void reconstructPlanes(SField c) {
    const I3 e = e_, n = n_;
    const int g = g_, gr = kPvHalf;
    const long sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField mx = mx_, my = my_, mz = mz_, al = alpha_;
    const double ieps = interfaceEps;
    const VofMetric gm = metric;  // `g` is the ghost width in this scope
    if (useWorklist) {
      // The same two things the dense kernel does, as two coherent kernels: a branchless zeroing
      // sweep (a pure store, bandwidth-bound) and the MYC reconstruction over the compacted list.
      // Every cell ends with exactly the value the dense kernel would have written. Design G §5.2:
      // the zero sweep also clears the slot map, and the list kernel builds each cell's polygon
      // cache entry from the plane it just wrote and records its slot.
      ensureCache_();
      IField slot = slot_;
      MomCache cache = cache_;
      Kokkos::parallel_for(
          "vof::curv::planes_zero",
          MDRange3<SExec>(SExec(), {g - gr, g - gr, g - gr},
                          {g + n.x + gr, g + n.y + gr, g + n.z + gr}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const long i = L3(x, y, z, e);
            mx(i) = 0.0;
            my(i) = 0.0;
            mz(i) = 0.0;
            al(i) = 0.0;
            slot(i) = -1;
          });
      LField list = listG_;
      Kokkos::parallel_for(
          "vof::curv::planes_list", Kokkos::RangePolicy<SExec>(SExec(), 0, nG_),
          KOKKOS_LAMBDA(long t) {
            const long i = list(t);
            wyReconstructCell(c, i, sy, sz, mx, my, mz, al);
            curvMomentsBuild(mx(i), my(i), mz(i), al(i), gm, cache(t));
            slot(i) = static_cast<int>(t);
          });
      return;
    }
    Kokkos::parallel_for(
        "vof::curv::planes",
        MDRange3<SExec>(SExec(), {g - gr, g - gr, g - gr},
                        {g + n.x + gr, g + n.y + gr, g + n.z + gr}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          if (!vofIsInterface(c(i), ieps)) {
            mx(i) = 0.0;
            my(i) = 0.0;
            mz(i) = 0.0;
            al(i) = 0.0;
            return;
          }
          wyReconstructCell(c, i, sy, sz, mx, my, mz, al);
        });
  }
 
  void heightPass(SField c) {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long st[3] = {1, e_.x, static_cast<long>(e_.x) * e_.y};
    SField mx = mx_, my = my_, mz = mz_, al = alpha_, kap = kappa_, br = branch_;
    const double mtol = monoTol, ptW = ptWeightWidth * metric.maxH(), ieps = interfaceEps;
    const double peps = pureEps;
    const VofMetric gm = metric;  // `g` is the ghost width in this scope
    const bool forceFb = debugForceFallback, oneDir = debugSingleDirection,
               useFit = useMixedHeightFit;
    const long s0 = st[0], s1 = st[1], s2 = st[2];
    if (useWorklist) {
      // Reset every inner cell (branchless: a pure pair of stores), then run the cascade over the
      // compacted interfacial list.  A non-interfacial cell ends at kappa = 0, branch = kCurvNone
      // and an interfacial one at whatever the shared body writes -- the dense kernel's outcome,
      // cell for cell.
      Kokkos::parallel_for(
          "vof::curv::hf_reset", 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);
            kap(i) = 0.0;
            br(i) = static_cast<double>(kCurvNone);
          });
      LField list = listI_;
      Kokkos::parallel_for(
          "vof::curv::hf_list", Kokkos::RangePolicy<SExec>(SExec(), 0, nI_), KOKKOS_LAMBDA(long t) {
            curvHeightCell(list(t), c, mx, my, mz, al, kap, br, s0, s1, s2, mtol, ptW, ieps,
                           forceFb, oneDir, useFit, gm, peps);
          });
      return;
    }
    Kokkos::parallel_for(
        "vof::curv::hf", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          curvHeightCell(L3(x, y, z, e), c, mx, my, mz, al, kap, br, s0, s1, s2, mtol, ptW, ieps,
                         forceFb, oneDir, useFit, gm, peps);
        });
  }
 
  void fallbackPass(SField c) {
    const I3 e = e_, n = n_;
    const int g = g_, gr = kPvHalf;
    SField mx = mx_, my = my_, mz = mz_, al = alpha_, kap = kappa_, br = branch_;
    const double dW = weightWidth * metric.maxH(), cmin = cosMin, ieps = interfaceEps;
    const VofMetric gm = metric;  // `g` is the ghost width in this scope
    const long sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    if (useWorklist) {
      // Tier 3 is a SUBSET of the interfacial cells (those tier 1/2 could not serve), so the
      // interfacial list already removes the whole empty domain; the branch guard inside the body
      // removes the rest.
      LField list = listI_;
      const VofPvCached acc{slot_.data(), cache_.data()};
      Kokkos::parallel_for(
          "vof::curv::pv_list", Kokkos::RangePolicy<SExec>(SExec(), 0, nI_), KOKKOS_LAMBDA(long t) {
            curvFallbackCell(list(t), acc, kap, br, sy, sz, gr, dW, cmin, gm);
          });
      return;
    }
    const VofPvOnTheFly<SField> acc{c, mx, my, mz, al, ieps, gm};
    Kokkos::parallel_for(
        "vof::curv::pv", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          curvFallbackCell(L3(x, y, z, e), acc, kap, br, sy, sz, gr, dW, cmin, gm);
        });
  }
 
  long clipPass() {
    const double km = kappaClipValue();
    if (!(km > 0.0))
      return 0;
    const I3 e = e_, n = n_;
    const int g = g_;
    SField kap = kappa_, br = branch_;
    long cnt = 0;
    if (useWorklist) {
      LField list = listI_;
      Kokkos::parallel_reduce(
          "vof::curv::clip_list", Kokkos::RangePolicy<SExec>(SExec(), 0, nI_),
          KOKKOS_LAMBDA(long t, long& acc) {
            const long i = list(t);
            // `!(|k| <= km)` rather than `|k| > km`: a NaN curvature is clipped (to +km, the sign
            // of a NaN is meaningless) and counted, never silently kept.
            if (csfKappaDefined(br(i)) && !(Kokkos::fabs(kap(i)) <= km)) {
              kap(i) = (kap(i) == kap(i)) ? Kokkos::copysign(km, kap(i)) : km;
              ++acc;
            }
          },
          cnt);
    } else {
      Kokkos::parallel_reduce(
          "vof::curv::clip", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
          KOKKOS_LAMBDA(int x, int y, int z, long& acc) {
            const long i = L3(x, y, z, e);
            // `!(|k| <= km)` rather than `|k| > km`: a NaN curvature is clipped (to +km, the sign
            // of a NaN is meaningless) and counted, never silently kept.
            if (csfKappaDefined(br(i)) && !(Kokkos::fabs(kap(i)) <= km)) {
              kap(i) = (kap(i) == kap(i)) ? Kokkos::copysign(km, kap(i)) : km;
              ++acc;
            }
          },
          cnt);
    }
    Kokkos::fence();
    return cnt;
  }
 
  Stats census() const {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField br = branch_;
    Stats s;
    Kokkos::parallel_reduce(
        "vof::curv::census", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z, long& ni, long& n1, long& n2, long& n2b, long& n3,
                      long& n4, long& n5) {
          const int b = static_cast<int>(br(L3(x, y, z, e)));
          if (b == kCurvNone)
            return;
          ++ni;
          if (b == kCurvHf)
            ++n1;
          else if (b == kCurvHfMixed)
            ++n2;
          else if (b == kCurvHfFit)
            ++n2b;
          else if (b == kCurvPv)
            ++n3;
          else if (b == kCurvPvReduced)
            ++n4;
          else
            ++n5;
        },
        s.interfacial, s.hf, s.hfMixed, s.hfFit, s.pv, s.pvReduced, s.noEstimate);
    Kokkos::fence();
    return s;
  }
 
 private:
  I3 n_{0, 0, 0}, e_{0, 0, 0};
  int g_ = 0;
  long len_ = 0;
  SField mx_, my_, mz_, alpha_, kappa_, branch_;
  LField listG_, listI_;  // WO-V9: the compacted interfacial-cell lists
  long nG_ = 0, nI_ = 0;  // ... and their lengths from the last compact()
  IField slot_;           // design G §5.2: grown-list position of each interfacial cell, else -1
  MomCache cache_;        // ... and each one's polygon, in grown-list order
 
  void ensureCache_() {
    const long cap = static_cast<long>(listG_.extent(0));
    if (cap >= (1L << 31))
      throw std::runtime_error(
          "peclet::flow::vof::VofCurvature: grown list capacity >= 2^31 (the slot map is int)");
    if (static_cast<long>(cache_.extent(0)) < cap)
      cache_ = MomCache(
          Kokkos::view_alloc(std::string("vof::curv::cache"), Kokkos::WithoutInitializing), cap);
  }
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_CURVATURE_FIELD_HPP