cut_cell.hpp

Go to the documentation of this file.

// core — Robust-Scaled cut-cell Dirichlet operator on a BlockOctree.
//
// A faithful port of flow's ξ-polynomial sub-cell boundary scheme
// (flow/src/cut_cell_ibm.hpp: poly_*, ibmFillEntry, ibmModifyStencil) onto the
// cell-centered octree. Where the openness/aperture scheme (poisson.hpp) imposes a
// *Neumann* wall (no flux through solid faces), this imposes a *Dirichlet* value
// u = u_bc on the immersed boundary located at the true sub-cell distance ξ·h
// (Shortley–Weller): for a cut cell whose neighbour in some direction is solid,
// the 7-point stencil is modified by the boundary-distance polynomials, with
// D_rescale row-scaling for the small-cell problem. This is the velocity-diffusion
// / scalar-Dirichlet half of the cut-cell IBM (the user's "ξ-polynomial sub-cell
// BC"); the cell fluid-volume fraction κ classifies solid/fluid/cut cells.
//
// Cut cells are assumed to have same-level face neighbours (the suite contract:
// resolve the immersed boundary in a uniformly-finest band, so cut cells never sit
// on a 2:1 interface — see docs/AMR.md). 3D (flow's 6-direction scheme).
#ifndef PECLET_AMR_CUT_CELL_HPP
#define PECLET_AMR_CUT_CELL_HPP
 
#include <algorithm>
#include <array>
#include <cmath>
#include <vector>
 
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/common.hpp"
#include "peclet/amr/face_csr.hpp"  // shared host+device assembled-operator row kernels
#include "peclet/amr/leaf_field.hpp"
#include "peclet/amr/poisson.hpp"
#include "peclet/core/common/host_parallel.hpp"
#include "peclet/core/common/types.hpp"
#include "peclet/core/scheme/cut_cell_closure.hpp"  // the shared poly_* (QUALITY_PLAN G.2 lift)
 
namespace peclet::amr {
 
// ---- boundary-distance polynomials (QUALITY_PLAN G.2: the templated copy shared with flow now
//      lives in core/include/peclet/core/scheme/cut_cell_closure.hpp; these are double-precision
//      forwards so every `cc::poly_*` call site below is unchanged) ----
namespace cc {
// MORTON_HD (from face_csr.hpp): KOKKOS_FUNCTION on a Kokkos build, empty otherwise — so
// buildCutStencil and these polynomials are device-callable under the AMR device assembler
// (assembly.hpp) yet compile unchanged in the pure-C++ host oracle build. poly_abs replaces
// std::fabs (host-only under a CUDA device pass); it is bit-identical to std::fabs for every value
// buildCutStencil feeds it (the only difference, fabs(-0.0)=+0.0 vs −0.0, never occurs and would
// not change the |·|< comparisons). NOT part of the closure family (core keeps no fabs
// substitute), so it stays local.
MORTON_HD inline double poly_abs(double x) {
  return x < 0.0 ? -x : x;
}
MORTON_HD inline double poly_D(double xi) {
  return peclet::core::scheme::poly_D<double>(xi);
}
MORTON_HD inline double poly_N_nb(double xi) {
  return peclet::core::scheme::poly_N_nb<double>(xi);
}
MORTON_HD inline double poly_Nc(double xi) {
  return peclet::core::scheme::poly_Nc<double>(xi);
}
MORTON_HD inline double poly_Nbc(double xi) {
  return peclet::core::scheme::poly_Nbc<double>(xi);
}
MORTON_HD inline double poly_D_sandwich(double xm, double xp) {
  return peclet::core::scheme::poly_D_sandwich<double>(xm, xp);
}
MORTON_HD inline double poly_N_c_sandwich(double xm, double xp) {
  return peclet::core::scheme::poly_N_c_sandwich<double>(xm, xp);
}
MORTON_HD inline double poly_Nbc_pp_sw(double xm, double xp) {
  return peclet::core::scheme::poly_Nbc_pp_sw<double>(xm, xp);
}
MORTON_HD inline double poly_Nbc_mp_sw(double xm, double xp) {
  return peclet::core::scheme::poly_Nbc_mp_sw<double>(xm, xp);
}
}  // namespace cc
 
template <unsigned Bits = 21u>
class AmrCutCell {
 public:
  static constexpr int Dim = 3;
  using Octree = BlockOctree<3, Bits>;
  using M = typename Octree::M;
  using Code = typename Octree::Code;
  using Coord = typename Octree::Coord;
 
  // direction k: 0=+x,1=-x,2=+y,3=-y,4=+z,5=-z (flow order); OPP swaps sides.
  static constexpr int OPP[6] = {1, 0, 3, 2, 5, 4};
 
