flow.hpp

Go to the documentation of this file.

// core — device (Kokkos) collocated incompressible Stokes step on a BlockOctree.
//
// The device counterpart of oracle::AmrFlow (flow.hpp): the whole cut-cell IBM projection step
// runs in Kokkos kernels instead of the host-serial gaussSeidel + host projection the drag
// study found to be the bottleneck. It reuses the host AmrCutCell / AmrPoisson to *build*
// the operators (geometry, openness, cut stencils — done once), then drives the time step
// entirely on the device:
//   * momentum predictor — MomentumOp (assembled cut-cell operator) solved with the
//     parallel BiCGStab of momentum.hpp;
//   * pressure projection — the openness Poisson on Multigrid / PCG;
//   * divergence, ABC gradient correction, rotational pressure update — face-CSR kernels
//     (FaceGeom) that mirror AmrPoisson::forEachFaceFull (same 2:1 sub-faces +
//     openness), so D / G / L stay consistent exactly as in the host collocated coupling.
//
// Stokes (advection off) and Navier–Stokes (setAdvection): implicit-FOU + explicit SOU/Koren
// deferred correction, advected by the divergence-free face field uf (built each projection;
// falls back to ½(u_i+u_j) until the first projection) — conservative (∇·uf = 0).
//
#ifndef PECLET_AMR_FLOW_HPP
#define PECLET_AMR_FLOW_HPP
 
#include <array>
#include <chrono>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
 
#include "peclet/amr/adapt.hpp"         // transferField (conservative remap for finishAdapt)
#include "peclet/amr/advect_recon.hpp"  // shared high-order face reconstruction (host+device)
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/cf_scheme.hpp"  // pluggable 2:1 C/F schemes (setCfScheme)
#include "peclet/amr/common.hpp"
#include "peclet/amr/cut_cell.hpp"
#include "peclet/amr/distributed_adapt.hpp"    // transferGradients (distributed finishAdapt)
#include "peclet/amr/distributed_flow_mg.hpp"  // distributed pressure MG (initMpi mode)
#include "peclet/amr/distributed_octree.hpp"
#include "peclet/amr/face_geom.hpp"          // FaceGeom (shared with the device assembler)
#include "peclet/amr/facegeom_assembly.hpp"  // assembleFaceGeom (D4/D6)
#include "peclet/amr/ghost_projection.hpp"   // directional ghost overlay (setGhostProjection)
#include "peclet/amr/ghost_projection_sampled.hpp"  // mixed-level sampled overlay (setGhostSampled)
#include "peclet/amr/leaf_halo.hpp"
#include "peclet/amr/momentum.hpp"
#include "peclet/amr/momentum_assembly.hpp"  // assembleMomentum (D3/D6)
#include "peclet/amr/multigrid.hpp"
#include "peclet/amr/pcg.hpp"
#include "peclet/amr/poisson.hpp"
#include "peclet/amr/seam_recon.hpp"  // seam-reconstruction tables (docs/amr_cf_convective.md)
#include "peclet/amr/velocity_mg.hpp"
#include "peclet/core/common/host_parallel.hpp"
#include "peclet/core/common/types.hpp"
#include "peclet/core/common/view.hpp"
 
namespace peclet::amr {
 
inline bool amrEnvFlag(const char* name) {
  const char* v = std::getenv(name);
  return v && v[0] && !(v[0] == '0' && v[1] == '\0');
}
 
// FaceGeom (the collocated projection's static face-geometry CSR) now lives in face_geom.hpp so the
// device assembler and this driver share the type without a circular include.
 
template <int Dim, unsigned Bits, class FluidFn>
FaceGeom buildFaceGeom(const AmrPoisson<Dim, Bits>& ap, FluidFn&& isFluid) {
  const Index n = ap.octree().numLeaves();
  std::vector<Index> start(static_cast<std::size_t>(n) + 1, 0);
  for (Index i = 0; i < n; ++i) {
    Index cnt = 0;
    ap.forEachFaceFull(i, [&](Index, int, int, double, double, double) { ++cnt; });
    start[static_cast<std::size_t>(i) + 1] = start[static_cast<std::size_t>(i)] + cnt;
  }
  const Index nf = start[static_cast<std::size_t>(n)];
  std::vector<Index> nbr(static_cast<std::size_t>(nf));
  std::vector<int> axis(static_cast<std::size_t>(nf)), dir(static_cast<std::size_t>(nf));
  std::vector<double> aArea(static_cast<std::size_t>(nf)), rArea(static_cast<std::size_t>(nf)),
      dist(static_cast<std::size_t>(nf)), alpha(static_cast<std::size_t>(nf));
  std::vector<Index> upupI(static_cast<std::size_t>(nf)), upupJ(static_cast<std::size_t>(nf));
  std::vector<double> invVol(static_cast<std::size_t>(n));
  std::vector<char> fluid(static_cast<std::size_t>(n));
  for (Index i = 0; i < n; ++i) {
    invVol[static_cast<std::size_t>(i)] = 1.0 / ap.cellVolume(i);
    fluid[static_cast<std::size_t>(i)] = isFluid(i) ? 1 : 0;
    Index k = start[static_cast<std::size_t>(i)];
    ap.forEachFaceFull(i, [&](Index j, int ax, int dr, double area, double d, double al) {
      nbr[static_cast<std::size_t>(k)] = j;
      axis[static_cast<std::size_t>(k)] = ax;
      dir[static_cast<std::size_t>(k)] = dr;
      aArea[static_cast<std::size_t>(k)] = al * area;
      rArea[static_cast<std::size_t>(k)] = area;
      dist[static_cast<std::size_t>(k)] = d;
      alpha[static_cast<std::size_t>(k)] = al;
      // SOU upstream-of-upwind probes (point neighbour one cell further upstream): if i is the
      // upwind cell the upstream is across i's −dir face; if j is upwind, across j's +dir face.
      upupI[static_cast<std::size_t>(k)] = ap.periodicNeighbor(i, ax, -dr);
      upupJ[static_cast<std::size_t>(k)] = ap.periodicNeighbor(j, ax, dr);
      ++k;
    });
  }
  FaceGeom g;
  g.n = n;
  g.start = toDevice(start, "fg_start");
  g.nbr = toDevice(nbr, "fg_nbr");
  g.axis = toDevice(axis, "fg_axis");
  g.dir = toDevice(dir, "fg_dir");
  g.alphaArea = toDevice(aArea, "fg_aarea");
  g.rawArea = toDevice(rArea, "fg_rarea");
  g.dist = toDevice(dist, "fg_dist");
  g.alpha = toDevice(alpha, "fg_alpha");
  g.upupI = toDevice(upupI, "fg_upupi");
  g.upupJ = toDevice(upupJ, "fg_upupj");
  g.invVol = toDevice(invVol, "fg_invvol");
  g.fluid = toDevice(fluid, "fg_fluid");
  return g;
}
 
inline void divergence(const FaceGeom& g, View<const double> u0, View<const double> u1,
                       View<const double> u2, View<double> div) {
  auto st = g.start;
  auto nb = g.nbr;
  auto ax = g.axis;
  auto dr = g.dir;
  auto aA = g.alphaArea;
  auto iv = g.invVol;
  auto fl = g.fluid;
  Kokkos::parallel_for(
      "amr::flow_div", g.n, KOKKOS_LAMBDA(const Index i) {
        if (!fl(i)) {
          div(i) = 0.0;
          return;
        }
        double d = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k) {
          const int a = ax(k);
          const double ui = (a == 0) ? u0(i) : (a == 1) ? u1(i) : u2(i);
          const Index j = nb(k);
          const double uj = (a == 0) ? u0(j) : (a == 1) ? u1(j) : u2(j);
          d += aA(k) * dr(k) * 0.5 * (ui + uj);
        }
        div(i) = d * iv(i);
      });
}
 
inline void buildFaceField(const FaceGeom& g, View<const double> u0, View<const double> u1,
                           View<const double> u2, View<const double> phi, View<double> uf) {
  auto st = g.start;
  auto nb = g.nbr;
  auto ax = g.axis;
  auto dr = g.dir;
  auto di = g.dist;
  Kokkos::parallel_for(
      "amr::flow_buildface", g.n, KOKKOS_LAMBDA(const Index i) {
        for (Index k = st(i); k < st(i + 1); ++k) {
          const int a = ax(k);
          const Index j = nb(k);
          const double ui = (a == 0) ? u0(i) : (a == 1) ? u1(i) : u2(i);
          const double uj = (a == 0) ? u0(j) : (a == 1) ? u1(j) : u2(j);
          const double gphi = (dr(k) > 0) ? (phi(j) - phi(i)) / di(k) : (phi(i) - phi(j)) / di(k);
          uf(k) = 0.5 * (ui + uj) - gphi;
        }
      });
}
 
inline double divFaceNorm(const FaceGeom& g, View<const double> uf) {
  auto st = g.start;
  auto dr = g.dir;
  auto aA = g.alphaArea;
  auto iv = g.invVol;
  auto fl = g.fluid;
  double s = 0.0;
  Kokkos::parallel_reduce(
      "amr::flow_divface", g.n,
      KOKKOS_LAMBDA(const Index i, double& acc) {
        if (!fl(i))
          return;
        double d = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k)
          d += aA(k) * dr(k) * uf(k);
        d *= iv(i);
        acc += d * d;
      },
      s);
  return std::sqrt(s);
}
 
inline void grad3(const FaceGeom& g, View<const double> f, View<double> gx, View<double> gy,
                  View<double> gz) {
  auto st = g.start;
  auto nb = g.nbr;
  auto ax = g.axis;
  auto dr = g.dir;
  auto di = g.dist;
  auto al = g.alpha;
  auto fl = g.fluid;
  Kokkos::parallel_for(
      "amr::flow_grad3", g.n, KOKKOS_LAMBDA(const Index i) {
        if (!fl(i)) {
          gx(i) = gy(i) = gz(i) = 0.0;
          return;
        }
        const double fi = f(i);
        double gp[3] = {0, 0, 0}, gm[3] = {0, 0, 0};
        int np[3] = {0, 0, 0}, nm[3] = {0, 0, 0};
        for (Index k = st(i); k < st(i + 1); ++k) {
          if (al(k) <= 1e-12)
            continue;
          const int a = ax(k);
          const double gg = (dr(k) > 0) ? (f(nb(k)) - fi) / di(k) : (fi - f(nb(k))) / di(k);
          if (dr(k) > 0) {
            gp[a] += gg;
            ++np[a];
          } else {
            gm[a] += gg;
            ++nm[a];
          }
        }
        double out[3];
        for (int a = 0; a < 3; ++a) {
          const double a1 = np[a] ? gp[a] / np[a] : 0.0;
          const double a2 = nm[a] ? gm[a] / nm[a] : 0.0;
          out[a] = 0.5 * (a1 + a2);
        }
        gx(i) = out[0];
        gy(i) = out[1];
        gz(i) = out[2];
      });
}
 
struct GhostGradOverlay {
  Index n = 0;       
  View<Index> cell;  
  View<Index> idx;   
  View<double> w;    
};
 
inline void applyGhostGrad(const GhostGradOverlay& ov, View<const double> f, View<double> gx,
                           View<double> gy, View<double> gz) {
  if (ov.n == 0)
    return;
  auto cell = ov.cell;
  auto idx = ov.idx;
  auto w = ov.w;
  Kokkos::parallel_for(
      "amr::flow_ghostgrad", 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 (int k = 0; k < 3; ++k)
            acc += w(s * 9 + a * 3 + k) * f(idx(s * 9 + a * 3 + k));
          g[a] = acc;
        }
        gx(i) = g[0];
        gy(i) = g[1];
        gz(i) = g[2];
      });
}
 
inline void momRhs(View<const double> uc, View<const double> gradP, View<const double> adv,
                   View<const double> rscale, View<const char> fluid, double idiag, double fc,
                   View<double> b, Index n) {
  Kokkos::parallel_for(
      "amr::flow_momrhs", n, KOKKOS_LAMBDA(const Index i) {
        b(i) = fluid(i) ? (idiag * uc(i) + fc - gradP(i) - adv(i)) * rscale(i) : 0.0;
      });
}
 
inline void buildFou(const FaceGeom& g, View<const double> u0, View<const double> u1,
                     View<const double> u2, double rho, View<const double> rscale,
                     View<double> advDiag, View<double> advCoef, View<const double> uf,
                     bool useFace) {
  auto st = g.start;
  auto nb = g.nbr;
  auto ax = g.axis;
  auto dr = g.dir;
  auto ra = g.rawArea;
  auto iv = g.invVol;
  auto fl = g.fluid;
  Kokkos::parallel_for(
      "amr::flow_buildfou", g.n, KOKKOS_LAMBDA(const Index i) {
        if (!fl(i)) {
          advDiag(i) = 0.0;
          for (Index k = st(i); k < st(i + 1); ++k)
            advCoef(k) = 0.0;
          return;
        }
        const double rs = rscale(i);
        double dsum = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k) {
          const Index j = nb(k);
          if (!fl(j)) {
            advCoef(k) = 0.0;
            continue;
          }
          const int a = ax(k);
          const double ui = (a == 0) ? u0(i) : (a == 1) ? u1(i) : u2(i);
          const double uj = (a == 0) ? u0(j) : (a == 1) ? u1(j) : u2(j);
          const double velOut = useFace ? dr(k) * uf(k) : dr(k) * 0.5 * (ui + uj);
          const double w = rs * rho * ra(k) * velOut * iv(i);
          if (velOut < 0.0) {
            advCoef(k) = w;  // inflow → off-diagonal toward upstream neighbour j
          } else {
            advCoef(k) = 0.0;
            dsum += w;  // outflow → diagonal
          }
        }
        advDiag(i) = dsum;
      });
}
 
template <class FldFn>
KOKKOS_INLINE_FUNCTION double advectedFaceValue(const SeamView& sv, Index k, bool upIsI,
                                                double phiUpUp, double phiUp, double phiDown,
                                                double distK, int advScheme, const FldFn& fld) {
  const Index s = sv.on ? sv.seam(k) : Index(-1);
  if (s < 0)
    return hoFaceValue(phiUpUp, phiUp, phiDown, advScheme);  // today's reconstruction, bit for bit
  const Index rs = upIsI ? sv.sampI(s) : sv.sampJ(s);        // φ_U* (case 1)
  const Index ru = upIsI ? sv.uuRecI(s) : sv.uuRecJ(s);      // φ_UU* (cases 2/3)
  const Index rd = upIsI ? sv.sampJ(s) : sv.sampI(s);        // φ_D*, Koren only
  const double d1 = upIsI ? sv.d1I(s) : sv.d1J(s);
  const double upStar = (rs >= 0) ? sv.gather(rs, fld) : phiUp;
  // No upstream record: the probe is at the upwind cell's own level (or was withheld by a gate),
  // which is the equal-spacing `dUU = 2·d1` today's formula already assumes.
  double uu = phiUpUp, dUU = 2.0 * d1;
  if (ru >= 0) {
    uu = sv.gather(ru, fld);
    dUU = sv.recDist(ru);
  }
  const double downStar = (advScheme != 0 && rd >= 0) ? sv.gather(rd, fld) : phiDown;
  return hoFaceValueSeam(uu, phiUp, upStar, downStar, d1, dUU, distK, advScheme);
}
 
