cf_scheme.hpp

Go to the documentation of this file.

// core — pluggable coarse/fine (2:1) interface schemes for the collocated AMR flow solver.
//
// The graded-accuracy limiter measured in docs/amr_collocated_projection.md (addendum): the
// standard two-point C/F flux is 1st-order at level boundaries, and BOTH the aperture and the
// ghost projection inherit it identically (+9% graded drag offset on the dilute sphere,
// scheme-independent). This header makes the C/F treatment a pluggable SCHEME and provides the
// Martin–Cartwright/Martin–Colella tangential-quadratic closure (CfScheme::quadratic — the same
// coarseStar substitution AmrPoisson::applyLaplacianQuad uses) for every operator the STEADY
// solution feels:
//
//   D  — the ½/½ face-average divergence constraint      (buildCfDivDelta,  vector → scalar)
//   G  — the ABC cell gradient (−∇pⁿ predictor + cell correction)  (buildCfGradDelta, scalar → 3)
//   ∇² — the momentum (velocity) diffusion at regular fluid rows    (buildCfLapDelta,  scalar →
//   scalar)
//
// Everything is the LINEAR substitution "coarse-side value → coarse*", where coarse* is the
// scheme's interpolation of the coarse cell at the fine cell's tangential position. So each
// operator's (quad − standard) difference is a precomputed CSR overlay over the level-boundary
// rows: built ONCE on the host (shared verbatim by the host oracle and the device — parity by
// construction), applied as one extra SpMV.
//
// STEADY-STATE PLACEMENT (the (1,2)-mixed philosophy, flow-validated): at the projection's fixed
// point φ→0, so the pressure MATRIX C/F order does not move the steady solution — the matrix (and
// the whole MG hierarchy / PCG / ghost BiCGStab) stays on the standard consistent operator, and
// the scheme enters through (a) the RHS divergence, (b) the pressure gradients, (c) the momentum
// operator (as a lagged deferred-correction RHS term, the same pattern as the implicit-FOU/SOU
// split and Multigrid::solveQuad). The per-step operator mismatch converges through the
// pseudo-transient stepping; the steady equations carry the 2nd-order C/F closure exactly.
//
// EXTENDING: add a value to CfScheme and a branch in cfAppendStencil — every operator delta picks
// it up. A scheme is fully described by the linear stencil it substitutes for the coarse-side
// value of one directed C/F sub-face.
//
// GATING — one predicate, per FACE, never per row (docs/amr_cf_flux_gate.md). The quadratic face
// value applies on a 2:1 sub-face IFF BOTH incident cells are REGULAR FLUID:
//
//     cfFace(i, j) := level(j) != level(i) && regular(i) && regular(j),   regular = fluid && !cut
//
// and that ONE predicate drives D, the G substitution and uf alike, through the one emitter
// cfAppendFaceValueDelta, so that D_std(Δuf) ≡ Δ_cfDiv holds entry by entry on every mesh. A face
// flux is a single number shared by two cells, so a gate that is a property of a CELL cannot be
// conservative: the old row gate let the regular cell book the correction in its constraint while
// the cut cell did not and uf carried it for both — a mass source proportional to the velocity,
// not to φ, so it did not decay at steady state. Withholding reverts the face to the standard
// two-point value (the cf=0 scheme), which is an O(h) face value on the set where a level
// boundary meets the wall — a curve, i.e. codimension 2, on the production graded meshes.
// `AmrFlow::numCfCutFaces()` counts that set, so the cost is measured rather than silent.
//
// The per-tangential-axis robustness fallbacks INSIDE the stencil are unchanged and independent
// of the above: the tangential coarse neighbours must exist, be same-level, be FLUID (the
// caller's `fluidOk` — never lean on a decoupled/held solid value), and the tangential faces must
// be sufficiently open (openness ≥ 0.5, matching applyLaplacianQuad's gate); one surviving side
// gives the linear closure, neither gives the raw coarse value. Reading a CUT cell as a
// tangential sample is fine — its velocity is a solved fluid value; the instability the row gate
// was added for is about which ROW owns a constraint, not which cells it reads.
#ifndef PECLET_AMR_CF_SCHEME_HPP
#define PECLET_AMR_CF_SCHEME_HPP
 
#include <array>
#include <cstdint>
#include <utility>
#include <vector>
 
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/common.hpp"
#include "peclet/amr/poisson.hpp"
#include "peclet/core/common/host_parallel.hpp"
#include "peclet/core/common/types.hpp"
 