  void init(const Octree& t, Real h0, Vec<3> origin = Vec<3>{}) {
    init(t, detail::filledVec<3>(h0), origin);
  }
  void init(const Octree& t, const Vec<3>& h0, Vec<3> origin = Vec<3>{}) {
    t_ = &t;
    h0_ = h0;
    origin_ = origin;
    extResolve_ = nullptr;  // distributed seam re-installed by the owner after every init
    ghostLo_.clear();
    ghostLv_.clear();
    frameShift_ = {};
    for (int d = 0; d < 3; ++d)
      fineExt_[d] = static_cast<Coord>(t.brick()[d] * (Index(1) << t.lmax()));
  }
 
  // ---- distributed seam (docs/amr_distributed_flow.md): mirrors AmrPoisson's. build() forwards
  // all three into the internal lap_ (which it re-inits), so the whole momentum geometry —
  // regular C/F rows and the ξ-overlay band — resolves cross-block probes through the LeafHalo
  // registry and evaluates the SDF in the GLOBAL frame (bit-identical samples on every rank).
  using ExtResolver = typename AmrPoisson<3, Bits>::ExtResolver;
  void setResolver(ExtResolver r) { extResolve_ = std::move(r); }
  void setGhosts(std::vector<std::array<long, 3>> lo, std::vector<unsigned> lv) {
    ghostLo_ = std::move(lo);
    ghostLv_ = std::move(lv);
  }
  Index numGhosts() const { return static_cast<Index>(ghostLv_.size()); }
  void setFrameShift(const std::array<long, 3>& s) { frameShift_ = s; }
 
  Index numLeaves() const { return t_->numLeaves(); }
  bool isFluid(Index i) const { return fluid_[static_cast<std::size_t>(i)]; }
  bool isCut(Index i) const { return cut_[static_cast<std::size_t>(i)] != 0; }
  Index neighborOf(Index i, int k) const { return nb_[static_cast<std::size_t>(i) * 6 + k]; }
  double kappa(Index i) const { return kappa_[static_cast<std::size_t>(i)]; }
  double rhsScale(Index i) const { return rscale_[static_cast<std::size_t>(i)]; }
 
  // ---- read-only views of the built geometry, for the device assembler (momentum_assembly.hpp)
  // to stage to the device and reproduce build()/assembleOperator there. Mirrors how AmrPoisson
  // exposes its openness to the device FV assembler.
  const AmrPoisson<3, Bits>& lap() const {
    return lap_;
  }  
  const std::vector<double>& sdfCRaw() const {
    return sdfC_;
  }  
  const std::vector<Index>& nbRaw() const { return nb_; }  
  const std::vector<char>& fluidRaw() const { return fluid_; }
  const std::vector<char>& cutRaw() const { return cut_; }
  const std::vector<double>& acRaw() const { return AC_; }
  const std::vector<double>& offRaw() const { return off_; }  
  const std::vector<double>& rscaleRaw() const { return rscale_; }
  double idiag() const { return idiag_; }
  double mu() const { return mu_; }
  Vec<3> beta() const {
    Vec<3> b{};
    for (int d = 0; d < 3; ++d)
      b[d] = mu_ / (h0_[d] * h0_[d]);
    return b;
  }
  bool hasAdv() const { return hasAdv_; }
  const std::vector<double>& advDiagRaw() const { return advDiag_; }
  const std::vector<double>& advCoefRaw() const { return advCoef_; }
  const std::vector<Index>& advStartRaw() const { return advStart_; }
  const std::vector<Index>& advNbrRaw() const { return advNbr_; }
 
