block_batch.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_BLOCK_BATCH_HPP
#define PECLET_FLOW_VOF_BLOCK_BATCH_HPP
 
#include <Kokkos_Core.hpp>
 
#include "mac_stencils.hpp"  // peclet::flow::SExec, SField, I3, L3
#include "vof/advect_wy.hpp"
#include "vof/curvature_field.hpp"  // VofRawField
 
namespace peclet::flow::vof {
 
struct VofRawMask {
  const unsigned char* p;
  KOKKOS_INLINE_FUNCTION unsigned char operator()(long i) const { return p[i]; }
};
 
struct VofBlockJob {
  double *c, *mx, *my, *mz, *al, *fl;
  double* u[3];                  
  unsigned char* cc;             
  const unsigned char* outside;  
  I3 e, n, o;                    
  int g;                         
  double dth;                    
  double weps;                   
  bool interfaceCfl;             
  bool span[3];                  
  unsigned char bcOwn[6];        
};
 
inline constexpr int kVofBlockBatch = 16;
 
struct VofBlockTable {
  VofBlockJob job[kVofBlockBatch];
  long off[kVofBlockBatch + 1];
  int nj;
  int base;
};
 
KOKKOS_INLINE_FUNCTION int vofJobOf(const long* off, long t) {
  int b = 0;
  while (t >= off[b + 1])
    ++b;
  return b;
}
 
template <class F>
inline long vofSetOffsets(VofBlockTable& T, F size) {
  T.off[0] = 0;
  for (int k = 0; k < T.nj; ++k)
    T.off[k + 1] = T.off[k] + size(T.job[k]);
  for (int k = T.nj; k < kVofBlockBatch; ++k)
    T.off[k + 1] = T.off[T.nj];
  return T.off[T.nj];
}
 
inline long vofJobLen(const VofBlockJob& J) {
  return static_cast<long>(J.e.x) * J.e.y * J.e.z;
}
inline long vofJobInner(const VofBlockJob& J) {
  return static_cast<long>(J.n.x) * J.n.y * J.n.z;
}
inline long vofJobGrown1(const VofBlockJob& J) {
  return static_cast<long>(J.n.x + 2) * (J.n.y + 2) * (J.n.z + 2);
}
 
// ---- H-4: the host launch forms (`doc/vof_step_performance_design.md` §14.3) --------------------
//
// On a host backend GCC compiles a static `RangePolicy` to `schedule(static)`: one contiguous block
// of the flat index per thread. The device forms above therefore serialise on a host -- a list
// kernel launched over the upper bound has each block's active entries at the head of its range,
// on one thread -- and every region kernel pays 64-bit div/mod per cell. The launchers below are
// the host branches; each batched launcher selects them at compile time on its execution space's
// memory space (as `ccdetail::ccRows3` does), so a device build never instantiates them and its
// launches are untouched. Every host form visits the same cells with the same per-cell body, and
// every body writes only its own cell(s) or adds integers atomically, so the state is bitwise the
// device form's on the same host.
 
template <class Exec>
inline constexpr bool kVofHostExec = std::is_same_v<typename Exec::memory_space, Kokkos::HostSpace>;
 
struct VofCountTable {
  long o[kVofBlockBatch + 1];
};
 
template <class Exec, class Cnt, class F>
inline void vofHostListFor(const char* name, int nj, Cnt cnt, F body) {
  Exec().fence();  // the counts were written by an earlier launch
  VofCountTable C;
  C.o[0] = 0;
  for (int k = 0; k < nj; ++k) {
    const long n = cnt(k);
    C.o[k + 1] = C.o[k] + (n > 0 ? n : 0);
  }
  for (int k = nj; k < kVofBlockBatch; ++k)
    C.o[k + 1] = C.o[nj];
  Kokkos::RangePolicy<Exec, Kokkos::Schedule<Kokkos::Dynamic>> pol(Exec(), 0, C.o[nj]);
  pol.set_chunk_size(16);
  Kokkos::parallel_for(name, pol, [=](const long t) {
    const int k = vofJobOf(C.o, t);
    body(k, t - C.o[k]);
  });
}
 
template <class Tab>
inline auto vofListCount(const Tab& T, LField start, LField end) {
  return [=](int k) { return end(T.base + k) - start(T.base + k); };
}
 
// The host x loops of the region kernels carry `omp simd` where `ccFor3`'s contract holds (an
// iteration writes only its own cell and reads no cell another iteration of the row writes); with
// the host's -ffp-contract=off that is bit-identical to the scalar loop. Same definition as
// `mac_cutcell.hpp`'s.
#ifndef PECLET_FLOW_OMP_SIMD
#if defined(_OPENMP)
#define PECLET_FLOW_OMP_SIMD _Pragma("omp simd")
#else
#define PECLET_FLOW_OMP_SIMD
#endif
#endif
 
struct VofRowTable {
  long off[kVofBlockBatch + 1];
};
 
template <class Exec, class Rows, class F>
inline void vofHostRowFor(const char* name, int nj, Rows rows, F row) {
  VofRowTable R;
  R.off[0] = 0;
  for (int k = 0; k < nj; ++k)
    R.off[k + 1] = R.off[k] + rows(k);
  for (int k = nj; k < kVofBlockBatch; ++k)
    R.off[k + 1] = R.off[nj];
  Kokkos::parallel_for(name, Kokkos::RangePolicy<Exec>(Exec(), 0, R.off[nj]), [=](const long t) {
    const int k = vofJobOf(R.off, t);
    row(k, t - R.off[k]);
  });
}
 
template <class Exec, class Rows, class RowOf, class Hit>
inline void vofHostRowScan(const char* name, int nj, int base, Rows rows, RowOf rowOf, Hit hit,
                           LField list, LField start, LField end) {
  VofRowTable R;
  R.off[0] = 0;
  for (int k = 0; k < nj; ++k)
    R.off[k + 1] = R.off[k] + rows(k);
  for (int k = nj; k < kVofBlockBatch; ++k)
    R.off[k + 1] = R.off[nj];
  Kokkos::parallel_scan(name, Kokkos::RangePolicy<Exec>(Exec(), 0, R.off[nj]),
                        [=](const long t, long& upd, const bool final) {
                          const int k = vofJobOf(R.off, t);
                          const long r = t - R.off[k];
                          int x0, x1;
                          const long i0 = rowOf(k, r, x0, x1);
                          long q = upd;
                          if (final) {
                            if (r == 0)
                              start(base + k) = q;
                            for (int x = x0; x < x1; ++x)
                              if (hit(k, i0 + x))
                                list(q++) = i0 + x;
                            if (t + 1 == R.off[k + 1])
                              end(base + k) = q;
                          } else {
                            for (int x = x0; x < x1; ++x)
                              if (hit(k, i0 + x))
                                ++q;
                          }
                          upd = q;
                        });
}
 
// ---- stage 1: the Courant numbers ---------------------------------------------------------------
 
inline void vofBatchCfl(const VofBlockTable& T, SField out) {
  using Team = Kokkos::TeamPolicy<SExec>;
  Kokkos::parallel_for(
      "vof::block::batch_cfl", Team(SExec(), T.nj, Kokkos::AUTO),
      KOKKOS_LAMBDA(const typename Team::member_type& tm) {
        const int k = tm.league_rank();
        const VofBlockJob& J = T.job[k];
        const I3 e = J.e, n = J.n;
        const int g = J.g;
        const long sx = 1, sy = e.x, sz = static_cast<long>(e.x) * e.y;
        const long region = static_cast<long>(n.x) * n.y * n.z;
        const VofRawField c{J.c}, u{J.u[0]}, v{J.u[1]}, w{J.u[2]};
        const double weps = J.weps;
        const bool band0 = J.interfaceCfl;
        double m = 0.0;
        Kokkos::parallel_reduce(
            Kokkos::TeamThreadRange(tm, region),
            [&](const long r, double& acc) {
              const int x = static_cast<int>(r % n.x) + g;
              const int y = static_cast<int>((r / n.x) % n.y) + g;
              const int z = static_cast<int>(r / (static_cast<long>(n.x) * n.y)) + g;
              const long i = L3(x, y, z, e);
              if (!band0) {  // WyAdvector::maxCourant
                acc = Kokkos::fmax(acc, Kokkos::fabs(u(i)));
                acc = Kokkos::fmax(acc, Kokkos::fabs(v(i)));
                acc = Kokkos::fmax(acc, Kokkos::fabs(w(i)));
                return;
              }
              // WyAdvector::maxCourantInterface
              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::single(Kokkos::PerTeam(tm), [&]() {
          // the two return expressions of maxCourant / maxCourantInterface
          out(T.base + k) = band0 ? Kokkos::fmax(m, 0.0) * J.dth : m * J.dth;
        });
      });
}
 
// ---- stage 2: the frozen dilation flag ----------------------------------------------------------
 
template <class Exec = SExec>
inline void vofBatchFreeze(const VofBlockTable& T) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_freeze", T.nj,
        [&](const int k) { return static_cast<long>(T.job[k].e.y) * T.job[k].e.z; },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const long i0 = r * J.e.x;  // extended row r = (r % e.y, r / e.y)
          const double* c = J.c;
          unsigned char* cc = J.cc;
          PECLET_FLOW_OMP_SIMD
          for (int x = 0; x < J.e.x; ++x)
            cc[i0 + x] = c[i0 + x] > 0.5 ? 1u : 0u;
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_freeze", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long i = t - T.off[k];
        J.cc[i] = J.c[i] > 0.5 ? 1u : 0u;
      });
}
 