namespace peclet::amr {
 
enum class CfScheme : int {
  standard = 0,   
  quadratic = 1,  
};
 
struct CfCsr {
  std::vector<Index> start;  
  std::vector<Index> slot;
  std::vector<double> coef;
};
 
struct CfCompCsr {
  std::vector<Index> start;
  std::vector<Index> slot;
  std::vector<double> coef;
  std::vector<int8_t> comp;
};
 
namespace detail {
 
template <unsigned Bits, class Ent, class FluidFn>
inline void cfAppendStencil(const AmrPoisson<3, Bits>& ap, std::vector<Ent>& out, Index coarse,
                            Index fine, int axis, double scale, FluidFn&& fluidOk, CfScheme scheme,
                            const Ent& proto = Ent{}) {
  if (scheme != CfScheme::quadratic)
    return;
  // DISTRIBUTED (2026-09-21): every geometric query here goes through `ap`, never through the
  // block octree. `coarse`/`fine` may be GHOST slots (loOf/levelOf serve those transparently) and
  // the tangential probes may leave the block (probeSlot routes them through the LeafHalo
  // resolver). Single-rank both spellings are the same lookup on the same coordinate — block ==
  // domain and probeSlot's periodic wrap is the one this function used to do by hand — so the
  // built CSRs are bit-identical there.
  const std::array<long, 3> bc = ap.loOf(coarse);
  const std::array<long, 3> bf = ap.loOf(fine);
  const unsigned Lc = ap.levelOf(coarse);
  const double sc = static_cast<double>(Index(1) << Lc);
  const double sf = static_cast<double>(Index(1) << ap.levelOf(fine));
  auto push = [&](Index cell, double w) {
    Ent e = proto;
    e.cell = cell;
    e.w = w;
    out.push_back(e);
  };
  for (int tt = 0; tt < 3; ++tt) {
    if (tt == axis)
      continue;
    // Phase 3: BOTH the tangential offset and the differencing width belong to the TANGENTIAL
    // axis `tt` — they were one `H` when the cells were cubes (`docs/amr_anisotropic.md` §3).
    const double H = ap.cellWidth(coarse, tt);
    // Phase 3: the tangential offset is a length along axis `tt` (`docs/amr_anisotropic.md` §3).
    const double dt =
        ((static_cast<double>(bf[tt]) + 0.5 * sf) - (static_cast<double>(bc[tt]) + 0.5 * sc)) *
        ap.h0()[tt];
    // Tangential samples at the coarse cell's ± neighbours. Same-level neighbours contribute
    // directly. A FINER neighbour (an island corner/edge: the region across the tangential face
    // is refined — by 2:1 exactly one level finer) is sampled by the volume average of the 2^Dim
    // leaves covering the coarse-size region: a 2nd-order sample of the region-center value, so
    // the corner rows keep the quadratic closure instead of falling back to the raw coarse value
    // (the old P5b-style skip, which left them locally 1st-order). A COARSER neighbour (or any
    // non-fluid / irregular cover) still falls back to the skip (robust).
    struct Samp {
      Index cell[1 << 3];
      double w[1 << 3];
      int n = 0;
      bool ok = false;
    };
    auto sampleTangential = [&](int dir) {
      Samp sm;
      const Index nb = ap.periodicNeighbor(coarse, tt, dir);
      if (nb < 0)
        return sm;
      const unsigned Lnb = ap.levelOf(nb);
      if (Lnb == Lc) {
        if (!fluidOk(nb))
          return sm;
        sm.cell[0] = nb;
        sm.w[0] = 1.0;
        sm.n = 1;
        sm.ok = true;
        return sm;
      }
      if (Lnb + 1 != Lc)
        return sm;  // coarser neighbour: keep the fallback
      // Finer neighbour: enumerate the 2^Dim children covering the coarse-size region. The child
      // probes go through ap.probeSlot — NOT a hand-wrapped BlockOctree::find, which is what this
      // did until 2026-09-21. probeSlot wraps in-block coords itself and hands out-of-block ones
      // to the LeafHalo resolver, so a child owned by another rank resolves to its ghost slot.
      // The hand wrap used the BLOCK period, which is the domain period only single-rank;
      // multi-rank it folded a neighbour rank's child back into this block and `find` returned a
      // real but unrelated leaf — either the guard below rejected it and the tangential side
      // dropped silently to the raw coarse value, or it passed and the stencil read the wrong
      // cell. Measured cost of that, np = 2 on tests/test_amr_distributed_cf_mpi: 131 % error.
      const long scL = 1L << Lc, shL = scL >> 1;
      std::array<long, 3> lo = bc;
      lo[tt] += (dir > 0) ? scL : -scL;
      for (int oct2 = 0; oct2 < (1 << 3); ++oct2) {
        std::array<long, 3> q = lo;
        for (int d = 0; d < 3; ++d)
          if ((oct2 >> d) & 1)
            q[d] += shL;
        const auto [ch, Lch] = ap.probeSlot(q);
        if (ch < 0 || Lch + 1 != Lc || !fluidOk(ch))
          return Samp{};  // irregular cover or solid child: fall back
        sm.cell[sm.n] = ch;
        sm.w[sm.n] = 1.0 / static_cast<double>(1 << 3);
        ++sm.n;
      }
      sm.ok = true;
      return sm;
    };
    const Samp sp = sampleTangential(+1);
    const Samp smi = sampleTangential(-1);
    const bool okP = sp.ok && ap.faceOpenness(coarse, tt, +1) >= 0.5;
    const bool okM = smi.ok && ap.faceOpenness(coarse, tt, -1) >= 0.5;
    if (okP && okM) {
      // Two-sided (the original path, bit-identical): tangential quadratic.
      const double cUp = dt / (2.0 * H) + 0.5 * dt * dt / (H * H);
      const double cUm = -dt / (2.0 * H) + 0.5 * dt * dt / (H * H);
      const double cUc = -dt * dt / (H * H);
      push(coarse, scale * cUc);
      for (int k2 = 0; k2 < sp.n; ++k2)
        push(sp.cell[k2], scale * cUp * sp.w[k2]);
      for (int k2 = 0; k2 < smi.n; ++k2)
        push(smi.cell[k2], scale * cUm * smi.w[k2]);
      continue;
    }
    // Wall-aware one-sided fallback (mixed-level cut-band plan §6.1): when exactly one
    // tangential side survives the gates (the other is solid / closed / irregular — the
    // near-wall seam case), substitute the one-sided LINEAR interpolation dt·u' instead of
    // falling all the way back to the raw coarse value: O(H²) tangential sample error vs the
    // raw fallback's O(H) offset. Both sides gated ⇒ raw fallback (unchanged).
    if (okP != okM) {
      const Samp& s1 = okP ? sp : smi;
      const double sgn = okP ? 1.0 : -1.0;
      push(coarse, -scale * sgn * dt / H);
      for (int k2 = 0; k2 < s1.n; ++k2)
        push(s1.cell[k2], scale * sgn * (dt / H) * s1.w[k2]);
    }
  }
}
 
struct ScalarEnt {
  Index cell = 0;
  double w = 0.0;
};
struct CompEnt {
  Index cell = 0;
  double w = 0.0;
  int8_t comp = 0;
};
 
template <unsigned Bits, class Ent, class FluidFn>
inline void cfAppendFaceValueDelta(const AmrPoisson<3, Bits>& ap, std::vector<Ent>& out,
                                   Index coarse, Index fine, int axis, double dist, double scale,
                                   FluidFn&& fluidOk, CfScheme scheme, const Ent& proto = Ent{}) {
  const double H = ap.cellWidth(coarse, axis), h = ap.cellWidth(fine, axis);
  const double wF = (0.5 * H) / dist, wC = (0.5 * h) / dist;
  Ent eF = proto, eC = proto;
  eF.cell = fine;
  eF.w = scale * (wF - 0.5);
  eC.cell = coarse;
  eC.w = scale * (wC - 0.5);
  out.push_back(eF);
  out.push_back(eC);
  cfAppendStencil(ap, out, coarse, fine, axis, scale * wC, fluidOk, scheme, proto);
}
 
template <class Ent, class Csr>
inline void compactCsr(const std::vector<std::vector<Ent>>& per, Csr& out, Index n) {
  out.start.assign(static_cast<std::size_t>(n) + 1, 0);
  for (Index i = 0; i < n; ++i)
    out.start[static_cast<std::size_t>(i) + 1] =
        out.start[static_cast<std::size_t>(i)] +
        static_cast<Index>(per[static_cast<std::size_t>(i)].size());
  const Index nz = out.start[static_cast<std::size_t>(n)];
  out.slot.resize(static_cast<std::size_t>(nz));
  out.coef.resize(static_cast<std::size_t>(nz));
  Index k = 0;
  for (Index i = 0; i < n; ++i)
    for (const auto& e : per[static_cast<std::size_t>(i)]) {
      out.slot[static_cast<std::size_t>(k)] = e.cell;
      out.coef[static_cast<std::size_t>(k)] = e.w;
      ++k;
    }
}
 
}  // namespace detail
 
template <unsigned Bits, class RowFn, class FluidFn>
inline CfCsr buildCfLapDelta(const AmrPoisson<3, Bits>& ap, const BlockOctree<3, Bits>& t,
                             double factor, RowFn&& rowOk, FluidFn&& fluidOk, CfScheme scheme) {
  const Index n = t.numLeaves();
  std::vector<std::vector<detail::ScalarEnt>> per(static_cast<std::size_t>(n));
  // Host-parallel (setup-parallel plan rung 1): row i writes ONLY per[i] and reads only the frozen
  // octree/openness — disjoint writes, so the staged rows are bitwise identical under any schedule
  // and compactCsr below still concatenates them in leaf order.
  hostParFor(n, [&](Index i) {
    if (!rowOk(i))
      return;
    const unsigned Li = t.level(i);
    const double invV = 1.0 / ap.cellVolume(i);
    ap.forEachFaceNeighbor(i, [&](Index j, Real c, int axis, double a) {
      const unsigned Lj = ap.levelOf(j);  // ghost-safe (j may be a ghost slot)
      if (Lj == Li)
        return;
      const Index coarse = (Lj > Li) ? j : i;
      const Index fine = (Lj > Li) ? i : j;
      const double scale = factor * invV * (a * c) * ((Lj > Li) ? 1.0 : -1.0);
      detail::cfAppendStencil(ap, per[static_cast<std::size_t>(i)], coarse, fine, axis, scale,
                              fluidOk, scheme);
    });
  });
  CfCsr csr;
  detail::compactCsr(per, csr, n);
  return csr;
}
 
template <unsigned Bits, class RegularFn, class FluidFn>
inline CfCompCsr buildCfDivDelta(const AmrPoisson<3, Bits>& ap, const BlockOctree<3, Bits>& t,
                                 RegularFn&& regularOk, FluidFn&& fluidOk, CfScheme scheme) {
  const Index n = t.numLeaves();
  std::vector<std::vector<detail::CompEnt>> per(static_cast<std::size_t>(n));
  if (scheme != CfScheme::standard) {
    hostParFor(n, [&](Index i) {  // per[i]-disjoint (rung 1)
      if (!regularOk(i))
        return;
      const unsigned Li = t.level(i);
      const double invV = 1.0 / ap.cellVolume(i);
      ap.forEachFaceFull(
          i, [&](Index j, int axis, int dir, double area, double dist, double alpha) {
            const unsigned Lj = ap.levelOf(j);  // ghost-safe
            // cfFace(i, j): regularOk(i) holds by the early-out above.
            if (Lj == Li || !regularOk(j))
              return;
            const Index coarse = (Lj > Li) ? j : i;
            const Index fine = (Lj > Li) ? i : j;
            const double scale = invV * alpha * area * static_cast<double>(dir);
            detail::CompEnt proto;
            proto.comp = static_cast<int8_t>(axis);
            detail::cfAppendFaceValueDelta(ap, per[static_cast<std::size_t>(i)], coarse, fine, axis,
                                           dist, scale, fluidOk, scheme, proto);
          });
    });
  }
  CfCompCsr csr;
  detail::compactCsr(per, csr, n);
  csr.comp.resize(csr.slot.size());
  Index k = 0;
  for (Index i = 0; i < n; ++i)
    for (const auto& e : per[static_cast<std::size_t>(i)])
      csr.comp[static_cast<std::size_t>(k++)] = e.comp;
  return csr;
}
 
template <unsigned Bits, class RegularFn, class FluidFn>
inline std::array<CfCsr, 3> buildCfGradDelta(const AmrPoisson<3, Bits>& ap,
                                             const BlockOctree<3, Bits>& t, RegularFn&& regularOk,
                                             FluidFn&& fluidOk, CfScheme scheme) {
  const Index n = t.numLeaves();
  std::array<std::vector<std::vector<detail::ScalarEnt>>, 3> per;
  for (int a = 0; a < 3; ++a)
    per[static_cast<std::size_t>(a)].resize(static_cast<std::size_t>(n));
  if (scheme != CfScheme::standard) {
    hostParFor(n, [&](Index i) {  // per[axis][i]-disjoint (rung 1)
      if (!regularOk(i))
        return;
      const unsigned Li = t.level(i);
      // Per axis/side: open-face count, the side's (uniform) center-to-center distance, and
      // whether the side has a level mismatch. slot [axis][side], side 0 = +, 1 = −.
      int cnt[3][2] = {};
      double sdist[3][2] = {};
      bool cf[3][2] = {};
      ap.forEachFaceFull(i, [&](Index j, int axis, int dir, double, double dist, double alpha) {
        if (alpha <= 1e-12)
          return;
        const int s = (dir > 0) ? 0 : 1;
        ++cnt[axis][s];
        sdist[axis][s] = dist;
        // D and G must be paired FACE BY FACE (docs/amr_cf_flux_gate.md §6.3, §8 point 3): a
        // face whose value D takes as standard must be one whose gradient G takes as standard, or
        // the constraint gains a component the gradient cannot see. So the side reweighting is
        // triggered by a PASSING C/F face, not by any level mismatch.
        if (ap.levelOf(j) != Li && regularOk(j))  // cfFace(i, j); regularOk(i) by the early-out
          cf[axis][s] = true;
      });
      for (int a = 0; a < 3; ++a) {
        if (!cf[a][0] && !cf[a][1])
          continue;  // no level mismatch on this axis: standard ½/½ is centered
        if (cnt[a][0] == 0 || cnt[a][1] == 0)
          continue;  // one-sided (closed) axis: leave the standard treatment
        const double dp = sdist[a][0], dm = sdist[a][1];
        const double wp = dm / (dp + dm), wm = dp / (dp + dm);
        auto& row = per[static_cast<std::size_t>(a)][static_cast<std::size_t>(i)];
        ap.forEachFaceFull(i, [&](Index j, int axis, int dir, double, double dist, double alpha) {
          if (axis != a || alpha <= 1e-12)
            return;
          const int s = (dir > 0) ? 0 : 1;
          const double w = (s == 0) ? wp : wm;
          const double inv = 1.0 / static_cast<double>(cnt[a][s]);
          const double gsgn = (dir > 0) ? 1.0 : -1.0;  // coefficient of f_j in g
          // (a) side reweighting of the STANDARD face gradient: (w − ½)·g/cnt.
          const double rw = (w - 0.5) * inv * gsgn / dist;
          if (rw != 0.0) {
            row.push_back({j, rw});
            row.push_back({i, -rw});
          }
          // (b) coarse* substitution inside the C/F face gradients, at the NEW side weight.
          const unsigned Lj = ap.levelOf(j);  // ghost-safe
          if (Lj == Li || !regularOk(j))      // cfFace(i, j)
            return;
          const Index coarse = (Lj > Li) ? j : i;
          const Index fine = (Lj > Li) ? i : j;
          const double ssgn = (Lj > Li) ? gsgn : -gsgn;  // sign of the substituted value in g
          detail::cfAppendStencil(ap, row, coarse, fine, a, w * inv * ssgn / dist, fluidOk, scheme);
        });
      }
    });
  }
  std::array<CfCsr, 3> out;
  for (int a = 0; a < 3; ++a)
    detail::compactCsr(per[static_cast<std::size_t>(a)], out[static_cast<std::size_t>(a)], n);
  return out;
}
 
struct CfUfDelta {
  CfCompCsr vel;  
};
 
template <unsigned Bits, class RegularFn, class FluidFn>
inline CfUfDelta buildCfUfDelta(const AmrPoisson<3, Bits>& ap, const BlockOctree<3, Bits>& t,
                                RegularFn&& regularOk, FluidFn&& fluidOk, CfScheme scheme) {
  const Index n = t.numLeaves();
  // Rung 1, two-pass. forEachFaceFull is CELL-MAJOR, so leaf i owns the contiguous slot range
  // [slotBase[i], slotBase[i+1]): count the slots per leaf in parallel, then an INTEGER exclusive
  // scan (exact under any combination order) reproduces today's serial slot numbering exactly.
  std::vector<Index> slotCnt(static_cast<std::size_t>(n), 0);
  hostParFor(n, [&](Index i) {
    Index c = 0;
    ap.forEachFaceFull(i, [&](Index, int, int, double, double, double) { ++c; });
    slotCnt[static_cast<std::size_t>(i)] = c;
  });
  std::vector<Index> slotBase(static_cast<std::size_t>(n) + 1, 0);
  hostParScan(n, [&](Index i, Index& acc, bool final) {
    if (final)
      slotBase[static_cast<std::size_t>(i)] = acc;
    acc += slotCnt[static_cast<std::size_t>(i)];
  });
  const Index nSlots =
      (n > 0) ? slotBase[static_cast<std::size_t>(n) - 1] + slotCnt[static_cast<std::size_t>(n) - 1]
              : 0;
  slotBase[static_cast<std::size_t>(n)] = nSlots;
  CfUfDelta d;
  d.vel.start.assign(static_cast<std::size_t>(nSlots) + 1, 0);
  if (scheme == CfScheme::standard)
    return d;
  // Pass 1 (parallel): each leaf stages its own entries in its own slot order, and writes its
  // slots' ROW WIDTHS into start[slot+1] (disjoint: one leaf owns each slot). The prefix sum is
  // taken serially below, exactly as the single-pass version accumulated it.
  std::vector<std::vector<detail::CompEnt>> velPer(static_cast<std::size_t>(n));
  hostParFor(n, [&](Index i) {
    const unsigned Li = t.level(i);
    Index slot = slotBase[static_cast<std::size_t>(i)];
    auto& vrow = velPer[static_cast<std::size_t>(i)];
    ap.forEachFaceFull(i, [&](Index j, int axis, int, double, double dist, double) {
      const unsigned Lj = ap.levelOf(j);               // ghost-safe
      if (Lj != Li && regularOk(i) && regularOk(j)) {  // cfFace(i, j); regular ⇒ fluid
        const Index coarse = (Lj > Li) ? j : i;
        const Index fine = (Lj > Li) ? i : j;
        std::vector<detail::CompEnt> ve;
        detail::CompEnt proto;
        proto.comp = static_cast<int8_t>(axis);
        detail::cfAppendFaceValueDelta(ap, ve, coarse, fine, axis, dist, 1.0, fluidOk, scheme,
                                       proto);
        for (const auto& e : ve)
          vrow.push_back(e);
        d.vel.start[static_cast<std::size_t>(slot) + 1] = static_cast<Index>(ve.size());
      }
      ++slot;
    });
  });
  // Pass 2 (serial, order-preserving): widths → offsets, then concatenate leaf-major. Identical to
  // the single pass's append order (leaf-major, slot-major within a leaf, entry order within slot).
  for (Index s = 0; s < nSlots; ++s)
    d.vel.start[static_cast<std::size_t>(s) + 1] += d.vel.start[static_cast<std::size_t>(s)];
  d.vel.slot.reserve(static_cast<std::size_t>(d.vel.start[static_cast<std::size_t>(nSlots)]));
  d.vel.coef.reserve(static_cast<std::size_t>(d.vel.start[static_cast<std::size_t>(nSlots)]));
  d.vel.comp.reserve(static_cast<std::size_t>(d.vel.start[static_cast<std::size_t>(nSlots)]));
  for (Index i = 0; i < n; ++i)
    for (const auto& e : velPer[static_cast<std::size_t>(i)]) {
      d.vel.slot.push_back(e.cell);
      d.vel.coef.push_back(e.w);
      d.vel.comp.push_back(e.comp);
    }
  return d;
}
 
template <unsigned Bits, class RegularFn, class FluidFn>
inline Index countCfCutFaces(const AmrPoisson<3, Bits>& ap, const BlockOctree<3, Bits>& t,
                             RegularFn&& regularOk, FluidFn&& fluidOk) {
  const Index n = t.numLeaves();
  // Per-leaf counts then a serial sum: disjoint writes (rung 1), and the total is an exact
  // integer under any schedule. One face traversal, the same shape as the CSR builders beside it.
  std::vector<Index> per(static_cast<std::size_t>(n), 0);
  hostParFor(n, [&](Index i) {
    if (!fluidOk(i))
      return;
    const unsigned Li = t.level(i);
    Index c = 0;
    ap.forEachFaceFull(i, [&](Index j, int, int, double, double, double) {
      if (ap.levelOf(j) == Li || !fluidOk(j))
        return;
      if (regularOk(i) && regularOk(j))
        return;
      ++c;
    });
    per[static_cast<std::size_t>(i)] = c;
  });
  Index cnt = 0;
  for (Index i = 0; i < n; ++i)
    cnt += per[static_cast<std::size_t>(i)];
  return cnt;
}
 
// ---- host applies (the oracle path; the device uses the same CSRs uploaded + SpMV kernels) ----
 
inline void cfApplyHost(const CfCsr& c, const std::vector<double>& f, std::vector<double>& out) {
  const Index n = static_cast<Index>(c.start.size()) - 1;
  for (Index i = 0; i < n; ++i) {
    double acc = 0.0;
    for (Index k = c.start[static_cast<std::size_t>(i)];
         k < c.start[static_cast<std::size_t>(i) + 1]; ++k)
      acc += c.coef[static_cast<std::size_t>(k)] *
             f[static_cast<std::size_t>(c.slot[static_cast<std::size_t>(k)])];
    out[static_cast<std::size_t>(i)] += acc;
  }
}
 
inline void cfApplyCompHost(const CfCompCsr& c, const std::array<std::vector<double>, 3>& u,
                            std::vector<double>& out) {
  const Index n = static_cast<Index>(c.start.size()) - 1;
  for (Index i = 0; i < n; ++i) {
    double acc = 0.0;
    for (Index k = c.start[static_cast<std::size_t>(i)];
         k < c.start[static_cast<std::size_t>(i) + 1]; ++k)
      acc += c.coef[static_cast<std::size_t>(k)] *
             u[static_cast<std::size_t>(c.comp[static_cast<std::size_t>(k)])]
              [static_cast<std::size_t>(c.slot[static_cast<std::size_t>(k)])];
    out[static_cast<std::size_t>(i)] += acc;
  }
}
 
// ---- device mirrors (Kokkos TUs only; include after a Kokkos-carrying header) ------------------
#ifdef KOKKOS_INLINE_FUNCTION
 
struct CfCsrDev {
  Index n = 0, nz = 0;
  View<Index> start, slot;
  View<double> coef;
};
struct CfCompCsrDev {
  Index n = 0, nz = 0;
  View<Index> start, slot;
  View<double> coef;
  View<int8_t> comp;
};
 
inline CfCsrDev uploadCfCsr(const CfCsr& h, const char* name) {
  CfCsrDev d;
  d.n = static_cast<Index>(h.start.size()) - 1;
  d.nz = static_cast<Index>(h.slot.size());
  if (d.nz == 0) {
    d.n = 0;
    return d;
  }
  d.start = toDevice(h.start, name);
  d.slot = toDevice(h.slot, name);
  d.coef = toDevice(h.coef, name);
  return d;
}
inline CfCompCsrDev uploadCfCompCsr(const CfCompCsr& h, const char* name) {
  CfCompCsrDev d;
  d.n = static_cast<Index>(h.start.size()) - 1;
  d.nz = static_cast<Index>(h.slot.size());
  if (d.nz == 0) {
    d.n = 0;
    return d;
  }
  d.start = toDevice(h.start, name);
  d.slot = toDevice(h.slot, name);
  d.coef = toDevice(h.coef, name);
  d.comp = toDevice(h.comp, name);
  return d;
}
 
inline void cfApply(const CfCsrDev& c, View<const double> f, View<double> out) {
  if (c.n == 0)
    return;
  auto st = c.start;
  auto sl = c.slot;
  auto w = c.coef;
  Kokkos::parallel_for(
      "amr::cf_apply", c.n, KOKKOS_LAMBDA(const Index i) {
        double acc = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k)
          acc += w(k) * f(sl(k));
        out(i) += acc;
      });
}
 
inline void cfApplyComp(const CfCompCsrDev& c, View<const double> u0, View<const double> u1,
                        View<const double> u2, View<double> out) {
  if (c.n == 0)
    return;
  auto st = c.start;
  auto sl = c.slot;
  auto w = c.coef;
  auto cp = c.comp;
  Kokkos::parallel_for(
      "amr::cf_apply_comp", c.n, KOKKOS_LAMBDA(const Index i) {
        double acc = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k) {
          const int a = cp(k);
          const double v = (a == 0) ? u0(sl(k)) : (a == 1) ? u1(sl(k)) : u2(sl(k));
          acc += w(k) * v;
        }
        out(i) += acc;
      });
}
 
#endif  // KOKKOS_INLINE_FUNCTION
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_CF_SCHEME_HPP