  template <class SdfFn>
  void build(SdfFn&& sdfFn, double idiag, double beta, int nsub = 4) {
    build(std::forward<SdfFn>(sdfFn), idiag, detail::filledVec<3>(beta), nsub);
  }
  template <class SdfFn>
  void build(SdfFn&& sdfFn, double idiag = 0.0, const Vec<3>& beta = detail::filledVec<3>(1.0),
             int nsub = 4) {
    const Index n = numLeaves();
    const Index ng = numGhosts();
    // C/F-aware Laplacian provider for regular (non-cut) fluid cells: an AmrPoisson
    // with NO openness (α=1) -> the plain ∇² with 2:1-interface coeff(si,sj). The
    // cut cells (finest, same-level) keep the ξ overlay below.
    lap_.init(*t_, h0_);
    if (extResolve_) {  // distributed: forward the seam into the freshly-init'd lap_
      lap_.setResolver(extResolve_);
      lap_.setGhosts(ghostLo_, ghostLv_);
      lap_.setFrameShift(frameShift_);
    }
    idiag_ = idiag;
    mu_ = beta[0] * h0_[0] * h0_[0];  // physical μ (operator A = idiag·I − μ∇²)
    sdfC_.assign(static_cast<std::size_t>(n + ng), 0.0);
    kappa_.assign(static_cast<std::size_t>(n), 0.0);
    fluid_.assign(static_cast<std::size_t>(n + ng), false);
    cut_.assign(static_cast<std::size_t>(n), 0);
    hasAdv_ = false;  // advection FOU is rebuilt per step via buildAdvectionFou
    AC_.assign(static_cast<std::size_t>(n), 1.0);
    rscale_.assign(static_cast<std::size_t>(n), 1.0);
    inhom_.assign(static_cast<std::size_t>(n), 0.0);
    off_.assign(static_cast<std::size_t>(n) * 6, 0.0);
    nb_.assign(static_cast<std::size_t>(n) * 6, -1);
 
    // Host-parallel leaf loop (setup-parallel plan rung 3; distributed included since the F1
    // resolution). `neighbor()` goes through AmrPoisson::probeSlot, which in a distributed
    // build calls `extResolve_` = LeafHalo::resolveGlobal; that path is now thread-safe — the
    // miss registration (the ONLY mutation reachable from a build pass) is mutex-guarded, and
    // the miss SET is order-canonical (coord-keyed std::map), so ghost numbering and the
    // fixpoint's rounds are bitwise-identical to the serial build at any thread count. See
    // LeafHalo::resolve and docs/amr_setup_parallel_plan.md ## Findings (F1 RESOLVED).
    auto forLeaves = [&](auto&& body) { hostParFor(n, body); };
 
    // Pass 1: cell-centre SDF, κ (subsampled), fluid flag, neighbour indices.
    // Single-rank: every write is to leaf i's OWN slots (sdfC_[i], fluid_[i], kappa_[i],
    // nb_[6i..6i+5]) and sdfFn is a pure geometry query — the disjoint-write pattern,
    // bitwise-invariant under any thread count. Pass 2 reads NEIGHBOURS' sdfC_, so it must be a
    // separate parallel region (the shim fences between them).
    forLeaves([&](Index i) {
      Vec<3> c = cellCenter(i);
      double sc = sdfFn(c);
      sdfC_[static_cast<std::size_t>(i)] = sc;
      fluid_[static_cast<std::size_t>(i)] = sc > 0.0;
      kappa_[static_cast<std::size_t>(i)] = volumeFraction(i, sdfFn, nsub);
      for (int k = 0; k < 6; ++k)
        nb_[static_cast<std::size_t>(i) * 6 + k] = neighbor(i, k);
    });
    // Ghost metadata (distributed): the SAME world SdfFn sampled at the ghost cell centre in
    // the GLOBAL frame — bit-identical to the owner's Pass-1 sample of that leaf, so the
    // ξ-classifications (and hence the CSR coefficients) agree across ranks exactly.
    for (Index g = 0; g < ng; ++g) {
      const auto& lo = ghostLo_[static_cast<std::size_t>(g)];
      const double s = static_cast<double>(1L << ghostLv_[static_cast<std::size_t>(g)]);
      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 sc = sdfFn(c);
      sdfC_[static_cast<std::size_t>(n + g)] = sc;
      fluid_[static_cast<std::size_t>(n + g)] = sc > 0.0;
    }
 
    // Pass 2: build per-leaf stencil.
    // Phase 3: the diagonal is 2*(beta_x + beta_y + beta_z), parenthesised as ONE sum so a cubic
    // octree reproduces the single correctly-rounded `6.0*beta` (Rule B).
    //
    // beta arrives as mu/h_finest^2 (AmrFlow::betaPerAxis), so it is the coefficient for a leaf at
    // the FINEST level. A leaf at level L has width 2^L * h_finest, so its viscous coefficient is
    // beta / 4^L. Scaling here, per leaf, rather than hoisting one beta/AC0 out of the loop: the
    // xi-overlay row is a pure multiple of the beta it is handed, so a level-L cut cell built with
    // the finest beta is 4^L too stiff -- an effective viscosity mu*4^L in the wall cells only.
    // That was measured as a silent constant velocity offset of -(3/8)(F/mu) h_leaf^2 (1 - 4^-L)
    // on plane Poiseuille (3.6-4.7% at L=1..3), with a perfect-looking parabola; see
    // docs/ROADMAP.md A5. Every other row in the operator was already level-aware: the regular
    // fluid rows use physical areas and volumes (assembleOperator below, -mu*invV*(a*c)), the
    // velocity-MG coarse operator computes mu/cellWidth^2 per row (velocity_mg.hpp), and the
    // mixed-level seam delta does exactly this scaling for seam rows only
    // (ghost_projection_sampled.hpp "the ROW-LOCAL beta per axis") -- which is why a uniformly
    // COARSE band, having no seam rows, kept the wrong coefficient.
    //
    // Exactly inert at the finest level: f == 1.0, so betaArr and AC0 are bit-for-bit the old beta
    // and AC0 expressions, and any run whose cut cells are all finest is byte-identical.
    forLeaves([&](Index i) {  // own-leaf slots only (rung 3, architect-pre-cleared: no
                              // neighbour-indexed stores in this body)
      const double f = static_cast<double>(Index(1) << t_->level(i));  // h_leaf / h_finest
      const double inv = 1.0 / (f * f);
      const double betaArr[3] = {beta[0] * inv, beta[1] * inv, beta[2] * inv};
      const double AC0 = idiag + ((2.0 * betaArr[0] + 2.0 * betaArr[1]) + 2.0 * betaArr[2]);
      if (!fluid_[static_cast<std::size_t>(i)]) {  // solid: identity row u=0
        AC_[static_cast<std::size_t>(i)] = 1.0;
        for (int k = 0; k < 6; ++k)
          off_[static_cast<std::size_t>(i) * 6 + k] = 0.0;
        return;
      }
      double sdf_n[6];
      bool anyGhost = false;
      for (int k = 0; k < 6; ++k) {
        Index j = nb_[static_cast<std::size_t>(i) * 6 + k];
        sdf_n[k] = (j >= 0) ? sdfC_[static_cast<std::size_t>(j)] : -1.0;  // missing => solid
        if (sdf_n[k] < 0.0)
          anyGhost = true;
      }
      cut_[static_cast<std::size_t>(i)] = anyGhost ? 1 : 0;
      double AC = AC0, off[6];
      for (int k = 0; k < 6; ++k)
        off[k] = -betaArr[k / 2];  // face k belongs to axis k/2
      double rscale = 1.0, inhomCoef = 0.0;
      if (anyGhost)
        buildCutStencil(sdfC_[static_cast<std::size_t>(i)], sdf_n, betaArr, AC0, AC, off, rscale,
                        inhomCoef);
      AC_[static_cast<std::size_t>(i)] = AC;
      for (int k = 0; k < 6; ++k)
        off_[static_cast<std::size_t>(i) * 6 + k] = off[k];
      rscale_[static_cast<std::size_t>(i)] = rscale;
      inhom_[static_cast<std::size_t>(i)] = inhomCoef;
    });
  }
 