// ---- stage 3a/b: the worklist and the PLIC reconstruction --------------------------------------
 
template <class Exec = SExec>
inline void vofBatchWorklist(const VofBlockTable& T, LField list, LField start, LField end) {
  if constexpr (kVofHostExec<Exec>) {  // H-4(d): the grown-by-one rows
    vofHostRowScan<Exec>(
        "vof::block::batch_worklist", T.nj, T.base,
        [&](const int k) { return static_cast<long>(T.job[k].n.y + 2) * (T.job[k].n.z + 2); },
        [=](const int k, const long r, int& x0, int& x1) {
          const VofBlockJob& J = T.job[k];
          const int ry = J.n.y + 2, g = J.g;
          x0 = g - 1;
          x1 = g + J.n.x + 1;
          return L3(0, g - 1 + static_cast<int>(r % ry), g - 1 + static_cast<int>(r / ry), J.e);
        },
        [=](const int k, const long i) { return wyIsMixed(T.job[k].c[i], T.job[k].weps); }, list,
        start, end);
    return;
  }
  Kokkos::parallel_scan(
      "vof::block::batch_worklist", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t, long& upd, const bool final) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long r = t - T.off[k];
        const int rx = J.n.x + 2, ry = J.n.y + 2;
        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 int g = J.g;
        const long i = L3(g - 1 + ix, g - 1 + iy, g - 1 + iz, J.e);
        const bool mixed = wyIsMixed(J.c[i], J.weps);
        if (final) {
          if (r == 0)
            start(T.base + k) = upd;
          if (mixed)
            list(upd) = i;
          if (t + 1 == T.off[k + 1])
            end(T.base + k) = upd + (mixed ? 1 : 0);
        }
        if (mixed)
          ++upd;
      });
}
 
template <class Exec = SExec>
inline void vofBatchPlic(const VofBlockTable& T, LField list, LField start, LField end) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostListFor<Exec>("vof::block::batch_plic", T.nj, vofListCount(T, start, end),
                         [=](const int k, const long t) {
                           const long q = start(T.base + k) + t;
                           const VofBlockJob& J = T.job[k];
                           const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
                           wyReconstructCell(VofRawField{J.c}, list(q), sy, sz, VofRawField{J.mx},
                                             VofRawField{J.my}, VofRawField{J.mz},
                                             VofRawField{J.al});
                         });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_plic", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const long s0 = start(T.base + k);
        const long q = s0 + (t - T.off[k]);
        if (q >= end(T.base + k))
          return;
        const VofBlockJob& J = T.job[k];
        const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
        wyReconstructCell(VofRawField{J.c}, list(q), sy, sz, VofRawField{J.mx}, VofRawField{J.my},
                          VofRawField{J.mz}, VofRawField{J.al});
      });
}
 
// ---- stage 3c/d: fluxes and the sweep update ----------------------------------------------------
 
inline long vofJobFaces(const VofBlockJob& J, int d) {
  return static_cast<long>(J.n.x + (d == 0)) * (J.n.y + (d == 1)) * (J.n.z + (d == 2));
}
 
KOKKOS_INLINE_FUNCTION void vofBatchFluxCell(const VofBlockJob& J, int d, long p) {
  const I3 e = J.e;
  const long sd = d == 0 ? 1 : (d == 1 ? static_cast<long>(e.x) : static_cast<long>(e.x) * e.y);
  const VofRawField c{J.c}, mx{J.mx}, my{J.my}, mz{J.mz}, al{J.al}, u{J.u[d]};
  if (J.outside != nullptr) {
    J.fl[p] =
        wyFaceFluxBc(u(p) * J.dth, p, sd, d, c, mx, my, mz, al, VofRawMask{J.outside}, J.weps);
    return;
  }
  J.fl[p] = wyFaceFlux(u(p) * J.dth, p, sd, d, c, mx, my, mz, al, J.weps);
}
 
template <class Exec = SExec>
inline void vofBatchFlux(const VofBlockTable& T, int d) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_flux", T.nj,
        [&](const int k) {
          const VofBlockJob& J = T.job[k];
          return static_cast<long>(J.n.y + (d == 1)) * (J.n.z + (d == 2));
        },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const int ay = J.n.y + (d == 1), g = J.g;
          const int y = static_cast<int>(r % ay) + g - (d == 1);
          const int z = static_cast<int>(r / ay) + g - (d == 2);
          const long i0 = L3(0, y, z, J.e);
          const int x0 = g - (d == 0), x1 = g + J.n.x;
          PECLET_FLOW_OMP_SIMD
          for (int x = x0; x < x1; ++x)
            vofBatchFluxCell(J, d, i0 + x);
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_flux", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long r = t - T.off[k];
        const int ax = J.n.x + (d == 0), ay = J.n.y + (d == 1);
        const int g = J.g;
        const int x = static_cast<int>(r % ax) + g - (d == 0);
        const int y = static_cast<int>((r / ax) % ay) + g - (d == 1);
        const int z = static_cast<int>(r / (static_cast<long>(ax) * ay)) + g - (d == 2);
        vofBatchFluxCell(J, d, L3(x, y, z, J.e));
      });
}
 
