ghost_projection_sampled.hpp

Go to the documentation of this file.

// core — SAMPLED ghost-projection overlay: cut cells at MULTIPLE octree levels.
//
// The D1 machinery of docs/amr_mixed_level_cut_band_plan.md: the ghost overlay's ±2
// same-level chain entries become SAMPLE SLOTS — precomputed linear functionals of nearby
// fluid leaves. A same-level entry is the identity (today's behaviour, bit-identical on a
// uniform band); an entry across a 2:1 level boundary is a VIRTUAL SAMPLE at the uniform-grid
// position, reconstructed by DEGREE-2 least squares over fluid leaf point values (the M2
// verdict, tests/study_amr_seam_sample_order.cpp: degree-1 is an O(1) non-decaying matrix
// perturbation at 53% of the physical scale; degree-2 decays ~O(h) to 0.4%). Fallback cascade
// (fluid-only, M1-measured 0.1–0.5%): LS2 -> LS1 -> covering-leaf identity.
//
// CLASSIFICATION CONSISTENCY (the invariant that prevents double-counted fluxes): a face is
// closed in the overlay iff it is closed in the binary operator. Both derive from the SAME
// primitive — a (sub)face is open iff BOTH adjacent ACTUAL leaf centers are fluid and the
// float mean of their center samples is >= 0 (`sampledFaceOpen`; use `makeBinaryOpenFnMixed`
// for AmrPoisson::buildOpenness / AmrMultigrid::setOpenness). A row's own-face state takes
// its SIGN from this canonical rule (any-sub-face-open for a finer-across face — the closed
// sub-faces of a mixed face are then Neumann-zero, a counted Design-A-like local defect) and
// its MAGNITUDE (theta, closure weights) from the row's virtual uniform-position samples (the
// plan's D2 route, level-independent — identical to the classic builder on a uniform band).
//
// Everything is a pure function of (octree, sdf): deterministic under adaptivity (D6).
//
// LAYOUT. The builder and the `*Host` appliers are plain host C++ — they are the oracle path
// (flow_oracle.hpp `setGhostSampled`) and they also BUILD the weights the device consumes, so
// oracle==device parity holds by construction rather than by a second implementation. The
// device mirror (`GhostOverlaySampledDev`, `ghostApplyDeltaSampled`, `ghostDivergDeltaSampled`,
// `GhostGradCsrDev`) lives in the `#ifdef KOKKOS_INLINE_FUNCTION` section at the bottom, so a
// non-Kokkos TU can include this header and get the host half only. Driven from the device
// solver by `AmrFlow::setGhostSampled` (flow.hpp), single-rank AND distributed.
//
// STATUS (2026-08-30, plans §Phase 2 + amr_march_perf_and_distributed_plan.md): Phase-1 rungs 1–4
// and Phase 2 are done — momentum ξ-row seam correction, wall-aware C/F gates, the device mirror,
// the graded refinement policy API (refine.hpp `refineToSdfGraded`), family-free on a seamed mesh,
// and a seam offset converging at ~1.7–1.9 so seams do not set the order. DISTRIBUTED since rungs
// D0–D2: the cloud candidates are a deterministic `probeSlot` set (never a search over the local
// leaf array — that would make a cloud depend on what this rank owns), the builder probes inside
// prepareDistributed's discovery fixpoint, and the clouds read GHOST slots. Two frames meet here
// and coincide only single-rank: geometry is sampled at GLOBAL world points, `probeSlot` takes
// BLOCK-LOCAL fine coords — `frameShift`/`globalFine` carry that, and every default reduces the
// expressions to their single-rank form bit for bit. Anything reached through a probe is an
// EXTENDED slot: use `pres.levelOf`/`loOf`, never `t.level`/`t.bounds`, or a ghost reads past the
// end of the local arrays.
// KNOWN GAPS, in the plan's risk register: sub-face closures are unimplemented (the closed
// sub-faces of a mixed face are Neumann-zero — counted by `nMixedFace`, measured non-vanishing)
// and pocket cells are not excluded from LS clouds.
#ifndef PECLET_AMR_GHOST_PROJECTION_SAMPLED_HPP
#define PECLET_AMR_GHOST_PROJECTION_SAMPLED_HPP
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <memory>
#include <stdexcept>
#include <utility>
#include <vector>
 
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/cf_scheme.hpp"  // CfCsr + detail::{ScalarEnt, compactCsr} (mom seam delta)
#include "peclet/amr/common.hpp"
#include "peclet/amr/ghost_projection.hpp"
#include "peclet/amr/poisson.hpp"
#include "peclet/core/common/host_parallel.hpp"
#include "peclet/core/common/types.hpp"
#include "peclet/core/scheme/ghost_closure.hpp"
 