  struct Assembled {
    std::vector<double> diag;  
    std::vector<Index> start;  
    std::vector<Index> nbr;    
    std::vector<double> coef;  
  };
  Assembled assembleOperator(bool scaleAdvByRscale = false) const {
    const Index n = numLeaves();
    Assembled A;
    A.diag.assign(static_cast<std::size_t>(n), 0.0);
    A.start.assign(static_cast<std::size_t>(n) + 1, 0);
    std::vector<std::vector<std::pair<Index, double>>> rows(static_cast<std::size_t>(n));
    for (Index i = 0; i < n; ++i) {
      const std::size_t s = static_cast<std::size_t>(i);
      if (!fluid_[s]) {  // solid: identity row (u = u_bc)
        A.diag[s] = 1.0;
        continue;
      }
      if (cut_[s]) {  // ξ-overlay stencil (same-level neighbours)
        A.diag[s] = AC_[s];
        for (int k = 0; k < 6; ++k) {
          double a = off_[s * 6 + k];
          if (a == 0.0)
            continue;
          Index j = nb_[s * 6 + k];
          if (j >= 0)
            rows[s].emplace_back(j, a);
        }
      } else {  // regular fluid: idiag·I − μ∇² with C/F-aware face coupling
        const double invV = 1.0 / lap_.cellVolume(i);
        double dsum = 0.0;
        lap_.forEachFaceNeighbor(i, [&](Index j, Real c, int, double a) {
          rows[s].emplace_back(j, -mu_ * invV * (a * c));
          dsum += a * c;
        });
        A.diag[s] = idiag_ + mu_ * invV * dsum;
      }
      if (hasAdv_) {  // implicit-FOU advection: diagonal (outflow) + CSR (inflow)
        const double as = scaleAdvByRscale ? rscale_[s] : 1.0;
        A.diag[s] += as * advDiag_[s];
        for (std::size_t p = static_cast<std::size_t>(advStart_[s]);
             p < static_cast<std::size_t>(advStart_[static_cast<std::size_t>(i) + 1]); ++p)
          rows[s].emplace_back(advNbr_[p], as * advCoef_[p]);
      }
    }
    for (Index i = 0; i < n; ++i)
      A.start[static_cast<std::size_t>(i) + 1] =
          A.start[static_cast<std::size_t>(i)] +
          static_cast<Index>(rows[static_cast<std::size_t>(i)].size());
    const Index nnz = A.start[static_cast<std::size_t>(n)];
    A.nbr.resize(static_cast<std::size_t>(nnz));
    A.coef.resize(static_cast<std::size_t>(nnz));
    for (Index i = 0; i < n; ++i) {
      Index k = A.start[static_cast<std::size_t>(i)];
      for (auto& e : rows[static_cast<std::size_t>(i)]) {
        A.nbr[static_cast<std::size_t>(k)] = e.first;
        A.coef[static_cast<std::size_t>(k)] = e.second;
        ++k;
      }
    }
    return A;
  }
 