inline void vofBatchBcLedger(const VofBlockTable& T, int d, SField out, int stride, int slot0) {
  Kokkos::parallel_for(
      "vof::block::batch_bcvol", Kokkos::RangePolicy<SExec>(SExec(), 0, 2 * T.nj),
      KOKKOS_LAMBDA(const long t) {
        const int k = static_cast<int>(t / 2), s = static_cast<int>(t % 2);
        const VofBlockJob& J = T.job[k];
        const int f = 2 * d + s;
        if (J.outside == nullptr || !J.bcOwn[f])
          return;
        const I3 e = J.e, n = J.n;
        const int g = J.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 nn[3] = {n.x, n.y, n.z};
        const int fa = (s == 0) ? (g - 1) : (g + nn[d] - 1);
        double acc = 0.0;
        for (int p1 = g; p1 < g + nn[c2]; ++p1)
          for (int p0 = g; p0 < g + nn[b]; ++p0)
            acc += J.fl[static_cast<long>(p0) * sb + static_cast<long>(p1) * sc +
                        static_cast<long>(fa) * sd];
        out(static_cast<long>(stride) * (T.base + k) + slot0 + f) += (s == 0) ? acc : -acc;
      });
}
 
KOKKOS_INLINE_FUNCTION void vofBatchUpdateCell(const VofBlockJob& J, int d, long i) {
  const I3 e = J.e;
  const long sd = d == 0 ? 1 : (d == 1 ? static_cast<long>(e.x) : static_cast<long>(e.x) * e.y);
  const double dth = J.dth;
  const VofRawField c{J.c}, fl{J.fl}, u{J.u[d]};
  const unsigned char* cc = J.cc;
  // The dilation term must scale the SAME uf by the SAME dt/h as the flux, or the exact
  // cancellation in full cells (advect_wy.hpp 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;
}
 
template <class Exec = SExec>
inline void vofBatchUpdate(const VofBlockTable& T, int d) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_update", T.nj,
        [&](const int k) { return static_cast<long>(T.job[k].n.y) * T.job[k].n.z; },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const int g = J.g;
          const long i0 =
              L3(0, static_cast<int>(r % J.n.y) + g, static_cast<int>(r / J.n.y) + g, J.e);
          PECLET_FLOW_OMP_SIMD
          for (int x = g; x < g + J.n.x; ++x)
            vofBatchUpdateCell(J, d, i0 + x);
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_update", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long r = t - T.off[k];
        const I3 e = J.e, n = J.n;
        const int g = J.g;
        const int x = static_cast<int>(r % n.x) + g;
        const int y = static_cast<int>((r / n.x) % n.y) + g;
        const int z = static_cast<int>(r / (static_cast<long>(n.x) * n.y)) + g;
        vofBatchUpdateCell(J, d, L3(x, y, z, e));
      });
}
 
// ---- stage 3e: the block ghost policy (`VofBlockSet::fillBlockGhosts`), pass by pass ------------
 
template <class Exec = SExec>
inline void vofBatchGhostZero(const VofBlockTable& T) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_ghost_zero", T.nj,
        [&](const int k) {
          const VofBlockJob& J = T.job[k];
          return static_cast<long>(J.e.z - J.n.z) * J.e.y +
                 static_cast<long>(J.n.z) * (J.e.y - J.n.y) + static_cast<long>(J.n.z) * J.n.y;
        },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const I3 e = J.e, n = J.n;
          const int g = J.g;
          const long rZ = static_cast<long>(e.z - n.z) * e.y,
                     rY = static_cast<long>(n.z) * (e.y - n.y);
          double* c = J.c;
          if (r < rZ + rY) {  // a whole ghost row
            int y, z;
            if (r < rZ) {  // z-ghost plane row
              const int l = static_cast<int>(r / e.y);
              y = static_cast<int>(r - static_cast<long>(l) * e.y);
              z = l < g ? l : n.z + l;
            } else {  // y-ghost row of the inner z range
              const long t2 = r - rZ;
              const int zi = static_cast<int>(t2 / (e.y - n.y));
              const int l = static_cast<int>(t2 - static_cast<long>(zi) * (e.y - n.y));
              y = l < g ? l : n.y + l;
              z = g + zi;
            }
            const long i0 = L3(0, y, z, e);
            PECLET_FLOW_OMP_SIMD
            for (int x = 0; x < e.x; ++x)
              c[i0 + x] = 0.0;
            return;
          }
          const long t3 = r - rZ - rY;  // the x-ghost runs of inner row (yi, zi)
          const int zi = static_cast<int>(t3 / n.y),
                    yi = static_cast<int>(t3 - static_cast<long>(zi) * n.y);
          const long i0 = L3(0, g + yi, g + zi, e);
          for (int x = 0; x < g; ++x)
            c[i0 + x] = 0.0;
          for (int x = g + n.x; x < e.x; ++x)
            c[i0 + x] = 0.0;
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_ghost_zero", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long i = t - T.off[k];
        const I3 e = J.e, n = J.n;
        const int g = J.g;
        const int x = static_cast<int>(i % e.x);
        const int y = static_cast<int>((i / e.x) % e.y);
        const int z = static_cast<int>(i / (static_cast<long>(e.x) * e.y));
        if (x >= g && x < g + n.x && y >= g && y < g + n.y && z >= g && z < g + n.z)
          return;
        J.c[i] = 0.0;
      });
}
 
template <class Exec = SExec>
inline void vofBatchPeriodic(const VofBlockTable& T, int a) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_periodic", T.nj,
        [&](const int k) {
          const I3 e = T.job[k].e;
          const int g = T.job[k].g;
          return a == 0 ? static_cast<long>(e.y) * e.z
                        : (a == 1 ? 2L * g * e.z : static_cast<long>(e.y) * 2 * g);
        },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const I3 e = J.e;
          const int g = J.g;
          double* f = J.c;
          if (a == 0) {  // extended row r: the g x-ghosts at each end
            const int N = e.x - 2 * g;
            const long base = r * e.x;
            for (int gl = 0; gl < g; ++gl) {
              f[base + gl] = f[base + gl + N];
              f[base + g + N + gl] = f[base + g + gl];
            }
            return;
          }
          long dst, src;
          if (a == 1) {  // ghost row (l, z): y = l (low) or g + N + (l - g) (high)
            const int N = e.y - 2 * g;
            const int z = static_cast<int>(r / (2 * g)), l = static_cast<int>(r % (2 * g));
            const int y = l < g ? l : N + l, sy = l < g ? l + N : l;  // high: y - N = g + (l - g)
            dst = L3(0, y, z, e);
            src = L3(0, sy, z, e);
          } else {  // ghost row (y, l)
            const int N = e.z - 2 * g;
            const int l = static_cast<int>(r / e.y), y = static_cast<int>(r % e.y);
            const int z = l < g ? l : N + l, sz = l < g ? l + N : l;
            dst = L3(0, y, z, e);
            src = L3(0, y, sz, e);
          }
          PECLET_FLOW_OMP_SIMD
          for (int x = 0; x < e.x; ++x)
            f[dst + x] = f[src + x];
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_periodic", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long r = t - T.off[k];
        const I3 e = J.e;
        const int g = J.g;
        const long st[3] = {1, e.x, static_cast<long>(e.x) * e.y};
        const int dims[3] = {e.x, e.y, e.z};
        const int b = (a + 1) % 3, c = (a + 2) % 3;
        const long sa = st[a], sb = st[b], sc = st[c];
        const int N = dims[a] - 2 * g;
        const int p0 = static_cast<int>(r % dims[b]);
        const int p1 = static_cast<int>(r / dims[b]);
        const VofRawField f{J.c};
        const long base = static_cast<long>(p0) * sb + static_cast<long>(p1) * sc;
        for (int gl = 0; gl < g; ++gl) {
          f(base + static_cast<long>(gl) * sa) = f(base + static_cast<long>(gl + N) * sa);
          f(base + static_cast<long>(g + N + gl) * sa) = f(base + static_cast<long>(g + gl) * sa);
        }
      });
}
 