namespace peclet::amr {
 
struct GhostOverlaySampled {
  GhostOverlay base;
  std::vector<Index> sampStart;  
  std::vector<Index> sampIdx;
  std::vector<double> sampW;
  std::vector<int8_t> sampFluid;  
  std::vector<Index> rowOf;       
  // build census (printed by the builder):
  long nIdentity = 0, nLS2 = 0, nLS1 = 0, nDegraded = 0, nSolidSlot = 0, nMixedFace = 0,
       nSignForced = 0;
};
 
namespace detail {
 
inline bool gpsSolveDense(int n, double* A, double* b) {
  for (int k = 0; k < n; ++k) {
    int piv = k;
    for (int r = k + 1; r < n; ++r)
      if (std::fabs(A[r * n + k]) > std::fabs(A[piv * n + k]))
        piv = r;
    if (std::fabs(A[piv * n + k]) < 1e-14)
      return false;
    if (piv != k) {
      for (int c = k; c < n; ++c) {
        const double tmp = A[k * n + c];
        A[k * n + c] = A[piv * n + c];
        A[piv * n + c] = tmp;
      }
      const double tb = b[k];
      b[k] = b[piv];
      b[piv] = tb;
    }
    for (int r = k + 1; r < n; ++r) {
      const double f = A[r * n + k] / A[k * n + k];
      for (int c = k; c < n; ++c)
        A[r * n + c] -= f * A[k * n + c];
      b[r] -= f * b[k];
    }
  }
  for (int k = n - 1; k >= 0; --k) {
    for (int c = k + 1; c < n; ++c)
      b[k] -= A[k * n + c] * b[c];
    b[k] /= A[k * n + k];
  }
  return true;
}
 
template <unsigned Bits, class WrapFn, class EmitFn>
inline void forEachCoveringSlot(const AmrPoisson<3, Bits>& pres, const long lo[3], const long hi[3],
                                WrapFn&& wrapLocal, EmitFn&& emit, bool descendPending = false) {
  struct Box {
    long lo[3], hi[3];
  };
  constexpr int kMaxDepth = 192;
  Box stack[kMaxDepth];
  int top = 0;
  for (int d = 0; d < 3; ++d) {
    stack[0].lo[d] = lo[d];
    stack[0].hi[d] = hi[d];
  }
  ++top;
  while (top > 0) {
    const Box r = stack[--top];
    const std::array<long, 3> q{r.lo[0], r.lo[1], r.lo[2]};
    const auto [slot, L] = pres.probeSlot(q);
    if (slot < 0) {
      emit(slot);
      if (descendPending) {
        // Discovery: probe an ALIGNED lattice inside the pending region instead of splitting it.
        // Stride = a power of two chosen so the lattice is at most ~9 points per axis, floored at
        // 4: a leaf of size >= stride is stride-aligned (octree) and so is caught this round;
        // finer leaves surface next round as pending sub-regions a stride smaller — the unknown
        // shrinks geometrically, so the fixpoint needs O(log span) CHEAP rounds. Both properties
        // were learnt the hard way on the lmax=6 seam mesh: a recursive blind split's terminal
        // corners inherit each box's own offset and defeat the miss map's dedup (560184 distinct
        // pending coords, 8.5 s round 1), and a FIXED stride of 4 turns a coarse background
        // row's cloud — rho = 2.2·h(L) is domain-sized for a level-6 row — into a lattice over
        // the whole remote domain (263876 coords). Power-of-two alignment makes coords SHARED
        // between overlapping boxes; blocks are root-aligned, so local alignment IS global.
        long ext = 0;
        for (int d = 0; d < 3; ++d)
          ext = std::max(ext, r.hi[d] - r.lo[d] + 1);
        int k = 2;
        while ((1L << (k + 3)) < ext)
          ++k;  // stride*8 >= ext ⇒ <= 9 lattice points per axis
        const long st = 1L << k;
        auto alignUp = [k, st](long v) { return ((v + st - 1) >> k) << k; };
        for (long z = alignUp(r.lo[2]); z <= r.hi[2]; z += st)
          for (long y = alignUp(r.lo[1]); y <= r.hi[1]; y += st)
            for (long x = alignUp(r.lo[0]); x <= r.hi[0]; x += st)
              (void)pres.probeSlot(std::array<long, 3>{x, y, z});
      }
      continue;
    } else {
      const std::array<long, 3> wq = wrapLocal(q);
      const std::array<long, 3> llo = pres.loOf(slot);
      const long s = static_cast<long>(1) << L;
      bool covers = true;
      for (int d = 0; d < 3; ++d) {
        const long u = llo[d] + (r.lo[d] - wq[d]);  // the leaf's lo in the probe's frame
        if (u + s - 1 < r.hi[d]) {
          covers = false;
          break;
        }
      }
      if (covers) {
        emit(slot);
        continue;
      }
    }
    int ax = 0;
    for (int d = 1; d < 3; ++d)
      if (r.hi[d] - r.lo[d] > r.hi[ax] - r.lo[ax])
        ax = d;
    // Coarsest 2^k boundary in (lo, hi]: the highest bit where lo and hi differ.
    unsigned long long x = static_cast<unsigned long long>(r.lo[ax] ^ r.hi[ax]);
    int k = 0;
    while (x >>= 1)
      ++k;
    const long cut = ((r.lo[ax] >> k) + 1) << k;
    if (top + 2 > kMaxDepth)
      throw std::runtime_error("amr::forEachCoveringSlot: descent stack overflow");
    stack[top] = r;
    stack[top].hi[ax] = cut - 1;
    ++top;
    stack[top] = r;
    stack[top].lo[ax] = cut;
    ++top;
  }
}
 
inline void gpsMonomials(const double d[3], int deg, double* m, int& nm) {
  m[0] = 1.0;
  m[1] = d[0];
  m[2] = d[1];
  m[3] = d[2];
  nm = 4;
  if (deg >= 2) {
    m[4] = d[0] * d[0];
    m[5] = d[1] * d[1];
    m[6] = d[2] * d[2];
    m[7] = d[0] * d[1];
    m[8] = d[0] * d[2];
    m[9] = d[1] * d[2];
    nm = 10;
  }
}
 
}  // namespace detail
 
template <unsigned Bits, class SdfFn>
inline auto makeBinaryOpenFnMixed(const BlockOctree<3, Bits>& t, const AmrPoisson<3, Bits>& pres,
                                  SdfFn sdfFn, const Vec<3>& h0, Vec<3> origin,
                                  std::array<long, 3> frameShift = {}) {
  return [&t, &pres, sdfFn, h0, origin, frameShift](const Vec<3>& fc, int axis) -> double {
    (void)t;
    auto centerSample = [&](int side) -> std::pair<bool, float> {
      Vec<3> probe = fc;
      probe[axis] += side * 0.25 * h0[axis];  // strictly inside the adjacent leaf (per axis)
      std::array<long, 3> q{};
      for (int d = 0; d < 3; ++d)
        q[d] = static_cast<long>(std::floor((probe[d] - origin[d]) / h0[d])) - frameShift[d];
      const Index j = pres.probeSlot(q).first;
      if (j < 0)
        return {false, -1.0f};
      const auto lo = pres.loOf(j);
      const double s = static_cast<double>(Index(1) << pres.levelOf(j));
      Vec<3> c{};
      for (int d = 0; d < 3; ++d)
        c[d] = origin[d] + (static_cast<double>(lo[d] + frameShift[d]) + 0.5 * s) * h0[d];
      const double sd = sdfFn(c);
      return {sd > 0.0, static_cast<float>(sd)};
    };
    const auto [flM, sM] = centerSample(-1);
    const auto [flP, sP] = centerSample(+1);
    return (flM && flP && 0.5f * (sM + sP) >= 0.0f) ? 1.0 : 0.0;
  };
}
 
template <unsigned Bits, class SdfFn>
inline GhostOverlaySampled buildGhostOverlaySampled(
    const BlockOctree<3, Bits>& t, const AmrPoisson<3, Bits>& pres, SdfFn&& sdf, int matrixOrder,
    int rhsOrder, Vec<3> origin = Vec<3>{}, const std::array<long, 3>* globalFine = nullptr,
    std::array<long, 3> frameShift = {}, bool discovery = false, double rhoFactor = 2.2,
    long maxSamples = 0) {
  GhostOverlaySampled ov;
  const Index n = t.numLeaves();
  // Phase 3 (`docs/amr_anisotropic.md` §5): the ROOT spacing per axis. `pres.h0()` already is
  // it; recovering it from `cellWidth(0)` (as this line did) only ever worked because the cells
  // were cubes.
  const Vec<3> h0 = pres.h0();
  const auto fe = pres.fineExt();
  const std::array<long, 3> shiftG = frameShift;
  const std::array<long, 3> gfine =
      globalFine ? *globalFine
                 : std::array<long, 3>{static_cast<long>(fe[0]), static_cast<long>(fe[1]),
                                       static_cast<long>(fe[2])};
  // World position -> block-local fine coord (the probe frame). Single-rank the shift is 0.
  auto probeCoord = [&](const Vec<3>& p, long q[3]) {
    for (int d = 0; d < 3; ++d)
      q[d] = static_cast<long>(std::floor((p[d] - origin[d]) / h0[d])) - shiftG[d];
  };
 
  // Leaf centers + fluid flags. Host-parallel (rung 4): per-leaf disjoint writes.
  std::vector<Vec<3>> cen(static_cast<std::size_t>(n));
  std::vector<char> fluid(static_cast<std::size_t>(n));
  hostParFor(n, [&](Index i) {
    auto b = t.bounds(i);
    const double s = static_cast<double>(Index(1) << t.level(i));
    Vec<3> c{};
    for (int d = 0; d < 3; ++d)
      c[d] = origin[d] +
             (static_cast<double>(b[0][d]) + static_cast<double>(shiftG[d]) + 0.5 * s) * h0[d];
    cen[static_cast<std::size_t>(i)] = c;
    fluid[static_cast<std::size_t>(i)] = sdf(c) > 0.0 ? 1 : 0;
  });
 
  // D2: the same two per-slot facts for the GHOST tail [n, n+nGhost). A ghost's geometry is
  // recomputed here rather than exchanged: `loOf` is the covering leaf's global anchor and the
  // level is the owner's, so `sdf(centre)` reproduces the owner's own value bit for bit (the same
  // argument buildOpenness's ghost rows rest on). Empty single-rank — there are no ghosts.
  const Index ng = pres.numGhosts();
  std::vector<Vec<3>> cenG(static_cast<std::size_t>(ng));
  std::vector<char> fluidG(static_cast<std::size_t>(ng));
  for (Index g = 0; g < ng; ++g) {
    const auto lo = pres.loOf(n + g);
    const double s = static_cast<double>(Index(1) << pres.levelOf(n + g));
    Vec<3> c{};
    for (int d = 0; d < 3; ++d)
      c[d] = origin[d] + (static_cast<double>(lo[d] + shiftG[d]) + 0.5 * s) * h0[d];
    cenG[static_cast<std::size_t>(g)] = c;
    fluidG[static_cast<std::size_t>(g)] = sdf(c) > 0.0 ? 1 : 0;
  }
  auto centerOf = [&](Index j) -> const Vec<3>& {
    return j < n ? cen[static_cast<std::size_t>(j)] : cenG[static_cast<std::size_t>(j - n)];
  };
  auto fluidOf = [&](Index j) -> char {
    return j < n ? fluid[static_cast<std::size_t>(j)] : fluidG[static_cast<std::size_t>(j - n)];
  };
 
  // ---- DD1 (docs/amr_march_perf_and_distributed_plan.md): the LS cloud candidates are a
  // DETERMINISTIC PROBE SET, not a search over the local leaf array.
  //
  // The set a cloud must contain is "every fluid leaf whose CENTER lies within rho of the sample
  // position" — a pure function of (octree, geometry). The old hash-bin search enumerated it by
  // walking THIS RANK's leaf array, which would make the cloud decomposition-DEPENDENT near a
  // block boundary (D6 forbids that). It is now enumerated by `detail::forEachCoveringSlot`, a
  // coordinate-space descent through `pres.probeSlot` — the route every other builder already uses
  // to reach across ranks.
  //
  // ORDER is contractual: it fixes the floating-point value of the normal-matrix sums and hence
  // the weights. The enumeration is therefore re-sorted into the OLD bin traversal order — bins of
  // 4*h0 in (bx, by, bz) loop order, and within a bin ascending Morton code of the leaf's global
  // lo corner. That is bit-for-bit the old order (the leaf array is Z-order sorted, so ascending
  // code == ascending leaf index) and, unlike the leaf index, it is the same under any
  // decomposition.
  //
  // The bins survive DD1 only as that order key, but `domain` below is also the minimum-image
  // PERIOD of the cloud metric, so it has to be the true one: F2 (§Findings, fixed on main before
  // this rung) derived it from the max leaf bound, which is INCLUSIVE, and the truncating /4 then
  // lost a whole bin — 4 fine cells short at every power-of-two domain. It is the octree's fine
  // extent, the period probeSlot and every other wrap already use. Single-rank the block IS the
  // domain; distributed (D1) it is `gfine`, the GLOBAL fine extent — a per-rank period would be
  // F2 all over again, in a decomposition-dependent coat.
  // AM4 (`docs/amr_anisotropic.md` §5): the ORDER-KEY bins and the minimum-image PERIOD both
  // become per-axis. `hb[d] = 4*h0[d]` keeps a bin four fine cells wide on every axis (so the
  // traversal order is the old one on a cubic grid, bit for bit), and the period is the true
  // world extent of each axis. NOTE the old single `domain` was the LONGEST axis's period, which
  // is wrong for the shorter axes of a NON-CUBIC brick whenever a cloud straddles the periodic
  // boundary there — pre-announced in the design note as the one place a bit can legitimately
  // move, and it moves toward the correct value.
  double hb[3], domain[3];
  long nb[3];
  for (int d = 0; d < 3; ++d) {
    hb[d] = 4.0 * h0[d];
    nb[d] = std::max<long>(1, gfine[d] / 4);
    domain[d] = static_cast<double>(gfine[d]) * h0[d];
  }
  const long nbx = std::max(std::max(nb[0], nb[1]), nb[2]);  // bin-index stride (order key only)
  auto wrapLocal = [gfine, shiftG](const std::array<long, 3>& p) {
    std::array<long, 3> q{};
    for (int d = 0; d < 3; ++d) {
      const long g = p[d] + shiftG[d];
      q[d] = ((g % gfine[d]) + gfine[d]) % gfine[d] - shiftG[d];
    }
    return q;
  };
  using Oct = BlockOctree<3, Bits>;
  using Code = typename Oct::Code;
  // The order key: the Morton code of the slot's GLOBAL lo corner. `loOf` serves local leaves and
  // ghosts alike (a ghost's lo IS its owner's, in this block's frame), so the key of a leaf is the
  // same number on every rank that can see it — which is what makes the cloud's accumulation order
  // decomposition-independent (D6) rather than merely deterministic.
  auto keyOf = [&](Index j) {
    const auto lo = pres.loOf(j);
    std::array<typename Oct::Coord, 3> q{};
    for (int d = 0; d < 3; ++d)
      q[d] = static_cast<typename Oct::Coord>(lo[d] + shiftG[d]);
    return Oct::M::encode(q).code();
  };
  struct Cand {
    long bin;
    Code key;
    Index slot;
    double d2;
  };
 
  // ---- M2a cloud-economy knobs (both default to the shipped behaviour exactly) -----------------
  // A degree-2 least squares needs 12 points; the shipped radius gathers 95-162 (M1). The two
  // arguments let the cloud-economy study vary the sampling without touching production defaults:
  //   rhoFactor   radius factor, default 2.2  (rho = factor * max(h, H))
  //   maxSamples  keep only the N nearest candidates, default 0 = no cap
  // (AmrFlow::setGhostSampled(on, rho, maxSamples) / Python set_ghost_sampled(on, rho=,
  // max_samples=); until 2026-09-10 these were the PECLET_CORE_GPS_RHO / _MAXN environment
  // variables — QUALITY_PLAN D3 makes every result-changing knob an explicit argument.)
  // The cap selects by (distance^2, then global Morton key) — a total order that is a pure
  // function of geometry, so the kept SET is decomposition-independent like everything else here;
  // the kept candidates are then emitted in the unchanged canonical (bin, Morton) order, so the
  // accumulation discipline is untouched. Capping changes the WEIGHTS, which is exactly what the
  // study measures — it is never on by default.
  const double gpsRhoFactor = rhoFactor > 0.0 ? rhoFactor : 2.2;
  const long gpsMaxN = maxSamples > 0 ? maxSamples : 0;
 
  // LS functional at world position p, degree deg, radius rho, scale H: returns the (idx, w)
  // list. Weight vector w_j = mono(d_j) . M^{-1} e0 with M the normal matrix.
  // AM3 (`docs/amr_anisotropic.md` §5): the cloud metric is the INDEX-SPACE one — the ball
  // `sum_d (del_d/rho_d)^2 <= 1` with `rho_d` the radius factor times the coarser of the two
  // widths ON AXIS d, and the monomials built on `del_d/H_d`. The degree-2 polynomial space is
  // closed under a per-axis scaling, so this is the SAME least-squares functional as the cubic
  // one up to round-off, and bitwise identical when the three spacings agree.
  auto lsFunctional = [&](const Vec<3>& p, const double rho[3], const double H[3], int deg,
                          std::vector<Index>& idx, std::vector<double>& w) -> bool {
    idx.clear();
    w.clear();
    // (1) every leaf overlapping the fine-cell box that contains the ball (a leaf whose center is
    //     in the ball contains that center, so it overlaps the box).
    long qlo[3], qhi[3];
    for (int d = 0; d < 3; ++d) {
      qlo[d] = static_cast<long>(std::floor((p[d] - origin[d] - rho[d]) / h0[d])) - shiftG[d];
      qhi[d] = static_cast<long>(std::floor((p[d] - origin[d] + rho[d]) / h0[d])) - shiftG[d];
    }
    std::vector<Index> slots;
    detail::forEachCoveringSlot(
        pres, qlo, qhi, wrapLocal,
        [&](Index s) {
          if (!discovery)
            slots.push_back(s);
        },
        /*descendPending=*/discovery);
    if (discovery)
      return false;  // probe-only: the misses are registered, the weights are never consumed
    std::sort(slots.begin(), slots.end());
    slots.erase(std::unique(slots.begin(), slots.end()), slots.end());
    // (2) keep the fluid ones whose center is inside the ball (the old predicate, verbatim).
    std::vector<Cand> cand;
    cand.reserve(slots.size());
    for (Index j : slots) {
      if (j < 0)
        continue;  // kPending: a discovery round; this pass's result is discarded anyway
      if (!fluidOf(j))
        continue;
      double r2 = 0, q2 = 0;
      for (int d = 0; d < 3; ++d) {
        double del = centerOf(j)[d] - p[d];
        if (del > 0.5 * domain[d])
          del -= domain[d];
        if (del < -0.5 * domain[d])
          del += domain[d];
        r2 += del * del;  // the (distance^2) tie-break key of the nearest-N cap
        const double u = del / rho[d];
        q2 += u * u;  // the axis-scaled ball predicate
      }
      if (q2 > 1.0)
        continue;
      const Vec<3>& cj = centerOf(j);
      const long bx = static_cast<long>((cj[0] - origin[0]) / hb[0]) % nbx;
      const long by = static_cast<long>((cj[1] - origin[1]) / hb[1]) % nbx;
      const long bz = static_cast<long>((cj[2] - origin[2]) / hb[2]) % nbx;
      cand.push_back(Cand{(bz * nbx + by) * nbx + bx, keyOf(j), j, r2});
    }
    if (gpsMaxN > 0 && static_cast<long>(cand.size()) > gpsMaxN) {
      // nearest-N by (distance^2, global Morton key): a strict, geometry-only total order.
      std::nth_element(
          cand.begin(), cand.begin() + gpsMaxN, cand.end(),
          [](const Cand& a, const Cand& b) { return a.d2 != b.d2 ? a.d2 < b.d2 : a.key < b.key; });
      cand.resize(static_cast<std::size_t>(gpsMaxN));
    }
    std::sort(cand.begin(), cand.end(), [](const Cand& a, const Cand& b) {
      return a.bin != b.bin ? a.bin < b.bin : a.key < b.key;
    });
    // (3) emit in the canonical bin traversal order.
    std::vector<Index> pts;
    pts.reserve(cand.size());
    long lo[3], hi[3];
    for (int d = 0; d < 3; ++d) {
      lo[d] = static_cast<long>(std::floor((p[d] - origin[d] - rho[d]) / hb[d]));
      hi[d] = static_cast<long>(std::floor((p[d] - origin[d] + rho[d]) / hb[d]));
    }
    for (long bx = lo[0]; bx <= hi[0]; ++bx)
      for (long by = lo[1]; by <= hi[1]; ++by)
        for (long bz = lo[2]; bz <= hi[2]; ++bz) {
          const long wx = (bx % nbx + nbx) % nbx, wy = (by % nbx + nbx) % nbx,
                     wz = (bz % nbx + nbx) % nbx;
          const long bl = (wz * nbx + wy) * nbx + wx;
          auto it = std::lower_bound(cand.begin(), cand.end(), bl,
                                     [](const Cand& c, long v) { return c.bin < v; });
          for (; it != cand.end() && it->bin == bl; ++it)
            pts.push_back(it->slot);
        }
    const int need = deg >= 2 ? 12 : 5;
    if (static_cast<long>(pts.size()) < need)
      return false;
    double A[100] = {}, e0[10] = {};
    int nm = 0;
    double mono[10];
    for (Index j : pts) {
      double d[3];
      for (int dd = 0; dd < 3; ++dd) {
        double del = centerOf(j)[dd] - p[dd];
        if (del > 0.5 * domain[dd])
          del -= domain[dd];
        if (del < -0.5 * domain[dd])
          del += domain[dd];
        d[dd] = del / H[dd];
      }
      detail::gpsMonomials(d, deg, mono, nm);
      for (int r = 0; r < nm; ++r)
        for (int c = 0; c < nm; ++c)
          A[r * nm + c] += mono[r] * mono[c];
    }
    e0[0] = 1.0;
    if (!detail::gpsSolveDense(nm, A, e0))
      return false;
    for (Index j : pts) {
      double d[3];
      for (int dd = 0; dd < 3; ++dd) {
        double del = centerOf(j)[dd] - p[dd];
        if (del > 0.5 * domain[dd])
          del -= domain[dd];
        if (del < -0.5 * domain[dd])
          del += domain[dd];
        d[dd] = del / H[dd];
      }
      detail::gpsMonomials(d, deg, mono, nm);
      double wj = 0;
      for (int k = 0; k < nm; ++k)
        wj += mono[k] * e0[k];
      idx.push_back(j);
      w.push_back(wj);
    }
    return true;
  };
 
  ov.rowOf.assign(static_cast<std::size_t>(n), -1);
  ov.sampStart.assign(1, 0);
  std::vector<Index> sIdx;
  std::vector<double> sW;
 
  // ---- Rung 4 (setup-parallel plan, D2's third pattern): the single sequential-append pass is
  // split in two. PASS 1 (below, host-parallel) does ALL the geometry per leaf and stages an
  // accepted row into that leaf's OWN heap-allocated RowStage — nothing is appended to `ov`, so
  // there is no `slot = ov.base.n` cursor to serialize on. PASS 2 (serial, after) walks the
  // leaves IN ORDER and appends the staged rows, which reproduces today's row order and CSR slot
  // order EXACTLY: row order = leaf order, and within a row the 15 chain slots in (a, q) order
  // with each slot's entries in the order lsFunctional produced them. Nothing is classified
  // twice — the plan permits it (row state is a pure function of geometry) but staging is
  // cheaper.
  //
  // Only ~13% of leaves become rows, so stages are allocated per accepted row (a null pointer for
  // every clean/solid leaf) rather than as a dense array.
  //
  // NOTE the cloud's candidate ORDER sets the order the LS normal-matrix sums are accumulated,
  // i.e. the floating-point value of the weights; DD1's canonical (bin, Morton) sort above pins it
  // independently of both the thread schedule and the decomposition.
  struct RowStage {
    Index cell = 0;
    float rescale = 0.0f;
    int8_t coupled = 0;
    int8_t state[6] = {};
    float th[6] = {}, w_bc[6] = {}, w_n1[6] = {}, w_n2[6] = {}, wm_n1[6] = {}, wm_n2[6] = {};
    double invh[3] = {0.0, 0.0, 0.0};  
    int8_t sampFluid[15] = {};
    Index sampCnt[15] = {};
    std::vector<Index> sIdx;
    std::vector<double> sW;
    long nIdentity = 0, nLS2 = 0, nLS1 = 0, nDegraded = 0, nSolidSlot = 0;
  };
  std::vector<std::unique_ptr<RowStage>> stages(static_cast<std::size_t>(n));
  // Per-leaf census slots for the two counters that are tallied over EVERY fluid leaf (clean rows
  // included), summed in leaf order in pass 2 — integers, so the total is order-independent.
  std::vector<uint8_t> sfPer(static_cast<std::size_t>(n), 0), mfPer(static_cast<std::size_t>(n), 0);
 
  hostParFor(n, [&](Index i) {
    if (!fluid[static_cast<std::size_t>(i)])
      return;
    const unsigned Li = t.level(i);
    const double h = pres.cellWidth(i);
    const Vec<3>& c = cen[static_cast<std::size_t>(i)];
 
    // Virtual uniform-position float samples (theta / closure magnitudes — the D2 route).
    // AM2 (`docs/amr_anisotropic.md` §5): a virtual sample sits `q` cells out along axis `a`, so
    // it steps by THAT axis's cell width. Identical to the single `h` on a cubic octree.
    float Cq[3][5], F[3][4];
    for (int a = 0; a < 3; ++a) {
      const double ha = pres.cellWidth(i, a);
      for (int q = -2; q <= 2; ++q) {
        Vec<3> p = c;
        p[a] += static_cast<double>(q) * ha;
        Cq[a][q + 2] = static_cast<float>(sdf(p));
      }
      for (int m = 0; m < 4; ++m)
        F[a][m] = 0.5f * (Cq[a][m] + Cq[a][m + 1]);
    }
    // Canonical own-face openness (any-sub-face-open across finer neighbours): force the SIGN
    // of the own-face F entries to the canonical rule so overlay-closed <=> binary-closed.
    bool anyOpen[6] = {};  // slot k = 2*axis + (0 plus, 1 minus)
    int nSub[6] = {};
    const float si0 = static_cast<float>(sdf(c));
    pres.forEachFaceFull(i, [&](Index j, int axis, int dir, double, double, double) {
      const int face = 2 * axis + (dir > 0 ? 0 : 1);
      ++nSub[face];
      if (j >= 0 && fluidOf(j)) {
        const float sj = static_cast<float>(sdf(centerOf(j)));
        if (0.5f * (si0 + sj) >= 0.0f)
          anyOpen[face] = true;
      }
    });
    for (int k = 0; k < 6; ++k) {
      const int a = k / 2;
      const int m = (k & 1) ? 1 : 2;  // minus own face = F[a][1], plus own face = F[a][2]
      const bool virtOpen = F[a][m] >= 0.0f && Cq[a][(k & 1) ? 1 : 3] > 0.0f && Cq[a][2] > 0.0f;
      if (virtOpen != anyOpen[k]) {
        // Marginal face where the virtual and canonical classifications disagree: force the
        // canonical sign (overlay-closed <=> binary-closed), keep the (small) magnitude.
        ++sfPer[static_cast<std::size_t>(i)];
        const float mag = std::fabs(F[a][m]) > 0.0f ? std::fabs(F[a][m]) : 1e-6f;
        F[a][m] = anyOpen[k] ? mag : -mag;
      }
      if (nSub[k] > 1 && !anyOpen[k])
        ++mfPer[static_cast<std::size_t>(i)];  // coarse face fully closed against finer neighbours
    }
    // Non-clean test (canonical own faces + virtual +/-1 centers), as the classic builder.
    bool clean = true;
    for (int a = 0; a < 3; ++a)
      clean = clean && Cq[a][1] >= 0.0f && Cq[a][3] >= 0.0f && F[a][1] >= 0.0f && F[a][2] >= 0.0f;
    if (clean)
      return;
 
    // Fill the row into a scratch ONE-row overlay (gpFillRow writes ov.field(slot), so it needs a
    // GhostOverlayRef; slot 0 of a private overlay gives the identical arithmetic).
    GhostOverlay one;
    one.cell.resize(1);
    one.rescale.resize(1);
    one.coupled.resize(1);
    one.state.resize(6);
    one.th.resize(6);
    one.w_bc.resize(6);
    one.w_n1.resize(6);
    one.w_n2.resize(6);
    one.wm_n1.resize(6);
    one.wm_n2.resize(6);
    detail::GhostOverlayRef ref{one};
    if (!scheme::gpFillRow(ref, 0, i, F, Cq, matrixOrder, rhsOrder))
      return;  // rejected row: nothing staged (the old code's resize(slot) rollback)
 
    auto st = std::make_unique<RowStage>();
    st->cell = one.cell[0];
    st->rescale = one.rescale[0];
    st->coupled = one.coupled[0];
    for (int k = 0; k < 6; ++k) {
      st->state[k] = one.state[static_cast<std::size_t>(k)];
      st->th[k] = one.th[static_cast<std::size_t>(k)];
      st->w_bc[k] = one.w_bc[static_cast<std::size_t>(k)];
      st->w_n1[k] = one.w_n1[static_cast<std::size_t>(k)];
      st->w_n2[k] = one.w_n2[static_cast<std::size_t>(k)];
      st->wm_n1[k] = one.wm_n1[static_cast<std::size_t>(k)];
      st->wm_n2[k] = one.wm_n2[static_cast<std::size_t>(k)];
    }
    for (int a = 0; a < 3; ++a)  // Phase 3: per axis
      st->invh[a] = 1.0 / pres.cellWidth(i, a);
 
    // Sample functionals for the 15 chain slots.
    for (int a = 0; a < 3; ++a)
      for (int q = -2; q <= 2; ++q) {
        const int sl = a * 5 + (q + 2);
        const Index nBefore = static_cast<Index>(st->sIdx.size());
        Vec<3> p = c;
        p[a] += static_cast<double>(q) * pres.cellWidth(i, a);  // AM2: per-axis step
        // Recover the covering leaf: floor the world position to fine units (a level-L cell
        // center is lo + 0.5*2^L in fine units, so the floor lands inside the cell), then into
        // the block-local probe frame.
        long lcq[3];
        probeCoord(p, lcq);
        const std::array<long, 3> lc{lcq[0], lcq[1], lcq[2]};
        const Index j = pres.probeSlot(lc).first;
        const bool fluidP = sdf(p) > 0.0;
        st->sampFluid[sl] = fluidP ? 1 : 0;
        // `pres.levelOf`, never `t.level`: D2 made `j` an EXTENDED slot, and a ghost index
        // read out of the local level array — silently reporting a wrong level, which sent
        // same-level covers down the least-squares branch instead of the identity one.
        if (q == 0 || (j >= 0 && pres.levelOf(j) == Li)) {
          // own cell or same-level cover: identity (the classic chain; bit-identical band).
          if (j >= 0 && fluidOf(j) == (fluidP ? 1 : 0)) {
            st->sIdx.push_back(j);
            st->sW.push_back(1.0);
            ++st->nIdentity;
          } else if (j >= 0) {
            st->sIdx.push_back(j);  // marginal float disagreement: still the covering value
            st->sW.push_back(1.0);
            ++st->nIdentity;
          }
          st->sampCnt[sl] = static_cast<Index>(st->sIdx.size()) - nBefore;
          continue;
        }
        if (!fluidP) {
          // solid virtual position: empty functional (reads the masked 0).
          ++st->nSolidSlot;
          st->sampCnt[sl] = 0;
          continue;
        }
        // AM3: the scale and the cloud radius are PER AXIS — the coarser of the two widths on
        // each axis. On a cubic octree all three are the same double, so the cloud, its
        // accumulation order and its weights are the pre-Phase-3 ones bit for bit.
        double H[3], rho[3];
        for (int d = 0; d < 3; ++d) {
          const double hd = pres.cellWidth(i, d);
          H[d] = (j >= 0) ? pres.cellWidth(j, d) : hd;
          rho[d] = gpsRhoFactor * std::max(hd, H[d]);
        }
        std::vector<Index> idx;
        std::vector<double> w;
        if (discovery) {
          // Probe-only: one descent registers the box's misses; the degree-1 retry would descend
          // the SAME box, so it is skipped along with everything below that reads the weights.
          (void)lsFunctional(p, rho, H, 2, idx, w);
          continue;
        }
        if (lsFunctional(p, rho, H, 2, idx, w)) {
          ++st->nLS2;
        } else if (lsFunctional(p, rho, H, 1, idx, w)) {
          ++st->nLS1;
        } else if (j >= 0 && fluidOf(j)) {
          idx.assign(1, j);
          w.assign(1, 1.0);
          ++st->nDegraded;
        } else {
          ++st->nDegraded;  // no support at all: reads 0 (fluid-only cascade floor)
        }
        for (std::size_t k = 0; k < idx.size(); ++k) {
          st->sIdx.push_back(idx[k]);
          st->sW.push_back(w[k]);
        }
        st->sampCnt[sl] = static_cast<Index>(st->sIdx.size()) - nBefore;
      }
    stages[static_cast<std::size_t>(i)] = std::move(st);
  });
 
  if (discovery) {
    // A discovery round consumes nothing staged: the misses are in the resolver, which is the
    // whole point. Return the empty overlay without pass 2 or the census line.
    return ov;
  }
 
  // PASS 2 (serial, order-preserving): append the staged rows in LEAF order — the row order, the
  // 15-slot order and the sampStart/sampIdx/sampW concatenation are exactly what the old single
  // pass produced.
  for (Index i = 0; i < n; ++i) {
    ov.nSignForced += sfPer[static_cast<std::size_t>(i)];
    ov.nMixedFace += mfPer[static_cast<std::size_t>(i)];
    const RowStage* st = stages[static_cast<std::size_t>(i)].get();
    if (!st)
      continue;
    ov.base.cell.push_back(st->cell);
    ov.base.rescale.push_back(st->rescale);
    ov.base.coupled.push_back(st->coupled);
    for (int k = 0; k < 6; ++k) {
      ov.base.state.push_back(st->state[k]);
      ov.base.th.push_back(st->th[k]);
      ov.base.w_bc.push_back(st->w_bc[k]);
      ov.base.w_n1.push_back(st->w_n1[k]);
      ov.base.w_n2.push_back(st->w_n2[k]);
      ov.base.wm_n1.push_back(st->wm_n1[k]);
      ov.base.wm_n2.push_back(st->wm_n2[k]);
    }
    Index acc = static_cast<Index>(sIdx.size());
    for (int sl = 0; sl < 15; ++sl) {
      ov.sampFluid.push_back(st->sampFluid[sl]);
      acc += st->sampCnt[sl];
      ov.sampStart.push_back(acc);
    }
    sIdx.insert(sIdx.end(), st->sIdx.begin(), st->sIdx.end());
    sW.insert(sW.end(), st->sW.begin(), st->sW.end());
    ov.nIdentity += st->nIdentity;
    ov.nLS2 += st->nLS2;
    ov.nLS1 += st->nLS1;
    ov.nDegraded += st->nDegraded;
    ov.nSolidSlot += st->nSolidSlot;
    for (int a = 0; a < 3; ++a)
      ov.base.invh.push_back(st->invh[a]);
    ov.rowOf[static_cast<std::size_t>(i)] = ov.base.n;
    ++ov.base.n;
  }
  ov.sampIdx = std::move(sIdx);
  ov.sampW = std::move(sW);
  std::fprintf(stderr,
               "[peclet.amr] sampled ghost overlay: %lld rows | slots identity %ld, LS2 %ld, "
               "LS1 %ld, degraded %ld, solid %ld | csr %lld | sign-forced faces %ld, closed mixed "
               "faces %ld\n",
               static_cast<long long>(ov.base.n), ov.nIdentity, ov.nLS2, ov.nLS1, ov.nDegraded,
               ov.nSolidSlot, static_cast<long long>(ov.sampIdx.size()), ov.nSignForced,
               ov.nMixedFace);
  return ov;
}
 
inline double gpsSample(const GhostOverlaySampled& ov, Index r, int a, int q,
                        const std::vector<double>& x) {
  const std::size_t s = static_cast<std::size_t>(r * 15 + a * 5 + (q + 2));
  double v = 0.0;
  for (Index k = ov.sampStart[s]; k < ov.sampStart[s + 1]; ++k)
    v += ov.sampW[static_cast<std::size_t>(k)] *
         x[static_cast<std::size_t>(ov.sampIdx[static_cast<std::size_t>(k)])];
  return v;
}
 
inline void ghostApplyDeltaSampledHost(const GhostOverlaySampled& ov, const std::vector<double>& x,
                                       std::vector<double>& y) {
  const GhostOverlay& g = ov.base;
  for (Index r = 0; r < g.n; ++r) {
    const std::size_t rr = static_cast<std::size_t>(r);
    const Index c = g.cell[rr];
    if (!g.coupled[rr]) {
      y[static_cast<std::size_t>(c)] = 0.0;
      continue;
    }
    auto X = [&](int a, int q) { return gpsSample(ov, r, a, q, x); };
    // Phase 3: the ratio form (see ghostApplyDeltaHost) — exactly 1.0 on a cubic octree.
    const double ih = g.invh[rr * 3];
    const double rx = g.invh[rr * 3 + 1] / ih, ry = g.invh[rr * 3 + 2] / ih;
    const double r2[3] = {1.0, rx * rx, ry * ry};
    double delta = 0.0;
    for (int k = 0; k < 6; ++k) {
      const int8_t st = g.state[rr * 6 + static_cast<std::size_t>(k)];
      if (st != scheme::GP_QUAD && st != scheme::GP_LIN)
        continue;
      const int a = k / 2;
      const int sgn = (k & 1) ? -1 : 1;
      const int mn = (k & 1) ? 1 : 0;
      const int mf = (k & 1) ? 2 : -1;
      const double w1 = g.wm_n1[rr * 6 + static_cast<std::size_t>(k)];
      const double w2 = g.wm_n2[rr * 6 + static_cast<std::size_t>(k)];
      delta += r2[a] * (sgn * w1 * (X(a, mn) - X(a, mn - 1)));
      if (st == scheme::GP_QUAD && w2 != 0.0)
        delta += r2[a] * (sgn * w2 * (X(a, mf) - X(a, mf - 1)));
    }
    y[static_cast<std::size_t>(c)] =
        g.rescale[rr] * (y[static_cast<std::size_t>(c)] + ih * ih * delta);
  }
}
 
inline void ghostDivergDeltaSampledHost(const GhostOverlaySampled& ov,
                                        const std::array<std::vector<double>, 3>& u,
                                        std::vector<double>& d) {
  const GhostOverlay& g = ov.base;
  for (Index r = 0; r < g.n; ++r) {
    const std::size_t rr = static_cast<std::size_t>(r);
    const Index c = g.cell[rr];
    if (!g.coupled[rr]) {
      d[static_cast<std::size_t>(c)] = 0.0;
      continue;
    }
    auto U = [&](int a, int m) {  // face-averaged value at face index i+m along axis a
      return 0.5 * (gpsSample(ov, r, a, m - 1, u[static_cast<std::size_t>(a)]) +
                    gpsSample(ov, r, a, m, u[static_cast<std::size_t>(a)]));
    };
    // Phase 3: the ratio form (see ghostDivergDeltaHost).
    const double ihd = g.invh[rr * 3];
    const double r1[3] = {1.0, g.invh[rr * 3 + 1] / ihd, g.invh[rr * 3 + 2] / ihd};
    double dd = 0.0;
    for (int k = 0; k < 6; ++k) {
      const int8_t st = g.state[rr * 6 + static_cast<std::size_t>(k)];
      if (st == scheme::GP_COUPLED)
        continue;
      const int a = k / 2;
      const int sgn = (k & 1) ? -1 : 1;
      const int mg = (k & 1) ? 0 : 1;
      const int mn = (k & 1) ? 1 : 0;
      const int mf = (k & 1) ? 2 : -1;
      if (st == scheme::GP_EXPLICIT) {
        dd += r1[a] * (sgn * U(a, mg));
        continue;
      }
      if (st == scheme::GP_BC_ONLY)
        continue;
      double val = g.w_n1[rr * 6 + static_cast<std::size_t>(k)] * U(a, mn);
      if (st == scheme::GP_QUAD)
        val += g.w_n2[rr * 6 + static_cast<std::size_t>(k)] * U(a, mf);
      dd += r1[a] * (sgn * val);
    }
    d[static_cast<std::size_t>(c)] = g.rescale[rr] * (d[static_cast<std::size_t>(c)] + ihd * dd);
  }
}
 
template <unsigned Bits, class Mom, class SdfFn>
inline CfCsr buildMomSeamDelta(const GhostOverlaySampled& ov, const BlockOctree<3, Bits>& t,
                               const AmrPoisson<3, Bits>& pres, const Mom& mom, SdfFn&& sdf,
                               Vec<3> origin, double idiag, double mu, long* nSkipped = nullptr,
                               std::array<long, 3> frameShift = {}) {
  const Index n = t.numLeaves();
  const double h0 = pres.cellWidth(0) / static_cast<double>(Index(1) << t.level(0));
  std::vector<std::vector<detail::ScalarEnt>> per(static_cast<std::size_t>(n));
  long skipped = 0;
  for (Index i = 0; i < n; ++i) {
    const Index r = ov.rowOf[static_cast<std::size_t>(i)];
    if (r < 0)
      continue;
    const unsigned Li = t.level(i);
    bool seam = false;
    for (int k = 0; k < 6; ++k) {
      const Index nb = mom.neighborOf(i, k);
      if (nb < 0 || pres.levelOf(nb) != Li)  // extended slot: ghost-safe (D2)
        seam = true;
    }
    if (!seam)
      continue;
    // Virtual classification + row-local ξ stencil.
    auto b = t.bounds(i);
    const double s = static_cast<double>(Index(1) << Li);
    Vec<3> c{};
    for (int d = 0; d < 3; ++d)
      c[d] = origin[d] +
             (static_cast<double>(b[0][d]) + static_cast<double>(frameShift[d]) + 0.5 * s) * h0;
    const double h = pres.cellWidth(i);
    const double sdfC = sdf(c);
    double sdfNv[6];
    bool anyGhost = false;
    for (int k = 0; k < 6; ++k) {
      Vec<3> p = c;
      p[k / 2] += ((k % 2 == 0) ? +1.0 : -1.0) * h;  // k even = plus side (AmrCutCell::neighbor)
      sdfNv[k] = sdf(p);
      if (sdfNv[k] < 0.0)
        anyGhost = true;
    }
    if (!mom.isCut(i) && !anyGhost)
      continue;  // regular both ways: the C/F-aware raw row is already the target
    if (!mom.isCut(i) && anyGhost) {
      ++skipped;  // raw-regular but virtually ghost (marginal): out of rung-1 scope
      continue;
    }
    auto push = [&](Index cell, double w) {
      if (w != 0.0)
        per[static_cast<std::size_t>(i)].push_back(detail::ScalarEnt{cell, w});
    };
    // + raw row / rscale_r (the assembled ξ row verbatim).
    const double rsr = mom.rhsScale(i);
    push(i, mom.acRaw()[static_cast<std::size_t>(i)] / rsr);
    for (int k = 0; k < 6; ++k) {
      const double a = mom.offRaw()[static_cast<std::size_t>(i) * 6 + k];
      const Index nb = mom.neighborOf(i, k);
      if (a != 0.0 && nb >= 0)
        push(nb, a / rsr);
    }
    // − virtual row / rscale_v.
    if (anyGhost) {
      // Phase 3: the ROW-LOCAL beta per axis (`docs/amr_anisotropic.md` §3). `h` is the row's
      // cellWidth on axis 0; the other two come from the same leaf.
      double beta[3];
      for (int d = 0; d < 3; ++d) {
        const double hd = pres.cellWidth(i, d);
        beta[d] = mu / (hd * hd);
      }
      const double AC0v = idiag + ((2.0 * beta[0] + 2.0 * beta[1]) + 2.0 * beta[2]);
      double ACv, offv[6], rsv = 1.0, inhomv = 0.0;
      Mom::buildCutStencil(sdfC, sdfNv, beta, AC0v, ACv, offv, rsv, inhomv);
      push(i, -ACv / rsv);
      for (int k = 0; k < 6; ++k) {
        if (offv[k] == 0.0)
          continue;
        const int q = (k % 2 == 0) ? +1 : -1;
        const std::size_t sl = static_cast<std::size_t>(r * 15 + (k / 2) * 5 + (q + 2));
        for (Index e = ov.sampStart[sl]; e < ov.sampStart[sl + 1]; ++e)
          push(ov.sampIdx[static_cast<std::size_t>(e)],
               -offv[k] / rsv * ov.sampW[static_cast<std::size_t>(e)]);
      }
    } else {
      // Virtually clean: target = the level-aware conservative regular row (assembleOperator's
      // regular branch, reproduced coefficient-for-coefficient).
      const double invV = 1.0 / mom.lap().cellVolume(i);
      double dsum = 0.0;
      mom.lap().forEachFaceNeighbor(i, [&](Index j, Real cf, int, double a) {
        push(j, -(-mu * invV * (a * cf)));  // minus the target's off-diagonal
        dsum += a * cf;
      });
      push(i, -(idiag + mu * invV * dsum));
    }
  }
  CfCsr out;
  detail::compactCsr(per, out, n);
  if (nSkipped)
    *nSkipped = skipped;
  return out;
}
 
inline double gpsDirGrad(const GhostOverlaySampled& ov, Index r, const std::vector<double>& fld,
                         int c, double invh) {
  auto fl = [&](int q) {
    return ov.sampFluid[static_cast<std::size_t>(r * 15 + c * 5 + (q + 2))] != 0;
  };
  auto F = [&](int q) { return gpsSample(ov, r, c, q, fld); };
  const bool ap = fl(+1), am = fl(-1);
  if (am && ap)
    return (F(+1) - F(-1)) * 0.5 * invh;
  if (ap)
    return fl(+2) ? (-3.0 * F(0) + 4.0 * F(+1) - F(+2)) * 0.5 * invh : (F(+1) - F(0)) * invh;
  if (am)
    return fl(-2) ? (3.0 * F(0) - 4.0 * F(-1) + F(-2)) * 0.5 * invh : (F(0) - F(-1)) * invh;
  return 0.0;
}
 
struct GhostGradCsrHost {
  std::vector<Index> cell;   
  std::vector<Index> start;  
  std::vector<Index> idx;
  std::vector<double> w;
};
 
template <unsigned Bits, class IsCutFn, class OkFn>
inline GhostGradCsrHost buildSampledGradCsr(const GhostOverlaySampled& ov,
                                            const BlockOctree<3, Bits>& t,
                                            const AmrPoisson<3, Bits>& pres, IsCutFn&& isCut,
                                            OkFn&& classicOk) {
  GhostGradCsrHost g;
  g.start.assign(1, 0);
  const Index n = t.numLeaves();
  auto pushSlot = [&](Index r, int a, int q, double wgt) {
    const std::size_t s = static_cast<std::size_t>(r * 15 + a * 5 + (q + 2));
    for (Index e = ov.sampStart[s]; e < ov.sampStart[s + 1]; ++e) {
      g.idx.push_back(ov.sampIdx[static_cast<std::size_t>(e)]);
      g.w.push_back(wgt * ov.sampW[static_cast<std::size_t>(e)]);
    }
  };
  auto pushCell = [&](Index j, double wgt) {
    g.idx.push_back(j);
    g.w.push_back(wgt);
  };
  for (Index i = 0; i < n; ++i) {
    if (!isCut(i))
      continue;
    g.cell.push_back(i);
    const double h = pres.cellWidth(i);
    const Index r = ov.rowOf[static_cast<std::size_t>(i)];
    for (int a = 0; a < 3; ++a) {
      if (r >= 0) {
        auto fl = [&](int q) {
          return ov.sampFluid[static_cast<std::size_t>(r * 15 + a * 5 + (q + 2))] != 0;
        };
        const bool ap = fl(+1), am = fl(-1);
        if (am && ap) {
          pushSlot(r, a, +1, 0.5 / h);
          pushSlot(r, a, -1, -0.5 / h);
        } else if (ap) {
          if (fl(+2)) {
            pushSlot(r, a, 0, -1.5 / h);
            pushSlot(r, a, +1, 2.0 / h);
            pushSlot(r, a, +2, -0.5 / h);
          } else {
            pushSlot(r, a, +1, 1.0 / h);
            pushSlot(r, a, 0, -1.0 / h);
          }
        } else if (am) {
          if (fl(-2)) {
            pushSlot(r, a, 0, 1.5 / h);
            pushSlot(r, a, -1, -2.0 / h);
            pushSlot(r, a, -2, 0.5 / h);
          } else {
            pushSlot(r, a, 0, 1.0 / h);
            pushSlot(r, a, -1, -1.0 / h);
          }
        }  // sandwiched: no entries (gradient 0)
      } else {
        // Classic same-level cascade (buildGhostGradOverlay's per-cell body).
        const Index jp = pres.periodicNeighbor(i, a, +1);
        const Index jm = pres.periodicNeighbor(i, a, -1);
        const bool ap = classicOk(jp, i), am = classicOk(jm, i);
        if (am && ap) {
          pushCell(jp, 0.5 / h);
          pushCell(jm, -0.5 / h);
        } else if (ap) {
          const Index jpp = pres.periodicNeighbor(jp, a, +1);
          if (classicOk(jpp, i)) {
            pushCell(i, -1.5 / h);
            pushCell(jp, 2.0 / h);
            pushCell(jpp, -0.5 / h);
          } else {
            pushCell(jp, 1.0 / h);
            pushCell(i, -1.0 / h);
          }
        } else if (am) {
          const Index jmm = pres.periodicNeighbor(jm, a, -1);
          if (classicOk(jmm, i)) {
            pushCell(i, 1.5 / h);
            pushCell(jm, -2.0 / h);
            pushCell(jmm, 0.5 / h);
          } else {
            pushCell(i, 1.0 / h);
            pushCell(jm, -1.0 / h);
          }
        }
      }
      g.start.push_back(static_cast<Index>(g.idx.size()));
    }
  }
  return g;
}
 
// ---- device mirror + kernels (Kokkos TUs only; include after a Kokkos-carrying header) ---------
#ifdef KOKKOS_INLINE_FUNCTION
 
struct GhostOverlaySampledDev {
  Index n = 0;
  View<Index> cell;
  View<float> rescale;
  View<int8_t> coupled, state;
  View<float> w_n1, w_n2, wm_n1, wm_n2;
  View<double> invh;
  View<Index> sampStart, sampIdx;
  View<double> sampW;
};
 
inline GhostOverlaySampledDev uploadGhostOverlaySampled(const GhostOverlaySampled& h) {
  GhostOverlaySampledDev d;
  d.n = h.base.n;
  if (h.base.n == 0)
    return d;
  d.cell = toDevice(h.base.cell, "gps_cell");
  d.rescale = toDevice(h.base.rescale, "gps_rescale");
  d.coupled = toDevice(h.base.coupled, "gps_coupled");
  d.state = toDevice(h.base.state, "gps_state");
  d.w_n1 = toDevice(h.base.w_n1, "gps_wn1");
  d.w_n2 = toDevice(h.base.w_n2, "gps_wn2");
  d.wm_n1 = toDevice(h.base.wm_n1, "gps_wmn1");
  d.wm_n2 = toDevice(h.base.wm_n2, "gps_wmn2");
  d.invh = toDevice(h.base.invh, "gps_invh");
  d.sampStart = toDevice(h.sampStart, "gps_sstart");
  d.sampIdx = toDevice(h.sampIdx, "gps_sidx");
  d.sampW = toDevice(h.sampW, "gps_sw");
  return d;
}
 
inline void ghostApplyDeltaSampled(const GhostOverlaySampledDev& ov, View<const double> x,
                                   View<double> y) {
  if (ov.n == 0)
    return;
  auto cell = ov.cell;
  auto resc = ov.rescale;
  auto cpl = ov.coupled;
  auto st = ov.state;
  auto wm1 = ov.wm_n1;
  auto wm2 = ov.wm_n2;
  auto invh = ov.invh;
  auto ss = ov.sampStart;
  auto si = ov.sampIdx;
  auto sw = ov.sampW;
  Kokkos::parallel_for(
      "amr::gps_apply_delta", ov.n, KOKKOS_LAMBDA(const Index r) {
        const Index c = cell(r);
        if (!cpl(r)) {
          y(c) = 0.0;
          return;
        }
        // M2c: the +/- faces of an axis read overlapping slots, so a row re-gathered the same
        // slot's CSR sum up to twice. k = 2a and 2a+1 are consecutive, so a FIVE-entry cache for
        // the current axis, invalidated when the axis changes, catches every repeat at negligible
        // register cost. Each slot is still summed in its own CSR order, exactly once, and delta
        // still accumulates over k in the same sequence ⇒ pure read-dedup, bit-identical.
        double Sv[5];
        bool Sh[5] = {false, false, false, false, false};
        int aCur = -1;
        auto X = [&](int a, int q) {
          if (a != aCur) {
            aCur = a;
            for (int m = 0; m < 5; ++m)
              Sh[m] = false;
          }
          const int idx = q + 2;
          if (!Sh[idx]) {
            const Index s = r * 15 + a * 5 + idx;
            double v = 0.0;
            for (Index e = ss(s); e < ss(s + 1); ++e)
              v += sw(e) * x(si(e));
            Sv[idx] = v;
            Sh[idx] = true;
          }
          return Sv[idx];
        };
        // Phase 3: the ratio form (per-axis weights, exactly 1.0 on a cubic octree).
        const double ih = invh(r * 3);
        const double rx = invh(r * 3 + 1) / ih, ry = invh(r * 3 + 2) / ih;
        const double r2[3] = {1.0, rx * rx, ry * ry};
        double delta = 0.0;
        for (int k = 0; k < 6; ++k) {
          const int8_t s = st(r * 6 + k);
          if (s != scheme::GP_QUAD && s != scheme::GP_LIN)
            continue;
          const int a = k / 2;
          const int sgn = (k & 1) ? -1 : 1;
          const int mn = (k & 1) ? 1 : 0;
          const int mf = (k & 1) ? 2 : -1;
          const double w1 = wm1(r * 6 + k), w2 = wm2(r * 6 + k);
          delta += r2[a] * (sgn * w1 * (X(a, mn) - X(a, mn - 1)));
          if (s == scheme::GP_QUAD && w2 != 0.0)
            delta += r2[a] * (sgn * w2 * (X(a, mf) - X(a, mf - 1)));
        }
        y(c) = resc(r) * (y(c) + ih * ih * delta);
      });
}
 
inline void ghostDivergDeltaSampled(const GhostOverlaySampledDev& ov, View<const double> u0,
                                    View<const double> u1, View<const double> u2, View<double> d) {
  if (ov.n == 0)
    return;
  auto cell = ov.cell;
  auto resc = ov.rescale;
  auto cpl = ov.coupled;
  auto st = ov.state;
  auto w1v = ov.w_n1;
  auto w2v = ov.w_n2;
  auto invh = ov.invh;
  auto ss = ov.sampStart;
  auto si = ov.sampIdx;
  auto sw = ov.sampW;
  Kokkos::parallel_for(
      "amr::gps_diverg_delta", ov.n, KOKKOS_LAMBDA(const Index r) {
        const Index c = cell(r);
        if (!cpl(r)) {
          d(c) = 0.0;
          return;
        }
        // M2c (see ghostApplyDeltaSampled): per-axis five-slot cache. Here the sharing is even
        // heavier — every U() straddles two slots and consecutive faces overlap — so the same
        // dedup applies verbatim, and each slot keeps its own CSR summation order.
        double Sv[5];
        bool Sh[5] = {false, false, false, false, false};
        int aCur = -1;
        auto S = [&](int a, int q) {
          if (a != aCur) {
            aCur = a;
            for (int m = 0; m < 5; ++m)
              Sh[m] = false;
          }
          const int idx = q + 2;
          if (!Sh[idx]) {
            const Index s = r * 15 + a * 5 + idx;
            double v = 0.0;
            for (Index e = ss(s); e < ss(s + 1); ++e) {
              const Index j = si(e);
              v += sw(e) * ((a == 0) ? u0(j) : (a == 1) ? u1(j) : u2(j));
            }
            Sv[idx] = v;
            Sh[idx] = true;
          }
          return Sv[idx];
        };
        auto U = [&](int a, int m) { return 0.5 * (S(a, m - 1) + S(a, m)); };
        // Phase 3: the ratio form (see ghostDivergDeltaSampledHost).
        const double ih = invh(r * 3);
        const double r1[3] = {1.0, invh(r * 3 + 1) / ih, invh(r * 3 + 2) / ih};
        double dd = 0.0;
        for (int k = 0; k < 6; ++k) {
          const int8_t s = st(r * 6 + k);
          if (s == scheme::GP_COUPLED)
            continue;
          const int a = k / 2;
          const int sgn = (k & 1) ? -1 : 1;
          const int mg = (k & 1) ? 0 : 1;
          const int mn = (k & 1) ? 1 : 0;
          const int mf = (k & 1) ? 2 : -1;
          if (s == scheme::GP_EXPLICIT) {
            dd += r1[a] * (sgn * U(a, mg));
            continue;
          }
          if (s == scheme::GP_BC_ONLY)
            continue;
          double val = w1v(r * 6 + k) * U(a, mn);
          if (s == scheme::GP_QUAD)
            val += w2v(r * 6 + k) * U(a, mf);
          dd += r1[a] * (sgn * val);
        }
        d(c) = resc(r) * (d(c) + ih * dd);
      });
}
 
struct GhostGradCsrDev {
  Index n = 0;
  View<Index> cell, start, idx;
  View<double> w;
};
 
inline GhostGradCsrDev uploadGhostGradCsr(const GhostGradCsrHost& h) {
  GhostGradCsrDev d;
  d.n = static_cast<Index>(h.cell.size());
  if (d.n == 0)
    return d;
  d.cell = toDevice(h.cell, "gcs_cell");
  d.start = toDevice(h.start, "gcs_start");
  d.idx = toDevice(h.idx, "gcs_idx");
  d.w = toDevice(h.w, "gcs_w");
  return d;
}
 
inline void applyGhostGradCsr(const GhostGradCsrDev& ov, View<const double> f, View<double> gx,
                              View<double> gy, View<double> gz) {
  if (ov.n == 0)
    return;
  auto cell = ov.cell;
  auto start = ov.start;
  auto idx = ov.idx;
  auto w = ov.w;
  Kokkos::parallel_for(
      "amr::flow_ghostgrad_csr", ov.n, KOKKOS_LAMBDA(const Index s) {
        const Index i = cell(s);
        double g[3];
        for (int a = 0; a < 3; ++a) {
          double acc = 0.0;
          for (Index e = start(s * 3 + a); e < start(s * 3 + a + 1); ++e)
            acc += w(e) * f(idx(e));
          g[a] = acc;
        }
        gx(i) = g[0];
        gy(i) = g[1];
        gz(i) = g[2];
      });
}
 
#endif  // KOKKOS_INLINE_FUNCTION
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_GHOST_PROJECTION_SAMPLED_HPP