  // ---- Runtime operator: the assembled CSR applied with the SHARED face_csr.hpp row kernels — the
  // exact same arithmetic the device runs (momentum.hpp), executed serially here. assembleOperator
  // folds the implicit-FOU advection into the single CSR, so the FaceCsrOpT view sets hasAdv=false.
  // The *Geometric variants above are the independent reference the oracle tests check this
  // against.
 
  FaceCsrOpT<HostArr<double>, HostArr<Index>> hostOp(const Assembled& A) const {
    FaceCsrOpT<HostArr<double>, HostArr<Index>> v;
    v.n = numLeaves();
    v.diag = HostArr<double>(A.diag.data());
    v.coef = HostArr<double>(A.coef.data());
    v.start = HostArr<Index>(A.start.data());
    v.nbr = HostArr<Index>(A.nbr.data());
    v.hasAdv = false;  // advection already folded into the single CSR by assembleOperator
    return v;
  }
 
  void applyOp(const std::vector<double>& u, std::vector<double>& out) const {
    const Assembled A = assembleOperator();
    const auto op = hostOp(A);
    const Index n = numLeaves();
    out.assign(static_cast<std::size_t>(n), 0.0);
    const HostArr<double> uacc(u.data());
    for (Index i = 0; i < n; ++i)
      out[static_cast<std::size_t>(i)] = faceCsrApplyRow(op, i, uacc);
  }
 
  double residual(const std::vector<double>& u, const std::vector<double>& b,
                  std::vector<double>& res) const {
    const Assembled A = assembleOperator();
    const auto op = hostOp(A);
    const Index n = numLeaves();
    res.assign(static_cast<std::size_t>(n), 0.0);
    const HostArr<double> uacc(u.data());
    double s = 0.0;
    for (Index i = 0; i < n; ++i) {
      double r = b[static_cast<std::size_t>(i)] - faceCsrApplyRow(op, i, uacc);
      res[static_cast<std::size_t>(i)] = r;
      if (fluid_[static_cast<std::size_t>(i)])
        s += r * r;
    }
    return std::sqrt(s);
  }
 
  void gaussSeidel(std::vector<double>& u, const std::vector<double>& b, int sweeps) const {
    const Assembled A = assembleOperator();
    const auto op = hostOp(A);
    const Index n = numLeaves();
    const HostArr<double> uacc(u.data());
    for (int s = 0; s < sweeps; ++s)
      for (Index i = 0; i < n; ++i) {
        double off, d;
        faceCsrOffDiag(op, i, uacc, off, d);
        u[static_cast<std::size_t>(i)] = faceCsrPointUpdate(b[static_cast<std::size_t>(i)], off, d,
                                                            u[static_cast<std::size_t>(i)], 1.0);
      }
  }
 
  void applyOpGeometric(const std::vector<double>& u, std::vector<double>& out) const {
    const Index n = numLeaves();
    out.assign(static_cast<std::size_t>(n), 0.0);
    std::vector<double> Lu;
    lap_.applyLaplacian(u, Lu);  // C/F-aware ∇² (α=1); used only for regular fluid cells
    for (Index i = 0; i < n; ++i) {
      const std::size_t s = static_cast<std::size_t>(i);
      if (!fluid_[s]) {
        out[s] = u[s];
      } else if (cut_[s]) {
        double acc = AC_[s] * u[s];
        for (int k = 0; k < 6; ++k) {
          double a = off_[s * 6 + k];
          if (a == 0.0)
            continue;
          Index j = nb_[s * 6 + k];
          if (j >= 0)
            acc += a * u[static_cast<std::size_t>(j)];
        }
        out[s] = acc;
      } else {
        out[s] = idiag_ * u[s] - mu_ * Lu[s];  // regular fluid (C/F-consistent)
      }
      if (hasAdv_ && fluid_[s])
        out[s] += advApply(i, u);  // implicit FOU advection
    }
  }
 