KOKKOS_INLINE_FUNCTION void vofBatchClampCell(const VofBlockJob& J, int x, int y, int z, int q0,
                                              int q1, int q2, bool p0, bool p1, bool p2) {
  const I3 e = J.e, o = J.o;
  const int g = J.g;
  const int gx = x - g + o.x, gy = y - g + o.y, gz = z - g + o.z;
  const int cx = p0 ? gx : (gx < 0 ? 0 : (gx >= q0 ? q0 - 1 : gx));
  const int cy = p1 ? gy : (gy < 0 ? 0 : (gy >= q1 ? q1 - 1 : gy));
  const int cz = p2 ? gz : (gz < 0 ? 0 : (gz >= q2 ? q2 - 1 : gz));
  if (cx != gx || cy != gy || cz != gz)
    J.c[L3(x, y, z, e)] = J.c[L3(cx - o.x + g, cy - o.y + g, cz - o.z + g, e)];
}
 
template <class Exec = SExec>
inline void vofBatchClamp(const VofBlockTable& T, I3 gs, bool px, bool py, bool pz) {
  const bool p0 = px, p1 = py, p2 = pz;
  const int q0 = gs.x, q1 = gs.y, q2 = gs.z;
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_clamp", T.nj,
        [&](const int k) { return static_cast<long>(T.job[k].e.y) * T.job[k].e.z; },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const I3 e = J.e, o = J.o;
          const int g = J.g;
          const int y = static_cast<int>(r % e.y), z = static_cast<int>(r / e.y);
          const int gy = y - g + o.y, gz = z - g + o.z, gx0 = o.x - g, gx1 = o.x - g + e.x;
          if ((p1 || (gy >= 0 && gy < q1)) && (p2 || (gz >= 0 && gz < q2)) &&
              (p0 || (gx0 >= 0 && gx1 <= q0)))
            return;  // inside on every non-periodic axis: nothing in this row writes
          for (int x = 0; x < e.x; ++x)
            vofBatchClampCell(J, x, y, z, q0, q1, q2, p0, p1, p2);
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_clamp", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long i = t - T.off[k];
        const I3 e = J.e;
        const int x = static_cast<int>(i % e.x);
        const int y = static_cast<int>((i / e.x) % e.y);
        const int z = static_cast<int>(i / (static_cast<long>(e.x) * e.y));
        vofBatchClampCell(J, x, y, z, q0, q1, q2, p0, p1, p2);
      });
}
 
// ---- stage 4: marker debris and sub-wispEps residue (`VofBlockSet::debrisPassBatchImpl`) --------
 
template <class SF>
KOKKOS_INLINE_FUNCTION bool vofDebrisCell(const SF& c, long i, long sy, long sz, double ieps,
                                          double cfull) {
  if (!vofIsInterface(c(i), ieps))
    return false;
  for (int oz = -2; oz <= 2; ++oz)
    for (int oy = -2; oy <= 2; ++oy)
      for (int ox = -2; ox <= 2; ++ox)
        if (c(i + ox + oy * sy + oz * sz) >= cfull)
          return false;  // attached: a cell at least half full of this marker within two
  return true;
}
 
template <class SF>
KOKKOS_INLINE_FUNCTION bool vofAttachedCell(const SF& c, long i, long sy, long sz, double ieps,
                                            double cfull) {
  if (!vofIsInterface(c(i), ieps))
    return false;
  bool attached = false;
  for (int oz = -2; oz <= 2 && !attached; ++oz)
    for (int oy = -2; oy <= 2 && !attached; ++oy)
      for (int ox = -2; ox <= 2 && !attached; ++ox)
        attached = c(i + ox + oy * sy + oz * sz) >= cfull;
  return attached;
}
 
KOKKOS_INLINE_FUNCTION bool vofResidueCell(double ci, double weps) {
  return ci != 0.0 && Kokkos::fabs(ci) <= weps;
}
 
template <class Exec = SExec>
inline void vofBatchMark(const VofBlockTable& T, int kind, double ieps, double cfull, double weps,
                         LField list, LField start, LField end) {
  if constexpr (kVofHostExec<Exec>) {  // H-4(d): the inner rows
    vofHostRowScan<Exec>(
        "vof::block::batch_debris_mark", T.nj, T.base,
        [&](const int k) { return static_cast<long>(T.job[k].n.y) * T.job[k].n.z; },
        [=](const int k, const long r, int& x0, int& x1) {
          const VofBlockJob& J = T.job[k];
          const int g = J.g;
          x0 = g;
          x1 = g + J.n.x;
          return L3(0, g + static_cast<int>(r % J.n.y), g + static_cast<int>(r / J.n.y), J.e);
        },
        [=](const int k, const long i) {
          const VofBlockJob& J = T.job[k];
          const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
          const VofRawField c{J.c};
          return kind == 0   ? vofDebrisCell(c, i, sy, sz, ieps, cfull)
                 : kind == 1 ? vofResidueCell(c(i), weps)
                             : vofAttachedCell(c, i, sy, sz, ieps, cfull);
        },
        list, start, end);
    return;
  }
  Kokkos::parallel_scan(
      "vof::block::batch_debris_mark", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t, long& upd, const bool final) {
        const int k = vofJobOf(T.off, t);
        const VofBlockJob& J = T.job[k];
        const long r = t - T.off[k];
        const int nx = J.n.x, ny = J.n.y, g = J.g;
        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, J.e);
        const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
        const VofRawField c{J.c};
        const bool hit = kind == 0   ? vofDebrisCell(c, i, sy, sz, ieps, cfull)
                         : kind == 1 ? vofResidueCell(c(i), weps)
                                     : vofAttachedCell(c, i, sy, sz, ieps, cfull);
        if (final) {
          if (r == 0)
            start(T.base + k) = upd;
          if (hit)
            list(upd) = i;
          if (t + 1 == T.off[k + 1])
            end(T.base + k) = upd + (hit ? 1 : 0);
        }
        if (hit)
          ++upd;
      });
}
 
enum VofDebrisSlot : int {
  kDbSumD = 0,  
  kDbSumR,      
  kDbW,         
  kDbVA,        
  kDbDV,        
  kDbActed,     
  kDbND,        
  kDbNR,        
  kDbPending,   
  kDbSlots
};
 
struct VofDebrisFlags {
  bool doDebris, doResidue, remove;
  double vaMin;
};
 
