advect_wy.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_ADVECT_WY_HPP
#define PECLET_FLOW_VOF_ADVECT_WY_HPP
 
#include <chrono>
#include <cstdio>
#include <functional>
#include <Kokkos_Core.hpp>
#include <stdexcept>
 
#include "mac_stencils.hpp"  // peclet::flow::SExec, SField, SMem, I3, L3
#include "policy.hpp"
#include "vof/cutcell.hpp"
#include "vof/plic.hpp"
#include "vof/wetting.hpp"
 
namespace peclet::flow::vof {
 
using UCField = Kokkos::View<unsigned char*, SMem>;
using LField = Kokkos::View<long*, SMem>;
 
KOKKOS_INLINE_FUNCTION bool wyIsMixed(double c) {
  return c > 0.0 && c < 1.0;
}
 
KOKKOS_INLINE_FUNCTION bool wyIsMixed(double c, double eps) {
  return c > eps && c < 1.0 - eps;
}
 
KOKKOS_INLINE_FUNCTION bool wyColourJump(double a, double b, double eps) {
  const double d = a - b;
  return (d > eps) || (d < -eps);
}
 
inline constexpr int kWySweepPerm[6][3] = {{0, 1, 2}, {1, 2, 0}, {2, 0, 1},
                                           {0, 2, 1}, {2, 1, 0}, {1, 0, 2}};
 
template <class SF>
KOKKOS_INLINE_FUNCTION void wyReconstructCell(const SF& c, long i, long sy, long sz, SF mx, SF my,
                                              SF mz, SF alpha) {
  double st[27];
  for (int kk = -1; kk <= 1; ++kk)
    for (int jj = -1; jj <= 1; ++jj)
      for (int ii = -1; ii <= 1; ++ii)
        st[plicSt(ii + 1, jj + 1, kk + 1)] = c(i + ii + jj * sy + kk * sz);
  double m[3];
  mycNormal(st, m);
  mx(i) = m[0];
  my(i) = m[1];
  mz(i) = m[2];
  alpha(i) = plicAlpha(m[0], m[1], m[2], c(i));
}
 
template <class SF>
KOKKOS_INLINE_FUNCTION double wyFaceFlux(double a, long p, long sd, int dir, const SF& c,
                                         const SF& mx, const SF& my, const SF& mz, const SF& alpha,
                                         double eps = 0.0) {
  if (a > 0.0) {  // donor is p; the outflow slab is the |a|-thick layer at its + face
    const double cd = c(p);
    return wyIsMixed(cd, eps) ? plicSlabVolume(mx(p), my(p), mz(p), alpha(p), dir, 1.0 - a, 1.0)
                              : cd * a;
  }
  if (a < 0.0) {  // donor is p + sd; the outflow slab is the |a|-thick layer at its - face
    const long q = p + sd;
    const double cd = c(q), aa = -a;
    return -(wyIsMixed(cd, eps) ? faceFluxVolume(mx(q), my(q), mz(q), alpha(q), dir, aa) : cd * aa);
  }
  return 0.0;
}
 
template <class SF, class MF>
KOKKOS_INLINE_FUNCTION double wyFaceFluxBc(double a, long p, long sd, int dir, const SF& c,
                                           const SF& mx, const SF& my, const SF& mz,
                                           const SF& alpha, const MF& outside, double eps = 0.0) {
  if (a > 0.0) {  // donor is p
    const double cd = c(p);
    return (wyIsMixed(cd, eps) && !outside(p))
               ? plicSlabVolume(mx(p), my(p), mz(p), alpha(p), dir, 1.0 - a, 1.0)
               : cd * a;
  }
  if (a < 0.0) {  // donor is p + sd
    const long q = p + sd;
    const double cd = c(q), aa = -a;
    return -((wyIsMixed(cd, eps) && !outside(q))
                 ? faceFluxVolume(mx(q), my(q), mz(q), alpha(q), dir, aa)
                 : cd * aa);
  }
  return 0.0;
}
 
class WyAdvector {
 public:
  struct Diagnostics {
    double sumC = 0.0;  
    double minC = 0.0;  
    double maxC = 0.0;  
    long mixed = 0;     
    long wisps = 0;     
    // --- rung V5a (WO-Q), all zero unless cut-cell geometry is attached -----------------------
    double volume = 0.0;     
    double rawVolume = 0.0;  
    double solidFillSum = 0.0;
    double solidSumC = 0.0;      
    double minCFluid = 0.0;      
    double maxCFluid = 0.0;      
    double clippedVolume = 0.0;  
    double clippedSigned = 0.0;  
    long cutCells = 0;           
    long solidCells = 0;         
    long clampedFaces = 0;       
  };
 
  void init(int nx, int ny, int nz, double h, int ghost = 3) {
    if (ghost < 2)
      throw std::invalid_argument("peclet::flow::vof::WyAdvector: ghost width must be >= 2");
    if (nx < 1 || ny < 1 || nz < 1)
      throw std::invalid_argument("peclet::flow::vof::WyAdvector: empty block");
    n_ = I3{nx, ny, nz};
    g_ = ghost;
    h_ = h;
    e_ = I3{nx + 2 * ghost, ny + 2 * ghost, nz + 2 * ghost};
    len_ = static_cast<long>(e_.x) * e_.y * e_.z;
    c_ = SField("vof::C", len_);
    mx_ = SField("vof::mx", len_);
    my_ = SField("vof::my", len_);
    mz_ = SField("vof::mz", len_);
    alpha_ = SField("vof::alpha", len_);
    flux_ = SField("vof::flux", len_);
    uf_ = SField("vof::uf", len_);
    vf_ = SField("vof::vf", len_);
    wf_ = SField("vof::wf", len_);
    cc_ = UCField("vof::cc", len_);
    listCap_ = static_cast<long>(n_.x + 2) * (n_.y + 2) * (n_.z + 2);
    list_ = LField("vof::worklist", listCap_);
  }
 
  // ---- geometry / storage ------------------------------------------------------------------
  I3 inner() const { return n_; }
  I3 extent() const { return e_; }
  int ghost() const { return g_; }
  double h() const { return h_; }
  long size() const { return len_; }
  long index(int x, int y, int z) const { return L3(x + g_, y + g_, z + g_, e_); }
 
  SField colour() const { return c_; }
  SField faceU() const { return uf_; }
  SField faceV() const { return vf_; }
  SField faceW() const { return wf_; }
  SField faceVel(int d) const { return d == 0 ? uf_ : (d == 1 ? vf_ : wf_); }
 