  void buildAdvectionFou(const std::array<std::vector<double>, 3>& uadv, double rho,
                         const std::vector<double>& uf, const std::vector<Index>& faceStart,
                         bool useFace) {
    const Index n = numLeaves();
    advDiag_.assign(static_cast<std::size_t>(n), 0.0);
    advStart_.assign(static_cast<std::size_t>(n) + 1, 0);
    hasAdv_ = true;
    auto velOutOf = [&](Index i, Index j, int axis, int dir, Index slot) {
      return useFace ? dir * uf[static_cast<std::size_t>(slot)]
                     : dir * 0.5 *
                           (uadv[axis][static_cast<std::size_t>(i)] +
                            uadv[axis][static_cast<std::size_t>(j)]);
    };
    // pass 1: count inflow (off-diagonal) fluid faces per cell.
    for (Index i = 0; i < n; ++i) {
      if (!fluid_[static_cast<std::size_t>(i)])
        continue;
      int cnt = 0;
      Index s = useFace ? faceStart[static_cast<std::size_t>(i)] : 0;
      lap_.forEachFaceFull(i, [&](Index j, int axis, int dir, double, double, double) {
        if (fluid_[static_cast<std::size_t>(j)] && velOutOf(i, j, axis, dir, s) < 0.0)
          ++cnt;
        ++s;
      });
      advStart_[static_cast<std::size_t>(i) + 1] = cnt;
    }
    for (Index i = 0; i < n; ++i)
      advStart_[static_cast<std::size_t>(i) + 1] += advStart_[static_cast<std::size_t>(i)];
    advNbr_.assign(advStart_[static_cast<std::size_t>(n)], -1);
    advCoef_.assign(advStart_[static_cast<std::size_t>(n)], 0.0);
    // pass 2: fill diagonal (outflow) + CSR off-diagonals (inflow).
    for (Index i = 0; i < n; ++i) {
      if (!fluid_[static_cast<std::size_t>(i)])
        continue;
      const double Vi = lap_.cellVolume(i);
      std::size_t pos = static_cast<std::size_t>(advStart_[static_cast<std::size_t>(i)]);
      Index s = useFace ? faceStart[static_cast<std::size_t>(i)] : 0;
      lap_.forEachFaceFull(i, [&](Index j, int axis, int dir, double area, double, double) {
        const Index slot = s++;
        if (!fluid_[static_cast<std::size_t>(j)])
          return;
        double velOut = velOutOf(i, j, axis, dir, slot);
        double w = rho * area * velOut / Vi;
        if (velOut < 0.0) {  // inflow → couple to upstream neighbour j (matches pass 1)
          advNbr_[pos] = j;
          advCoef_[pos] = w;
          ++pos;
        } else {
          advDiag_[static_cast<std::size_t>(i)] += w;  // outflow → diagonal
        }
      });
    }
  }
 
  std::vector<double> makeRhs(const std::vector<double>& src, double u_bc) const {
    const Index n = numLeaves();
    std::vector<double> b(static_cast<std::size_t>(n), 0.0);
    for (Index i = 0; i < n; ++i) {
      if (!fluid_[static_cast<std::size_t>(i)]) {
        b[static_cast<std::size_t>(i)] = u_bc;  // solid held at u_bc
        continue;
      }
      b[static_cast<std::size_t>(i)] =
          src[static_cast<std::size_t>(i)] * rscale_[static_cast<std::size_t>(i)] +
          inhom_[static_cast<std::size_t>(i)] * u_bc;
    }
    return b;
  }
 
  double residualGeometric(const std::vector<double>& u, const std::vector<double>& b,
                           std::vector<double>& res) const {
    applyOpGeometric(u, res);
    double s = 0.0;
    const Index n = numLeaves();
    for (Index i = 0; i < n; ++i) {
      double r = b[static_cast<std::size_t>(i)] - res[static_cast<std::size_t>(i)];
      res[static_cast<std::size_t>(i)] = r;
      if (fluid_[static_cast<std::size_t>(i)])
        s += r * r;
    }
    return std::sqrt(s);
  }
 