inline void vofBatchDebrisSums(const VofBlockTable& T, VofDebrisFlags F, LField listD, LField listR,
                               LField listA, LField sD0, LField eD0, LField sR0, LField eR0,
                               LField sA0, LField eA0, SField out, int stride, int slot0) {
  Kokkos::parallel_for(
      "vof::block::batch_debris_sums", Kokkos::RangePolicy<SExec>(SExec(), 0, T.nj),
      KOKKOS_LAMBDA(const long q) {
        const int j = T.base + static_cast<int>(q);
        const VofRawField c{T.job[q].c};
        const long nD = F.doDebris ? eD0(j) - sD0(j) : 0;
        const long nR = F.doResidue ? eR0(j) - sR0(j) : 0;
        const bool removeD = F.doDebris && F.remove && nD > 0;
        const bool pending = removeD || nR > 0;
        const long ob = static_cast<long>(stride) * j + slot0;
        double sD = 0.0, sR = 0.0, W = 0.0, VA = 0.0, dV = 0.0;
        const long d0 = F.doDebris ? sD0(j) : 0;
        for (long t = 0; t < nD; ++t)
          sD += c(listD(d0 + t));
        double acted = 0.0;
        if (pending) {
          const long r0 = F.doResidue ? sR0(j) : 0, a0 = sA0(j), nA = eA0(j) - sA0(j);
          for (long t = 0; t < nR; ++t)
            sR += c(listR(r0 + t));
          for (long t = 0; t < nA; ++t) {
            const double ca = c(listA(a0 + t));
            W += ca * (1.0 - ca);
            VA += ca;
          }
          acted = (VA < F.vaMin || !(W > 0.0)) ? 0.0 : 1.0;
          // the removed volume, in the fixed order D (if removing) then R -- the serial order
          const long nDr = removeD ? nD : 0;
          for (long t = 0; t < nDr; ++t)
            dV += c(listD(d0 + t));
          for (long t = 0; t < nR; ++t)
            dV += c(listR(r0 + t));
        }
        out(ob + kDbSumD) = sD;
        out(ob + kDbSumR) = sR;
        out(ob + kDbW) = W;
        out(ob + kDbVA) = VA;
        out(ob + kDbDV) = dV;
        out(ob + kDbActed) = acted;
        out(ob + kDbND) = static_cast<double>(nD);
        out(ob + kDbNR) = static_cast<double>(nR);
        out(ob + kDbPending) = pending ? 1.0 : 0.0;
      });
}
 
KOKKOS_INLINE_FUNCTION long vofDebrisActCount(const VofDebrisFlags& F, LField sA0, LField eA0,
                                              SField out, long ob, int j) {
  if (out(ob + kDbActed) == 0.0)
    return 0;
  const long nD = static_cast<long>(out(ob + kDbND)), nR = static_cast<long>(out(ob + kDbNR));
  const bool removeD = F.doDebris && F.remove && nD > 0;
  const long nDr = removeD ? nD : 0, nA = eA0(j) - sA0(j);
  return nDr + nR + nA;
}
 
KOKKOS_INLINE_FUNCTION void vofDebrisActEntry(const VofBlockTable& T, const VofDebrisFlags& F,
                                              int k, long t, LField listD, LField listR,
                                              LField listA, LField sD0, LField sR0, LField sA0,
                                              LField eA0, SField out, int stride, int slot0,
                                              SField ex, Kokkos::View<int*, SMem> capf) {
  const int j = T.base + k;
  const long ob = static_cast<long>(stride) * j + slot0;
  if (out(ob + kDbActed) == 0.0)
    return;
  const long nD = static_cast<long>(out(ob + kDbND)), nR = static_cast<long>(out(ob + kDbNR));
  const bool removeD = F.doDebris && F.remove && nD > 0;
  const long nDr = removeD ? nD : 0, nA = eA0(j) - sA0(j);
  const long nZ = nDr + nR, nTot = nZ + nA;
  if (t >= nTot)
    return;
  const VofRawField c{T.job[k].c};
  if (t < nDr) {
    c(listD(sD0(j) + t)) = 0.0;
    return;
  }
  if (t < nZ) {
    c(listR(sR0(j) + t - nDr)) = 0.0;
    return;
  }
  const double dV = out(ob + kDbDV), W = out(ob + kDbW);
  const long ta = t - nZ;
  const long i = listA(sA0(j) + 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;
    Kokkos::atomic_max(&capf(j), 1);
  }
  ex(T.off[k] + ta) = e;
  c(i) = cn;
}
 
template <class Exec = SExec>
inline void vofBatchDebrisWrite(const VofBlockTable& T, VofDebrisFlags F, LField listD,
                                LField listR, LField listA, LField sD0, LField sR0, LField sA0,
                                LField eA0, SField out, int stride, int slot0, SField ex,
                                Kokkos::View<int*, SMem> capf) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostListFor<Exec>(
        "vof::block::batch_debris_write", T.nj,
        [&](const int k) {
          const int j = T.base + k;
          return vofDebrisActCount(F, sA0, eA0, out, static_cast<long>(stride) * j + slot0, j);
        },
        [=](const int k, const long t) {
          vofDebrisActEntry(T, F, k, t, listD, listR, listA, sD0, sR0, sA0, eA0, out, stride, slot0,
                            ex, capf);
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_debris_write", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long tt) {
        const int k = vofJobOf(T.off, tt);
        vofDebrisActEntry(T, F, k, tt - T.off[k], listD, listR, listA, sD0, sR0, sA0, eA0, out,
                          stride, slot0, ex, capf);
      });
}
 
inline void vofBatchDebrisLost(const VofBlockTable& T, LField sA0, LField eA0, SField out,
                               int stride, int slot0, SField ex, Kokkos::View<int*, SMem> capf) {
  Kokkos::parallel_for(
      "vof::block::batch_debris_lost", Kokkos::RangePolicy<SExec>(SExec(), 0, T.nj),
      KOKKOS_LAMBDA(const long q) {
        const int j = T.base + static_cast<int>(q);
        const long ob = static_cast<long>(stride) * j + slot0;
        if (out(ob + kDbActed) == 0.0)
          return;
        double lost = 0.0;
        if (capf(j) > 0) {
          const long nA = eA0(j) - sA0(j), e0 = T.off[q];
          double a = 0.0;  // list order, as the serial loop accumulated it
          for (long t = 0; t < nA; ++t)
            a += ex(e0 + t);
          lost = a;
        }
        out(ob + kDbDV) = lost;
      });
}
 
// ---- stage 5: the bubble boxes (`VofBlockSet::bubbleBox`) ------------------------------------
 