inline void deferredSou(const FaceGeom& g, View<const double> u0, View<const double> u1,
                        View<const double> u2, int comp, double rho, int advScheme,
                        View<double> defc, View<const double> uf, bool useFace, bool seamOn) {
  auto st = g.start;
  auto nb = g.nbr;
  auto ax = g.axis;
  auto dr = g.dir;
  auto ra = g.rawArea;
  auto ds = g.dist;
  auto iv = g.invVol;
  auto fl = g.fluid;
  auto uiP = g.upupI;
  auto ujP = g.upupJ;
  const SeamView sv = makeSeamView(g, seamOn);
  Kokkos::parallel_for(
      "amr::flow_defsou", g.n, KOKKOS_LAMBDA(const Index i) {
        if (!fl(i)) {
          defc(i) = 0.0;
          return;
        }
        auto fld = [&](Index c) { return (comp == 0) ? u0(c) : (comp == 1) ? u1(c) : u2(c); };
        double sou = 0.0, fou = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k) {
          const Index j = nb(k);
          if (!fl(j))
            continue;
          const int a = ax(k);
          const double uai = (a == 0) ? u0(i) : (a == 1) ? u1(i) : u2(i);
          const double uaj = (a == 0) ? u0(j) : (a == 1) ? u1(j) : u2(j);
          const double velOut = useFace ? dr(k) * uf(k) : dr(k) * 0.5 * (uai + uaj);
          const bool upIsI = velOut > 0.0;
          const Index up = upIsI ? i : j;
          const Index down = upIsI ? j : i;
          const Index upup = upIsI ? uiP(k) : ujP(k);
          const double phiUp = fld(up);
          const double phiUpUp = (upup >= 0 && fl(upup)) ? fld(upup) : phiUp;
          const double phiDown = fld(down);
          const double phiFace = advectedFaceValue(sv, k, upIsI, phiUpUp, phiUp, phiDown, ds(k),
                                                   advScheme, fld);  // shared recon
          sou += ra(k) * velOut * phiFace;
          fou += ra(k) * velOut * fld(up);  // FOU flux = velOut · upwind value
        }
        defc(i) = rho * iv(i) * (sou - fou);  // ρ·(SOU − FOU), unscaled
      });
}
 
inline void advectExplicit(const FaceGeom& g, View<const double> u0, View<const double> u1,
                           View<const double> u2, int comp, double rho, int advScheme,
                           View<double> defc, View<const double> uf, bool useFace, bool seamOn) {
  auto st = g.start;
  auto nb = g.nbr;
  auto ax = g.axis;
  auto dr = g.dir;
  auto ra = g.rawArea;
  auto ds = g.dist;
  auto iv = g.invVol;
  auto fl = g.fluid;
  auto uiP = g.upupI;
  auto ujP = g.upupJ;
  const SeamView sv = makeSeamView(g, seamOn);
  Kokkos::parallel_for(
      "amr::flow_advexpl", g.n, KOKKOS_LAMBDA(const Index i) {
        if (!fl(i)) {
          defc(i) = 0.0;
          return;
        }
        auto fld = [&](Index c) { return (comp == 0) ? u0(c) : (comp == 1) ? u1(c) : u2(c); };
        double sou = 0.0;
        for (Index k = st(i); k < st(i + 1); ++k) {
          const Index j = nb(k);
          if (!fl(j))
            continue;
          const int a = ax(k);
          const double uai = (a == 0) ? u0(i) : (a == 1) ? u1(i) : u2(i);
          const double uaj = (a == 0) ? u0(j) : (a == 1) ? u1(j) : u2(j);
          const double velOut = useFace ? dr(k) * uf(k) : dr(k) * 0.5 * (uai + uaj);
          const bool upIsI = velOut > 0.0;
          const Index up = upIsI ? i : j;
          const Index down = upIsI ? j : i;
          const Index upup = upIsI ? uiP(k) : ujP(k);
          const double phiUp = fld(up);
          const double phiUpUp = (upup >= 0 && fl(upup)) ? fld(upup) : phiUp;
          const double phiDown = fld(down);
          const double phiFace = advectedFaceValue(sv, k, upIsI, phiUpUp, phiUp, phiDown, ds(k),
                                                   advScheme, fld);  // shared recon
          sou += ra(k) * velOut * phiFace;
        }
        defc(i) = rho * sou * iv(i);  // ρ·SOU (fully explicit)
      });
}
 
inline void correct(View<double> uc, View<const double> gphi, View<const char> fluid, Index n) {
  Kokkos::parallel_for(
      "amr::flow_correct", n, KOKKOS_LAMBDA(const Index i) {
        if (fluid(i))
          uc(i) -= gphi(i);
      });
}
 
inline void presUpdate(View<double> p, View<const double> phi, View<const double> div,
                       View<const char> fluid, double rho_dt, double mu, Index n) {
  Kokkos::parallel_for(
      "amr::flow_presupd", n, KOKKOS_LAMBDA(const Index i) {
        if (fluid(i))
          p(i) += rho_dt * phi(i) - mu * div(i);
      });
}
 