  void gaussSeidelGeometric(std::vector<double>& u, const std::vector<double>& b,
                            int sweeps) const {
    const Index n = numLeaves();
    for (int s = 0; s < sweeps; ++s)
      for (Index i = 0; i < n; ++i) {
        const std::size_t si = static_cast<std::size_t>(i);
        if (!fluid_[si]) {
          u[si] = b[si];
        } else if (cut_[si]) {  // ξ-overlay stencil (same-level neighbours)
          double sum = b[si];
          for (int k = 0; k < 6; ++k) {
            double a = off_[si * 6 + k];
            if (a == 0.0)
              continue;
            Index j = nb_[si * 6 + k];
            if (j >= 0)
              sum -= a * u[static_cast<std::size_t>(j)];
          }
          double d = AC_[si];
          if (hasAdv_) {
            sum -= advOffSum(i, u);
            d += advDiag_[si];
          }
          if (d != 0.0)
            u[si] = sum / d;
        } else {  // regular fluid: idiag·u − μ∇² with C/F-aware face coupling
          double offsum = 0.0, dsum = 0.0;
          lap_.forEachFaceNeighbor(i, [&](Index j, Real c, int, double) {
            offsum += c * u[static_cast<std::size_t>(j)];
            dsum += c;
          });
          double Vi = lap_.cellVolume(i);
          double diagA = idiag_ + mu_ * dsum / Vi;
          double sum = b[si] + mu_ * offsum / Vi;
          if (hasAdv_) {
            sum -= advOffSum(i, u);
            diagA += advDiag_[si];
          }
          u[si] = sum / diagA;
        }
      }
  }
 
 public:
  double fouApply(Index i, const std::vector<double>& field) const {
    return advDiag_[static_cast<std::size_t>(i)] * field[static_cast<std::size_t>(i)] +
           advOffSum(i, field);
  }
  bool hasAdvection() const { return hasAdv_; }
 
 private:
  // Σ_csr coef·field[nbr] — the implicit-FOU advection off-diagonal coupling.
  double advOffSum(Index i, const std::vector<double>& field) const {
    double s = 0.0;
    for (std::size_t p = static_cast<std::size_t>(advStart_[static_cast<std::size_t>(i)]);
         p < static_cast<std::size_t>(advStart_[static_cast<std::size_t>(i) + 1]); ++p)
      s += advCoef_[p] * field[static_cast<std::size_t>(advNbr_[p])];
    return s;
  }
  double advApply(Index i, const std::vector<double>& u) const { return fouApply(i, u); }
 
 public:
  // Port of ibmFillEntry<0> + ibmModifyStencil for one cut cell (Dirichlet). Public + MORTON_HD so
  // the device assembler (momentum_assembly.hpp) runs the SAME per-cell stencil build on device.
  MORTON_HD static void buildCutStencil(double sdf_c, const double sdf_n[6], const double beta[3],
                                        double AC0, double& ACout, double off[6], double& rscaleOut,
                                        double& inhomOut) {
    bool ghost[6];
    double xi[6], D[6];
    for (int k = 0; k < 6; ++k) {
      if (sdf_n[k] < 0.0) {
        ghost[k] = true;
        double th = sdf_c / (sdf_c - sdf_n[k]);
        th = th < 1e-4 ? 1e-4 : (th > 1.0 ? 1.0 : th);
        xi[k] = th;
        D[k] = cc::poly_D(th);
      } else {
        ghost[k] = false;
        xi[k] = 1.0;
        D[k] = 1e9;
      }
    }
    bool sand[3] = {ghost[0] && ghost[1], ghost[2] && ghost[3], ghost[4] && ghost[5]};
    double Dsand[3] = {0, 0, 0};
    for (int a = 0; a < 3; ++a)
      if (sand[a])
        Dsand[a] = cc::poly_D_sandwich(xi[2 * a + 1], xi[2 * a]);
    double minAbs = 1e30, descale = 1.0;
    auto upd = [&](double v) {
      if (cc::poly_abs(v) < minAbs) {
        minAbs = cc::poly_abs(v);
        descale = v;
      }
    };
    for (int a = 0; a < 3; ++a) {
      if (sand[a])
        upd(Dsand[a]);
      else {
        if (ghost[2 * a])
          upd(D[2 * a]);
        if (ghost[2 * a + 1])
          upd(D[2 * a + 1]);
      }
    }
    double K[6] = {0}, Mf[6] = {1, 1, 1, 1, 1, 1}, X[6] = {0}, Nbc[6] = {0},
           R[6] = {1, 1, 1, 1, 1, 1};
    for (int a = 0; a < 3; ++a) {
      int km = 2 * a + 1, kp = 2 * a;
      double Daxis = sand[a] ? Dsand[a] : (ghost[kp] ? D[kp] : (ghost[km] ? D[km] : descale));
      double r = descale / Daxis;
      if (cc::poly_abs(Daxis) < 1e-9)
        r = 1.0;
      R[kp] = R[km] = r;
      if (sand[a]) {
        K[kp] = cc::poly_N_c_sandwich(xi[km], xi[kp]) * r;
        K[km] = cc::poly_N_c_sandwich(xi[kp], xi[km]) * r;
        Nbc[kp] = (cc::poly_Nbc_pp_sw(xi[km], xi[kp]) + cc::poly_Nbc_mp_sw(xi[km], xi[kp])) * r;
        Nbc[km] = (cc::poly_Nbc_pp_sw(xi[kp], xi[km]) + cc::poly_Nbc_mp_sw(xi[kp], xi[km])) * r;
        Mf[kp] = Mf[km] = 0.0;
      } else {
        for (int side = 0; side < 2; ++side) {
          int kk = side == 0 ? kp : km;
          if (ghost[kk]) {
            K[kk] = cc::poly_Nc(xi[kk]) * r;
            X[kk] = cc::poly_N_nb(xi[kk]) * r;
            Nbc[kk] = cc::poly_Nbc(xi[kk]) * r;
            Mf[kk] = 0.0;
          } else {
            K[kk] = 0.0;
            Mf[kk] = 1.0;
            X[kk] = 0.0;
            Nbc[kk] = 0.0;
          }
        }
      }
    }
    // ibmModifyStencil (orig off-diagonal = -beta for every direction). OPP is a function-local
    // constexpr (not the static member) so the runtime index mod[OPP[k]] is device-safe under CUDA.
    constexpr int OPP_[6] = {1, 0, 3, 2, 5, 4};
    double aC = AC0 * descale, mod[6] = {0, 0, 0, 0, 0, 0}, inhom = 0.0;
    for (int k = 0; k < 6; ++k) {
      double vnb = -beta[k / 2];
      aC += vnb * K[k];
      inhom += Nbc[k] * vnb;
      mod[k] += vnb * (descale * Mf[k] - 1.0);
      mod[OPP_[k]] += vnb * X[k];
    }
    ACout = aC;
    for (int k = 0; k < 6; ++k)
      off[k] = -beta[k / 2] + mod[k];
    rscaleOut = descale;
    inhomOut = inhom;
  }
 