template <class Exec = SExec>
inline void vofBatchBox(const VofBlockTable& T, double eps, SField out, int stride, int slot0) {
  if constexpr (kVofHostExec<Exec>) {
    int box[6 * kVofBlockBatch];
    for (int k = 0; k < T.nj; ++k)
      for (int a = 0; a < 6; ++a)
        box[6 * k + a] =
            a < 3 ? Kokkos::reduction_identity<int>::min() : Kokkos::reduction_identity<int>::max();
    int* bx = box;
    vofHostRowFor<Exec>(
        "vof::block::batch_bbox", T.nj,
        [&](const int k) { return static_cast<long>(T.job[k].n.y) * T.job[k].n.z; },
        [=](const int k, const long r) {
          const VofBlockJob& J = T.job[k];
          const I3 e = J.e, n = J.n;
          const int g = J.g;
          const int py = static_cast<int>(r % n.y), pz = static_cast<int>(r / n.y);
          const long i0 = L3(g, py + g, pz + g, e);
          const VofRawField c{J.c};
          int lo = -1, hi = -1;
          for (int px = 0; px < n.x; ++px)
            if (Kokkos::fabs(c(i0 + px)) > eps) {
              if (lo < 0)
                lo = px;
              hi = px;
            }
          if (lo < 0)
            return;  // no cell of the row is in the box
          int* s = bx + 6 * k;
          Kokkos::atomic_min(&s[0], lo);
          Kokkos::atomic_min(&s[1], py);
          Kokkos::atomic_min(&s[2], pz);
          Kokkos::atomic_max(&s[3], hi);
          Kokkos::atomic_max(&s[4], py);
          Kokkos::atomic_max(&s[5], pz);
        });
    Exec().fence();
    for (int k = 0; k < T.nj; ++k) {
      const long ob = static_cast<long>(stride) * (T.base + k) + slot0;
      for (int a = 0; a < 6; ++a)
        out(ob + a) = static_cast<double>(box[6 * k + a]);
    }
    return;
  }
  using Team = Kokkos::TeamPolicy<SExec>;
  Kokkos::parallel_for(
      "vof::block::batch_bbox", Team(SExec(), T.nj, Kokkos::AUTO),
      KOKKOS_LAMBDA(const typename Team::member_type& tm) {
        const int k = tm.league_rank();
        const VofBlockJob& J = T.job[k];
        const I3 e = J.e, n = J.n;
        const int g = J.g;
        const long region = static_cast<long>(n.x) * n.y * n.z;
        const VofRawField c{J.c};
        int v[6];
        for (int a = 0; a < 6; ++a) {
          int r = 0;
          auto body = [&](const long q, int& acc) {
            const int px = static_cast<int>(q % n.x);
            const int py = static_cast<int>((q / n.x) % n.y);
            const int pz = static_cast<int>(q / (static_cast<long>(n.x) * n.y));
            if (!(Kokkos::fabs(c(L3(px + g, py + g, pz + g, e))) > eps))
              return;
            const int p = (a % 3 == 0) ? px : (a % 3 == 1 ? py : pz);
            if (a < 3)
              acc = p < acc ? p : acc;
            else
              acc = p > acc ? p : acc;
          };
          if (a < 3)
            Kokkos::parallel_reduce(Kokkos::TeamThreadRange(tm, region), body, Kokkos::Min<int>(r));
          else
            Kokkos::parallel_reduce(Kokkos::TeamThreadRange(tm, region), body, Kokkos::Max<int>(r));
          v[a] = r;
        }
        Kokkos::single(Kokkos::PerTeam(tm), [&]() {
          const long ob = static_cast<long>(stride) * (T.base + k) + slot0;
          for (int a = 0; a < 6; ++a)
            out(ob + a) = static_cast<double>(v[a]);
        });
      });
}
 
// ---- stage 6: the curvature cascade and the CSF face force (`VofBlockSet::computeCsf`) ---------
 
struct VofCurvJob {
  double *c, *mx, *my, *mz, *al, *kap, *br;
  int* slot;     
  double* f[3];  
  I3 e, n;
  int g;
  double ieps, mtol, ptW, peps, dW, cmin, km;
  bool forceFb, oneDir, useFit;
  VofMetric gm;
};
 
struct VofCurvTable {
  VofCurvJob job[kVofBlockBatch];
  long off[kVofBlockBatch + 1];
  int nj;
  int base;
  PvMoments* cache;  
};
 
template <class F>
inline long vofSetOffsets(VofCurvTable& T, F size) {
  T.off[0] = 0;
  for (int k = 0; k < T.nj; ++k)
    T.off[k + 1] = T.off[k] + size(T.job[k]);
  for (int k = T.nj; k < kVofBlockBatch; ++k)
    T.off[k + 1] = T.off[T.nj];
  return T.off[T.nj];
}
inline long vofCurvInner(const VofCurvJob& J) {
  return static_cast<long>(J.n.x) * J.n.y * J.n.z;
}
inline long vofCurvGrown(const VofCurvJob& J) {
  return static_cast<long>(J.n.x + 2 * kPvHalf) * (J.n.y + 2 * kPvHalf) * (J.n.z + 2 * kPvHalf);
}
 
KOKKOS_INLINE_FUNCTION long vofCurvCell(const VofCurvJob& J, long r, bool grown) {
  const int gr = grown ? kPvHalf : 0;
  const int rx = J.n.x + 2 * gr, ry = J.n.y + 2 * gr;
  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));
  return L3(J.g - gr + ix, J.g - gr + iy, J.g - gr + iz, J.e);
}
 
inline long vofCurvRows(const VofCurvJob& J, bool grown) {
  const int gr = grown ? kPvHalf : 0;
  return static_cast<long>(J.n.y + 2 * gr) * (J.n.z + 2 * gr);
}
 
inline long vofCurvRow(const VofCurvJob& J, long r, bool grown, int& x0, int& x1) {
  const int gr = grown ? kPvHalf : 0;
  const int ry = J.n.y + 2 * gr;
  x0 = J.g - gr;
  x1 = J.g + J.n.x + gr;
  return L3(0, J.g - gr + static_cast<int>(r % ry), J.g - gr + static_cast<int>(r / ry), J.e);
}
 
template <class Exec = SExec>
inline void vofCurvCompact(const VofCurvTable& T, bool grown, LField list, LField start,
                           LField end) {
  if constexpr (kVofHostExec<Exec>) {  // H-4(d): the grown or inner rows
    vofHostRowScan<Exec>(
        "vof::block::batch_curv_compact", T.nj, T.base,
        [&](const int k) { return vofCurvRows(T.job[k], grown); },
        [=](const int k, const long r, int& x0, int& x1) {
          return vofCurvRow(T.job[k], r, grown, x0, x1);
        },
        [=](const int k, const long i) { return vofIsInterface(T.job[k].c[i], T.job[k].ieps); },
        list, start, end);
    return;
  }
  Kokkos::parallel_scan(
      "vof::block::batch_curv_compact", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t, long& upd, const bool final) {
        const int k = vofJobOf(T.off, t);
        const VofCurvJob& J = T.job[k];
        const long r = t - T.off[k];
        const long i = vofCurvCell(J, r, grown);
        const bool hit = vofIsInterface(J.c[i], J.ieps);
        if (final) {
          if (r == 0)
            start(T.base + k) = upd;
          if (hit)
            list(upd) = i;
          if (t + 1 == T.off[k + 1])
            end(T.base + k) = upd + (hit ? 1 : 0);
        }
        if (hit)
          ++upd;
      });
}
 
template <class Exec = SExec>
inline void vofCurvPlanesZero(const VofCurvTable& T) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_curv_planes_zero", T.nj,
        [&](const int k) { return vofCurvRows(T.job[k], true); },
        [=](const int k, const long r) {
          const VofCurvJob& J = T.job[k];
          int x0, x1;
          const long i0 = vofCurvRow(J, r, true, x0, x1);
          PECLET_FLOW_OMP_SIMD
          for (int x = x0; x < x1; ++x) {
            J.mx[i0 + x] = 0.0;
            J.my[i0 + x] = 0.0;
            J.mz[i0 + x] = 0.0;
            J.al[i0 + x] = 0.0;
            J.slot[i0 + x] = -1;
          }
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_curv_planes_zero", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofCurvJob& J = T.job[k];
        const long i = vofCurvCell(J, t - T.off[k], true);
        J.mx[i] = 0.0;
        J.my[i] = 0.0;
        J.mz[i] = 0.0;
        J.al[i] = 0.0;
        J.slot[i] = -1;
      });
}
 