  // ---- rung V5a (WO-Q): optional cut-cell geometry ------------------------------------------
  //
  // OFF by default and EVERY kernel branches on `hasGeometry()` OUTSIDE its lambda, so with no
  // geometry attached the arithmetic executed is literally the V1 expression (not a multiplication
  // by 1.0) and the whole V1/V2a/V3/V4 battery stays byte-identical. See `vof/cutcell.hpp` for the
  // three rules the geometry branch adds and why the conserved functional is `sum eps_eff C`.
  //
  // USAGE (the solver does exactly this in `buildVofBlock`):
  //   a.enableGeometry();
  //   <fill faceOpenness(d) with the +d face openness and epsFraction() with buildCellFraction,
  //    on the INNER region, then run the block's own ghost policy on each of them>
  //   a.classifyGeometry();            // writes kindDouble() over every index except the low plane
  //   <run the block's ghost policy on kindDouble()>   // makes the low plane / ghosts the owner's
  //   a.finalizeGeometry();            // kindDouble() -> the unsigned-char classification
  bool hasGeometry() const { return hasGeom_; }
  void enableGeometry() {
    if (len_ <= 0)
      throw std::runtime_error("peclet::flow::vof::WyAdvector: enableGeometry before init()");
    if (!hasGeom_) {
      for (int d = 0; d < 3; ++d)
        of_[d] = SField("vof::openness", len_);
      eps_ = SField("vof::eps", len_);
      kindD_ = SField("vof::kindD", len_);
      kind_ = UCField("vof::kind", len_);
      fill_ = UCField("vof::fillState", len_);
      mark_ = UCField("vof::fillMark", len_);
      hasGeom_ = true;
    }
  }
  void disableGeometry() { hasGeom_ = false; }
  SField faceOpenness(int d) const { return of_[d]; }
  SField epsFraction() const { return eps_; }
  SField kindDouble() const { return kindD_; }
  UCField cellKind() const { return kind_; }
  bool cutFluxClamp = true;
  long clampedFaces() const { return clampedFaces_; }
  double clippedVolume() const { return clippedVolume_; }
  double clippedSigned() const { return clippedSigned_; }
 
  void classifyGeometry() {
    if (!hasGeom_)
      return;
    const I3 e = e_;
    const long sx = 1, sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField ox = of_[0], oy = of_[1], oz = of_[2], ep = eps_, kd = kindD_;
    Kokkos::parallel_for(
        "vof::wy::classify", MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          if (x == 0 || y == 0 || z == 0) {
            kd(i) = 0.0;  // the exchange owns this plane
            return;
          }
          const double oLo[3] = {ox(i - sx), oy(i - sy), oz(i - sz)};
          const double oHi[3] = {ox(i), oy(i), oz(i)};
          kd(i) = vofIsSolidCell(ep(i), oLo, oHi) ? 1.0 : 0.0;
        });
    Kokkos::fence();
  }
  void finalizeGeometry() {
    if (!hasGeom_)
      return;
    SField kd = kindD_;
    UCField kk = kind_;
    Kokkos::parallel_for(
        "vof::wy::finalize_geom", Kokkos::RangePolicy<SExec>(SExec(), 0, len_),
        KOKKOS_LAMBDA(long i) { kk(i) = kd(i) > 0.5 ? kVofSolid : kVofFluid; });
    Kokkos::fence();
  }
 
  // ---- rung V5b (WO-S): the theta-consistent band fill ---------------------------------------
  //
  // OFF by default: `hasWetting()` false runs WO-Q's neutral pass 1 verbatim, so the whole V5a
  // battery stays byte-identical. When on, PASS 1 ONLY is replaced (passes 2-3 and the shrinking
  // depth budget are untouched) by the theta-plane of `vof/wetting.hpp`.
  //
  // USAGE (the solver does exactly this):
  //   a.enableWetting();
  //   <fill wallSdf() and contactAngle() on the INNER region, then run the block's ghost policy>
  //   ... per fill:  a.buildWettingNormals();  <exchange wettingNormal(0..2)>;  a.solidBandFill();
  //
  // Every read the theta pass makes is owner-correct out to ghost depth 3 (the exchanged sdf, the
  // exchanged colour, the exchanged fluid-only normals), which is what keeps the INNER result
  // decomposition-independent — the same argument WO-Q finding 5 makes for the classification.
  bool hasWetting() const { return hasWet_; }
  void enableWetting() {
    if (!hasGeom_)
      throw std::runtime_error("peclet::flow::vof::WyAdvector: enableWetting needs geometry");
    if (!hasWet_) {
      sdfB_ = SField("vof::wallSdf", len_);
      thetaB_ = SField("vof::contactAngle", len_);
      for (int d = 0; d < 3; ++d)
        mfl_[d] = SField("vof::fluidNormal", len_);
      appB_ = SField("vof::apparentAngle", len_);
      wetB_ = UCField("vof::wetBranch", len_);
      hasWet_ = true;
    }
  }
  void disableWetting() { hasWet_ = false; }
  SField wallSdf() const { return sdfB_; }
  SField contactAngle() const { return thetaB_; }
  SField wettingNormal(int d) const { return mfl_[d]; }
 
  VofMetric metric;
  int wettingPivot = kVofPivotVolume;
  double wettingTangentEps = 1e-6;
  double wettingPureEps = 1e-8;
 