 private:
  Vec<3> cellCenter(Index i) const {
    auto b = t_->bounds(i);
    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>(static_cast<long>(b[0][d]) + frameShift_[d]) + 0.5 * s) * h0_[d];
    return c;
  }
 
  template <class SdfFn>
  double volumeFraction(Index i, SdfFn&& sdfFn, int nsub) const {
    auto b = t_->bounds(i);
    double s = static_cast<double>(Index(1) << t_->level(i));
    // Phase 3: the subsample grid spans the BOX leaf, so the step is per axis.
    Vec<3> w{};
    int inside = 0, total = nsub * nsub * nsub;
    Vec<3> base{};
    for (int d = 0; d < 3; ++d) {
      w[d] = s * h0_[d];
      base[d] =
          origin_[d] + static_cast<double>(static_cast<long>(b[0][d]) + frameShift_[d]) * h0_[d];
    }
    for (int a = 0; a < nsub; ++a)
      for (int bb = 0; bb < nsub; ++bb)
        for (int cc2 = 0; cc2 < nsub; ++cc2) {
          Vec<3> p{base[0] + (a + 0.5) / nsub * w[0], base[1] + (bb + 0.5) / nsub * w[1],
                   base[2] + (cc2 + 0.5) / nsub * w[2]};
          if (sdfFn(p) > 0.0)
            ++inside;
        }
    return static_cast<double>(inside) / total;
  }
 
  Index neighbor(Index i, int k) const {
    return lap_.periodicNeighbor(i, k / 2, (k % 2 == 0) ? +1 : -1);
  }
 
  const Octree* t_ = nullptr;
  Vec<3> h0_ = detail::filledVec<3>(1.0);
  Vec<3> origin_{};
  std::array<Coord, 3> fineExt_{};
  ExtResolver extResolve_;                    // distributed seam (forwarded into lap_ by build)
  std::vector<std::array<long, 3>> ghostLo_;  // ghost slot → block-local lo (longs)
  std::vector<unsigned> ghostLv_;             // ghost slot → covering-leaf level
  std::array<long, 3> frameShift_{};          // block global fine origin (0 single-rank)
  std::vector<double> sdfC_, kappa_, AC_, rscale_, inhom_, off_;
  std::vector<Index> nb_;
  std::vector<char> fluid_, cut_;
  AmrPoisson<3, Bits> lap_;  // C/F-aware ∇² provider for regular fluid cells (α=1)
  double idiag_ = 0.0, mu_ = 1.0;
  std::vector<double> advDiag_, advCoef_;  // implicit-FOU advection (rebuilt per step)
  std::vector<Index> advStart_, advNbr_;   // CSR off-diagonals (C/F-conservative)
  bool hasAdv_ = false;
};
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_CUT_CELL_HPP