KOKKOS_INLINE_FUNCTION void vofCurvPlaneEntry(const VofCurvTable& T, int k, long q, LField list) {
  const VofCurvJob& J = T.job[k];
  const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
  const long i = list(q);
  wyReconstructCell(VofRawField{J.c}, i, sy, sz, VofRawField{J.mx}, VofRawField{J.my},
                    VofRawField{J.mz}, VofRawField{J.al});
  curvMomentsBuild(J.mx[i], J.my[i], J.mz[i], J.al[i], J.gm, T.cache[q]);
  J.slot[i] = static_cast<int>(q);
}
 
template <class Exec = SExec>
inline void vofCurvPlanesList(const VofCurvTable& T, LField list, LField start, LField end) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostListFor<Exec>(
        "vof::block::batch_curv_planes", T.nj, vofListCount(T, start, end),
        [=](const int k, const long t) { vofCurvPlaneEntry(T, k, start(T.base + k) + t, list); });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_curv_planes", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const long q = start(T.base + k) + (t - T.off[k]);
        if (q >= end(T.base + k))
          return;
        vofCurvPlaneEntry(T, k, q, list);
      });
}
 
template <class Exec = SExec>
inline void vofCurvHfReset(const VofCurvTable& T) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_curv_hf_reset", T.nj,
        [&](const int k) { return vofCurvRows(T.job[k], false); },
        [=](const int k, const long r) {
          const VofCurvJob& J = T.job[k];
          int x0, x1;
          const long i0 = vofCurvRow(J, r, false, x0, x1);
          PECLET_FLOW_OMP_SIMD
          for (int x = x0; x < x1; ++x) {
            J.kap[i0 + x] = 0.0;
            J.br[i0 + x] = static_cast<double>(kCurvNone);
          }
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_curv_hf_reset", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const VofCurvJob& J = T.job[k];
        const long i = vofCurvCell(J, t - T.off[k], false);
        J.kap[i] = 0.0;
        J.br[i] = static_cast<double>(kCurvNone);
      });
}
 
inline constexpr int kVofTier3MaxTeams = 1 << 16;
inline void vofCurvFallbackTeams(const VofCurvTable& T, LField list, LField start, LField end) {
  using Policy = Kokkos::TeamPolicy<SExec>;
  constexpr int side = 2 * kPvHalf + 1, nk = side * side * side;
  const int league =
      std::max(1, std::min(kVofTier3MaxTeams, static_cast<int>(SExec().concurrency() / kVofWarp)));
  Policy pol(SExec(), league, kVofWarp);
  pol.set_scratch_size(0, Kokkos::PerTeam(VofPvTeamScratch::shmem(nk)));
  Kokkos::parallel_for(
      "vof::block::batch_curv_pv_team", pol, KOKKOS_LAMBDA(const typename Policy::member_type& tm) {
        const VofPvTeamScratch sc(tm, nk);
        long cum[kVofBlockBatch + 1];  // device prefix of the jobs' entry counts
        cum[0] = 0;
        for (int k = 0; k < T.nj; ++k)
          cum[k + 1] = cum[k] + (end(T.base + k) - start(T.base + k));
        for (long g = tm.league_rank(); g < cum[T.nj]; g += tm.league_size()) {
          int k = 0;
          while (g >= cum[k + 1])
            ++k;
          const VofCurvJob& J = T.job[k];
          const long i = list(start(T.base + k) + (g - cum[k]));
          const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
          curvFallbackTeam(tm, sc, i, VofPvCached{J.slot, T.cache}, VofRawField{J.kap},
                           VofRawField{J.br}, sy, sz, kPvHalf, J.dW, J.cmin, J.gm);
          tm.team_barrier();  // the single lane is done with the scratch before the next entry's
                              // scan
        }
      });
}
 
KOKKOS_INLINE_FUNCTION void vofCurvListEntry(const VofCurvTable& T, int pass, int k, long q,
                                             LField list, LField cnt) {
  const VofCurvJob& J = T.job[k];
  const long i = list(q);
  const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
  const VofRawField c{J.c}, mx{J.mx}, my{J.my}, mz{J.mz}, al{J.al}, kap{J.kap}, br{J.br};
  const bool all = pass == 4;
  if (pass == 0 || all) {
    curvHeightCell(i, c, mx, my, mz, al, kap, br, 1L, sy, sz, J.mtol, J.ptW, J.ieps, J.forceFb,
                   J.oneDir, J.useFit, J.gm, J.peps);
    if (!all)
      return;
  }
  if (pass == 1 || all) {
    curvFallbackCell(i, VofPvCached{J.slot, T.cache}, kap, br, sy, sz, kPvHalf, J.dW, J.cmin, J.gm);
    if (!all)
      return;
  }
  long* ct = &cnt(8 * (T.base + k));
  if (pass == 2 || all) {
    const double km = J.km;
    // km <= 0: this cascade's clip is off (VofCurvature::clipPass returns 0).
    // `!(|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 (km > 0.0 && csfKappaDefined(br(i)) && !(Kokkos::fabs(kap(i)) <= km)) {
      kap(i) = (kap(i) == kap(i)) ? Kokkos::copysign(km, kap(i)) : km;
      Kokkos::atomic_add(&ct[7], 1L);
    }
    if (!all)
      return;
  }
  const int b = static_cast<int>(br(i));
  if (b == kCurvNone)
    return;
  Kokkos::atomic_add(&ct[0], 1L);
  if (b == kCurvHf)
    Kokkos::atomic_add(&ct[1], 1L);
  else if (b == kCurvHfMixed)
    Kokkos::atomic_add(&ct[2], 1L);
  else if (b == kCurvHfFit)
    Kokkos::atomic_add(&ct[3], 1L);
  else if (b == kCurvPv)
    Kokkos::atomic_add(&ct[4], 1L);
  else if (b == kCurvPvReduced)
    Kokkos::atomic_add(&ct[5], 1L);
  else
    Kokkos::atomic_add(&ct[6], 1L);
}
 
template <class Exec = SExec>
inline void vofCurvListPass(const VofCurvTable& T, int pass, LField list, LField start, LField end,
                            LField cnt) {
  if constexpr (kVofHostExec<Exec>) {  // tier 3 included: the host keeps the loop
    vofHostListFor<Exec>("vof::block::batch_curv_list", T.nj, vofListCount(T, start, end),
                         [=](const int k, const long t) {
                           vofCurvListEntry(T, pass, k, start(T.base + k) + t, list, cnt);
                         });
    return;
  }
  if constexpr (!Kokkos::SpaceAccessibility<Kokkos::HostSpace, SField::memory_space>::accessible) {
    if (pass == 1) {  // tier 3 on a device: persistent warp-teams (C3); the host keeps the loop
      vofCurvFallbackTeams(T, list, start, end);
      return;
    }
  }
  Kokkos::parallel_for(
      "vof::block::batch_curv_list", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        const long q = start(T.base + k) + (t - T.off[k]);
        if (q >= end(T.base + k))
          return;
        vofCurvListEntry(T, pass, k, q, list, cnt);
      });
}
 