  void buildWettingNormals() {
    if (!hasWet_)
      return;
    const I3 e = e_, n = n_;
    const int g = g_;
    SField c = c_, mx = mfl_[0], my = mfl_[1], mz = mfl_[2];
    UCField kk = kind_;
    Kokkos::parallel_for(
        "vof::wy::wetting_normals",
        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);
          mx(i) = 0.0;
          my(i) = 0.0;
          mz(i) = 0.0;
          if (kk(i) != kVofFluid)
            return;
          double c27[27];
          unsigned char fl[27];
          for (int kz = -1; kz <= 1; ++kz)
            for (int ky = -1; ky <= 1; ++ky)
              for (int kx = -1; kx <= 1; ++kx) {
                const int q = plicSt(kx + 1, ky + 1, kz + 1);
                const long j = L3(x + kx, y + ky, z + kz, e);
                c27[q] = c(j);
                fl[q] = kk(j) == kVofFluid ? 1u : 0u;
              }
          double m[3];
          if (!youngsNormalFluidOnly(c27, fl, m))
            return;
          mx(i) = m[0];
          my(i) = m[1];
          mz(i) = m[2];
        });
    Kokkos::fence();
  }
 
  void wettingCensus(long counts[kVofWetCount], double& meanAppDeg, long& nApp) const {
    wettingCensusGhost(counts, meanAppDeg, nApp, 0);
  }
  void wettingCensusGhost(long counts[kVofWetCount], double& meanAppDeg, long& nApp,
                          int ghost) const {
    for (int b = 0; b < kVofWetCount; ++b)
      counts[b] = 0;
    meanAppDeg = 0.0;
    nApp = 0;
    if (!hasWet_)
      return;
    const I3 e = e_, n = n_;
    const int g = g_;
    UCField wb = wetB_, kk = kind_;
    SField ap = appB_;
    using MD = MDRange3<SExec>;
    const int gw = ghost < 0 ? 0 : (ghost > g ? g : ghost);
    MD pol(SExec(), {g - gw, g - gw, g - gw}, {g + n.x + gw, g + n.y + gw, g + n.z + gw});
    for (int b = 0; b < kVofWetCount; ++b) {
      long acc = 0;
      const unsigned char bb = static_cast<unsigned char>(b);
      Kokkos::parallel_reduce(
          "vof::wy::wet_census", pol,
          KOKKOS_LAMBDA(int x, int y, int z, long& a) {
            const long i = L3(x, y, z, e);
            if (kk(i) == kVofSolid && wb(i) == bb)  // SOLID cells only: the band is the write set
              ++a;
          },
          acc);
      counts[b] = acc;
    }
    double sum = 0.0;
    long cnt = 0;
    Kokkos::parallel_reduce(
        "vof::wy::wet_apparent", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& a, long& c2) {
          const long i = L3(x, y, z, e);
          if (wb(i) == static_cast<unsigned char>(kVofWetTheta)) {
            a += ap(i);
            ++c2;
          }
        },
        sum, cnt);
    Kokkos::fence();
    nApp = cnt;
    meanAppDeg = cnt ? sum / static_cast<double>(cnt) * (180.0 / 3.14159265358979323846) : 0.0;
  }
 
  void solidBandFill() {
    if (!hasGeom_)
      return;
    {
      UCField fs = fill_, mk = mark_, kk = kind_;
      Kokkos::parallel_for(
          "vof::wy::fill_reset", Kokkos::RangePolicy<SExec>(SExec(), 0, len_),
          KOKKOS_LAMBDA(long i) {
            fs(i) = kk(i) == kVofSolid ? kVofFillNone : kVofFillFluid;
            mk(i) = 0u;
          });
    }
    if (hasWet_) {
      UCField wb = wetB_;
      Kokkos::parallel_for(
          "vof::wy::wet_reset", Kokkos::RangePolicy<SExec>(SExec(), 0, len_),
          KOKKOS_LAMBDA(long i) { wb(i) = static_cast<unsigned char>(kVofWetNone); });
    }
    for (int k = 1; k <= 3; ++k) {
      if (hasWet_ && k == 1)
        solidBandFillPassWetting();
      else
        solidBandFillPass(k);
    }
    Kokkos::fence();
  }
 
  // ---- hooks -------------------------------------------------------------------------------
  std::function<void(SField)> exchange;
  std::function<double(double)> globalMax;
 
  // ---- rung V-BC (WO-R) hooks ----------------------------------------------------------------
  void setOutsideMask(UCField m) { outside_ = m; }
  UCField outsideMask() const { return outside_; }
  bool hasOutsideMask() const { return outside_.extent(0) == static_cast<std::size_t>(len_); }
 
  double bcFaceVolume(int f) const { return bcVol_[f]; }
  void resetBcFaceVolume() {
    for (int f = 0; f < 6; ++f)
      bcVol_[f] = 0.0;
  }
  void setBcFaceOwned(int f, bool on) { bcOwn_[f] = on ? 1u : 0u; }
 
  // ---- WO-V9: per-kernel timing (OFF by default, and inert when off) ------------------------
  //
  // A phase boundary on a device backend is meaningless without a fence — queued kernels would be
  // billed to whichever phase happens to call the clock next — so every boundary here fences,
  // exactly as `IbmSolver::phaseTick` already does for the step's three coarse phases. That fence
  // is a synchronisation, not a computation: it changes WHEN work happens, never WHAT it computes,
  // and the gate on this instrument is that a run with `timingOn` true and one with it false agree
  // bit for bit. When it is false, `tick_()` returns 0.0 without fencing and `addT_` returns
  // without reading the clock, so the instrumented build is the uninstrumented one plus a
  // predictable branch.
  struct Timing {
    double freeze = 0.0;       
    double reconstruct = 0.0;  
    double fluxes = 0.0;       
    double sweep = 0.0;        
    double clip = 0.0;         
    double exchange = 0.0;     
    long sweeps = 0;           
  };
  bool timingOn = false;
  Timing tm;
  void resetTiming() { tm = Timing(); }
  const Timing& timing() const { return tm; }
  static double wallSeconds() {
    return std::chrono::duration<double>(std::chrono::steady_clock::now().time_since_epoch())
        .count();
  }
  double tick_() const {
    if (!timingOn)
      return 0.0;
    Kokkos::fence();
    return wallSeconds();
  }
  void addT_(double& acc, double t0) {
    if (!timingOn)
      return;
    Kokkos::fence();
    acc += wallSeconds() - t0;
  }
  template <class F>
  void timedExchange(F&& fn) {
    const double t0 = tick_();
    fn();
    addT_(tm.exchange, t0);
  }
  void exchangeTimed(SField f) {
    timedExchange([&] { exchange(f); });
  }
 
  bool useWorklist = true;
  double cflLimit = 0.25;
  bool interfaceLocalCfl = false;
  bool debugRecomputeDilation = false;
  double wispEps = 0.0;
 
  void syncGhosts() {
    requireExchange();
    exchange(c_);
  }
 
  void advect(double dt, long step) {
    requireExchange();
    const double dth = dt / h_;
    const double cflLocal =
        interfaceLocalCfl ? (hasGeom_ ? maxCourantInterfaceCut(dth) : maxCourantInterface(dth))
                          : maxCourant(dth);
    enforceCfl(globalMax ? globalMax(cflLocal) : cflLocal, dt);
 
    // (1) THE dilation flag: frozen ONCE from C^n, used unchanged by all three sweeps.
    freezeDilationFlag();
    if (hasGeom_) {  // rung V5a: the clip / clamp census is per STEP, not per sweep
      clippedVolume_ = 0.0;
      clippedSigned_ = 0.0;
      clampedFaces_ = 0;
    }
 
    // (2) three directional sweeps
    const int* perm = kWySweepPerm[static_cast<int>(step % 6)];
    for (int s = 0; s < 3; ++s) {
      const int d = perm[s];
      if (debugRecomputeDilation && s > 0)
        freezeDilationFlag();  // the trap, on purpose — see `debugRecomputeDilation`
      reconstruct();
      computeFluxes(d, dth);
      applySweep(d, dth);
      if (hasGeom_)
        clipCutCells();  // rung V5a rule 3 (diagnostic clip, cut cells only)
      exchangeTimed(c_);
    }
    ++steps_;
  }
 
  double lastCfl() const { return lastCfl_; }
  long lastMixedCount() const { return mixedCount_; }
 
  // ---- rung V2b (WO-K) hooks -----------------------------------------------------------------
  // The momentum-consistent transport of `rho^c u_c` must be driven by the SAME PLIC planes, the
  // SAME face Courant numbers and the SAME sweep order as this colour advection. The planes are
  // re-reconstructed before EVERY sweep and overwritten in place, so a sibling advector that ran
  // after `advect()` returned would see only the last sweep's planes — the momentum sweeps have to
  // be INTERLEAVED with the colour sweeps. `MomentumConsistentAdvector` therefore drives the step
  // itself out of the public pieces below (`freezeDilationFlag` / `reconstruct` / `computeFluxes` /
  // `applySweep` / `exchange`) instead of calling `advect()`, and reads the planes through these
  // accessors between `computeFluxes` and `applySweep`, i.e. while they still describe the field
  // the sweep is about to move. Nothing here changes `advect()`.
 
  SField planeM(int d) const { return d == 0 ? mx_ : (d == 1 ? my_ : mz_); }
  SField planeAlpha() const { return alpha_; }
  UCField dilationFlag() const { return cc_; }
  SField faceFlux() const { return flux_; }
 
  double checkCourant(double dt) {
    const double dth = dt / h_;
    const double cflLocal =
        interfaceLocalCfl ? (hasGeom_ ? maxCourantInterfaceCut(dth) : maxCourantInterface(dth))
                          : maxCourant(dth);
    const double cfl = globalMax ? globalMax(cflLocal) : cflLocal;
    enforceCfl(cfl, dt);
    return cfl;
  }
 
  void enforceCfl(double cfl, double dt) {
    lastCfl_ = cfl;
    if (!(cfl <= cflLimit)) {
      char msg[256];
      std::snprintf(msg, sizeof(msg),
                    "peclet::flow::vof::WyAdvector: CFL = max|uf| dt/h = %.6g exceeds the "
                    "Weymouth-Yue boundedness cap %.6g (dt = %.6g, h = %.6g) - reduce dt",
                    cfl, cflLimit, dt, h_);
      throw std::runtime_error(msg);
    }
  }
 
  void noteBatchedStep() { ++steps_; }
  bool bcFaceOwned(int f) const { return bcOwn_[f] != 0; }
  void addBcFaceVolume(int f, double v) { bcVol_[f] += v; }
 
  Diagnostics diagnostics() const {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField c = c_;
    using MD = MDRange3<SExec>;
    double sum = 0.0, mn = 0.0, mx = 0.0;
    long mixed = 0, wisps = 0;
    MD pol(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z});
    Kokkos::parallel_reduce(
        "vof::wy::diag_sum", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& acc) { acc += c(L3(x, y, z, e)); }, sum);
    Kokkos::parallel_reduce(
        "vof::wy::diag_min", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
          acc = Kokkos::fmin(acc, c(L3(x, y, z, e)));
        },
        Kokkos::Min<double>(mn));
    Kokkos::parallel_reduce(
        "vof::wy::diag_max", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
          acc = Kokkos::fmax(acc, c(L3(x, y, z, e)));
        },
        Kokkos::Max<double>(mx));
    Kokkos::parallel_reduce(
        "vof::wy::diag_counts", pol,
        KOKKOS_LAMBDA(int x, int y, int z, long& nm, long& nw) {
          const double v = c(L3(x, y, z, e));
          if (wyIsMixed(v)) {
            ++nm;
            if (v < 1e-8 || v > 1.0 - 1e-8)
              ++nw;
          }
        },
        mixed, wisps);
    Kokkos::fence();
    Diagnostics dg;
    dg.sumC = sum;
    dg.minC = mn;
    dg.maxC = mx;
    dg.mixed = mixed;
    dg.wisps = wisps;
    if (hasGeom_)
      cutDiagnostics(dg);
    return dg;
  }
 
  double maxCourant(double dth) const {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField u = uf_, v = vf_, w = wf_;
    double m = 0.0;
    Kokkos::parallel_reduce(
        "vof::wy::cfl", 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);
          acc = Kokkos::fmax(acc, Kokkos::fabs(u(i)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(v(i)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(w(i)));
        },
        Kokkos::Max<double>(m));
    Kokkos::fence();
    return m * dth;
  }
 
  double maxCourantInterface(double dth) const {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long sx = 1, sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField c = c_, u = uf_, v = vf_, w = wf_;
    const double weps = wispEps;
    double m = 0.0;
    Kokkos::parallel_reduce(
        "vof::wy::cfl_interface", 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);
          const bool band = wyIsMixed(ci, weps) || wyColourJump(c(i - sx), ci, weps) ||
                            wyColourJump(c(i + sx), ci, weps) ||
                            wyColourJump(c(i - sy), ci, weps) ||
                            wyColourJump(c(i + sy), ci, weps) ||
                            wyColourJump(c(i - sz), ci, weps) || wyColourJump(c(i + sz), ci, weps);
          if (!band)
            return;
          acc = Kokkos::fmax(acc, Kokkos::fabs(u(i)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(u(i - sx)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(v(i)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(v(i - sy)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(w(i)));
          acc = Kokkos::fmax(acc, Kokkos::fabs(w(i - sz)));
        },
        Kokkos::Max<double>(m));
    Kokkos::fence();
    // Max's identity is -inf; an empty band (a uniform colour field) means no constraint, not an
    // unconstrained step. Clamp so the reported CFL is a number and the cap test is meaningful.
    return Kokkos::fmax(m, 0.0) * dth;
  }
 
  double maxCourantInterfaceAuto(double dth) const {
    return hasGeom_ ? maxCourantInterfaceCut(dth) : maxCourantInterface(dth);
  }
 
  double maxDiscreteDivergence(double dth) const {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField u = uf_, v = vf_, w = wf_;
    double m = 0.0;
    Kokkos::parallel_reduce(
        "vof::wy::divmax", 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 d = (u(i) - u(i - 1)) + (v(i) - v(i - sy)) + (w(i) - w(i - sz));
          acc = Kokkos::fmax(acc, Kokkos::fabs(d) * dth);
        },
        Kokkos::Max<double>(m));
    Kokkos::fence();
    return m;
  }
 
 public:
  // ---- implementation detail -----------------------------------------------------------------
  // These stay PUBLIC only because nvcc rejects an extended __host__ __device__ lambda whose
  // enclosing member function has private or protected access. Treat them as private.
 
  void requireExchange() const {
    if (!exchange)
      throw std::runtime_error("peclet::flow::vof::WyAdvector: the `exchange` hook is not set");
  }
 
  void freezeDilationFlag() {
    const double _t0 = tick_();
    SField c = c_;
    UCField cc = cc_;
    Kokkos::parallel_for(
        "vof::wy::freeze", Kokkos::RangePolicy<SExec>(SExec(), 0, len_),
        // thesis eq. A.28: g = 1 if C^n > 1/2, else 0. Strict `>` (the tie is measure zero and the
        // proof only needs the flag to be a constant of the step).
        KOKKOS_LAMBDA(long i) { cc(i) = c(i) > 0.5 ? 1u : 0u; });
    addT_(tm.freeze, _t0);
  }
 
  void reconstruct() {
    const double _t0 = tick_();
    reconstructImpl();
    addT_(tm.reconstruct, _t0);
  }
  void reconstructImpl() {
    const I3 e = e_;
    const int g = g_;
    const int rx = n_.x + 2, ry = n_.y + 2, rz = n_.z + 2;
    const long region = static_cast<long>(rx) * ry * rz;
    const long sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField c = c_, mx = mx_, my = my_, mz = mz_, al = alpha_;
    const double weps = wispEps;  // WO-R2 item 4: the reconstruction and the flux share it
 
    if (useWorklist) {
      LField list = list_;
      long cnt = 0;
      Kokkos::parallel_scan(
          "vof::wy::worklist", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
          KOKKOS_LAMBDA(const long r, long& upd, const bool final) {
            const int ix = static_cast<int>(r % rx);
            const int iy = static_cast<int>((r / rx) % ry);
            const int iz = static_cast<int>(r / (static_cast<long>(rx) * ry));
            const long i = L3(g - 1 + ix, g - 1 + iy, g - 1 + iz, e);
            if (wyIsMixed(c(i), weps)) {
              if (final)
                list(upd) = i;
              ++upd;
            }
          },
          cnt);
      Kokkos::fence();
      mixedCount_ = cnt;
      Kokkos::parallel_for(
          "vof::wy::plic", Kokkos::RangePolicy<SExec>(SExec(), 0, cnt),
          KOKKOS_LAMBDA(long t) { wyReconstructCell(c, list(t), sy, sz, mx, my, mz, al); });
    } else {
      Kokkos::parallel_for(
          "vof::wy::plic_dense", Kokkos::RangePolicy<SExec>(SExec(), 0, region),
          KOKKOS_LAMBDA(const long r) {
            const int ix = static_cast<int>(r % rx);
            const int iy = static_cast<int>((r / rx) % ry);
            const int iz = static_cast<int>(r / (static_cast<long>(rx) * ry));
            const long i = L3(g - 1 + ix, g - 1 + iy, g - 1 + iz, e);
            if (wyIsMixed(c(i), weps))
              wyReconstructCell(c, i, sy, sz, mx, my, mz, al);
          });
    }
  }
 
  void computeFluxes(int d, double dth) {
    const double _t0 = tick_();
    computeFluxesImpl(d, dth);
    addT_(tm.fluxes, _t0);
  }
  void computeFluxesImpl(int d, double dth) {
    if (hasGeom_) {  // rung V5a: openness-weighted flux (`cutcell.hpp` eq. 1)
      // COMPOSED with the V-BC boundary rule since WO-R2 item 2: `computeFluxesCut` carries the
      // out-of-domain mask and the per-face ledger itself, so `F = o_f * C_datum * a` at a domain
      // face of a geometry-carrying block, including a domain face a solid cuts (o_f < 1).
      computeFluxesCut(d, dth);
      return;
    }
    const I3 e = e_;
    const int g = g_;
    const long sd =
        d == 0 ? 1 : (d == 1 ? static_cast<long>(e_.x) : static_cast<long>(e_.x) * e_.y);
    // d-index [g-1, g+n_d): the `-` face of the first inner cell through the `+` face of the last.
    // Transverse indices stay on the inner region.
    int lo[3] = {g, g, g};
    const int hi[3] = {g + n_.x, g + n_.y, g + n_.z};
    lo[d] -= 1;
    SField c = c_, mx = mx_, my = my_, mz = mz_, al = alpha_, fl = flux_, u = faceVel(d);
    const double weps = wispEps;
    if (hasOutsideMask()) {  // WO-R: boundary donors are DATA, not a reconstructable interface
      UCField ob = outside_;
      Kokkos::parallel_for(
          "vof::wy::flux_bc",
          MDRange3<SExec>(SExec(), {lo[0], lo[1], lo[2]}, {hi[0], hi[1], hi[2]}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const long p = L3(x, y, z, e);
            fl(p) = wyFaceFluxBc(u(p) * dth, p, sd, d, c, mx, my, mz, al, ob, weps);
          });
      accumulateBcFaceVolume(d);
      return;
    }
    Kokkos::parallel_for(
        "vof::wy::flux", MDRange3<SExec>(SExec(), {lo[0], lo[1], lo[2]}, {hi[0], hi[1], hi[2]}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long p = L3(x, y, z, e);
          fl(p) = wyFaceFlux(u(p) * dth, p, sd, d, c, mx, my, mz, al, weps);
        });
  }
 
  void accumulateBcFaceVolume(int d) {
    const I3 e = e_, n = n_;
    const int g = g_;
    const int b = (d + 1) % 3, c2 = (d + 2) % 3;
    const long st[3] = {1, e_.x, static_cast<long>(e_.x) * e_.y};
    const long sd = st[d], sb = st[b], sc = st[c2];
    const int nb = (b == 0) ? n.x : (b == 1) ? n.y : n.z;
    const int nc = (c2 == 0) ? n.x : (c2 == 1) ? n.y : n.z;
    const int nd = (d == 0) ? n.x : (d == 1) ? n.y : n.z;
    SField fl = flux_;
    for (int s = 0; s < 2; ++s) {
      if (!bcOwn_[2 * d + s])
        continue;  // not a global domain face this block owns: nothing to account for
      const int fa = (s == 0) ? (g - 1) : (g + nd - 1);
      double acc = 0.0;
      Kokkos::parallel_reduce(
          "vof::wy::bcvol", MDRange2<SExec>(SExec(), {g, g}, {g + nb, g + nc}),
          KOKKOS_LAMBDA(int p0, int p1, double& r) {
            r += fl(static_cast<long>(p0) * sb + static_cast<long>(p1) * sc +
                    static_cast<long>(fa) * sd);
          },
          acc);
      Kokkos::fence();
      if (bcOwn_[2 * d + s])
        bcVol_[2 * d + s] += (s == 0) ? acc : -acc;
    }
  }
 
  void applySweep(int d, double dth) {
    const double _t0 = tick_();
    applySweepImpl(d, dth);
    addT_(tm.sweep, _t0);
    ++tm.sweeps;
  }
  void applySweepImpl(int d, double dth) {
    if (hasGeom_) {  // rung V5a: the update in FLUID-VOLUME units
      applySweepCut(d, dth);
      return;
    }
    const I3 e = e_, n = n_;
    const int g = g_;
    const long sd =
        d == 0 ? 1 : (d == 1 ? static_cast<long>(e_.x) : static_cast<long>(e_.x) * e_.y);
    SField c = c_, fl = flux_, u = faceVel(d);
    UCField cc = cc_;
    Kokkos::parallel_for(
        "vof::wy::update", 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);
          // The dilation term must scale the SAME uf by the SAME dt/h as the flux, or the exact
          // cancellation in full cells (file header) is lost to rounding.
          const double aP = u(i) * dth, aM = u(i - sd) * dth;
          const double dil = cc(i) ? (aP - aM) : 0.0;
          c(i) = c(i) + (fl(i - sd) - fl(i)) + dil;
        });
  }
 
  // ---- rung V5a (WO-Q) kernel bodies ---------------------------------------------------------
  // Siblings of the V1 kernels above, reached only through the `hasGeom_` branch at the top of
  // each. The V1 bodies are untouched, so "no geometry" is byte-identical by construction rather
  // than by a tolerance.
 
  void computeFluxesCut(int d, double dth) {
    const I3 e = e_;
    const int g = g_;
    const long sd =
        d == 0 ? 1 : (d == 1 ? static_cast<long>(e_.x) : static_cast<long>(e_.x) * e_.y);
    int lo[3] = {g, g, g};
    const int hi[3] = {g + n_.x, g + n_.y, g + n_.z};
    lo[d] -= 1;
    SField c = c_, mx = mx_, my = my_, mz = mz_, al = alpha_, fl = flux_, u = faceVel(d),
           o = of_[d], ep = eps_;
    UCField kd = kind_;
    const bool clamp = cutFluxClamp;
    const double weps = wispEps;
    // WO-R2 item 2 — V5a x V-BC. With BOTH rungs armed a domain-face flux in a geometry-carrying
    // block is `o_f * C_datum * a`: the out-of-domain donor is fluxed ALGEBRAICALLY (it is a
    // prescribed datum with no reconstructable interface — `wyFaceFluxBc`), and the openness
    // weight `o_f` still applies, because a solid may cut the domain face itself (o_f < 1 there,
    // 1 on the clear part). The per-face boundary ledger then accumulates the SAME `o_f`-weighted
    // number the update telescopes on, so the WO-R budget identity holds unchanged in cut cells.
    // The clamp keeps its contract: it applies only to a MIXED donor, and an out-of-domain donor
    // has taken the algebraic branch, whose flux is already exactly bounded.
    const bool bc = hasOutsideMask();
    UCField ob = bc ? outside_ : UCField();
    long nclamp = 0;
    Kokkos::parallel_reduce(
        "vof::wy::flux_cut", MDRange3<SExec>(SExec(), {lo[0], lo[1], lo[2]}, {hi[0], hi[1], hi[2]}),
        KOKKOS_LAMBDA(int x, int y, int z, long& acc) {
          const long p = L3(x, y, z, e);
          const double op = o(p);
          if (op == 0.0) {
            fl(p) = 0.0;
            return;
          }
          const double a = u(p) * dth;
          const double raw =
              vofCutFlux(op, bc ? wyFaceFluxBc(a, p, sd, d, c, mx, my, mz, al, ob, weps)
                                : wyFaceFlux(a, p, sd, d, c, mx, my, mz, al, weps));
          if (!clamp || a == 0.0) {
            fl(p) = raw;
            return;
          }
          const long q = a > 0.0 ? p : p + sd;
          const double cq = c(q);
          // algebraic branch: already exactly bounded, and clamping it would break the exact
          // full-cell cancellation. An out-of-domain donor is on that branch by construction.
          if (!wyIsMixed(cq, weps) || (bc && ob(q))) {
            fl(p) = raw;
            return;
          }
          const double aa = a > 0.0 ? a : -a;
          const double mag = raw > 0.0 ? raw : -raw;
          const double cl = vofCutFluxClamp(mag, op * aa, vofEpsEff(ep(q), kd(q) == kVofSolid), cq);
          if (cl != mag)
            ++acc;
          fl(p) = a > 0.0 ? cl : -cl;
        },
        nclamp);
    Kokkos::fence();
    clampedFaces_ += nclamp;
    if (bc)
      accumulateBcFaceVolume(d);
  }
 
  void applySweepCut(int d, double dth) {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long sd =
        d == 0 ? 1 : (d == 1 ? static_cast<long>(e_.x) : static_cast<long>(e_.x) * e_.y);
    SField c = c_, fl = flux_, u = faceVel(d), o = of_[d], ep = eps_;
    UCField cc = cc_, kd = kind_;
    Kokkos::parallel_for(
        "vof::wy::update_cut", 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);
          if (kd(i) == kVofSolid)
            return;
          const double epsEff = vofEpsEff(ep(i), false);
          // Same `o_f a_f` in the flux and in the dilation term — that shared factor is what makes
          // the two contributions telescope against the projection's openness-weighted divergence.
          // `o * (u * dth)` and NOT `(o * u) * dth`: the flux forms the Courant number FIRST
          // (`wyFaceFlux(u(p) * dth, ...)`), so only this association makes the two cancel to the
          // last bit in a full cell — the property that keeps the conservation floor on the
          // interface area rather than on the domain volume (see the file header).
          const double aP = o(i) * (u(i) * dth), aM = o(i - sd) * (u(i - sd) * dth);
          const double dil = cc(i) ? (aP - aM) : 0.0;
          c(i) = c(i) + ((fl(i - sd) - fl(i)) + dil) / epsEff;
        });
  }
 
  void clipCutCells() {
    const double _t0 = tick_();
    clipCutCellsImpl();
    addT_(tm.clip, _t0);
  }
  void clipCutCellsImpl() {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField c = c_, ep = eps_;
    UCField kd = kind_;
    double acc = 0.0, sgn = 0.0;
    Kokkos::parallel_reduce(
        "vof::wy::clip_cut", 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, double& b) {
          const long i = L3(x, y, z, e);
          if (kd(i) == kVofSolid)
            return;
          const double epsEff = vofEpsEff(ep(i), false);
          if (!(epsEff < 1.0))
            return;  // whole fluid cell: V1 boundedness applies verbatim, do not touch it
          const double cv = c(i);
          if (cv >= 0.0 && cv <= 1.0)
            return;
          double moved = 0.0;
          c(i) = vofClipCut(cv, epsEff, moved);
          a += Kokkos::fabs(moved);
          b += moved;
        },
        acc, sgn);
    Kokkos::fence();
    clippedVolume_ += acc;
    clippedSigned_ += sgn;
  }
 
  double maxCourantInterfaceCut(double dth) const {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long sx = 1, sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField c = c_, u = uf_, v = vf_, w = wf_, ox = of_[0], oy = of_[1], oz = of_[2], ep = eps_;
    UCField kd = kind_;
    const double weps = wispEps;
    double m = 0.0;
    Kokkos::parallel_reduce(
        "vof::wy::cfl_interface_cut",
        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);
          if (kd(i) == kVofSolid)
            return;
          const double ci = c(i);
          const bool band = wyIsMixed(ci, weps) || wyColourJump(c(i - sx), ci, weps) ||
                            wyColourJump(c(i + sx), ci, weps) ||
                            wyColourJump(c(i - sy), ci, weps) ||
                            wyColourJump(c(i + sy), ci, weps) ||
                            wyColourJump(c(i - sz), ci, weps) || wyColourJump(c(i + sz), ci, weps);
          if (!band)
            return;
          const double ei = ep(i);
          acc = Kokkos::fmax(acc, vofCutCourant(ox(i), u(i) * dth, ei));
          acc = Kokkos::fmax(acc, vofCutCourant(ox(i - sx), u(i - sx) * dth, ei));
          acc = Kokkos::fmax(acc, vofCutCourant(oy(i), v(i) * dth, ei));
          acc = Kokkos::fmax(acc, vofCutCourant(oy(i - sy), v(i - sy) * dth, ei));
          acc = Kokkos::fmax(acc, vofCutCourant(oz(i), w(i) * dth, ei));
          acc = Kokkos::fmax(acc, vofCutCourant(oz(i - sz), w(i - sz) * dth, ei));
        },
        Kokkos::Max<double>(m));
    Kokkos::fence();
    return Kokkos::fmax(m, 0.0);
  }
 
  void solidBandFillPass(int k) {
    const I3 e = e_, n = n_;
    const int g = g_, maxDepth = 3 - k;
    const long sx = 1, sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField c = c_;
    UCField fs = fill_, mk = mark_;
    Kokkos::parallel_for(
        "vof::wy::band_fill", MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          if (vofGhostDepth(x, y, z, g, n.x, n.y, n.z) > maxDepth)
            return;
          const long i = L3(x, y, z, e);
          if (fs(i) != kVofFillNone)
            return;
          const long nb[6] = {i - sx, i + sx, i - sy, i + sy, i - sz, i + sz};
          double acc = 0.0;
          int cnt = 0;
          for (int q = 0; q < 6; ++q)
            if (vofFillReadable(fs(nb[q]), k)) {
              acc += c(nb[q]);
              ++cnt;
            }
          if (cnt == 0)
            return;  // no donor yet — a later pass may reach it, or it stays as it was
          c(i) = acc / cnt;
          mk(i) = 1u;
        });
    Kokkos::fence();
    Kokkos::parallel_for(
        "vof::wy::band_fill_commit", Kokkos::RangePolicy<SExec>(SExec(), 0, len_),
        KOKKOS_LAMBDA(long i) {
          if (mk(i)) {
            fs(i) = static_cast<unsigned char>(kVofFillPass0 + k);
            mk(i) = 0u;
          }
        });
    Kokkos::fence();
  }
 
  void solidBandFillPassWetting() {
    const I3 e = e_, n = n_;
    const int g = g_, maxDepth = 2;  // pass k = 1
    const long sx = 1, sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField c = c_, sdf = sdfB_, th = thetaB_, ap = appB_;
    SField mx = mfl_[0], my = mfl_[1], mz = mfl_[2];
    UCField fs = fill_, mk = mark_, kk = kind_, wb = wetB_;
    const int pivot = wettingPivot;
    const double tEps = wettingTangentEps, pureEps = wettingPureEps;
    const VofMetric gme = metric;  // `g` is the ghost width in this scope
    Kokkos::parallel_for(
        "vof::wy::band_fill_theta", MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          if (vofGhostDepth(x, y, z, g, n.x, n.y, n.z) > maxDepth)
            return;
          const long i = L3(x, y, z, e);
          if (fs(i) != kVofFillNone)
            return;  // fluid, or already carrying data
          // The wall normal, solid -> fluid, from the exchanged SDF. Depth <= 2 guarantees the
          // central difference stays inside the g = 3 block. The difference is taken on the INDEX
          // lattice, so `gw` is an index-space gradient; `vofWettingPlane` pulls it back to the
          // physical direction itself (V4.1).
          double nw[3] = {0.5 * (sdf(i + sx) - sdf(i - sx)), 0.5 * (sdf(i + sy) - sdf(i - sy)),
                          0.5 * (sdf(i + sz) - sdf(i - sz))};
          // Phase 3 (V4.2): the WALK is an index walk along the PHYSICAL wall normal, i.e. along
          // `H^-1 n_w ~ H^-2 grad_xi(sdf)`. Written as one division per axis followed by ONE L2
          // normalisation, so at `metric.h = {1,1,1}` it is `gw/|gw|` — today's vector, bit for
          // bit (a second normalisation of an already-unit vector would NOT be).
          double wk[3] = {nw[0] / (gme.h[0] * gme.h[0]), nw[1] / (gme.h[1] * gme.h[1]),
                          nw[2] / (gme.h[2] * gme.h[2])};
          const double nn = Kokkos::sqrt(wk[0] * wk[0] + wk[1] * wk[1] + wk[2] * wk[2]);
          if (nn > 1e-12) {
            wk[0] /= nn;
            wk[1] /= nn;
            wk[2] /= nn;
            for (int step = 1; step <= 4; ++step) {
              const int fx = x + static_cast<int>(Kokkos::round(step * wk[0]));
              const int fy = y + static_cast<int>(Kokkos::round(step * wk[1]));
              const int fz = z + static_cast<int>(Kokkos::round(step * wk[2]));
              if (fx < 0 || fy < 0 || fz < 0 || fx >= e.x || fy >= e.y || fz >= e.z)
                break;
              const long fi = L3(fx, fy, fz, e);
              if (kk(fi) != kVofFluid)
                continue;  // still inside the solid: keep walking outward
              const double cf = c(fi);
              if (cf <= pureEps || cf >= 1.0 - pureEps) {
                // The anchor column carries no interface, but the CONTINUED interface may still
                // reach this band cell from the column next door (the wetting case: the contact
                // line crosses the first fluid row one column in). Average the theta-planes of the
                // anchor's MIXED fluid neighbours; with none, the pure continuation stands.
                double acc2 = 0.0;
                int cnt2 = 0;
                const bool inner =
                    fx >= 1 && fy >= 1 && fz >= 1 && fx + 1 < e.x && fy + 1 < e.y && fz + 1 < e.z;
                if (inner)
                  for (int kz = -1; kz <= 1; ++kz)
                    for (int ky = -1; ky <= 1; ++ky)
                      for (int kx = -1; kx <= 1; ++kx) {
                        const long gi = L3(fx + kx, fy + ky, fz + kz, e);
                        if (kk(gi) != kVofFluid)
                          continue;
                        const double cg = c(gi);
                        if (cg <= pureEps || cg >= 1.0 - pureEps)
                          continue;
                        const double mg[3] = {mx(gi), my(gi), mz(gi)};
                        if (Kokkos::fabs(mg[0]) + Kokkos::fabs(mg[1]) + Kokkos::fabs(mg[2]) <= 0.0)
                          continue;
                        double mt2[3], al2, ca2;
                        const double t2 = th(i);
                        vofWettingPlane(mg, cg, nw, Kokkos::cos(t2), Kokkos::sin(t2), sdf(gi),
                                        pivot, tEps, mt2, al2, ca2, gme);
                        const int ds2[3] = {x - (fx + kx), y - (fy + ky), z - (fz + kz)};
                        acc2 += vofWettingFraction(mt2, al2, ds2);
                        ++cnt2;
                      }
                c(i) = cnt2 ? acc2 / cnt2 : cf;
                mk(i) = 1u;
                wb(i) = static_cast<unsigned char>(cnt2 ? kVofWetNeighbour : kVofWetPure);
                return;
              }
              const double mf[3] = {mx(fi), my(fi), mz(fi)};
              if (Kokkos::fabs(mf[0]) + Kokkos::fabs(mf[1]) + Kokkos::fabs(mf[2]) <= 0.0)
                break;  // no usable fluid-only normal -> the neutral fallback below
              double mth[3], alphaTh, cosApp;
              const double t0 = th(i);
              const int br = vofWettingPlane(mf, cf, nw, Kokkos::cos(t0), Kokkos::sin(t0), sdf(fi),
                                             pivot, tEps, mth, alphaTh, cosApp, gme);
              const int ds[3] = {x - fx, y - fy, z - fz};
              c(i) = vofWettingFraction(mth, alphaTh, ds);
              ap(i) = Kokkos::acos(cosApp < -1.0 ? -1.0 : (cosApp > 1.0 ? 1.0 : cosApp));
              mk(i) = 1u;
              wb(i) = static_cast<unsigned char>(br);
              return;
            }
          }
          // Fallback: WO-Q's neutral pass-1 rule, verbatim.
          const long nb[6] = {i - sx, i + sx, i - sy, i + sy, i - sz, i + sz};
          double acc = 0.0;
          int cnt = 0;
          for (int q = 0; q < 6; ++q)
            if (vofFillReadable(fs(nb[q]), 1)) {
              acc += c(nb[q]);
              ++cnt;
            }
          if (cnt == 0)
            return;
          c(i) = acc / cnt;
          mk(i) = 1u;
          wb(i) = static_cast<unsigned char>(kVofWetNeutral);
        });
    Kokkos::fence();
    Kokkos::parallel_for(
        "vof::wy::band_fill_theta_commit", Kokkos::RangePolicy<SExec>(SExec(), 0, len_),
        KOKKOS_LAMBDA(long i) {
          if (mk(i)) {
            fs(i) = static_cast<unsigned char>(kVofFillPass0 + 1);
            mk(i) = 0u;
          }
        });
    Kokkos::fence();
  }
 
  void cutDiagnostics(Diagnostics& dg) const {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField c = c_, ep = eps_;
    UCField kd = kind_;
    using MD = MDRange3<SExec>;
    MD pol(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z});
    double vol = 0.0, raw = 0.0, solidC = 0.0;
    long ncut = 0, nsolid = 0;
    Kokkos::parallel_reduce(
        "vof::wy::cut_sums", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& v, double& r, double& sc) {
          const long i = L3(x, y, z, e);
          if (kd(i) == kVofSolid) {
            sc += c(i);
            return;
          }
          v += vofEpsEff(ep(i), false) * c(i);
          r += ep(i) * c(i);
        },
        vol, raw, solidC);
    Kokkos::parallel_reduce(
        "vof::wy::cut_counts", pol,
        KOKKOS_LAMBDA(int x, int y, int z, long& nc, long& ns) {
          const long i = L3(x, y, z, e);
          if (kd(i) == kVofSolid) {
            ++ns;
            return;
          }
          const double ee = vofEpsEff(ep(i), false);
          if (ee < 1.0)
            ++nc;
        },
        ncut, nsolid);
    double mn = 1e300, mx = -1e300;
    Kokkos::parallel_reduce(
        "vof::wy::cut_minfluid", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& a) {
          const long i = L3(x, y, z, e);
          if (kd(i) == kVofSolid || ep(i) < 1.0)
            return;
          a = Kokkos::fmin(a, c(i));
        },
        Kokkos::Min<double>(mn));
    Kokkos::parallel_reduce(
        "vof::wy::cut_maxfluid", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& a) {
          const long i = L3(x, y, z, e);
          if (kd(i) == kVofSolid || ep(i) < 1.0)
            return;
          a = Kokkos::fmax(a, c(i));
        },
        Kokkos::Max<double>(mx));
    Kokkos::fence();
    dg.volume = vol;
    dg.rawVolume = raw;
    dg.solidFillSum = solidC;
    dg.cutCells = ncut;
    dg.solidCells = nsolid;
    dg.minCFluid = mn > 1e299 ? 0.0 : mn;
    dg.maxCFluid = mx < -1e299 ? 0.0 : mx;
    dg.clippedVolume = clippedVolume_;
    dg.clippedSigned = clippedSigned_;
    dg.clampedFaces = clampedFaces_;
  }
 
 private:
  I3 n_{0, 0, 0}, e_{0, 0, 0};
  int g_ = 3;
  double h_ = 1.0;
  long len_ = 0, listCap_ = 0;
  SField c_, mx_, my_, mz_, alpha_, flux_, uf_, vf_, wf_;
  UCField cc_, outside_;  // outside_: the WO-R out-of-domain mask (empty = the V1 path)
  double bcVol_[6] = {0, 0, 0, 0, 0, 0};
  unsigned char bcOwn_[6] = {0, 0, 0, 0, 0, 0};
  LField list_;
  long mixedCount_ = 0, steps_ = 0;
  double lastCfl_ = 0.0;
  // rung V5a (WO-Q) cut-cell geometry; `hasGeom_` false => every V1 kernel runs unchanged.
  bool hasGeom_ = false;
  SField of_[3], eps_, kindD_;
  UCField kind_, fill_, mark_;
  // rung V5b (WO-S) wetting; `hasWet_` false => pass 1 is WO-Q's neutral rule, byte-identically.
  bool hasWet_ = false;
  SField sdfB_, thetaB_, mfl_[3], appB_;
  UCField wetB_;
  double clippedVolume_ = 0.0, clippedSigned_ = 0.0;
  long clampedFaces_ = 0;
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_ADVECT_WY_HPP