// ===========================================================================
// AmrFlow — collocated Stokes projection step, fully on device.
// ===========================================================================
template <unsigned Bits = 21u>
class AmrFlow {
 public:
  using Octree = BlockOctree<3, Bits>;
 
  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;
    dist_ = nullptr;  // single-rank unless initMpi is used
  }
 
  void initMpi(DistributedOctree<3, Bits>& d) {
    init(d.local(), d.h0(), d.globalGeometry().origin);
    dist_ = &d;
  }
  Vec<3> betaPerAxis() const {
    Vec<3> b{};
    for (int d = 0; d < 3; ++d)
      b[d] = mu_ / (h0_[d] * h0_[d]);
    return b;
  }
  void setDensity(double rho) { rho_ = rho; }
  void setViscosity(double mu) { mu_ = mu; }
  void setDt(double dt) { dt_ = dt; }
  void setBodyForce(double fx, double fy, double fz) { f_ = {fx, fy, fz}; }
  void setPressurePCG(bool on) { presPCG_ = on; }
  void setGhostGradient(bool on) { ghostGrad_ = on; }
  void setApertureOrder(int order) { apertureOrder_ = order; }
  void setGhostProjection(bool on, int matrixOrder = 2, int rhsOrder = 2) {
    ghostProjReq_ = on ? 1 : 0;  // explicit selection disables the AUTO default
    gpMatrixOrder_ = matrixOrder;
    gpRhsOrder_ = rhsOrder;
  }
  void setGhostSampled(bool on, double rho = 2.2, long maxSamples = 0) {
    ghostSampledReq_ = on ? 1 : 0;
    gpsRho_ = rho;
    gpsMaxN_ = maxSamples;
  }
  void setCfScheme(int scheme) { cfScheme_ = static_cast<CfScheme>(scheme); }
  void setAdvection(bool on) { advect_ = on; }
  void setImplicitAdvection(bool on) { implicitFou_ = on; }
  void setAdvectionScheme(int s) { advScheme_ = s; }
  void setMomentumTol(double tol) { momTol_ = tol; }
  void setPressureTol(double tol) { presTol_ = tol; }
  void setUfAdvection(bool on) { ufAdvect_ = on; }
  void setSeamReconstruction(bool on) { seamRecon_ = on; }
  void setMomentumMG(bool on) { momMGon_ = on; }
 
  void setVelocityMGStaircase(bool on) { useStaircaseMG_ = on; }
  void setVelocityMGMinCoarse(Index m) { mgMinCoarse_ = m; }
  void setMomentumGS(bool on) { momGS_ = on; }
 
  void setMomentumMGSolver(bool on) { momMGSolver_ = on; }
 
  void setOuterIterations(int n, double tol = 1e-6) {
    outerIters_ = (n < 1) ? 1 : n;
    outerTol_ = tol;
  }
 
  template <class SdfFn>
  void setSolid(SdfFn&& sdfFn) {
    const Index n = t_->numLeaves();
    // Phase profiler (PECLET_AMR_PROFILE_SETUP=1): setSolid is the whole setup cost at bed
    // scale (measured 127 us/leaf single-threaded, [[performance-sota-yardstick]]) — the
    // per-phase breakdown is what any optimization must start from.
    const bool prof_ = [] {
      const char* e = std::getenv("PECLET_AMR_PROFILE_SETUP");
      return e && e[0] == '1';
    }();
    auto profT0_ = std::chrono::steady_clock::now();
    auto profPhase = [&](const char* name) {
      if (!prof_)
        return;
      Kokkos::fence();
      auto now = std::chrono::steady_clock::now();
      double ms = std::chrono::duration<double, std::milli>(now - profT0_).count();
      std::fprintf(stderr, "[setSolid] %-28s %9.1f ms  (%6.1f us/leaf)\n", name, ms,
                   1e3 * ms / static_cast<double>(n));
      profT0_ = std::chrono::steady_clock::now();
    };
    // Host operator build (geometry, openness, cut stencils) — same as oracle::AmrFlow::setSolid.
    mom_.init(*t_, h0_, origin_);
    pres_.init(*t_, h0_);
    pres_.setOrigin(origin_);
    profPhase("mom.init + pres.init");
    // Resolve the projection mode. DEFAULT SWITCH (2026-08-25, user decision): AUTO = the
    // GHOST (fluid-only) projection — family-free/stable/unique (see the setter doc) — falling
    // back to the aperture projection with a stderr notice when the finest band is too thin for
    // the overlay (probed below; the July AUTO arm restored, now unconditional). Explicit
    // setGhostProjection(true/false) pins the scheme (band violations then THROW as before).
    // The overlay is built below (it needs only mom_'s sdf samples + pres_'s topology walk,
    // no openness).
    const bool wantGhost = (ghostProjReq_ != 0);  // explicit on (1) or AUTO (-1)
    ghostProj_ = wantGhost;                       // provisional; AUTO may fall back below
    ghostSampled_ = (ghostSampledReq_ == 1) && wantGhost;
    // D1/D2 (docs/amr_march_perf_and_distributed_plan.md): the mixed-level cut band is
    // DISTRIBUTED. The sampled builders probe through the LeafHalo seam inside
    // prepareDistributed's fixpoint (D1) and their least-squares clouds read ghost slots like any
    // other CSR the step consumes (D2), so there is no single-rank guard left. Acceptance:
    // np=1 bitwise, np=2,4 within the established decomposition-independence class
    // (tests/test_amr_distributed_seam_mpi.cpp).
    if (dist_) {
      // Distributed: install the resolver seams and run every prober to the miss-collect
      // fixpoint (docs/amr_distributed_flow.md). Freezes the ±2 halo, leaves mom_ FULLY built
      // (its final round ran with every ghost resolved) and the solver hooks installed.
      prepareDistributed(sdfFn);
    } else {
      nExt_ = n;
      allred_ = {};
      momSolver_.setDistributed({}, {}, 0);
      pcg_.setDistributed({}, {}, 0);
      mom_.build(sdfFn, /*idiag=*/rho_ / dt_, /*beta=*/betaPerAxis());
    }
    profPhase("mom.build (SDF+cut stencils)");
    GhostOverlay hov;
    GhostOverlaySampled hovS;
    if (ghostSampled_) {  // mixed-level cut band: sample-slot overlay, no band-margin probe
      const auto gfine = globalFineExtent();
      hovS = buildGhostOverlaySampled(*t_, pres_, sdfFn, gpMatrixOrder_, gpRhsOrder_, origin_,
                                      &gfine, shiftD_, /*discovery=*/false, gpsRho_, gpsMaxN_);
      profPhase("buildGhostOverlaySampled");
    } else if (ghostProj_) {  // probe the band margin; explicit request throws on violation, AUTO
                              // falls
      bool viol = false;
      hov = buildGhostOverlay(*t_, pres_, mom_.sdfCRaw(), gpMatrixOrder_, gpRhsOrder_, &viol);
      if (dist_) {
        // COLLECTIVE band-margin decision: the flag is agreed across ranks before anyone
        // commits (one rank throwing/falling back alone would deadlock the MG collectives).
        int lv = viol ? 1 : 0, gv = 0;
        MPI_Allreduce(&lv, &gv, 1, MPI_INT, MPI_LOR, dist_->comm());
        viol = gv != 0;
      }
      if (viol) {
        if (ghostProjReq_ == 1)
          throw std::runtime_error(
              "amr ghost projection: an overlay row's ±2 closure reach crosses a 2:1 level "
              "boundary — widen the refineToSdf band margin");
        ghostProj_ = false;  // AUTO: fall back to the aperture projection
        fprintf(stderr,
                "peclet::core AmrFlow: AUTO scheme fell back to the aperture projection (the "
                "finest band is too thin for the ghost overlay). Select explicitly with "
                "setGhostProjection to silence this notice.\n");
      }
    }
    if (ghostProj_) {
      // Ghost projection: the pressure geometry is the BINARY openness on the unchanged MG
      // rails; the closure physics lives in the overlay (built above). Sampled mode: the
      // level-aware canonical rule (the sampled overlay's face states are FORCED to it — the
      // overlay-closed <=> binary-closed invariant).
      if (ghostSampled_) {
        auto binFn = makeBinaryOpenFnMixed(
            *t_, pres_, [&sdfFn](const Vec<3>& p) { return sdfFn(p); }, h0_, origin_, shiftD_);
        pres_.buildOpenness(binFn);
        profPhase("buildOpenness (mixed)");
        // D1: the pressure hierarchy has to be the DISTRIBUTED one when a seam is installed —
        // gpOp0() reads presMGD_ whenever dist_ is set, so building presMG_ here left the ghost
        // solver pointing at an empty operator.
        if (dist_)
          presMGD_.build(*dist_, h0_, binFn, &dhalo_, /*liftRoot=*/true, presBottomExtent_);
        else
          presMG_.build(*t_, h0_, binFn, /*periodic=*/true, /*immersedWall=*/false,
                        /*liftRoot=*/true, presBottomExtent_);
        profPhase("presMG.build");
      } else {
        auto binFn = makeBinaryOpenFn([&sdfFn](const Vec<3>& p) { return sdfFn(p); }, h0_);
        pres_.buildOpenness(binFn);
        if (dist_)
          presMGD_.build(*dist_, h0_, binFn, &dhalo_, /*liftRoot=*/true, presBottomExtent_);
        else
          presMG_.build(*t_, h0_, binFn, /*periodic=*/true, /*immersedWall=*/false,
                        /*liftRoot=*/true, presBottomExtent_);
      }
      ghostGrad_ = true;  // the directional gradient is part of the scheme
      // Fragmentation guard: pockets outside the main binary component are decoupled (see
      // findPocketCells) — folded into maskC_ below and hidden from the directional gradients.
      // Single-rank host BFS only: the DISTRIBUTED label-propagation guard is a rung-6 item
      // (a rank-local BFS would mislabel components that span ranks), so multi-rank runs are
      // unprotected on fragmenting geometries until it lands.
      gpPocket_ = dist_ ? std::vector<char>{} : findPocketCells(*t_, pres_, mom_.sdfCRaw());
      profPhase("findPocketCells");
    } else {
      gpPocket_.clear();
      auto openFn = [&](const Vec<3>& fc, int axis) { return faceFrac(sdfFn, fc, axis); };
      pres_.buildOpenness(openFn);
      if (dist_)
        presMGD_.build(*dist_, h0_, openFn, &dhalo_, /*liftRoot=*/true, presBottomExtent_);
      else
        presMG_.build(*t_, h0_, openFn, /*periodic=*/true, /*immersedWall=*/false,
                      /*liftRoot=*/true, presBottomExtent_);
    }
    // What solves the coarsest level (docs/amr_mg_depth.md §6.6): `auto` engages the agglomerated
    // GraphAMG-PCG bottom only where the ladder ran out above `presBottomExtent_`.
    if (dist_)
      presMGD_.setBottom(presBottom_);
    else
      presMG_.setBottom(presBottom_);
    // Singular periodic pressure: per-level nullspace projection.
    if (dist_) {
      presMGD_.setRemoveMean(true);
      pcg_.setDistributed([this](View<double> v) { presMGD_.sync(0, v); }, allred_,
                          presMGD_.extendedSize(0));
    } else {
      presMG_.setRemoveMean(true);
    }
    profPhase("nullspace setup");
 
    // Device assembly (D6): the static cut-cell momentum operator CSR and the collocated face
    // geometry are assembled ON THE DEVICE (assembleMomentum / assembleFaceGeom), and the
    // pressure MG operators are device-assembled per level (Multigrid D5) — so no host CSR walk and
    // no operator round-trip. Each is bit-for-bit identical to the host assembler on OpenMP (locked
    // in test_amr_momentum / test_amr_facegeom), so the flow result is unchanged. A shared device
    // octree view backs both assemblers.
    // DISTRIBUTED: assemble on the HOST through the resolver seam instead (the device walkers
    // cannot resolve cross-block probes) and upload — same parity-locked builders, ghost
    // columns included.
    BlockOctreeView<3, Bits> ov;
    if (!dist_)
      ov.upload(*t_);
    if (dist_) {
      auto A = mom_.assembleOperator();
      momOp_ = MomentumOp{};
      momOp_.n = n;
      momOp_.diag = toDevice(A.diag, "mo_diag");
      momOp_.faceStart = toDevice(A.start, "mo_start");
      momOp_.faceNbr = toDevice(A.nbr, "mo_nbr");
      momOp_.faceCoef = toDevice(A.coef, "mo_coef");
    } else {
      momOp_ = assembleMomentum<Bits>(mom_, ov);  // static Stokes operator (hasAdv stays false)
      profPhase("assembleMomentum (device)");
    }
    // Velocity multigrid (momentum preconditioner): the Galerkin hierarchy A_c = R·A·P built
    // directly from the exact assembled momentum CSR. Consistent with the fine cut-cell
    // operator by construction (inherits the ξ-overlay + D_rescale row scaling; a coarse cell
    // of all-solid children stays an identity row). It only changes the preconditioner (the
    // BiCGStab matvec is the exact operator) ⇒ same converged solution, but the iteration
    // count stays ~flat with N instead of growing like the Jacobi-preconditioned BiCGStab.
    // Build the chosen momentum-MG: Galerkin (MomentumMG) by default, or the rediscretized
    // staircase (VelocityMG) — both from the static Stokes operator, once.
    // DISTRIBUTED: the Galerkin MG becomes RANK-LOCAL (additive-Schwarz preconditioner): the
    // local rows of the exact operator with the ghost COLUMNS dropped. At np=1 nothing is
    // dropped, so the preconditioner — and hence the whole step — stays bit-identical to the
    // single-rank path; at np>1 it is a block preconditioner (iterations may grow with np,
    // the converged step is unchanged — the preconditioner never moves the solution). The
    // EXACT cross-rank Galerkin RAP is the documented follow-up
    // (docs/amr_distributed_flow.md §4c); the staircase variant is single-rank-only for now
    // (distributed falls back to Jacobi-with-halo preconditioning).
    if (momMGon_ && dist_ && !useStaircaseMG_) {
      auto A = mom_.assembleOperator();
      std::vector<Index> lstart(1, 0);
      std::vector<Index> lnbr;
      std::vector<double> lcoef;
      for (Index i = 0; i < n; ++i) {
        for (Index k = A.start[static_cast<std::size_t>(i)];
             k < A.start[static_cast<std::size_t>(i) + 1]; ++k)
          if (A.nbr[static_cast<std::size_t>(k)] < n) {  // drop ghost columns (Schwarz cut)
            lnbr.push_back(A.nbr[static_cast<std::size_t>(k)]);
            lcoef.push_back(A.coef[static_cast<std::size_t>(k)]);
          }
        lstart.push_back(static_cast<Index>(lnbr.size()));
      }
      momMG_.build(*t_, A.diag, lstart, lnbr, lcoef);
      momMG_.setGaussSeidel(momGS_);
    } else if (momMGon_ && !dist_) {
      if (useStaircaseMG_) {
        std::vector<double> kap(static_cast<std::size_t>(n));
        std::vector<char> fl(static_cast<std::size_t>(n)), cu(static_cast<std::size_t>(n));
        for (Index i = 0; i < n; ++i) {
          kap[static_cast<std::size_t>(i)] = mom_.kappa(i);
          fl[static_cast<std::size_t>(i)] = mom_.isFluid(i) ? 1 : 0;
          cu[static_cast<std::size_t>(i)] = mom_.isCut(i) ? 1 : 0;
        }
        velMG_.build(*t_, h0_, rho_ / dt_, mu_, momOp_, kap, fl, cu, mgMinCoarse_);
        velMG_.setGaussSeidel(momGS_);
      } else {
        // The Galerkin RAP hierarchy is a host triple-product over the fine CSR, so it needs the
        // operator on the host; assemble it there for the MG build only (bit-identical to momOp_).
        auto A = mom_.assembleOperator();
        momMG_.build(*t_, A.diag, A.start, A.nbr, A.coef);
        momMG_.setGaussSeidel(momGS_);
      }
    }
    profPhase("velocity MG build");
    std::vector<char> fluidVec(static_cast<std::size_t>(n));
    for (Index i = 0; i < n; ++i)
      fluidVec[static_cast<std::size_t>(i)] = mom_.isFluid(i) ? 1 : 0;
    if (dist_) {
      // Host walker through the resolver seam (nbr/upup may be ghost slots), then EXTEND the
      // fluid flags over the ghost tail — the advection kernels read fl(j)/fl(upup) at ghost
      // slots (mom_.fluidRaw() carries the ghost metadata filled from the world SdfFn).
      geom_ = buildFaceGeom(pres_, [&](Index i) { return mom_.isFluid(i); });
      geom_.fluid = toDevice(mom_.fluidRaw(), "fg_fluid_ext");
    } else {
      geom_ = assembleFaceGeom<Bits>(pres_, fluidVec, ov);
    }
    profPhase("assembleFaceGeom");
    std::vector<double> rs(static_cast<std::size_t>(n));
    for (Index i = 0; i < n; ++i)
      rs[static_cast<std::size_t>(i)] = mom_.rhsScale(i);
    rscale_ = toDevice(rs, "df_rscale");
    fluid_ = geom_.fluid;
    if (ghostGrad_ && !ghostSampled_)
      buildGhostGradOverlay();
    else
      gc_ = GhostGradOverlay{};  // sampled mode: the CSR overlay gcS_ owns all cut cells
    // C/F interface scheme overlays (cf_scheme.hpp): the same host builders the oracle uses
    // (parity by construction), uploaded once. Momentum delta = ×μ on the α=1 velocity geometry.
    // The per-FACE C/F gate (docs/amr_cf_flux_gate.md §6.1, §6.4) needs `regular = fluid && !cut`
    // at GHOST slots too: the gate is a property of the FACE, and a 2:1 sub-face at a block seam
    // has one of its two cells on another rank. `cut(g)` cannot be recomputed here — it is the SDF
    // sign at g's SIX face-neighbour centres, and the discovery fixpoint registers the C/F
    // tangential reach and the face sweep's ±2 upstream reach, not those six — so the OWNER's flag
    // is exchanged instead. One `syncScalar` on an nExt-sized View, exact by definition, and
    // decomposition-independent because the owner's own SDF samples are.
    std::vector<char> cfRegular(static_cast<std::size_t>(nExt_), 0);
    for (Index i = 0; i < n; ++i)
      cfRegular[static_cast<std::size_t>(i)] = (mom_.isFluid(i) && !mom_.isCut(i)) ? 1 : 0;
    if (dist_) {
      std::vector<double> rx(static_cast<std::size_t>(nExt_), 0.0);
      for (Index i = 0; i < n; ++i)
        rx[static_cast<std::size_t>(i)] = cfRegular[static_cast<std::size_t>(i)] ? 1.0 : 0.0;
      View<double> rv = toDevice(rx, "cf_regular");
      syncScalar(rv);
      const std::vector<double> rh = toVector(rv);
      for (Index s = n; s < nExt_; ++s)
        cfRegular[static_cast<std::size_t>(s)] = (rh[static_cast<std::size_t>(s)] > 0.5) ? 1 : 0;
    }
    // Programming-error guard (§6.4): a builder that asks for a slot the flag was never filled for
    // is a bug in the discovery, not a configuration the solver should tolerate silently.
    auto regularOk = [&](Index s) {
      if (s < 0 || s >= static_cast<Index>(cfRegular.size()))
        throw std::logic_error("peclet::amr: C/F regular flag queried outside [0, nExt)");
      return cfRegular[static_cast<std::size_t>(s)] != 0;
    };
    cfCutFaces_ = countCfCutFaces(pres_, *t_, regularOk, [&](Index s) { return mom_.isFluid(s); });
    if (cfScheme_ != CfScheme::standard) {
      auto fluidOk = [&](Index i) { return mom_.isFluid(i); };
      auto rowRegular = [&](Index i) { return mom_.isFluid(i) && !mom_.isCut(i); };
      auto lapD = buildCfLapDelta(mom_.lap(), *t_, mu_, rowRegular, fluidOk, cfScheme_);
      // The projection family (D, the ABC gradient substitution, and the advecting face field)
      // is gated PER FACE, never per row: `cfFace(i, j) = level mismatch && regular(i) &&
      // regular(j)`, one predicate through one emitter (docs/amr_cf_flux_gate.md). A face flux is
      // ONE number shared by two cells, so a gate that is a property of a cell cannot be
      // conservative — under the old row gate the regular cell booked the C/F correction in its
      // constraint, the cut cell did not, and uf carried it for both: a permanent mass source at
      // the cut cell, ‖D(Δvel) − Δ_cfDiv‖ = 2.018e-01 against ‖div(uf)‖ = 2.016e-01 on a mesh
      // whose cut band meets the level boundary. Cut rows still receive no delta (the 2026-08-27
      // stability fix — 2 of 12 throat-graded meshes marched to k~1e12 with cut rows in cfDiv —
      // is retained in full by the same predicate, since a cut row has no passing face); what
      // changed is that the REGULAR row across a cut-adjacent 2:1 face no longer receives it
      // either, and that face reverts to the standard two-point value. Inert by geometry on any
      // uniform or finest-band mesh (cut cells there have no C/F face); `numCfCutFaces()` counts
      // the sub-faces where it is not.
      auto divD = buildCfDivDelta(pres_, *t_, regularOk, fluidOk, cfScheme_);
      auto gd = buildCfGradDelta(pres_, *t_, regularOk, fluidOk, cfScheme_);
      auto ufd = buildCfUfDelta(pres_, *t_, regularOk, fluidOk, cfScheme_);
      // Distributed diagnostic (numCfGhostColumns): how many overlay entries read a GHOST slot.
      // It is the observable that says the seam is actually exercised — a build that silently
      // fell back to the raw coarse value at every block boundary reports zero.
      cfGhostCols_ = 0;
      auto countGhost = [&](const std::vector<Index>& col) {
        for (Index c : col)
          if (c >= n)
            ++cfGhostCols_;
      };
      countGhost(lapD.slot);
      countGhost(divD.slot);
      for (int a = 0; a < 3; ++a)
        countGhost(gd[static_cast<std::size_t>(a)].slot);
      countGhost(ufd.vel.slot);
      cfMom_ = uploadCfCsr(lapD, "cf_mom");
      cfDiv_ = uploadCfCompCsr(divD, "cf_div");
      for (int a = 0; a < 3; ++a)
        cfGrad_[static_cast<std::size_t>(a)] =
            uploadCfCsr(gd[static_cast<std::size_t>(a)], "cf_grad");
      cfUfVel_ = uploadCfCompCsr(ufd.vel, "cf_ufvel");
      // The convective flux at a 2:1 seam (ROADMAP B5, docs/amr_cf_convective.md §5.2): the
      // descriptors + records the two advective kernels reconstruct the ADVECTED value from.
      // It reuses the SAME per-face gate (`regularOk`) and the SAME tangential sample operator
      // (`cfAppendStencil`, quadratic — §5.1 Q-A: one such operator suite-wide) as the four
      // overlays above, which is why it is built here rather than beside the FaceGeom.
      //
      // It RIDES THE C/F SCHEME (the design note pins the stencil — "cfAppendStencil ...
      // scheme = quadratic" §5.1 — but not this gating): with CfScheme::standard the whole
      // projection family takes the raw coarse value at a 2:1 face, and a quadratic tangential
      // sample in the advection alone would be a closure the rest of the operator does not
      // share. It is also what keeps the distributed probe story exact — probeCfScheme() and
      // probeSeamLayer() run under the same condition, so every coordinate this builder asks
      // for is already in the halo registry.
      uploadSeamRecon(buildSeamRecon(pres_, regularOk, fluidOk, CfScheme::quadratic));
    } else {
      cfGhostCols_ = 0;
      cfMom_ = CfCsrDev{};
      cfDiv_ = CfCompCsrDev{};
      for (int a = 0; a < 3; ++a)
        cfGrad_[static_cast<std::size_t>(a)] = CfCsrDev{};
      cfUfVel_ = CfCompCsrDev{};
      clearSeamRecon();
    }
    profPhase("cf overlays");
    if (ghostProj_) {
      // Closure overlay (pre-built in the mode-resolve step above) + the coupled-subspace mask
      // for the BiCGStab projection. Sampled mode uploads the sample-slot overlay instead (the
      // classic gpOv_ stays empty — every delta site calls both, empties no-op).
      gpOv_ = ghostSampled_ ? GhostOverlayDev{} : uploadGhostOverlay(hov);
      gpOvS_ = ghostSampled_ ? uploadGhostOverlaySampled(hovS) : GhostOverlaySampledDev{};
      const GhostOverlay& rows = ghostSampled_ ? hovS.base : hov;
      std::vector<double> mc(static_cast<std::size_t>(n), 0.0);
      for (Index i = 0; i < n; ++i)
        mc[static_cast<std::size_t>(i)] =
            (mom_.isFluid(i) && !(!gpPocket_.empty() && gpPocket_[static_cast<std::size_t>(i)]))
                ? 1.0
                : 0.0;
      for (Index r = 0; r < rows.n; ++r)
        if (!rows.coupled[static_cast<std::size_t>(r)])
          mc[static_cast<std::size_t>(rows.cell[static_cast<std::size_t>(r)])] = 0.0;
      maskC_ = toDevice(mc, "gp_maskc");
      // Sampled mode: momentum ξ-row seam correction (host CSR folds 1/rscale for the oracle's
      // pre-rscale source; the device adds POST-rscale via cfApply, so re-fold rscale per row)
      // + the CSR directional-gradient overlay (cascade over sample functionals; classic
      // same-level fallback on cut cells without a row — the oracle's gradP split).
      if (ghostSampled_) {
        CfCsr msd = buildMomSeamDelta(hovS, *t_, pres_, mom_, sdfFn, origin_, rho_ / dt_, mu_,
                                      /*nSkipped=*/nullptr, shiftD_);
        for (Index i = 0; i < n; ++i) {
          const double rs = mom_.rhsScale(i);
          for (Index k = msd.start[static_cast<std::size_t>(i)];
               k < msd.start[static_cast<std::size_t>(i) + 1]; ++k)
            msd.coef[static_cast<std::size_t>(k)] *= rs;
        }
        gpsMomDelta_ = uploadCfCsr(msd, "gps_momdelta");
        auto classicOk = [&](Index j, Index i) {
          return j >= 0 && mom_.isFluid(j) && pres_.levelOf(j) == pres_.levelOf(i) &&
                 !(!gpPocket_.empty() && j < n && gpPocket_[static_cast<std::size_t>(j)]);
        };
        gcS_ = uploadGhostGradCsr(buildSampledGradCsr(
            hovS, *t_, pres_, [&](Index i) { return mom_.isCut(i); }, classicOk));
      } else {
        gpsMomDelta_ = CfCsrDev{};
        gcS_ = GhostGradCsrDev{};
      }
      // Krylov scratch carries the ghost tail (matvec inputs) in distributed mode.
      auto mk = [&](const char* l) { return View<double>(l, static_cast<std::size_t>(nExt_)); };
      gpr_ = mk("gp_r");
      gprh_ = mk("gp_rhat");
      gpp_ = mk("gp_p");
      gpph_ = mk("gp_phat");
      gpv_ = mk("gp_v");
      gps_ = mk("gp_s");
      gpsh_ = mk("gp_shat");
      gpt_ = mk("gp_t");
    }
 
    // Device state. Fields whose ghost entries are READ through the CSRs / overlays (u, p, φ)
    // — plus div (its View is deep_copied into the nExt-sized distributed MG rhs) and the u
    // snapshots (full-extent deep_copies of u) — carry the ghost tail [n, nExt); nExt == n
    // single-rank, so those allocations are unchanged there.
    for (int c = 0; c < 3; ++c) {
      u_[c] = View<double>("df_u", static_cast<std::size_t>(nExt_));
      gx_[c] = View<double>("df_g", static_cast<std::size_t>(n));
      Kokkos::deep_copy(u_[c], 0.0);
    }
    p_ = View<double>("df_p", static_cast<std::size_t>(nExt_));
    phi_ = View<double>("df_phi", static_cast<std::size_t>(nExt_));
    div_ = View<double>("df_div", static_cast<std::size_t>(nExt_));
    uf_ = View<double>("df_uf",
                       geom_.nbr.extent(0));  // ABC divergence-free face field (per CSR face)
    faceFieldBuilt_ = false;
    bmom_ = View<double>("df_bmom", static_cast<std::size_t>(n));
    Kokkos::deep_copy(p_, 0.0);
    // Implicit-FOU advection state. The momentum operator + its velocity-MG are rebuilt each
    // step from the *full* operator (viscous + FOU) so the MG is advection-aware (the viscous-
    // only MG diverges on the advection operator at cut cells); the FOU is baked into the CSR
    // (momOp_.hasAdv stays false). defc holds the device-computed explicit ρ(SOU−FOU)
    // deferred correction. uadvHost_ caches u^n on the host for the per-step operator rebuild.
    for (int c = 0; c < 3; ++c) {
      defc_[c] = View<double>("df_defc", static_cast<std::size_t>(n));
      Kokkos::deep_copy(defc_[c], 0.0);
      // uⁿ snapshot for the backward-Euler mass term (frozen across Picard outer iters) + the
      // previous outer iterate for the outer-loop convergence test.
      u0_[c] = View<double>("df_u0", static_cast<std::size_t>(nExt_));
      uprev_[c] = View<double>("df_uprev", static_cast<std::size_t>(nExt_));
    }
    // Device-resident implicit-FOU advection: the FOU operator (advDiag + per-face advCoef over the
    // face-geometry CSR) is rebuilt on device each step from uⁿ and added to the static Stokes
    // operator in the matvec (no host round-trip). The static operator + Galerkin velocity-MG are
    // built once (above); the advection is a perturbation the static MG still preconditions.
    advDiag_ = View<double>("df_advdiag", static_cast<std::size_t>(n));
    advCoef_ = View<double>("df_advcoef", geom_.nbr.extent(0));
    momOp_.advStart = geom_.start;
    momOp_.advNbr = geom_.nbr;
    momOp_.advDiag = advDiag_;
    momOp_.advCoef = advCoef_;
    momOp_.hasAdv = false;  // set per-step when advection is on
    momSolver_.setJacobi(2, 0.7);
    // Generic MG preconditioner: dispatch the chosen hierarchy's V-cycle (z = M⁻¹ r), decoupling
    // the solver from the MG type so Galerkin and staircase are interchangeable. Distributed:
    // the rank-local Galerkin hierarchy built above (staircase unsupported there ⇒ Jacobi).
    if (momMGon_ && !(dist_ && useStaircaseMG_)) {
      if (useStaircaseMG_)
        momSolver_.setPreconditioner(
            [this](View<const double> r, View<double> z) { runMgVcycle(velMG_, r, z); });
      else
        momSolver_.setPreconditioner(
            [this](View<const double> r, View<double> z) { runMgVcycle(momMG_, r, z); });
    }
    pcg_.setVcycle(2, 2, 60, 0.8);
    pcg_.setSingular(true);
    n_ = n;
    ready_ = true;
    profPhase("overlays upload + state alloc");
  }
 
  // ---- M0 (docs/amr_march_perf_and_distributed_plan.md): the per-phase STEP profiler -----------
  // `PECLET_AMR_PROFILE_STEP=1` accumulates FENCED per-phase timings over a window of steps
  // (`PECLET_AMR_PROFILE_STEP_WINDOW`, default 50) and prints ms/step and µs/leaf per phase next
  // to the iteration counters — the instrument Phase M's attribution matrix reads. The four
  // pressure-solve rows are NESTED inside `pressure solve` (the residual is Krylov vector work and
  // the reductions), so they are printed indented and must not be added to the total.
  //
  // OFF (the default) every hook is one bool test: no fence, no clock, no change to the kernel
  // sequence — so a profiler-off march is bitwise what it was before the profiler existed.
  enum StepPhase {
    SP_GLUE,
    SP_GRAD,
    SP_CF,
    SP_OVERLAY,
    SP_ADV,
    SP_MOMRHS,
    SP_MOMSOLVE,
    SP_DIV,
    SP_PRES,
    SP_PRES_MV,
    SP_PRES_MVOV,
    SP_PRES_PROJ,
    SP_PRES_PC,
    SP_FINISH,
    SP_N
  };
  using SpClock = std::chrono::steady_clock;
  SpClock::time_point spMark() const {
    if (!stepProf_)
      return SpClock::time_point{};
    Kokkos::fence();
    return SpClock::now();
  }
  void spAdd(int ph, SpClock::time_point t0) {
    if (!stepProf_)
      return;
    Kokkos::fence();
    spAcc_[ph] += std::chrono::duration<double, std::milli>(SpClock::now() - t0).count();
  }
  void spEndStep() {
    if (!stepProf_)
      return;
    ++spSteps_;
    spMom_ += lastMomIters_;
    spPres_ += lastPresIters_;
    spOuter_ += lastOuterIters_;
    if (!spHeader_) {  // the static hierarchy shape, on the first profiled step (H-mg's
                       // denominator, and the one line a contention-exposed box can still trust)
      const std::vector<Index> lad = pressureMgLevels();
      std::fprintf(stderr, "[step-prof] leaves %lld | pressure MG levels %zu:", (long long)n_,
                   lad.size());
      for (Index n : lad)
        std::fprintf(stderr, " %lld", (long long)n);
      std::fprintf(stderr, " | bottom %s\n", pressureMgBottom().c_str());
      spHeader_ = true;
    }
    if (spSteps_ < spWindow_)
      return;
    static const char* kName[SP_N] = {
        "glue (copies + syncs)", "grad3 (binary)",    "cf delta kernels",   "overlay delta kernels",
        "advection build",       "momentum rhs",      "momentum solve",     "divergence",
        "pressure solve",        "  . binary matvec", "  . overlay matvec", "  . projection",
        "  . MG preconditioner", "finish projection"};
    const double w = static_cast<double>(spSteps_);
    const double leaves = static_cast<double>(n_);
    double tot = 0.0;
    for (int k = 0; k < SP_N; ++k)
      if (k < SP_PRES_MV || k > SP_PRES_PC)
        tot += spAcc_[k];
    std::fprintf(stderr,
                 "[step-prof] %d steps | %.3f ms/step (%.3f us/leaf) | mom %.1f it, "
                 "pres %.1f it, outer %.2f\n",
                 spSteps_, tot / w, 1e3 * tot / w / leaves, spMom_ / w, spPres_ / w, spOuter_ / w);
    for (int k = 0; k < SP_N; ++k) {
      const bool nested = (k >= SP_PRES_MV && k <= SP_PRES_PC);
      std::fprintf(stderr, "[step-prof]   %-22s %8.3f ms/step  %8.4f us/leaf  %5.1f%%%s\n",
                   kName[k], spAcc_[k] / w, 1e3 * spAcc_[k] / w / leaves,
                   100.0 * spAcc_[k] / (tot > 0 ? tot : 1.0), nested ? "  (nested)" : "");
    }
    for (int k = 0; k < SP_N; ++k)
      spAcc_[k] = 0.0;
    spSteps_ = 0;
    spMom_ = spPres_ = spOuter_ = 0.0;
  }
 
  void requireOperator(const char* what) const {
    if (!ready_)
      throw std::runtime_error(std::string("amr::AmrFlow::") + what +
                               ": no cut-cell operator -- call setSolid() first (Python: "
                               "Flow.set_solid / set_solid_spheres; after begin_adapt, "
                               "finish_adapt)");
  }
 
  void step(int momIters = 100, int presIters = 60) {
    requireOperator("step");
    const Index n = n_;
    const double idiag = rho_ / dt_;
    lastMomIters_ = 0;
    lastOuterIters_ = 1;
    const auto spT0 = spMark();
    // Freeze the time-level uⁿ for the backward-Euler mass term + warm start; it stays anchored
    // across the Picard outer iterations (only the advecting velocity re-lags). For outerIters_==1
    // this is just a copy of uⁿ ⇒ bit-identical to the single lagged step.
    for (int c = 0; c < 3; ++c)
      Kokkos::deep_copy(u0_[c], View<const double>(u_[c]));
    // −∇p^n is constant across the outer iterations (pressure is projected once, after the loop,
    // like flow's single per-step projection) ⇒ hoist it out.
    syncScalar(p_);  // ghost tail of pⁿ before the gradient reads (no-op single-rank)
    spAdd(SP_GLUE, spT0);
    auto spT = spMark();
    grad3(geom_, View<const double>(p_), gx_[0], gx_[1], gx_[2]);
    spAdd(SP_GRAD, spT);
    spT = spMark();
    for (int a = 0; a < 3; ++a)  // 2nd-order C/F face gradients (level-boundary rows)
      cfApply(cfGrad_[static_cast<std::size_t>(a)], View<const double>(p_), gx_[a]);
    spAdd(SP_CF, spT);
    spT = spMark();
    applyGhostGrad(gc_, View<const double>(p_), gx_[0], gx_[1], gx_[2]);
    applyGhostGradCsr(gcS_, View<const double>(p_), gx_[0], gx_[1], gx_[2]);
    spAdd(SP_OVERLAY, spT);
    // Picard outer loop over the lagged advection only (the momentum nonlinearity); for
    // outerIters_==1 this is the single lagged predictor, then one projection — bit-identical to
    // before.
    for (int outer = 0; outer < outerIters_; ++outer) {
      // --- advection lagged to the *current* predictor iterate (uⁿ on the first pass):
      // implicit-FOU operator + explicit ρ(SOU−FOU) deferred correction. The matvec stays linear
      // during each solve (the advecting velocity is frozen in advDiag_/advCoef_). With advection
      // OFF the operator and RHS are identical every pass, so a second pass reproduces the first ⇒
      // instant early-stop. ---
      spT = spMark();
      if (advect_) {
        momOp_.hasAdv = implicitFou_;
        syncVel();  // ghost tails of the (lagged) advecting velocity for buildFou/deferredSou
        // Advect with the divergence-free face field uf (from the previous projection); fall back
        // to ½(u_i+u_j) before the first projection has built it. The implicit FOU and the explicit
        // SOU/FOU deferred correction use the SAME velocity ⇒ the FOU cancels at steady state
        // (host-parity).
        const View<const double> ufv(uf_);
        const bool useUf = faceFieldBuilt_ && ufAdvect_;
        if (implicitFou_)
          buildFou(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
                   View<const double>(u_[2]), rho_, View<const double>(rscale_), advDiag_, advCoef_,
                   ufv, useUf);
        for (int c = 0; c < 3; ++c) {
          if (implicitFou_)
            deferredSou(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
                        View<const double>(u_[2]), c, rho_, advScheme_, defc_[c], ufv, useUf,
                        seamRecon_);
          else
            advectExplicit(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
                           View<const double>(u_[2]), c, rho_, advScheme_, defc_[c], ufv, useUf,
                           seamRecon_);
        }
        // The staircase MG's fine level mirrors the sharp operator; refresh it so it picks up the
        // current advection state (hasAdv). (The Galerkin MG is the static viscous operator.)
        if (momMGon_ && useStaircaseMG_)
          velMG_.setFineOp(momOp_);
      }
      spAdd(SP_ADV, spT);
      // --- predictor: incremental BE viscous (+ implicit-FOU) solve per component, RHS carries
      // −∇p^n and −ρ(SOU−FOU); the mass term is anchored at uⁿ (u0_), the solve warm-starts at the
      // current iterate. ---
      // D2: the momentum ξ-seam delta is a CSR over the overlay's SAMPLE functionals, so on a
      // mixed-level band it reads ghost slots — and with advection off nothing else has refreshed
      // uⁿ's ghost tail since the previous step's projection. The C/F momentum delta (cfMom_) is
      // in exactly the same position: it reads u at the coarse cell's tangential neighbours,
      // which are ghost slots at a block seam. (Zero ghosts single-rank / np=1, so this is a
      // no-op exchange there and the path stays bit-identical.)
      if (dist_ && (ghostSampled_ || cfScheme_ != CfScheme::standard))
        syncVel();
      for (int c = 0; c < 3; ++c) {
        spT = spMark();
        momRhs(View<const double>(u0_[c]), View<const double>(gx_[c]), View<const double>(defc_[c]),
               View<const double>(rscale_), View<const char>(fluid_), idiag, f_[c], bmom_, n);
        spAdd(SP_MOMRHS, spT);
        spT = spMark();
        // C/F-scheme deferred correction on the velocity diffusion: b += μ(∇²_scheme − ∇²_std)
        // of the lagged component (regular rows only, rscale = 1 there).
        cfApply(cfMom_, View<const double>(u_[c]), bmom_);
        // Momentum ξ-row seam correction (sampled mode; rscale re-folded at upload).
        cfApply(gpsMomDelta_, View<const double>(u_[c]), bmom_);
        spAdd(SP_CF, spT);
        spT = spMark();
        // P4 (opt-in): the velocity-MG used as the *solver* — MG-preconditioned defect correction,
        // no Krylov (the flow RB-GS/velocity-MG mirror; cannot break down on the non-symmetric
        // operator) — vs the default MG-preconditioned BiCGStab. Both reach the same solution (the
        // matvec is the exact operator); the choice only trades robustness for convergence rate.
        lastMomIters_ +=
            (momMGSolver_ ? momSolver_.solveDefectCorrection(
                                momOp_, u_[c], View<const double>(bmom_), momIters, momTol_)
                          : momSolver_.solveBiCGStab(momOp_, u_[c], View<const double>(bmom_),
                                                     momIters, momTol_))
                .iters;
        spAdd(SP_MOMSOLVE, spT);
      }
      // Outer-loop convergence on the predictor velocity (skipped for the default outerIters_==1,
      // so that path is untouched). With advection off the second iterate equals the first ⇒ stops
      // at 2.
      if (outerIters_ > 1) {
        lastOuterIters_ = outer + 1;
        if (outer > 0) {
          double dmax = 0.0;
          for (int c = 0; c < 3; ++c)
            dmax = std::max(
                dmax, maxAbsDiff(View<const double>(u_[c]), View<const double>(uprev_[c]), n));
          if (dmax < outerTol_)
            break;
        }
        for (int c = 0; c < 3; ++c)
          Kokkos::deep_copy(uprev_[c], View<const double>(u_[c]));
      }
    }
    project(presIters);  // single pressure projection per step (flow structure)
    spEndStep();
  }
 
  void project(int presIters = 60) {
    requireOperator("project");
    const Index n = n_;
    auto spT = spMark();
    syncVel();  // ghost tails of u* for the divergence (+ the ghost-closed overlay delta)
    divergence(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
               View<const double>(u_[2]), div_);
    spAdd(SP_DIV, spT);
    spT = spMark();
    cfApplyComp(cfDiv_, View<const double>(u_[0]), View<const double>(u_[1]),
                View<const double>(u_[2]), div_);  // 2nd-order C/F face averages (setCfScheme)
    spAdd(SP_CF, spT);
    spT = spMark();
    if (ghostProj_) {  // ghost-closed constraint: binary div (geom_ carries binary α) + overlay
      ghostDivergDelta(gpOv_, View<const double>(u_[0]), View<const double>(u_[1]),
                       View<const double>(u_[2]), div_);
      ghostDivergDeltaSampled(gpOvS_, View<const double>(u_[0]), View<const double>(u_[1]),
                              View<const double>(u_[2]), div_);
    }
    spAdd(SP_OVERLAY, spT);
    Kokkos::deep_copy(phi_, 0.0);
    if (ghostProj_) {
      spT = spMark();
      lastPresIters_ = solveGhostBiCGStab(phi_, View<const double>(div_), presIters, presTol_);
      spAdd(SP_PRES, spT);
      spT = spMark();
      finishProjection(n);
      spAdd(SP_FINISH, spT);
      return;
    }
    // Two selectable pressure drivers, like flow's CutcellMG: MG-PCG (default, presPCG_) and the
    // bounded stationary V-cycle (setPressurePCG(false)). MG-PCG covers ADVECTION too — the
    // historic exclusion ("transient near-nullspace issue") was characterised 2026-08-19 and was
    // NOT a property of the operator (which is geometry-only, advection-independent, SPD in the
    // volume-weighted inner product): the aperture RHS is INCOMPATIBLE by a fluid-mean component
    // (div_ is zeroed at solid-centered cells whose faces are partially open, breaking the
    // telescoping that would make Σ V·div = 0 over the operator's DOF set; the defect grows with
    // the developed flow, ~3e-3 relative at steady state on the Z&H sphere). The un-deflated
    // V-cycle's residual therefore STALLS at exactly |mean|·sqrt(V_fluid) — the old "60 bounded
    // cycles" was a stagnation cap, not a convergence count. The PCG's per-iteration fluid-range
    // projection (maskSolid + volume-weighted fluid-mean removal) deflates exactly that component,
    // so CG is valid and healthy: measured flat 15–17 iters/step (tol 1e-10) across the whole
    // impulsive N=32 transient, steady K identical to the V-cycle path to 4+ digits.
    // Debug knob (env, default-off): PECLET_AMR_PRES_DEBUG=1 — per-cycle/-solve residual +
    // RHS-compatibility trace to stderr (the characterisation instrumentation, kept).
    const bool dbg = amrEnvFlag("PECLET_AMR_PRES_DEBUG");
    if (dbg && !dist_) {
      // RHS compatibility: the operator's left null vector is the constant over fluid cells in the
      // volume-weighted inner product, so a solvable RHS needs Σ V_i·div_i ≈ 0 over fluid cells.
      const FvOp& op0 = presMG_.op(0);
      auto st = op0.faceStart;
      auto w = op0.faceW;
      auto bcD = op0.bcDiag;
      auto iv = op0.invVol;
      auto dv = div_;
      double su = 0.0, sv = 0.0, sn = 0.0;
      Kokkos::parallel_reduce(
          "amr::dbg_compat", n,
          KOKKOS_LAMBDA(const Index i, double& a, double& b, double& c) {
            double d = bcD(i);
            for (Index k = st(i); k < st(i + 1); ++k)
              d += w(k);
            const double m = (d > 1e-30) ? 1.0 : 0.0;
            a += m * dv(i) / iv(i);
            b += m / iv(i);
            c += m * dv(i) * dv(i) / iv(i);
          },
          su, sv, sn);
      std::fprintf(stderr, "[amr pres] rhs fluid-mean=%.3e |rhs|_D=%.3e (rel mean %.3e)\n", su / sv,
                   std::sqrt(sn), (su / sv) / (std::sqrt(sn) + 1e-300));
      if (!presDbgSpdDone_) {
        // One-shot direct SPD probe of the assembled pressure operator (advection cannot enter
        // the assembly — this measures it): symmetry <y,Lx>_D vs <x,Ly>_D on deterministic
        // pseudo-random vectors, and the Rayleigh quotient sign (L negative-semidefinite).
        presDbgSpdDone_ = true;
        View<double> xr("dbg_x", static_cast<std::size_t>(n)),
            yr("dbg_y", static_cast<std::size_t>(n)), Ax("dbg_Ax", static_cast<std::size_t>(n)),
            Ay("dbg_Ay", static_cast<std::size_t>(n));
        Kokkos::parallel_for(
            "amr::dbg_fill", n, KOKKOS_LAMBDA(const Index i) {
              xr(i) = std::sin(0.7 * static_cast<double>(i) + 0.3);
              yr(i) = std::cos(1.3 * static_cast<double>(i) + 1.1);
            });
        applyFv(op0, View<const double>(xr), Ax);
        applyFv(op0, View<const double>(yr), Ay);
        auto dotD = [&](View<const double> a, View<const double> b) {
          double s = 0.0;
          Kokkos::parallel_reduce(
              "amr::dbg_dot", n,
              KOKKOS_LAMBDA(const Index i, double& acc) { acc += a(i) * b(i) / iv(i); }, s);
          return s;
        };
        const double yLx = dotD(View<const double>(yr), View<const double>(Ax));
        const double xLy = dotD(View<const double>(xr), View<const double>(Ay));
        const double xLx = dotD(View<const double>(xr), View<const double>(Ax));
        const double yLy = dotD(View<const double>(yr), View<const double>(Ay));
        std::fprintf(stderr,
                     "[amr pres] SPD probe: <y,Lx>_D=%.15e <x,Ly>_D=%.15e rel asym=%.2e; "
                     "Rayleigh <x,Lx>_D=%.3e <y,Ly>_D=%.3e (must be <=0)\n",
                     yLx, xLy, std::fabs(yLx - xLy) / (std::fabs(yLx) + 1e-300), xLx, yLy);
      }
    }
    spT = spMark();
    if (presPCG_) {
      const auto R = dist_
                         ? pcg_.solve(presMGD_, phi_, View<const double>(div_), presIters, presTol_)
                         : pcg_.solve(presMG_, phi_, View<const double>(div_), presIters, presTol_);
      lastPresIters_ = R.iters;
      if (dbg)
        std::fprintf(stderr, "[amr pres] pcg iters=%d res0=%.3e res=%.3e rel=%.3e\n", R.iters,
                     R.res0, R.res, R.res0 > 0 ? R.res / R.res0 : 0.0);
    } else if (dist_) {
      Kokkos::deep_copy(presMGD_.b(0), div_);
      Kokkos::deep_copy(presMGD_.x(0), 0.0);
      for (int it = 0; it < presIters; ++it)
        presMGD_.vcycle(2, 2, 60, 0.8);
      Kokkos::deep_copy(phi_, presMGD_.x(0));
      lastPresIters_ = presIters;
    } else {
      Kokkos::deep_copy(presMG_.b(0), div_);
      Kokkos::deep_copy(presMG_.x(0), 0.0);
      View<double> rdbg;
      if (dbg)
        rdbg = View<double>("pres_dbg_res", static_cast<std::size_t>(n));
      for (int it = 0; it < presIters; ++it) {
        presMG_.vcycle(2, 2, 60, 0.8);
        if (dbg) {
          residualFv(presMG_.op(0), View<const double>(presMG_.x(0)),
                     View<const double>(presMG_.b(0)), rdbg);
          const double rn =
              std::sqrt(dotPlain(View<const double>(rdbg), View<const double>(rdbg), n));
          std::fprintf(stderr, "[amr pres] vcycle %2d |r|=%.6e\n", it + 1, rn);
        }
      }
      Kokkos::deep_copy(phi_, presMG_.x(0));
      lastPresIters_ = presIters;
    }
    spAdd(SP_PRES, spT);
    spT = spMark();
    finishProjection(n);
    spAdd(SP_FINISH, spT);
  }
 
  void finishProjection(Index n) {
    // u*'s ghost tail is current (project() synced it and u is untouched since); φ's is not.
    syncScalar(phi_);
    buildFaceField(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
                   View<const double>(u_[2]), View<const double>(phi_), uf_);
    // 2nd-order C/F face values (setCfScheme), VELOCITY part only: the distance-weighted
    // {u_fine, coarse*} interpolation on the 2:1 sub-faces. Its divergence IS the constraint's
    // Δ_cfDiv — measured identical to 8.7e-17 — so the projection has already balanced it and uf
    // stays divergence-free with it in.
    cfApplyComp(cfUfVel_, View<const double>(u_[0]), View<const double>(u_[1]),
                View<const double>(u_[2]), uf_);
    // The φ part is NOT applied, and must not be. It substitutes the quadratic coarse* into uf's
    // FACE GRADIENT, but the pressure matrix is L = D_std·G_std (a recorded decision: at the fixed
    // point φ→0, so the matrix's C/F order cannot move the steady solution). A gradient the solve
    // never inverted cannot be balanced by it, so adding it breaks exactly the property uf exists
    // to have: Σ_faces α·A·dir·uf = 0, i.e. that advecting a scalar with uf conserves it and
    // leaves a uniform field uniform.
    //
    // Measured on a graded sphere at cf=1 (N=32, lmax=2, band=3, 25 steps), decomposing
    //     d = (rhs − Lφ) − Δ_cfDiv + D(Δvel) + D(Δφ)      [identity holds to 2.3e-15]
    //     solve residual ‖rhs − Lφ‖      1.34e-12
    //     ‖D(Δvel) − Δ_cfDiv‖            8.71e-17   ← velocity part: balanced, keep it
    //     ‖D(Δφ)‖                        6.675e-03  ← THE ENTIRE flux imbalance
    // It is proportional to φ (0.240 at step 1, 6.7e-3 at step 25, 3.1e-4 at step 60) so it
    // vanishes at steady state and never showed up in a steady drag gate — but with advection on
    // it sat at 4.0e-2, the same order as the cell field's O(h²) divergence, i.e. uf had stopped
    // being a conservative flux for the whole of the transient.
    //
    // If a quadratic face gradient in uf is wanted, the pressure matrix has to invert the same
    // operator (L = D·G_quad). Mixing them is what this removes. The φ CSR is no longer built at
    // all (docs/amr_cf_flux_gate.md §6.7) — a CSR that must never be applied is a trap; §6.8
    // records the deferred-correction path that would buy a quadratic uf face gradient legally.
    faceFieldBuilt_ = true;
    grad3(geom_, View<const double>(phi_), gx_[0], gx_[1], gx_[2]);
    for (int a = 0; a < 3; ++a)  // 2nd-order C/F face gradients (level-boundary rows)
      cfApply(cfGrad_[static_cast<std::size_t>(a)], View<const double>(phi_), gx_[a]);
    applyGhostGrad(gc_, View<const double>(phi_), gx_[0], gx_[1], gx_[2]);
    applyGhostGradCsr(gcS_, View<const double>(phi_), gx_[0], gx_[1], gx_[2]);
    for (int c = 0; c < 3; ++c)
      correct(u_[c], View<const double>(gx_[c]), View<const char>(fluid_), n);
    presUpdate(p_, View<const double>(phi_), View<const double>(div_), View<const char>(fluid_),
               rho_ / dt_, mu_, n);
  }
 
  // ---- adaptivity during a run (ladder step 5) -------------------------------------------------
  // The octree is borrowed by pointer and mutated EXTERNALLY (refine/coarsen/balance/adapt on the
  // same object). beginAdapt snapshots the current topology + fields; after the mutation,
  // finishAdapt conservatively remaps u and the accumulated rotational p onto the new mesh
  // (minmod-limited linear transferField) and rebuilds every solver structure via setSolid. The
  // caller must keep the cut band at the finest level on the new mesh (refineToSdf the geometry
  // band again after a solution-driven adapt): the ghost overlay build throws / auto-falls-back
  // exactly as in setSolid. uf restarts from the ½-average fallback for one step.
 
  void beginAdapt() {
    requireOperator("beginAdapt");
    ready_ = false;  // no operator until finishAdapt rebuilds it on the new topology
    adaptOldT_ = std::make_unique<Octree>(*t_);
    for (int c = 0; c < 3; ++c)
      adaptU_[static_cast<std::size_t>(c)] = velocity(c);
    adaptP_ = pressure();
    if (dist_) {
      // Halo-completed prolongation gradients on the OLD mesh (while dist_ still holds it):
      // the block-local transfer stencil zeroes gradients at interior block boundaries,
      // which would make the remap np-dependent there (measured ~5% field divergence).
      for (int c = 0; c < 3; ++c)
        adaptGradU_[static_cast<std::size_t>(c)] =
            transferGradients(*dist_, adaptU_[static_cast<std::size_t>(c)]);
      adaptGradP_ = transferGradients(*dist_, adaptP_);
    }
  }
 
  template <class SdfFn>
  void finishAdapt(SdfFn&& sdfFn) {
    if (!adaptOldT_)
      throw std::runtime_error("amr::AmrFlow::finishAdapt called without beginAdapt");
    std::array<std::vector<double>, 3> nu;
    for (int c = 0; c < 3; ++c)
      nu[static_cast<std::size_t>(c)] =
          transferField(*adaptOldT_, adaptU_[static_cast<std::size_t>(c)], *t_, /*linear=*/true,
                        dist_ ? &adaptGradU_[static_cast<std::size_t>(c)] : nullptr);
    std::vector<double> np =
        transferField(*adaptOldT_, adaptP_, *t_, /*linear=*/true, dist_ ? &adaptGradP_ : nullptr);
    setSolid(sdfFn);  // full operator/overlay rebuild on the new topology (zeroes the fields)
    for (int c = 0; c < 3; ++c) {
      setVelocity(c, nu[static_cast<std::size_t>(c)]);
      zeroSolid(u_[c]);  // cells that became solid on the new mesh hold 0 (no-slip state)
    }
    setPressure(np);
    zeroSolid(p_);  // solid p is pinned/decoupled
    adaptOldT_.reset();
    for (int c = 0; c < 3; ++c)
      adaptU_[static_cast<std::size_t>(c)].clear();
    adaptP_.clear();
  }
 
  template <class SdfFn>
  void rebalanceMpi(SdfFn&& sdfFn) {
    if (!dist_)
      throw std::runtime_error("amr::AmrFlow::rebalanceMpi requires initMpi");
    std::vector<std::vector<double>> cols(4);
    for (int c = 0; c < 3; ++c)
      cols[static_cast<std::size_t>(c)] = velocity(c);
    cols[3] = pressure();
    dist_->rebalance(cols);  // t_ still points at dist_->local(), now the new block's octree
    setSolid(sdfFn);         // full rebuild (registry, halos, operators) on the new block
    for (int c = 0; c < 3; ++c) {
      setVelocity(c, cols[static_cast<std::size_t>(c)]);
      zeroSolid(u_[c]);
    }
    setPressure(cols[3]);
    zeroSolid(p_);
  }
 
  void setPressure(const std::vector<double>& h) {
    auto m = Kokkos::create_mirror_view(p_);
    for (Index i = 0; i < n_; ++i)
      m(i) = h[static_cast<std::size_t>(i)];
    Kokkos::deep_copy(p_, m);
  }
 
  void zeroSolid(View<double> v) {
    auto fl = fluid_;
    Kokkos::parallel_for(
        "amr::flow_zerosolid", n_, KOKKOS_LAMBDA(const Index i) {
          if (!fl(i))
            v(i) = 0.0;
        });
  }
 
  std::vector<double> debugSou(int comp) {
    syncVel();
    View<double> s("dbg_sou", static_cast<std::size_t>(n_));
    advectExplicit(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
                   View<const double>(u_[2]), comp, 1.0, advScheme_, s, View<const double>(uf_),
                   false, seamRecon_);
    std::vector<double> h(static_cast<std::size_t>(n_));
    auto m = Kokkos::create_mirror_view(s);
    Kokkos::deep_copy(m, s);
    for (Index i = 0; i < n_; ++i)
      h[static_cast<std::size_t>(i)] = m(i);
    return h;
  }
  template <class MG>
  void runMgVcycle(MG& mg, View<const double> r, View<double> z) {
    // Element copies over the LOCAL rows (not deep_copy): in distributed mode the Krylov
    // scratch carries a ghost tail while the (rank-local) MG levels are local-sized; the
    // solver refreshes z's ghost tail before the next matvec. Same values single-rank.
    const Index nl = mg.numLeaves(0);
    {
      auto b0 = mg.b(0);
      Kokkos::parallel_for("amr::mgv_b", nl, KOKKOS_LAMBDA(const Index i) { b0(i) = r(i); });
    }
    Kokkos::deep_copy(mg.x(0), 0.0);
    mg.vcycle(mgVcPre_, mgVcPre_, mgVcBottom_, 0.7);
    {
      auto x0 = mg.x(0);
      Kokkos::parallel_for("amr::mgv_z", nl, KOKKOS_LAMBDA(const Index i) { z(i) = x0(i); });
    }
  }
  static double maxAbsDiff(View<const double> a, View<const double> b, Index n) {
    double m = 0.0;
    Kokkos::parallel_reduce(
        "amr::flow_maxdiff", n,
        KOKKOS_LAMBDA(const Index i, double& lm) {
          double d = a(i) - b(i);
          if (d < 0.0)
            d = -d;
          if (d > lm)
            lm = d;
        },
        Kokkos::Max<double>(m));
    return m;
  }
  void copyToHost(const View<double>& d, std::vector<double>& h) const {
    auto m = Kokkos::create_mirror_view(d);
    Kokkos::deep_copy(m, d);
    for (Index i = 0; i < n_; ++i)
      h[static_cast<std::size_t>(i)] = m(i);
  }
  void setVelocity(int c, const std::vector<double>& h) {
    auto m = Kokkos::create_mirror_view(u_[c]);
    for (Index i = 0; i < n_; ++i)
      m(i) = h[static_cast<std::size_t>(i)];
    Kokkos::deep_copy(u_[c], m);
  }
  std::vector<double> velocity(int c) const { return localVector(u_[c]); }
  std::vector<double> pressure() const { return localVector(p_); }
 
  std::vector<double> localVector(const View<double>& v) const {
    if (v.extent(0) == static_cast<std::size_t>(n_))
      return peclet::core::toVector(v);
    View<double> packed(Kokkos::view_alloc("amr::local_packed", Kokkos::WithoutInitializing),
                        static_cast<std::size_t>(n_));
    Kokkos::parallel_for("amr::pack_local", n_, KOKKOS_LAMBDA(const Index i) { packed(i) = v(i); });
    return peclet::core::toVector(packed);
  }
 
  std::vector<double> velocities() const {
    View<double> packed(Kokkos::view_alloc("amr::vel_packed", Kokkos::WithoutInitializing),
                        static_cast<std::size_t>(n_) * 3);
    for (int c = 0; c < 3; ++c) {
      auto uc = u_[c];
      auto p = packed;
      const int cc = c;
      Kokkos::parallel_for(
          "amr::pack_vel", n_, KOKKOS_LAMBDA(const Index i) { p(i * 3 + cc) = uc(i); });
    }
    return peclet::core::toVector(packed);
  }
  double divNormL2() {
    syncVel();
    divergence(geom_, View<const double>(u_[0]), View<const double>(u_[1]),
               View<const double>(u_[2]), div_);
    cfApplyComp(cfDiv_, View<const double>(u_[0]), View<const double>(u_[1]),
                View<const double>(u_[2]), div_);
    if (ghostProj_) {
      ghostDivergDelta(gpOv_, View<const double>(u_[0]), View<const double>(u_[1]),
                       View<const double>(u_[2]), div_);
      ghostDivergDeltaSampled(gpOvS_, View<const double>(u_[0]), View<const double>(u_[1]),
                              View<const double>(u_[2]), div_);
    }
    return std::sqrt(allSum(dotPlain(View<const double>(div_), View<const double>(div_), n_)));
  }
  double divNormFace() {
    const double l = divFaceNorm(geom_, View<const double>(uf_));
    return std::sqrt(allSum(l * l));
  }
  std::vector<Index> pressureMgLevels() const {
    std::vector<Index> out;
    if (!dist_) {
      for (std::size_t L = 0; L < presMG_.numLevels(); ++L)
        out.push_back(presMG_.numLeaves(L));
      return out;
    }
    for (std::size_t L = 0; L < presMGD_.numLevels(); ++L)
      out.push_back(presMGD_.numLeaves(L));
    for (std::size_t L = 0; L < presMGD_.numStageLevels(); ++L)
      out.push_back(presMGD_.stageLeaves(L));  // §6.5: the stage's level 0 IS the moved level
    return out;
  }
  std::string pressureMgBottom() const {
    return dist_ ? presMGD_.bottomName() : presMG_.bottomName();
  }
 
  void setPressureBottom(const std::string& kind) {
    using B = typename Multigrid<3, Bits>::Bottom;
    if (kind == "auto")
      presBottom_ = B::Auto;
    else if (kind == "smoother")
      presBottom_ = B::Smoother;
    else if (kind == "agglomerated")
      presBottom_ = B::Agglomerated;
    else
      throw std::runtime_error(
          "amr::AmrFlow::setPressureBottom: expected 'auto', 'smoother' or 'agglomerated'");
    if (dist_)
      presMGD_.setBottom(presBottom_);
    else
      presMG_.setBottom(presBottom_);
  }
 
  void setPressureBottomExtent(Index e) {
    if (e < 1)
      throw std::runtime_error("amr::AmrFlow::setPressureBottomExtent: must be >= 1");
    presBottomExtent_ = e;
  }
  Index pressureBottomExtent() const { return presBottomExtent_; }
 
  std::vector<double> faceField() const { return peclet::core::toVector(uf_); }
  struct FaceTopologyHost {
    std::vector<Index> start, nbr, upupI, upupJ;
    std::vector<int> axis, dir;
    std::vector<double> rawArea, dist, alpha;
    std::vector<Index> seam, sampI, sampJ, uuRecI, uuRecJ, recStart, recCell;
    std::vector<double> d1I, d1J, recDist, recW;
    std::vector<double> center;
  };
  FaceTopologyHost faceTopology() const {
    FaceTopologyHost t;
    t.start = peclet::core::toVector(geom_.start);
    t.nbr = peclet::core::toVector(geom_.nbr);
    t.upupI = peclet::core::toVector(geom_.upupI);
    t.upupJ = peclet::core::toVector(geom_.upupJ);
    t.axis = peclet::core::toVector(geom_.axis);
    t.dir = peclet::core::toVector(geom_.dir);
    t.rawArea = peclet::core::toVector(geom_.rawArea);
    t.dist = peclet::core::toVector(geom_.dist);
    t.alpha = peclet::core::toVector(geom_.alpha);
    t.seam = peclet::core::toVector(geom_.seam);
    t.sampI = peclet::core::toVector(geom_.sampI);
    t.sampJ = peclet::core::toVector(geom_.sampJ);
    t.uuRecI = peclet::core::toVector(geom_.uuRecI);
    t.uuRecJ = peclet::core::toVector(geom_.uuRecJ);
    t.recStart = peclet::core::toVector(geom_.recStart);
    t.recCell = peclet::core::toVector(geom_.recCell);
    t.d1I = peclet::core::toVector(geom_.d1I);
    t.d1J = peclet::core::toVector(geom_.d1J);
    t.recDist = peclet::core::toVector(geom_.recDist);
    t.recW = peclet::core::toVector(geom_.recW);
    t.center = slotCenters();
    return t;
  }
  std::vector<double> slotCenters() const {
    std::vector<double> c(static_cast<std::size_t>(nExt_) * 3, 0.0);
    for (Index s = 0; s < nExt_; ++s) {
      const std::array<long, 3> lo = pres_.loOf(s);
      const double w = static_cast<double>(1L << pres_.levelOf(s));
      for (int a = 0; a < 3; ++a)
        c[static_cast<std::size_t>(s) * 3 + static_cast<std::size_t>(a)] =
            origin_[a] + h0_[a] * (static_cast<double>(lo[a] + shiftD_[a]) + 0.5 * w);
    }
    return c;
  }
  Index numSeamSampleRecords() const { return seamSampleRecords_; }
  Index numSeamLayerRecords() const { return seamLayerRecords_; }
  Index numLeaves() const { return n_; }
  Index numGhostCells() const { return nExt_ - n_; }
  Index numCfGhostColumns() const { return cfGhostCols_; }
  Index numCfCutFaces() const { return cfCutFaces_; }
  bool isFluid(Index i) const { return mom_.isFluid(i); }
  int lastMomIters() const { return lastMomIters_; }
  int lastPresIters() const { return lastPresIters_; }
  int lastOuterIters() const { return lastOuterIters_; }
 
  // ---- distributed plumbing (docs/amr_distributed_flow.md, rung 4) ----------------------------
 
  void syncVel() {
    if (dist_)
      dhex_.exchange3(u_[0], u_[1], u_[2]);
  }
  void syncScalar(View<double> v) {
    if (dist_)
      dhex_.exchange(v);
  }
  double allSum(double s) const {
    if (!dist_)
      return s;
    double g = 0.0;
    MPI_Allreduce(&s, &g, 1, MPI_DOUBLE, MPI_SUM, dist_->comm());
    return g;
  }
  std::array<long, 3> globalFineExtent() const {
    std::array<long, 3> g{};
    for (int a = 0; a < 3; ++a)
      g[a] = dist_ ? static_cast<long>(dist_->globalFineSize()[a])
                   : static_cast<long>(pres_.fineExt()[a]);
    return g;
  }
 
  const FvOp& gpOp0() { return dist_ ? presMGD_.op(0) : presMG_.op(0); }
 
  template <class SdfFn>
  void prepareDistributed(SdfFn&& sdfFn) {
    dhalo_.init(*dist_);
    for (int a = 0; a < 3; ++a)
      shiftD_[a] = dist_->blockFineOrigin()[a];
    auto resv = [hp = &dhalo_, sh = shiftD_](const std::array<long, 3>& p) -> Index {
      std::array<long, 3> g = p;
      for (int a = 0; a < 3; ++a)
        g[a] += sh[a];
      return hp->resolveGlobal(g);
    };
    mom_.setFrameShift(shiftD_);
    mom_.setResolver(resv);
    pres_.setFrameShift(shiftD_);
    pres_.setResolver(resv);
    const Index n = t_->numLeaves();
    const Vec<3> beta = betaPerAxis();
    int rounds = 0;
    const bool profSetup = amrEnvFlag("PECLET_AMR_PROFILE_SETUP");
    for (;;) {
      ++rounds;
      const auto rT0 = std::chrono::steady_clock::now();
      installGhostMeta();  // metadata for every ghost known so far (same-round hits read it)
      // Classic path: mom_.build stays in the loop because its sdfC feeds the classic overlay
      // prober below. Sampled path (D3 setup-cost fix): its ±1 probes are a subset of the face
      // sweep's explicit periodicNeighbor probes, so the full build runs ONCE after the fixpoint
      // instead of once per round.
      if (!ghostSampled_)
        mom_.build(sdfFn, rho_ / dt_, beta);  // ±1 probes (also fills the extended sdfC/fluid)
      // FaceGeom probes: the full face enumeration + the SOU upstream-of-upwind ±2 reach.
      // Host-parallel since the F1 resolution (miss registration is mutex-guarded and the miss
      // SET is order-canonical — see LeafHalo::resolve); pure discovery, results discarded.
      hostParFor(n, [&](Index i) {
        pres_.forEachFaceFull(i, [&](Index j, int ax, int dr, double, double, double) {
          (void)pres_.periodicNeighbor(i, ax, -dr);
          if (j >= 0)
            (void)pres_.periodicNeighbor(j, ax, dr);
        });
      });
      // Overlay / directional-gradient ±2 chains (every cut cell is a non-clean overlay row,
      // so this covers buildGhostGradOverlay's probes too). Discovery only — result discarded.
      if (ghostSampled_) {
        // D1 (DD2): the SAMPLED builders are probers too, and they reach further than the classic
        // ±2 chain — an LS cloud descends a whole box of radius 2.2·max(h,H) around each virtual
        // position, and the mixed openness rule probes ±¼h0 off every face centroid. The overlay
        // runs in PROBE-ONLY discovery mode (D3 setup-cost fix): identical control flow and probe
        // pattern, pending regions blind-descended at 4-cell granularity so their misses register
        // in one collective round instead of one octree level per round, and none of the
        // least-squares work a discovery round would discard is done. The REAL build (setSolid,
        // after this fixpoint) runs with every ghost resolved and never sets the flag.
        const auto gfine = globalFineExtent();
        (void)buildGhostOverlaySampled(*t_, pres_, sdfFn, gpMatrixOrder_, gpRhsOrder_, origin_,
                                       &gfine, shiftD_, /*discovery=*/true, gpsRho_, gpsMaxN_);
        auto binFn = makeBinaryOpenFnMixed(
            *t_, pres_, [&sdfFn](const Vec<3>& p) { return sdfFn(p); }, h0_, origin_, shiftD_);
        pres_.buildOpenness(binFn);  // discovery only; setSolid rebuilds it after the fixpoint
      } else if (ghostProj_ || ghostGrad_) {
        bool viol = false;
        (void)buildGhostOverlay(*t_, pres_, mom_.sdfCRaw(), gpMatrixOrder_, gpRhsOrder_, &viol);
      }
      if (cfScheme_ != CfScheme::standard) {
        probeCfScheme();   // the four buildCf*Delta builders' tangential reach
        probeSeamLayer();  // the seam builder's far-side reach across a GHOST neighbour
      }
      const double rBuildMs =
          std::chrono::duration<double, std::milli>(std::chrono::steady_clock::now() - rT0).count();
      const long pend = dhalo_.resolveMisses();
      if (profSetup)
        std::fprintf(stderr, "[setSolid]   round %d: build %.1f ms, %lld ghosts, %ld pending\n",
                     rounds, rBuildMs, (long long)dhalo_.numGhosts(), pend);
      if (pend == 0)
        break;
    }
    if (profSetup)
      // The discovery fixpoint's round count. The classic +-2 chain converges in <=3; the SAMPLED
      // builders reach a whole LS box, and D0's descent deliberately learns one octree level per
      // round rather than registering a miss per fine cell, so the sampled path needs more (it
      // grows with the level span the clouds cover). Each round is a full rebuild + a collective.
      std::fprintf(stderr, "[setSolid] distributed discovery fixpoint: %d rounds (%lld ghosts)\n",
                   rounds, (long long)dhalo_.numGhosts());
    installGhostMeta();
    if (ghostSampled_)
      mom_.build(sdfFn, rho_ / dt_, beta);  // once, every coord it probes already resolved
    dhalo_.finalize();
    nExt_ = dhalo_.extendedSize();
    dhex_.init(dhalo_);
    allred_ = [this](double s) { return allSum(s); };
    momSolver_.setDistributed([this](View<double> v) { dhex_.exchange(v); }, allred_, nExt_);
  }
 
  void probeCfScheme() {
    const Index n = t_->numLeaves();
    hostParFor(n, [&](Index i) {
      const unsigned Li = pres_.levelOf(i);
      pres_.forEachFaceFull(i, [&](Index j, int axis, int, double, double, double) {
        const unsigned Lj = pres_.levelOf(j);
        if (Lj == Li)
          return;
        probeCfTangential((Lj > Li) ? j : i, axis);
      });
    });
  }
 
  void probeCfTangential(Index coarse, int axis) {
    const unsigned Lc = pres_.levelOf(coarse);
    const std::array<long, 3> bc = pres_.loOf(coarse);
    const long sc = 1L << Lc, sh = sc >> 1;
    for (int tt = 0; tt < 3; ++tt) {
      if (tt == axis)
        continue;
      for (int dir = -1; dir <= 1; dir += 2) {
        (void)pres_.periodicNeighbor(coarse, tt, dir);
        std::array<long, 3> lo = bc;
        lo[tt] += (dir > 0) ? sc : -sc;
        for (int oct = 0; oct < (1 << 3); ++oct) {
          std::array<long, 3> q = lo;
          for (int d = 0; d < 3; ++d)
            if ((oct >> d) & 1)
              q[d] += sh;
          (void)pres_.probeSlot(q);
        }
      }
    }
  }
 
  void probeSeamLayer() {
    const Index n = t_->numLeaves();
    hostParFor(n, [&](Index i) {
      pres_.forEachFaceFull(i, [&](Index j, int axis, int dir, double, double, double) {
        if (j < 0)
          return;
        const unsigned Lj = pres_.levelOf(j);
        if (Lj > 0) {  // the four face-layer corners beyond j (a layer record with C == j)
          const std::array<long, 3> lo = pres_.loOf(j);
          const long sj = 1L << Lj, sh = sj >> 1;
          std::array<long, 3> p{};
          for (int d = 0; d < 3; ++d)
            p[d] = lo[d];
          p[axis] = (dir > 0) ? lo[axis] + sj : lo[axis] - 1;
          for (int k = 0; k < 4; ++k) {
            std::array<long, 3> q = p;
            int bit = 0;
            for (int t = 0; t < 3; ++t) {
              if (t == axis)
                continue;
              q[t] = lo[t] + (((k >> bit) & 1) ? sh : 0L);
              ++bit;
            }
            (void)pres_.probeSlot(q);
          }
        }
        const Index uJ = pres_.periodicNeighbor(j, axis, dir);
        if (uJ >= 0 && pres_.levelOf(uJ) > Lj)
          probeCfTangential(uJ, axis);  // a sample record with C == uJ, F == j (case 2)
      });
    });
  }
 
  void installGhostMeta() {
    const Index ng = dhalo_.numGhosts();
    std::vector<std::array<long, 3>> glo(static_cast<std::size_t>(ng));
    std::vector<unsigned> glv(static_cast<std::size_t>(ng));
    for (Index g = 0; g < ng; ++g) {
      for (int a = 0; a < 3; ++a)
        glo[static_cast<std::size_t>(g)][a] =
            static_cast<long>(dhalo_.ghostCoord(g)[a]) - shiftD_[a];
      glv[static_cast<std::size_t>(g)] = static_cast<unsigned>(dhalo_.level(dhalo_.numLocal() + g));
    }
    mom_.setGhosts(glo, glv);  // copies — pres_ takes the originals
    pres_.setGhosts(std::move(glo), std::move(glv));
  }
 
  // (public like runMgVcycle: nvcc rejects extended device lambdas in private member functions)
  // Project onto the coupled subspace: pin decoupled rows (solid-centered + no-phi-coupling
  // overlay rows) to 0 and remove the volume-weighted mean over the coupled cells (the constant
  // null mode of the connected fluid region) — removeMeanVol with the coupled mask (the mean
  // reduced globally in distributed mode; allred_ is empty single-rank ⇒ bit-identical).
  void gpProject(View<double> v) { removeMeanVolReduced(v, gpOp0().invVol, maskC_, n_, allred_); }
 
  // Nonsymmetric ghost pressure matvec: y = P[rho·(L_bin x + Delta x)]. The caller keeps x's
  // ghost tail current (syncScalar before every call in distributed mode).
  void ghostMatvec(View<const double> x, View<double> y) {
    auto spT = spMark();
    applyFv(gpOp0(), x, y);
    spAdd(SP_PRES_MV, spT);
    spT = spMark();
    ghostApplyDelta(gpOv_, x, y);
    ghostApplyDeltaSampled(gpOvS_, x, y);
    spAdd(SP_PRES_MVOV, spT);
    spT = spMark();
    gpProject(y);
    spAdd(SP_PRES_PROJ, spT);
  }
 
  // Preconditioner: two binary-openness V-cycles (the unchanged MG hierarchy) + projection.
  void ghostPrec(View<const double> r, View<double> z) {
    const auto spT = spMark();
    if (dist_) {
      Kokkos::deep_copy(presMGD_.b(0), r);
      Kokkos::deep_copy(presMGD_.x(0), 0.0);
      presMGD_.vcycle(2, 2, 60, 0.8);
      presMGD_.vcycle(2, 2, 60, 0.8);
      Kokkos::deep_copy(z, presMGD_.x(0));
    } else {
      Kokkos::deep_copy(presMG_.b(0), r);
      Kokkos::deep_copy(presMG_.x(0), 0.0);
      presMG_.vcycle(2, 2, 60, 0.8);
      presMG_.vcycle(2, 2, 60, 0.8);
      Kokkos::deep_copy(z, presMG_.x(0));
    }
    gpProject(z);
    spAdd(SP_PRES_PC, spT);
  }
 
  // MG-preconditioned BiCGStab on the ghost pressure operator (device mirror of the oracle's
  // solveGhostBiCGStab: same projection, same stagnation guard against the small
  // attainable-residual floor of the slightly incompatible ghost system). Returns iterations.
  int solveGhostBiCGStab(View<double> x, View<const double> b, int maxIters, double tol = 1e-10) {
    const Index n = n_;
    // Dots are locally summed then globally reduced (allSum = identity single-rank); every
    // matvec input's ghost tail is refreshed first. Stagnation/early-break branches depend on
    // reduced scalars only ⇒ every rank takes the same branch.
    syncScalar(x);
    ghostMatvec(View<const double>(x), gpr_);
    {
      auto r = gpr_;
      auto bb = b;
      Kokkos::parallel_for("amr::gp_r0", n, KOKKOS_LAMBDA(const Index i) { r(i) = bb(i) - r(i); });
    }
    gpProject(gpr_);
    Kokkos::deep_copy(gprh_, gpr_);
    const double res0 =
        std::sqrt(allSum(dotPlain(View<const double>(gpr_), View<const double>(gpr_), n)));
    if (res0 == 0.0)
      return 0;
    double rho = 1, alpha = 1, omega = 1, best = res0;
    int noImprove = 0;
    Kokkos::deep_copy(gpv_, 0.0);
    Kokkos::deep_copy(gpp_, 0.0);
    int it = 0;
    for (; it < maxIters; ++it) {
      const double rhoNew =
          allSum(dotPlain(View<const double>(gprh_), View<const double>(gpr_), n));
      if (rhoNew == 0.0)
        break;
      const double beta = (rhoNew / rho) * (alpha / omega);
      bicgPUpdate(gpp_, View<const double>(gpr_), View<const double>(gpv_), beta, omega, n);
      ghostPrec(View<const double>(gpp_), gpph_);
      syncScalar(gpph_);
      ghostMatvec(View<const double>(gpph_), gpv_);
      const double rhatV = allSum(dotPlain(View<const double>(gprh_), View<const double>(gpv_), n));
      if (rhatV == 0.0)
        break;
      alpha = rhoNew / rhatV;
      Kokkos::deep_copy(gps_, gpr_);
      axpy(gps_, -alpha, View<const double>(gpv_), n);
      const double snorm =
          std::sqrt(allSum(dotPlain(View<const double>(gps_), View<const double>(gps_), n)));
      if (snorm <= tol * res0) {
        axpy(x, alpha, View<const double>(gpph_), n);
        ++it;
        break;
      }
      ghostPrec(View<const double>(gps_), gpsh_);
      syncScalar(gpsh_);
      ghostMatvec(View<const double>(gpsh_), gpt_);
      const double tt = allSum(dotPlain(View<const double>(gpt_), View<const double>(gpt_), n));
      omega = (tt != 0.0)
                  ? allSum(dotPlain(View<const double>(gpt_), View<const double>(gps_), n)) / tt
                  : 0.0;
      axpy(x, alpha, View<const double>(gpph_), n);
      axpy(x, omega, View<const double>(gpsh_), n);
      Kokkos::deep_copy(gpr_, gps_);
      axpy(gpr_, -omega, View<const double>(gpt_), n);
      const double rnorm =
          std::sqrt(allSum(dotPlain(View<const double>(gpr_), View<const double>(gpr_), n)));
      if (rnorm <= tol * res0) {
        ++it;
        break;
      }
      if (rnorm < 0.999 * best) {
        best = rnorm;
        noImprove = 0;
      } else if (++noImprove >= 6) {
        ++it;
        break;  // attainable-residual floor (compatibility gap) — stagnation guard
      }
      rho = rhoNew;
      if (omega == 0.0)
        break;
    }
    gpProject(x);
    return it;
  }
 
 private:
  void uploadSeamRecon(const SeamReconHost& h) {
    geom_.seam = toDevice(h.seam, "seam_desc");
    geom_.sampI = toDevice(h.sampI, "seam_sampi");
    geom_.sampJ = toDevice(h.sampJ, "seam_sampj");
    geom_.uuRecI = toDevice(h.uuRecI, "seam_uureci");
    geom_.uuRecJ = toDevice(h.uuRecJ, "seam_uurecj");
    geom_.d1I = toDevice(h.d1I, "seam_d1i");
    geom_.d1J = toDevice(h.d1J, "seam_d1j");
    geom_.recStart = toDevice(h.recStart, "seam_recstart");
    geom_.recDist = toDevice(h.recDist, "seam_recdist");
    geom_.recCell = toDevice(h.recCell, "seam_reccell");
    geom_.recW = toDevice(h.recW, "seam_recw");
    seamSampleRecords_ = h.numSampleRecords;
    seamLayerRecords_ = h.numLayerRecords;
  }
  void clearSeamRecon() {
    geom_.seam = View<Index>{};
    geom_.sampI = View<Index>{};
    geom_.sampJ = View<Index>{};
    geom_.uuRecI = View<Index>{};
    geom_.uuRecJ = View<Index>{};
    geom_.d1I = View<double>{};
    geom_.d1J = View<double>{};
    geom_.recStart = View<Index>{};
    geom_.recDist = View<double>{};
    geom_.recCell = View<Index>{};
    geom_.recW = View<double>{};
    seamSampleRecords_ = 0;
    seamLayerRecords_ = 0;
  }
 
  // Host build of the ghost-gradient overlay (setGhostGradient): one row per cut cell (fluid
  // with a solid face neighbour — the cells where the ABC grad3 is gauge-dependent O(1/h)),
  // holding a 3-point directional FD stencil per axis. Mirrors oracle::AmrFlow::gradOfDir: cut
  // cells have same-level face neighbours by the finest-band contract; the ±2 probe falls back
  // to the 2-point one-sided closure when that cell is solid or not same-level.
  void buildGhostGradOverlay() {
    const Index n = t_->numLeaves();
    std::vector<Index> cells;
    for (Index i = 0; i < n; ++i)
      if (mom_.isCut(i))
        cells.push_back(i);
    const Index m = static_cast<Index>(cells.size());
    std::vector<Index> idx(static_cast<std::size_t>(m) * 9, 0);
    std::vector<double> w(static_cast<std::size_t>(m) * 9, 0.0);
    // Pocket cells (fragmentation guard) count as solid for the directional gradient: their φ is
    // pinned/decoupled, and reading the pinned 0 is the gauge-dependent defect the ghost
    // gradient exists to avoid. Levels via pres_.levelOf: ghost-slot-safe in distributed
    // builds (identical to t_->level for local leaves; the pocket guard is single-rank-only —
    // gpPocket_ stays empty distributed until the label-propagation guard lands).
    auto ok = [&](Index j, Index i) {
      return j >= 0 && mom_.isFluid(j) && pres_.levelOf(j) == pres_.levelOf(i) &&
             !(!gpPocket_.empty() && j < n && gpPocket_[static_cast<std::size_t>(j)]);
    };
    for (Index s = 0; s < m; ++s) {
      const Index i = cells[static_cast<std::size_t>(s)];
      const double h = pres_.cellWidth(i);
      for (int a = 0; a < 3; ++a) {
        const std::size_t o = static_cast<std::size_t>(s) * 9 + static_cast<std::size_t>(a) * 3;
        for (int k = 0; k < 3; ++k)
          idx[o + static_cast<std::size_t>(k)] = i;  // safe defaults (w = 0)
        const Index jp = pres_.periodicNeighbor(i, a, +1);
        const Index jm = pres_.periodicNeighbor(i, a, -1);
        const bool ap = ok(jp, i), am = ok(jm, i);
        if (am && ap) {
          idx[o] = jp;
          w[o] = 0.5 / h;
          idx[o + 1] = jm;
          w[o + 1] = -0.5 / h;
        } else if (ap) {
          const Index jpp = pres_.periodicNeighbor(jp, a, +1);
          if (ok(jpp, i)) {
            idx[o] = i;
            w[o] = -1.5 / h;
            idx[o + 1] = jp;
            w[o + 1] = 2.0 / h;
            idx[o + 2] = jpp;
            w[o + 2] = -0.5 / h;
          } else {
            idx[o] = jp;
            w[o] = 1.0 / h;
            idx[o + 1] = i;
            w[o + 1] = -1.0 / h;
          }
        } else if (am) {
          const Index jmm = pres_.periodicNeighbor(jm, a, -1);
          if (ok(jmm, i)) {
            idx[o] = i;
            w[o] = 1.5 / h;
            idx[o + 1] = jm;
            w[o + 1] = -2.0 / h;
            idx[o + 2] = jmm;
            w[o + 2] = 0.5 / h;
          } else {
            idx[o] = i;
            w[o] = 1.0 / h;
            idx[o + 1] = jm;
            w[o + 1] = -1.0 / h;
          }
        }  // sandwiched: all weights stay 0
      }
    }
    gc_.n = m;
    gc_.cell = toDevice(cells, "gc_cell");
    gc_.idx = toDevice(idx, "gc_idx");
    gc_.w = toDevice(w, "gc_w");
  }
 
  // Fluid area fraction of a face. DEFAULT (order 2, 2026-08-26 user decision): triangle-fan
  // marching squares on FIVE ANALYTIC samples (4 corners + centre) -- exact linear fraction per
  // triangle, O(h^2), no saddle ambiguity, and the sub-resolution floor 1e-6 (alpha ~ 1e-12 rows
  // destroy the pressure conditioning -- flow tracker rows 51/52). Because the AMR solver holds
  // the analytic sdfFn, this samples the TRUE geometry (no trilinear ceiling). order 1 = the
  // legacy one-sample gradient-normalised linear model (kept for A/B), which carries a SIGNED
  // convexity bias measured at +0.59%/+0.27% bed permeability at R=8/12 (decay ~h^2).
  static double triFrac(double a, double b, double c) {
    const bool pa = a >= 0.0, pb = b >= 0.0, pc = c >= 0.0;
    const int np = (pa ? 1 : 0) + (pb ? 1 : 0) + (pc ? 1 : 0);
    if (np == 3)
      return 1.0;
    if (np == 0)
      return 0.0;
    double x, y, z;
    if (np == 1) {
      if (pa) {
        x = a;
        y = b;
        z = c;
      } else if (pb) {
        x = b;
        y = c;
        z = a;
      } else {
        x = c;
        y = a;
        z = b;
      }
      const double den = (x - y) * (x - z);
      return den > 1e-300 ? (x * x) / den : 1.0;
    }
    if (!pa) {
      x = a;
      y = b;
      z = c;
    } else if (!pb) {
      x = b;
      y = c;
      z = a;
    } else {
      x = c;
      y = a;
      z = b;
    }
    const double den = (x - y) * (x - z);
    return 1.0 - (den > 1e-300 ? (x * x) / den : 1.0);
  }
  template <class SdfFn>
  double faceFrac(SdfFn&& sdfFn, const Vec<3>& fc, int axis) const {
    if (apertureOrder_ >= 2) {
      if (sdfFn(fc) <= 0.0)
        return 0.0;  // center gate (kept from order 1; see flow ccFaceOpenMS)
      const int t1 = (axis + 1) % 3, t2 = (axis + 2) % 3;
      // Phase 3: the four corner probes sit half a cell out on EACH tangential axis.
      const double e1 = 0.5 * h0_[t1], e2 = 0.5 * h0_[t2];
      auto at = [&](double d1, double d2) {
        Vec<3> p = fc;
        p[t1] += d1;
        p[t2] += d2;
        return sdfFn(p);
      };
      const double c00 = at(-e1, -e2), c10 = at(e1, -e2), c11 = at(e1, e2), c01 = at(-e1, e2);
      const double cc = sdfFn(fc);
      const double frac = 0.25 * (triFrac(c00, c10, cc) + triFrac(c10, c11, cc) +
                                  triFrac(c11, c01, cc) + triFrac(c01, c00, cc));
      return frac < 1e-3 ? 0.0 : (frac > 1.0 - 1e-12 ? 1.0 : frac);  // floor: see flow ccFaceOpenMS
    }
    double sd = sdfFn(fc);
    if (sd <= 0.0)
      return 0.0;
    Vec<3> g{};
    for (int d = 0; d < 3; ++d) {
      Vec<3> pp = fc, pm = fc;
      pp[d] += h0_[d];
      pm[d] -= h0_[d];
      g[d] = (sdfFn(pp) - sdfFn(pm)) / (2.0 * h0_[d]);
    }
    double gmag = std::sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
    if (gmag < 1e-6)
      gmag = 1e-6;
    int t1 = (axis + 1) % 3, t2 = (axis + 2) % 3;
    // Phase 3: the projected in-plane extent of the face, per tangential axis.
    double denom = (std::fabs(g[t1]) * h0_[t1] + std::fabs(g[t2]) * h0_[t2]) / gmag;
    if (denom < 1e-9)
      denom = 1e-9;
    double frac = 0.5 + sd / denom;
    return frac < 0.0 ? 0.0 : (frac > 1.0 ? 1.0 : frac);
  }
 
  const Octree* t_ = nullptr;
  Vec<3> h0_ = detail::filledVec<3>(1.0);
  Vec<3> origin_{};
  double rho_ = 1.0, mu_ = 1.0, dt_ = 1e6;
  Vec<3> f_{};
  bool presPCG_ = true;
  bool presDbgSpdDone_ = false;  // one-shot debug SPD probe (PECLET_AMR_PRES_DEBUG)
  bool momMGon_ = true;  // velocity-MG momentum preconditioner (scalable; see setMomentumMG)
  bool useStaircaseMG_ = false;  // false = Galerkin (MomentumMG), true = staircase (VelocityMG)
  int mgVcPre_ = 2, mgVcBottom_ = 30;  // momentum-MG V-cycle pre/post sweeps + bottom sweeps
  Index mgMinCoarse_ = 256;            // staircase velocity-MG pore-scale cap (coarsest cell count)
  bool momGS_ = false;  // opt-in: multicolour Gauss–Seidel smoother in the momentum MG
  bool momMGSolver_ =
      false;  // opt-in (P4): velocity-MG as the solver (defect correction), not BiCGStab
  bool ghostGrad_ = true;     // gauge-exact directional cell gradient (setGhostGradient) — the
                              // DEFAULT since 2026-08-18, mirroring flow's collocated
                              // set_collocated_scheme("gauge-exact")
  bool ghostProj_ = false;    // RESOLVED projection mode (set by setSolid from the request)
  int8_t ghostProjReq_ = -1;  // -1 = AUTO (DEFAULT since 2026-08-25: ghost, aperture fallback
                              // on thin band), 0 = explicit aperture, 1 = explicit ghost
  int apertureOrder_ = 2;     // aperture estimator order (setApertureOrder; default 2)
  int gpMatrixOrder_ = 2, gpRhsOrder_ = 2;  // closure orders (2,2 = the production pair; the
                                            // (1,2) mixed form is march-unstable at scale)
  // DEFAULT since 2026-09-21 (user decision): the Martin-Cartwright tangential quadratic, NOT the
  // standard two-point flux. The register had already settled that "cf=1 is not optional on graded
  // meshes" and rejected cf=0 there, but the shipped default stayed standard, so a user who built a
  // graded mesh silently got the configuration the project had decided against. Measured on a
  // self-similar ladder (docs/amr_graded_convergence.md §3), as orders:
  //
  //                          standard        quadratic
  //      interface NORMAL to the variation      2.00    2.00 (inert: no tangential variation to
  //      interface TANGENTIAL                   0.41    1.60  correct, so the delta is empty)
  //
  // with the absolute gap on the tangential arm widening 3.8x at n=32 to 8.8x at n=64. INERT BY
  // GEOMETRY on any uniform or finest-band mesh: no C/F faces means an empty delta CSR and the
  // identical code path, so only graded runs move.
  CfScheme cfScheme_ = CfScheme::quadratic;  // 2:1 C/F interface scheme (setCfScheme)
  int outerIters_ = 1;       // Picard outer iterations over the lagged advection (default 1)
  double outerTol_ = 1e-6;   // outer-loop early-stop tolerance on max|Δu|
  double momTol_ = 1e-8;     // per-step momentum BiCGStab relative tolerance (Phase-0 knob)
  double presTol_ = 1e-10;   // pressure solve relative tolerance (MG-PCG and ghost BiCGStab)
  bool advect_ = false;      // momentum advection ∇·(u u) (off ⇒ Stokes)
  bool implicitFou_ = true;  // implicit-FOU deferred-correction (stable) vs fully-explicit
  int advScheme_ = 0;        // high-order flux: 0 = SOU (default), 1 = Koren TVD
  Index n_ = 0;
  bool ready_ = false;  // an operator exists (setSolid ran; cleared by beginAdapt)
  int lastMomIters_ = 0, lastPresIters_ = 0, lastOuterIters_ = 1;
  // M0 step profiler (PECLET_AMR_PROFILE_STEP): all inert unless stepProf_.
  bool stepProf_ = amrEnvFlag("PECLET_AMR_PROFILE_STEP");
  bool spHeader_ = false;
  int spWindow_ = [] {
    const char* e = std::getenv("PECLET_AMR_PROFILE_STEP_WINDOW");
    const int v = e ? std::atoi(e) : 0;
    return v > 0 ? v : 50;
  }();
  int spSteps_ = 0;
  double spAcc_[SP_N] = {};
  double spMom_ = 0.0, spPres_ = 0.0, spOuter_ = 0.0;
 
  AmrCutCell<Bits> mom_;
  AmrPoisson<3, Bits> pres_;
  Multigrid<3, Bits> presMG_;
  typename Multigrid<3, Bits>::Bottom presBottom_ = Multigrid<3, Bits>::Bottom::Auto;
  Index presBottomExtent_ = 4;
  MomentumMG<Bits> momMG_;  // Galerkin velocity multigrid (momentum preconditioner)
  VelocityMG<Bits> velMG_;  // rediscretized staircase velocity multigrid (alternative)
  MomentumOp momOp_;
  MomentumSolver<Bits> momSolver_;
  PCG<3, Bits> pcg_;
  std::array<View<double>, 3> defc_;  // explicit ρ(SOU−FOU) deferred correction per component
  View<double> advDiag_, advCoef_;    // device-resident implicit-FOU operator (rebuilt each step)
  FaceGeom geom_;
  GhostGradOverlay gc_;        // directional ghost-gradient overlay (empty unless setGhostGradient)
  GhostOverlayDev gpOv_;       // closure overlay (empty unless setGhostProjection)
  bool ghostSampled_ = false;  // RESOLVED sampled mode (set by setSolid)
  int8_t ghostSampledReq_ = 0;      // setGhostSampled request (mixed-level cut band)
  double gpsRho_ = 2.2;             // sampled band: LS cloud radius factor (setGhostSampled)
  long gpsMaxN_ = 0;                // sampled band: nearest-N candidate cap, 0 = uncapped
  GhostOverlaySampledDev gpOvS_;    // sample-slot overlay (empty unless sampled)
  GhostGradCsrDev gcS_;             // sampled CSR directional-gradient overlay
  CfCsrDev gpsMomDelta_;            // momentum ξ-row seam correction (rscale-folded)
  std::vector<char> gpPocket_;      // fragmentation guard: 1 = decoupled pocket cell (ghost mode)
  CfCsrDev cfMom_;                  // +μ(∇²_scheme − ∇²_std) momentum RHS overlay
  CfCompCsrDev cfDiv_;              // (D_scheme − D_std) divergence overlay
  std::array<CfCsrDev, 3> cfGrad_;  // (G_scheme − G_std) per gradient axis
  CfCompCsrDev cfUfVel_;            // (uf_scheme − uf_std) face-field overlay: velocity part
  Index cfGhostCols_ = 0;           // overlay entries reading a ghost slot (diagnostic)
  Index cfCutFaces_ = 0;            // C/F slots where the face gate withholds the quadratic
  bool seamRecon_ = true;           // seam reconstruction of the advected value (§5.5, default on)
  Index seamSampleRecords_ = 0;     // seam tables: sample records built (docs/amr_cf_convective.md)
  Index seamLayerRecords_ = 0;      // seam tables: face-layer records built
  View<double> maskC_;              // 1 = coupled row (Krylov subspace), 0 = pinned
  View<double> gpr_, gprh_, gpp_, gpph_, gpv_, gps_, gpsh_, gpt_;  // ghost BiCGStab scratch
  View<double> rscale_;
  View<char> fluid_;
  std::array<View<double>, 3> u_, gx_;
  std::array<View<double>, 3> u0_,
      uprev_;  // frozen uⁿ (BE mass term) + previous Picard outer iterate
  View<double> p_, phi_, div_, bmom_;
  View<double> uf_;  // ABC/Basilisk divergence-free face field (one per CSR (sub)face)
  bool faceFieldBuilt_ =
      false;  // uf_ populated by a projection (else advection falls back to ½(u_i+u_j))
  bool ufAdvect_ = true;  // ABLATION (diagnostics): advect with uf_ (ON, the shipped scheme) or
                          // with the un-projected ½(u_i+u_j) — see setUfAdvection
  std::unique_ptr<Octree> adaptOldT_;          // beginAdapt topology snapshot
  std::array<std::vector<double>, 3> adaptU_;  // beginAdapt field snapshots
  std::vector<double> adaptP_;
  std::array<std::vector<std::array<double, 3>>, 3> adaptGradU_;  // distributed transfer grads
  std::vector<std::array<double, 3>> adaptGradP_;
 
  // ---- distributed context (initMpi; all null/empty single-rank) -----------------------------
  DistributedOctree<3, Bits>* dist_ = nullptr;  // the ORB block + communicator
  LeafHalo<3, Bits> dhalo_;                    // the flow's ±2 ghost registry (frozen in setSolid)
  LeafHaloExchange dhex_;                      // device value refresh over dhalo_
  DistributedFlowMultigrid<3, Bits> presMGD_;  // distributed pressure MG (level 0 on dhalo_)
  std::array<long, 3> shiftD_{};               // block global fine origin
  Index nExt_ = 0;                             // n_ + ghosts (== n_ single-rank)
  std::function<double(double)> allred_;       // Allreduce hook (empty single-rank)
};
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_FLOW_HPP