template <class Exec = SExec>
inline void vofCurvListPasses(const VofCurvTable& T, bool anyClip, LField list, LField start,
                              LField end, LField cnt) {
  if constexpr (kVofHostExec<Exec>) {
    vofCurvListPass<Exec>(T, 4, list, start, end, cnt);
    return;
  }
  vofCurvListPass<Exec>(T, 0, list, start, end, cnt);  // tiers 1-2
  vofCurvListPass<Exec>(T, 1, list, start, end, cnt);  // tier 3
  if (anyClip)
    vofCurvListPass<Exec>(T, 2, list, start, end, cnt);  // the admissibility clip
  vofCurvListPass<Exec>(T, 3, list, start, end, cnt);    // the census
}
 
KOKKOS_INLINE_FUNCTION void vofCsfForceCell(const VofCurvJob& J, long i, double sig, double w0,
                                            double w1, double w2) {
  const long sy = J.e.x, sz = static_cast<long>(J.e.x) * J.e.y;
  const VofRawField cv{J.c}, kp{J.kap}, kb{J.br};
  for (int cc = 0; cc < 3; ++cc) {
    const double wc = cc == 0 ? w0 : (cc == 1 ? w1 : w2);
    const long strd = (cc == 0) ? 1 : (cc == 1 ? sy : sz);
    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);
    }
    J.f[cc][i] = f;
  }
}
 
template <class Exec = SExec>
inline void vofCsfForceBatch(const VofCurvTable& T, double sig, double w0, double w1, double w2) {
  if constexpr (kVofHostExec<Exec>) {
    vofHostRowFor<Exec>(
        "vof::block::batch_csf_force", T.nj,
        [&](const int k) { return vofCurvRows(T.job[k], false); },
        [=](const int k, const long r) {
          const VofCurvJob& J = T.job[k];
          int x0, x1;
          const long i0 = vofCurvRow(J, r, false, x0, x1);
          PECLET_FLOW_OMP_SIMD
          for (int x = x0; x < x1; ++x)
            vofCsfForceCell(J, i0 + x, sig, w0, w1, w2);
        });
    return;
  }
  Kokkos::parallel_for(
      "vof::block::batch_csf_force", Kokkos::RangePolicy<SExec>(SExec(), 0, T.off[T.nj]),
      KOKKOS_LAMBDA(const long t) {
        const int k = vofJobOf(T.off, t);
        vofCsfForceCell(T.job[k], vofCurvCell(T.job[k], t - T.off[k], false), sig, w0, w1, w2);
      });
}
 
// ---- C2: the block statistics (§5.10) ----------------------------------------------------------
 
template <int N>
struct VofSumN {
  double v[N];
  KOKKOS_INLINE_FUNCTION VofSumN() {
    for (int i = 0; i < N; ++i)
      v[i] = 0.0;
  }
  KOKKOS_INLINE_FUNCTION VofSumN& operator+=(const VofSumN& o) {
    for (int i = 0; i < N; ++i)
      v[i] += o.v[i];
    return *this;
  }
};
 
}  // namespace peclet::flow::vof
 
namespace Kokkos {
template <int N>
struct reduction_identity<peclet::flow::vof::VofSumN<N>> {
  KOKKOS_FORCEINLINE_FUNCTION static peclet::flow::vof::VofSumN<N> sum() {
    return peclet::flow::vof::VofSumN<N>();
  }
};
}  // namespace Kokkos
 
namespace peclet::flow::vof {
 
struct VofStatsPrev {
  double c[kVofBlockBatch][3];
  unsigned char valid[kVofBlockBatch];
};
 
inline constexpr int kVofStats = 14;
 
inline void vofBatchStats(const VofBlockTable& T, const VofStatsPrev& P, double hh, double aeps,
                          double dt, SField out, long off) {
  using Team = Kokkos::TeamPolicy<SExec>;
  Kokkos::parallel_for(
      "vof::block::batch_stats", Team(SExec(), T.nj, Kokkos::AUTO),
      KOKKOS_LAMBDA(const typename Team::member_type& tm) {
        const int k = tm.league_rank();
        const VofBlockJob& J = T.job[k];
        const I3 e = J.e, n = J.n, o = J.o;
        const int g = J.g;
        const long sy = e.x, sz = static_cast<long>(e.x) * e.y;
        const long region = static_cast<long>(n.x) * n.y * n.z;
        const VofRawField c{J.c};
        double* q0 = out.data() + off + static_cast<long>(kVofStats) * (T.base + k);
        VofSumN<4> m1;
        Kokkos::parallel_reduce(
            Kokkos::TeamThreadRange(tm, region),
            [&](const long r, VofSumN<4>& a) {
              const int x = static_cast<int>(r % n.x) + g;
              const int y = static_cast<int>((r / n.x) % n.y) + g;
              const int z = static_cast<int>(r / (static_cast<long>(n.x) * n.y)) + g;
              const double q = c(L3(x, y, z, e));
              a.v[0] += q;
              a.v[1] += q * (x - g + o.x + 0.5) * hh;
              a.v[2] += q * (y - g + o.y + 0.5) * hh;
              a.v[3] += q * (z - g + o.z + 0.5) * hh;
            },
            Kokkos::Sum<VofSumN<4>>(m1));
        const double v = m1.v[0];
        if (v <= 0.0) {  // measure()'s empty-marker branch: everything but the volume is zero
          Kokkos::single(Kokkos::PerTeam(tm), [&]() {
            q0[0] = v;
            for (int d = 1; d < kVofStats; ++d)
              q0[d] = 0.0;
          });
          return;
        }
        const double cx = m1.v[1] / v, cy = m1.v[2] / v, cz = m1.v[3] / v;
        VofSumN<7> m2;
        Kokkos::parallel_reduce(
            Kokkos::TeamThreadRange(tm, region),
            [&](const long r, VofSumN<7>& a) {
              const int x = static_cast<int>(r % n.x) + g;
              const int y = static_cast<int>((r / n.x) % n.y) + g;
              const int z = static_cast<int>(r / (static_cast<long>(n.x) * n.y)) + g;
              const long i = L3(x, y, z, e);
              const double q = c(i);
              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;
              a.v[0] += q * px * px;
              a.v[1] += q * py * py;
              a.v[2] += q * pz * pz;
              a.v[3] += q * px * py;
              a.v[4] += q * px * pz;
              a.v[5] += q * py * pz;
              if (!(q > aeps) || !(q < 1.0 - aeps))  // interfaceArea: mixed cells only
                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], q);
              double pv[8][3], ctr[3], ar = 0.0;
              const int nv = vof::plicPolygon(m[0], m[1], m[2], al, pv);
              vof::polygonAreaCentroid(pv, nv, ctr, ar);
              a.v[6] += ar;
            },
            Kokkos::Sum<VofSumN<7>>(m2));
        Kokkos::single(Kokkos::PerTeam(tm), [&]() {
          const double nc[3] = {cx, cy, cz};
          q0[0] = v;
          for (int d = 0; d < 3; ++d) {
            q0[1 + d] = nc[d];
            q0[4 + d] = (P.valid[k] && dt > 0.0) ? (nc[d] - P.c[k][d]) / dt : 0.0;
          }
          q0[7] = m2.v[6];
          const double iv = 1.0 / v;
          for (int d = 0; d < 6; ++d)
            q0[8 + d] = m2.v[d] * iv;
        });
      });
}
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_BLOCK_BATCH_HPP