mac_approx_projection.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_MAC_APPROX_PROJECTION_HPP
#define PECLET_FLOW_MAC_APPROX_PROJECTION_HPP
 
#include <Kokkos_Core.hpp>
 
#include "mac_cutcell.hpp"  // peclet::flow::C3, CCField, CCConst, CCExec, ccSampleExt
 
namespace peclet::flow {
 
// const view of the (float) IBM stencil coefficients (same alias mac_ibm.hpp defines; this header
// is included first, so declare it here for the mode-4 stencilMatvec).
using MConst = Kokkos::View<const float*, CCMem>;
 
// Average cell-centered velocities onto the staggered face layout: uf(i) is the velocity at the low
// (-x) face of cell i (located at i-1/2) = ½(U(i-1)+U(i)); likewise vf/wf along y/z. This
// reproduces the exact layout the staggered solver stores directly, so divergOpen / projectCorrect
// (mac_pressure.hpp) act on uf/vf/wf unchanged. Computed over the block where the axis neighbour
// exists ([1,e) per axis); the cell velocity ghosts must be filled first.
inline void centerToFace(CCField uf, CCField vf, CCField wf, CCConst U, CCConst V, CCConst W, C3 e,
                         int g) {
  (void)g;
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::center_to_face", MD(space, {1, 1, 1}, {e.x, e.y, e.z}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        uf(i) = 0.5 * (U(i) + U(i - sx));
        vf(i) = 0.5 * (V(i) + V(i - sy));
        wf(i) = 0.5 * (W(i) + W(i - sz));
      });
}
 
// Wall-aware variant of centerToFace (opt-in: Solver::setFaceInterp(1); collocated only). At a face
// whose adjacent cell centers straddle the immersed boundary (one SDF negative), the plain average
// ½(U_fluid + U_solid=0) implicitly places the wall AT THE SOLID CELL CENTER (theta = 1) regardless
// of the true axis intercept theta = sdf_c/(sdf_c - sdf_n) — an O(1) relative flux error at every
// cut face that makes the collocated drag first-order (doc/collocated_first_order_analysis.md).
// Here the face value is reconstructed from the fluid side by a WALL-ANCHORED quadratic in
// (theta - s) — exactly zero at the wall — least-squares-fitted to the fluid-side samples
// {U_c, U_far, U_far2} with the near-wall cell SHED smoothly as theta -> 0 (weight
// theta²/(theta²+0.15²)): a sliver cell's value is slaved to the wall by the momentum IBM and
// carries no independent information, so downweighting it bounds the noise amplification without
// dropping the order — the theta->0 limit is the (still 2nd-order) quadratic through
// {wall, far, far2}. A small Tikhonov term on the curvature keeps the 2x2 solve well-posed when
// far2 is unavailable; a fluid cell sandwiched along the axis falls back to the wall-anchored
// linear through U_c (bounded: xe/theta <= 1/2). Faces whose center lies behind the wall
// (theta <= 1/2) get 0, continuously (the parabola vanishes at the wall); their open area is a
// corner sliver. Interior fluid faces reproduce the plain average bit-exactly. Validated a priori
// against analytic Stokes flow past a sphere: face-value error O(h²) (the plain average is O(h)).
// Per-face stencil of the wall-aware map: up to 3 (cell index, weight) pairs such that
// uf(face) = sum_k w[k]*U(idx[k]); returns the entry count (0 = closed face or corner sliver:
// uf = 0). Shared by the forward map (centerToFaceWallAware) and its TRANSPOSE
// (transposeGradWallAware) so the two stay exact adjoints — the accuracy of the constrained steady
// state hinges on that pairing, not on the forward truncation alone.
KOKKOS_INLINE_FUNCTION int wallAwareFaceStencil(CCConst sdf, long i, long sa, int cFace, int eAxis,
                                                double xcen, long idx[3], double w[3]) {
  const double sm = sdf(i - sa), sp = sdf(i);  // low cell (coord cFace-1) / high cell (cFace)
  const bool fm = sm >= 0.0, fp = sp >= 0.0;
  if (fm && fp) {  // both fluid: the plain 1/2-1/2 average, unchanged
    idx[0] = i;
    w[0] = 0.5;
    idx[1] = i - sa;
    w[1] = 0.5;
    return 2;
  }
  if (!fm && !fp)
    return 0;                       // fully solid face (openness 0; value irrelevant)
  const long ic = fp ? i : i - sa;  // the fluid cell
  const double sc = fp ? sp : sm, ss = fp ? sm : sp;
  double th = sc / (sc - ss);  // axis intercept from the fluid center (ibmFillEntry convention)
  th = th < 1e-4 ? 1e-4 : (th > 1.0 ? 1.0 : th);
  // Evaluation abscissa in wall coordinates (theta - s): the face-CENTER distance theta - 1/2 by
  // default, or (xcen >= 0, mode 3) the OPEN-FACE-CENTROID wall distance from the static geometry
  // build. The projection's flux model is o_f * uf, so the correct uf is the open-area MEAN of the
  // velocity — the midpoint-at-centroid quadrature — not the face-center point value: the open area
  // lies on the fluid side, farther from the wall, so the (accurate) face-center value UNDERCOUNTS
  // the flux at O(h). Measured a priori (Stokes-past-sphere): centroid evaluation drops the
  // systematic flux bias by ~2 orders of magnitude at N=128.
  const double xe = (xcen >= 0.0) ? xcen : th - 0.5;
  if (xe <= 0.0)
    return 0;                             // face (or its open centroid) behind the wall
  const int cc = fp ? cFace : cFace - 1;  // fluid cell's axis coordinate
  const int up = fp ? 1 : -1;             // coordinate step AWAY from the solid
  const long dir = fp ? -sa : sa;         // index step TOWARD the solid
  const bool ok1 = (cc + up >= 0) && (cc + up < eAxis) && sdf(ic - dir) >= 0.0;
  if (!ok1) {  // sandwiched along this axis: wall-anchored linear through U_c
    idx[0] = ic;
    w[0] = Kokkos::fmin(xe / th, 2.0);  // capped (xe can exceed theta at centroid evaluation)
    return 1;
  }
  const bool ok2 = (cc + 2 * up >= 0) && (cc + 2 * up < eAxis) && sdf(ic - 2 * dir) >= 0.0;
  const double wc = th * th / (th * th + 0.0225);  // shed the sliver cell (theta* = 0.15)
  const double w2 = ok2 ? 0.25 : 0.0;
  const double xc = th, xf = th + 1.0, xg = th + 2.0;
  const double S11 = wc * xc * xc + xf * xf + w2 * xg * xg;
  const double S12 = wc * xc * xc * xc + xf * xf * xf + w2 * xg * xg * xg;
  const double S22 = wc * xc * xc * xc * xc + xf * xf * xf * xf + w2 * xg * xg * xg * xg + 0.01;
  const double det = S11 * S22 - S12 * S12;
  // coefficient of datum k (abscissa x_k, weight w_k) in T = a*xe + b*xe^2:
  //   c_k = w_k * [ x_k*(S22*xe - S12*xe^2) + x_k^2*(S11*xe^2 - S12*xe) ] / det
  const double A = (S22 * xe - S12 * xe * xe) / det, B = (S11 * xe * xe - S12 * xe) / det;
  idx[0] = ic;
  w[0] = wc * (xc * A + xc * xc * B);
  idx[1] = ic - dir;
  w[1] = xf * A + xf * xf * B;
  if (!ok2)
    return 2;
  idx[2] = ic - 2 * dir;
  w[2] = w2 * (xg * A + xg * xg * B);
  return 3;
}
 
KOKKOS_INLINE_FUNCTION double wallAwareFaceValue(CCConst U, CCConst sdf, long i, long sa, int cFace,
                                                 int eAxis, double xcen) {
  long idx[3];
  double w[3];
  const int n = wallAwareFaceStencil(sdf, i, sa, cFace, eAxis, xcen, idx, w);
  double v = 0.0;
  for (int k = 0; k < n; ++k)
    v += w[k] * U(idx[k]);
  return v;
}
 
// Static per-face geometry for the mode-3 centroid quadrature: for each face whose adjacent cell
// centers straddle the boundary, subsample the face plane (4x4 trilinear SDF), locate the OPEN-area
// centroid, and store its axis wall-distance (in cells, measured toward the solid side; two-point
// line intercept from the trilinear SDF, clamped to [0,2]). 0 => no open samples (no flux). Faces
// with both centers on one side store -1 (sentinel: stencil falls back to theta-1/2, which the
// both-fluid branch never reads anyway). One-time cost at setSolid; the geometry is static.
inline void buildFaceCentroidDist(CCField xcx, CCField xcy, CCField xcz, CCConst sdf, C3 e) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::face_centroid_dist", MD(space, {1, 1, 1}, {e.x, e.y, e.z}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        // face center in the sampling frame where CELL CENTERS sit at integer coordinates
        // (buildOpenness convention): -x face of cell (x,y,z) is at (x-1/2, y, z), etc.
        const double fc[3][3] = {{x - 0.5, (double)y, (double)z},
                                 {(double)x, y - 0.5, (double)z},
                                 {(double)x, (double)y, z - 0.5}};
        const long st[3] = {sx, sy, sz};
        CCField out[3] = {xcx, xcy, xcz};
        for (int a = 0; a < 3; ++a) {
          const double sm = sdf(i - st[a]), sp = sdf(i);
          if ((sm >= 0.0) == (sp >= 0.0)) {
            out[a](i) = -1.0;  // not a solid-bordering face: stencil never reads it
            continue;
          }
          const double dir = (sm < 0.0) ? -1.0 : 1.0;  // axis step TOWARD the solid side
          const int t1 = (a + 1) % 3, t2 = (a + 2) % 3;
          double c1 = 0.0, c2 = 0.0;
          int nOpen = 0;
          for (int j1 = 0; j1 < 4; ++j1)
            for (int j2 = 0; j2 < 4; ++j2) {
              const double o1 = -0.375 + 0.25 * j1, o2 = -0.375 + 0.25 * j2;
              double p[3] = {fc[a][0], fc[a][1], fc[a][2]};
              p[t1] += o1;
              p[t2] += o2;
              if (ccSampleExt(sdf, e, p[0], p[1], p[2]) > 0.0) {
                c1 += o1;
                c2 += o2;
                ++nOpen;
              }
            }
          if (nOpen == 0) {
            out[a](i) = 0.0;  // fully covered face: zero flux
            continue;
          }
          double p0[3] = {fc[a][0], fc[a][1], fc[a][2]};
          p0[t1] += c1 / nOpen;
          p0[t2] += c2 / nOpen;  // open-area centroid on the face
          double p1[3] = {p0[0], p0[1], p0[2]};
          p1[a] += dir;  // one cell toward the solid
          const double s0 = ccSampleExt(sdf, e, p0[0], p0[1], p0[2]);
          const double s1 = ccSampleExt(sdf, e, p1[0], p1[1], p1[2]);
          double d = (s0 > 0.0 && s0 - s1 > 1e-12) ? s0 / (s0 - s1) : 0.0;
          out[a](i) = d < 0.0 ? 0.0 : (d > 2.0 ? 2.0 : d);  // cap at 2 cells: never leaves the
          // fitted data range [theta, theta+2], so the stencil weights stay interpolation-bounded
        }
      });
}
 
// centerToFace with the wall-aware reconstruction at solid-bordering faces (see above). The SDF is
// the cell-centered field (ghosts filled); U/V/W ghosts must be filled first, as for centerToFace.
// useCen (mode 3): evaluate at the stored open-centroid wall distance xcx/xcy/xcz instead of the
// face center (the flux-consistent quadrature).
inline void centerToFaceWallAware(CCField uf, CCField vf, CCField wf, CCConst U, CCConst V,
                                  CCConst W, CCConst sdf, CCConst xcx, CCConst xcy, CCConst xcz,
                                  bool useCen, C3 e, int g) {
  (void)g;
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::center_to_face_wall", MD(space, {1, 1, 1}, {e.x, e.y, e.z}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        uf(i) = wallAwareFaceValue(U, sdf, i, sx, x, e.x, useCen ? xcx(i) : -1.0);
        vf(i) = wallAwareFaceValue(V, sdf, i, sy, y, e.y, useCen ? xcy(i) : -1.0);
        wf(i) = wallAwareFaceValue(W, sdf, i, sz, z, e.z, useCen ? xcz(i) : -1.0);
      });
}
 
// The O-WEIGHTED EXACT ADJOINT of the wall-aware constraint, applied to a cell field p:
// out(i) = sum over the (up to 6) axis faces f whose stencil references cell i of
// c_i(f)·o(f)·(p(f) - p(f-sa)) — i.e. Cᵀp for the constraint C = D·diag(o)·T. For interior fluid
// cells (o = 1, ½/½ weights) this is exactly the central difference ½(p(i+sa) - p(i-sa)).
//
// WHY the o-weighting is essential (learned the hard way): the wall-anchored T deliberately has
// row sums != 1 at cut faces (the no-slip anchor pulls them down), so the UNWEIGHTED transpose
// force Tᵀ·G·P no longer telescopes to zero over a periodic box — the fluid feels a net spurious
// pressure force proportional to the near-wall gradients, which feeds back through the flow and
// DIVERGES (the mode-3 runaway; capping the abscissa does not help because the weights were never
// the problem). With the exact adjoint the pressure does NO work on the constraint manifold
// ((u, Cᵀp) = (Cu, p) = 0), cutting the feedback unconditionally, and o·T·Tᵀ stays
// interpolation-bounded so the per-step projection remains contractive. Used (setFaceInterp(2/3))
// for BOTH the incremental predictor's -grad(P) and the projection's cell correction: at steady
// state the predictor gradient IS the pressure force, so it must be the adjoint of the constraint —
// upgrading T alone (mode 1) measurably WORSENS the drag.
inline void transposeGradWallAware(CCField out, CCConst p, CCConst sdf, CCConst o, CCConst xc,
                                   bool useCen, int axis, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::transpose_grad_wall", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const long sa = (axis == 0) ? sx : (axis == 1) ? sy : sz;
        const int eA = (axis == 0) ? e.x : (axis == 1) ? e.y : e.z;
        const int ca = (axis == 0) ? x : (axis == 1) ? y : z;
        double acc = 0.0;
        long idx[3];
        double w[3];
        for (int m = -2; m <= 3; ++m) {  // all faces whose stencil can reference this cell
          const int cf = ca + m;
          if (cf < 1 || cf >= eA)
            continue;
          const long f = i + (long)m * sa;
          const int n = wallAwareFaceStencil(sdf, f, sa, cf, eA, useCen ? xc(f) : -1.0, idx, w);
          for (int k = 0; k < n; ++k)
            if (idx[k] == i)
              acc += w[k] * o(f) * (p(f) - p(f - sa));
        }
        out(i) = acc;
      });
}
 
// Cell fluid volume fraction cs (mode 4): 1 in clear fluid, 0 in solid, subsampled (4x4x4 trilinear
// SDF) at cut cells. Static, one-time at setSolid. Used to weight the FV momentum time term and the
// body/pressure forcing so a cut cell is driven over its FLUID volume only (omitting cs is an O(h)
// surface error — part of what keeps modes 0-3 first-order).
inline void buildCellFraction(CCField cs, CCConst sdf, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::cell_fraction", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        const double sc = sdf(i);
        // fully clear one cell away on all axes -> not a cut cell; cheap exact classification
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const bool f0 = sc >= 0.0;
        const bool cut = (sdf(i + sx) >= 0.0) != f0 || (sdf(i - sx) >= 0.0) != f0 ||
                         (sdf(i + sy) >= 0.0) != f0 || (sdf(i - sy) >= 0.0) != f0 ||
                         (sdf(i + sz) >= 0.0) != f0 || (sdf(i - sz) >= 0.0) != f0;
        if (!cut) {
          cs(i) = sc >= 0.0 ? 1.0 : 0.0;
          return;
        }
        int nf = 0;
        for (int j0 = 0; j0 < 4; ++j0)
          for (int j1 = 0; j1 < 4; ++j1)
            for (int j2 = 0; j2 < 4; ++j2) {
              const double o0 = -0.375 + 0.25 * j0, o1 = -0.375 + 0.25 * j1,
                           o2 = -0.375 + 0.25 * j2;
              if (ccSampleExt(sdf, e, x + o0, y + o1, z + o2) > 0.0)
                ++nf;
            }
        cs(i) = nf / 64.0;
      });
}
 
// FINITE-VOLUME MOMENTUM OPERATOR (mode 4): applies L_FV(U) over the fluid control volume of each
// cell (the centroid wall-gradient is validated a priori in tests/study/fv_wallflux_apriori.py):
//   L_FV(U)_i = idt·cs_i·U_i + mu·[ Σ_f o_f·(U_i − U_nbr) + Σ_a W_a·g_a^centroid(U) ]
// The o_f-weighted two-point face fluxes + the fragment-centroid wall drag = the same finite-volume
// boundary geometry the projection uses. In the interior (cs=1, o_f=1, W=0) this is exactly the
// backward-Euler diffusion operator idt·U − mu·Lap(U) — identical to the IBM matrix M there, so the
// defect correction (M·u − L_FV·u) vanishes and interior cells stay byte-identical to mode 0.
inline void fvViscousApply(CCField Lu, CCConst U, CCConst sdf, CCConst cs, CCConst ox, CCConst oy,
                           CCConst oz, double mu, double idt, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::fv_viscous_apply", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long st[3] = {sx, sy, sz};
        const long i = (long)x + (long)y * sy + (long)z * sz;
        CCConst oa[3] = {ox, oy, oz};
        double W[3], of[3];  // fragment normal per axis; low/high face openness
        double diag = 0.0, offs = 0.0, aw = 0.0;
        for (int a = 0; a < 3; ++a) {
          const double om = oa[a](i), op = oa[a](i + st[a]);
          W[a] = om - op;
          of[a] = om + op;  // low + high face openness
          diag += om + op;
          offs += om * U(i - st[a]) + op * U(i + st[a]);
          aw += (W[a] < 0.0 ? -W[a] : W[a]);
        }
        double wall = 0.0;
        if (aw > 1e-12) {  // centroid-anchored wall drag mu·Σ_a W_a g_a (a-priori-validated:
                           // fv_wallflux_apriori.py)
          double nx = 0.5 * (sdf(i + sx) - sdf(i - sx));
          double ny = 0.5 * (sdf(i + sy) - sdf(i - sy));
          double nz = 0.5 * (sdf(i + sz) - sdf(i - sz));
          double nn = Kokkos::sqrt(nx * nx + ny * ny + nz * nz);
          if (nn > 1e-12) {
            nx /= nn;
            ny /= nn;
            nz /= nn;
            const double nv[3] = {nx, ny, nz};
            const double sdi = sdf(i);
            // foot point p* = x − sdi·n̂, clamped into the sampleable block [1, e−2] so the two
            // off-cell trilinear taps (p*+2σ) never index outside the padded field (guards a NaN/
            // large sdi too — ccSampleExt only clamps integer indices, not a NaN coordinate).
            auto cl = [](double v, double hi) { return v < 1.0 ? 1.0 : (v > hi ? hi : v); };
            const double px = cl(x - sdi * nx, e.x - 2.0);
            const double py = cl(y - sdi * ny, e.y - 2.0);
            const double pz = cl(z - sdi * nz, e.z - 2.0);
            for (int a = 0; a < 3; ++a) {
              const double sg = nv[a] >= 0.0 ? 1.0 : -1.0;
              const double u1 = ccSampleExt(U, e, px + (a == 0 ? sg : 0.0),
                                            py + (a == 1 ? sg : 0.0), pz + (a == 2 ? sg : 0.0));
              const double u2 =
                  ccSampleExt(U, e, px + (a == 0 ? 2.0 * sg : 0.0), py + (a == 1 ? 2.0 * sg : 0.0),
                              pz + (a == 2 ? 2.0 * sg : 0.0));
              wall += W[a] * sg * (2.0 * u1 - 0.5 * u2);
            }
          }
        }
        // FV viscous operator = idt·cs·U + μ·[Σo_f(U_i−U_nbr) − Σ_a W_a g_a].  The wall term sign
        // is −μ Σ W_a g_a: from ∫_CV −μ∇²u = −μ[flux_out − flux_in], the wall flux g_a enters with
        // a minus, so a resistive wall (∂u/∂n<0) ADDS to the operator (dissipative), matching the
        // interior −μLap. (+ would be anti-dissipative and blows the solve up.)
        Lu(i) = idt * cs(i) * U(i) + mu * (diag * U(i) - offs - wall);
        (void)of;
      });
}
 
// Lagrange quadratic through (-1,a1),(0,a2),(+1,a3) evaluated at xx (Basilisk embed.h `quadratic`).
KOKKOS_INLINE_FUNCTION double eQuad(double xx, double a1, double a2, double a3) {
  return (a1 * (xx - 1.0) + a3 * (xx + 1.0)) * xx * 0.5 - a2 * (xx - 1.0) * (xx + 1.0);
}
 
// TRUE-NORMAL wall gradient d(U)/dn at the embedded boundary of a cut cell — the Basilisk embed.h
// `dirichlet_gradient` (no-slip U_wall = 0), in cell units (h = 1). n̂ = unit inward-to-FLUID normal
// (∇sdf), p = boundary-point offset from the cell centre (cell units). Along the dominant-|n̂| axis
// it places two image points 1 and 2 cells into the fluid, interpolates U there by TRANSVERSE
// bi-quadratic interpolation of the cell-centred values (`eQuad`×`eQuad`), and fits the quadratic
// {0 at wall, v0 at d0, v1 at d1} for a 2nd-order derivative. Fallbacks (Basilisk): if the near
// image point's 3×3 transverse stencil straddles the wall, a biased-linear `embed_interpolate`
// anchored at the fluid home cell (uses only fluid cells); if even the home cell is solid, the
// degenerate 1-point estimate through the cell centre. A-priori-validated O(h²) with 0% degenerate
// on the Stokes sphere (tests/study/fv_wallflux_apriori.py, variant C). Interior/transverse reads
// of U and sdf must be halo-filled (the velocity block carries G=2 ghosts).
KOKKOS_INLINE_FUNCTION double embedDirichletGradient(CCConst U, CCConst sdf, C3 e, int x, int y,
                                                     int z, double nx, double ny, double nz,
                                                     double px, double py, double pz) {
  const double nv[3] = {nx, ny, nz}, pv[3] = {px, py, pz};
  const int ci[3] = {x, y, z}, ext[3] = {e.x, e.y, e.z};
  const long st[3] = {1, (long)e.x, (long)e.x * e.y};
  int da = 0;
  double best = Kokkos::fabs(nv[0]);
  if (Kokkos::fabs(nv[1]) > best) {
    best = Kokkos::fabs(nv[1]);
    da = 1;
  }
  if (Kokkos::fabs(nv[2]) > best)
    da = 2;
  const int t1 = (da + 1) % 3, t2 = (da + 2) % 3;
  const double sgn = nv[da] >= 0.0 ? 1.0 : -1.0;
  auto clampi = [](int v, int n) { return v < 0 ? 0 : (v >= n ? n - 1 : v); };
  auto rd = [](double v) {  // round to nearest integer, clamped to {-1,0,1} (Basilisk j/k)
    double r = v >= 0.0 ? Kokkos::floor(v + 0.5) : Kokkos::ceil(v - 0.5);
    return r > 1.0 ? 1.0 : (r < -1.0 ? -1.0 : r);
  };
  double vL[2], dL[2];
  bool defd[2];
  for (int l = 0; l < 2; ++l) {
    const int io = (l + 1) * (int)sgn;  // ±1, ±2 cells into the fluid along the dominant axis
    const double dl = ((double)io - pv[da]) / nv[da];  // distance (cells) from wall to image plane
    const double y1 = pv[t1] + dl * nv[t1], z1 = pv[t2] + dl * nv[t2];
    const int jj = (int)rd(y1), kk = (int)rd(z1);
    const double ly = y1 - jj, lz = z1 - kk;
    double vv[3][3];
    bool sf[3][3];
    for (int dk = -1; dk <= 1; ++dk)
      for (int dj = -1; dj <= 1; ++dj) {
        int cc[3];
        cc[da] = ci[da] + io;
        cc[t1] = ci[t1] + jj + dj;
        cc[t2] = ci[t2] + kk + dk;
        const long idx = (long)clampi(cc[0], ext[0]) * st[0] + (long)clampi(cc[1], ext[1]) * st[1] +
                         (long)clampi(cc[2], ext[2]) * st[2];
        vv[dj + 1][dk + 1] = U(idx);
        sf[dj + 1][dk + 1] = sdf(idx) >= 0.0;
      }
    bool full = true;
    for (int a2 = 0; a2 < 3; ++a2)
      for (int b2 = 0; b2 < 3; ++b2)
        full = full && sf[a2][b2];
    const bool home = sf[1][1];
    const double vbq =
        eQuad(lz, eQuad(ly, vv[0][0], vv[1][0], vv[2][0]), eQuad(ly, vv[0][1], vv[1][1], vv[2][1]),
              eQuad(ly, vv[0][2], vv[1][2], vv[2][2]));
    double vbl = vv[1][1];  // biased-linear embed_interpolate (fluid-side only), anchored at home
    {
      const int fp = ly >= 0.0 ? 2 : 0, fm = ly >= 0.0 ? 0 : 2;
      if (sf[fp][1])
        vbl += Kokkos::fabs(ly) * (vv[fp][1] - vv[1][1]);
      else if (sf[fm][1])
        vbl += Kokkos::fabs(ly) * (vv[1][1] - vv[fm][1]);
    }
    {
      const int fp = lz >= 0.0 ? 2 : 0, fm = lz >= 0.0 ? 0 : 2;
      if (sf[1][fp])
        vbl += Kokkos::fabs(lz) * (vv[1][fp] - vv[1][1]);
      else if (sf[1][fm])
        vbl += Kokkos::fabs(lz) * (vv[1][1] - vv[1][fm]);
    }
    vL[l] = full ? vbq : vbl;
    dL[l] = dl;
    defd[l] = full || home;
  }
  if (defd[0] && defd[1]) {  // 2-point quadratic fit (3rd-order value -> 2nd-order gradient)
    const double d0 = dL[0], d1 = dL[1], v0 = vL[0], v1 = vL[1];
    return (v0 * d1 / d0 - v1 * d0 / d1) / (d1 - d0);
  }
  if (defd[0])  // near image point only -> 1-point linear
    return vL[0] / dL[0];
  double d0 = Kokkos::fabs(pv[da] / nv[da]);  // degenerate sliver: 1-point through the cell centre
  // Floor d0 at the resolution scale. The degenerate estimate U/d0 enters the momentum march as an
  // EXPLICIT lagged wall flux (defect correction), whose per-step gain ~ mu*area/d0 against the row
  // diagonal mu*sum(o) + rho/dt must stay < 1: with area <= sqrt(3), d0 >= 0.5 bounds it by ~0.6 for
  // any (mu, rho, dt). The previous 1e-3 floor let a cell centre grazing the surface (|sdf|~1e-3,
  // guaranteed on multi-sphere beds) drive gains of O(10^2) -- the dt-independent x~100/step blowup
  // behind the mode-6/7 "CutcellMG preconditioner non-finite z" abort. (Basilisk's identical
  // dirichlet_gradient stencil is safe there because its home coefficient goes into the IMPLICIT
  // diagonal via *coef; this port applies it explicitly, so the floor carries the stability.)
  if (d0 < 0.5)
    d0 = 0.5;
  const long ii = (long)ci[0] * st[0] + (long)ci[1] * st[1] + (long)ci[2] * st[2];
  return U(ii) / d0;
}
 
// EMBED viscous operator (setFaceInterp(5)): identical to fvViscousApply except the wall drag is
// the TRUE-NORMAL Basilisk gradient (embedDirichletGradient) rather than the axis-by-axis
// fragment-normal estimate. Per cut cell the fragment area is |W| = |o_{a−}−o_{a+}|
// (divergence-theorem normal) and the wall flux is +μ·area·d(U)/dn — algebraically the −μ·(W·∇U) of
// fvViscousApply but with the O(h²) true-normal derivative in place of the O(h) axis
// reconstruction. Interior cells (area→0) reduce to idt·cs·U + μ(diag·U − offs) exactly, so the
// mode-0 defect vanishes there.
inline void embedViscousApply(CCField Lu, CCConst U, CCConst sdf, CCConst cs, CCConst ox,
                              CCConst oy, CCConst oz, double mu, double idt, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::embed_viscous_apply", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long st[3] = {sx, sy, sz};
        const long i = (long)x + (long)y * sy + (long)z * sz;
        CCConst oa[3] = {ox, oy, oz};
        double Wv[3];
        double diag = 0.0, offs = 0.0;
        for (int a = 0; a < 3; ++a) {
          const double om = oa[a](i), op = oa[a](i + st[a]);
          Wv[a] = om - op;
          diag += om + op;
          offs += om * U(i - st[a]) + op * U(i + st[a]);
        }
        const double area = Kokkos::sqrt(Wv[0] * Wv[0] + Wv[1] * Wv[1] + Wv[2] * Wv[2]);
        double wall = 0.0;
        if (area > 1e-12) {
          double nx = 0.5 * (sdf(i + sx) - sdf(i - sx));
          double ny = 0.5 * (sdf(i + sy) - sdf(i - sy));
          double nz = 0.5 * (sdf(i + sz) - sdf(i - sz));
          const double nn = Kokkos::sqrt(nx * nx + ny * ny + nz * nz);
          if (nn > 1e-12) {
            nx /= nn;
            ny /= nn;
            nz /= nn;
            const double sdi = sdf(i);  // foot-point (boundary centroid proxy) offset, cell units
            const double dudn = embedDirichletGradient(U, sdf, e, x, y, z, nx, ny, nz, -sdi * nx,
                                                       -sdi * ny, -sdi * nz);
            wall = area * dudn;  // +μ·area·dudn == −μ·(W·∇U) with the true-normal derivative
          }
        }
        Lu(i) = idt * cs(i) * U(i) + mu * (diag * U(i) - offs + wall);
      });
}
 
// 7-point stencil matvec y = M·u using the stored (rs-scaled) IBM operator coefficients. Used to
// form the mode-4 defect-correction RHS  b = M·u^k − rs·L_FV(u^k) + rs·b_FV, whose fixed point
// satisfies L_FV·u* = b_FV exactly, with M only the (stable, small-cell-safe) preconditioner.
inline void stencilMatvec(CCField y, CCConst u, MConst AC, MConst AW, MConst AE, MConst AS,
                          MConst AN, MConst AB, MConst AT, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::stencil_matvec", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y2, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y2 * sy + (long)z * sz;
        y(i) = (double)AC(i) * u(i) + (double)AW(i) * u(i - sx) + (double)AE(i) * u(i + sx) +
               (double)AS(i) * u(i - sy) + (double)AN(i) * u(i + sy) + (double)AB(i) * u(i - sz) +
               (double)AT(i) * u(i + sz);
      });
}
 
// u -= d over the inner cells (the mode-2 correction applies the transposeGradWallAware field).
inline void subtractField(CCField u, CCConst d, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::subtract_field", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        u(i) -= d(i);
      });
}
 
// OPENNESS-WEIGHTED cell pressure gradient along one axis (Basilisk centered_gradient / embed cell
// correction): out(i) = (o(i)·(p(i)−p(i−sa)) + o(i+sa)·(p(i+sa)−p(i))) / (o(i)+o(i+sa)). In the
// interior (o=1) this is the central difference ½(p(i+sa)−p(i−sa)); at a CUT cell with one closed
// face (o=0) it returns the OPEN face gradient at FULL weight (not the plain projectCorrectCenter's
// ½), which is the flux-consistent pressure force the approximate projection needs at the embedded
// boundary — the O(h) under-correction the plain map makes there (analysis defect (b)). Used for
// both the incremental −grad(P^n) predictor and the projection correction of the EMBED path (mode
// 6), so the two stay the openness-adjoint of the fs-weighted divergence constraint.
inline void centerGradOpen(CCField out, CCConst p, CCConst o, int axis, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::center_grad_open", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const long sa = (axis == 0) ? sx : (axis == 1) ? sy : sz;
        const double om = o(i), op = o(i + sa);
        out(i) = (om * (p(i) - p(i - sa)) + op * (p(i + sa) - p(i))) / (om + op + 1e-12);
      });
}
 
// ADJOINT-APERTURE cell pressure gradient along one axis (setFaceInterp(11)):
// out(i) = 1/2·(o(i)·(p(i)−p(i−sa)) + o(i+sa)·(p(i+sa)−p(i))) — centerGradOpen WITHOUT the
// normalization, which makes it the exact TRANSPOSE of the aperture divergence of the 1/2-1/2
// face average: G = −(D_α Π)^T. Two structural consequences (collocated_invisible_subspace.md):
// (i) support-consistent — it reads a solid-centred φ wherever the constraint couples it
// (α_f > 0), so the invisible pressure subspace of the gauge-exact/plain gradients collapses;
// (ii) adjoint — the dt→∞ Uzawa pressure map has an SPSD Schur complement, the stabilizable
// case (the normalized embed pair (modes 6/7) is non-adjoint and measured unconditionally
// unstable on beds, dt-free doubling ~75 steps). The price is the 1/2·α under-weighting of the
// pressure force at nearly-closed cut cells (accuracy measured on the ladder, not assumed).
inline void centerGradAperture(CCField out, CCConst p, CCConst o, int axis, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::center_grad_aperture", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const long sa = (axis == 0) ? sx : (axis == 1) ? sy : sz;
        out(i) = 0.5 * (o(i) * (p(i) - p(i - sa)) + o(i + sa) * (p(i + sa) - p(i)));
      });
}
 
// PER-CELL-RESCALED adjoint-aperture gradient (setFaceInterp(12)): centerGradAperture times the
// scalar S(i) = 6 / max(sum_a(o_a(i) + o_a(i+sa)), 0.5), one weight per CELL (not per axis).
// S == 1 in the bulk (sum = 6), so interior cells reproduce mode 11 / mode 0 exactly; at cut
// cells it restores (on average) the full-weight pressure force the 1/2*alpha under-weighting of
// the pure adjoint removes -- the mode-11 accuracy price -- while keeping the corrector a
// positive-DIAGONAL rescaling S*G of the adjoint pair (support unchanged; contrast mode 6's
// per-axis normalization, which is NOT such a rescaling and is measured unstable). The cap
// sum >= 0.5 bounds S <= 12 at nearly-closed cells.
inline void centerGradApertureScaled(CCField out, CCConst p, CCConst ox, CCConst oy, CCConst oz,
                                     int axis, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::center_grad_aperture_scaled", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const long sa = (axis == 0) ? sx : (axis == 1) ? sy : sz;
        CCConst o = (axis == 0) ? ox : (axis == 1) ? oy : oz;
        double osum = ox(i) + ox(i + sx) + oy(i) + oy(i + sy) + oz(i) + oz(i + sz);
        if (osum < 0.5)
          osum = 0.5;
        out(i) = (6.0 / osum) * 0.5 *
                 (o(i) * (p(i) - p(i - sa)) + o(i + sa) * (p(i + sa) - p(i)));
      });
}
 
// CAPPED per-axis normalized aperture gradient (setFaceInterp(13)): the mode-6 embed form
// (om*d- + op*d+)/(om+op) with the denominator FLOORED at omin -- equivalently the adjoint
// gradient centerGradAperture times the per-cell-per-COMPONENT diagonal S = 2/max(om+op, omin).
// Rationale: mode 6's floor is 1e-12, an unbounded diagonal gain (up to ~1e12) at nearly-closed
// axes, and the mode-6 pressure loop is measured unconditionally unstable on beds; mode 11
// (S = 1) is unconditionally clean but under-weights the cut-cell pressure force by ~(om+op)/2
// (k gap -11% at R=8). This kernel keeps the full-weight embed gradient wherever om+op >= omin
// (every ordinary cut cell) and only caps the gain where an axis is nearly closed.
inline void centerGradOpenCapped(CCField out, CCConst p, CCConst o, int axis, double omin, C3 e,
                                 int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::center_grad_open_capped", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const long sa = (axis == 0) ? sx : (axis == 1) ? sy : sz;
        const double om = o(i), op = o(i + sa);
        double den = om + op;
        if (den < omin)
          den = omin;
        out(i) = (om * (p(i) - p(i - sa)) + op * (p(i + sa) - p(i))) / den;
      });
}
 
// u,v,w -= openness-weighted cell pressure gradient of phi (centerGradOpen per axis). The embed
// analogue of projectCorrectCenter: the cut-cell cell velocity gets the full open-face pressure
// force.
inline void projectCorrectCenterOpen(CCField u, CCField v, CCField w, CCConst phi, CCConst ox,
                                     CCConst oy, CCConst oz, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::correct_center_open", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const double omx = ox(i), opx = ox(i + sx);
        const double omy = oy(i), opy = oy(i + sy);
        const double omz = oz(i), opz = oz(i + sz);
        u(i) -= (omx * (phi(i) - phi(i - sx)) + opx * (phi(i + sx) - phi(i))) / (omx + opx + 1e-12);
        v(i) -= (omy * (phi(i) - phi(i - sy)) + opy * (phi(i + sy) - phi(i))) / (omy + opy + 1e-12);
        w(i) -= (omz * (phi(i) - phi(i - sz)) + opz * (phi(i + sz) - phi(i))) / (omz + opz + 1e-12);
      });
}
 
// Cell-centered velocity correction u_c -= grad_c(phi), grad_c per axis = ½·(g⁻ + g⁺) of the two
// adjacent FACE pressure-gradients, but a face that is fully CLOSED (openness 0 — a solid
// neighbour) contributes a ZERO gradient instead of reading that neighbour's φ. Rationale:
//   * interior fluid cell (both faces open): ½(φᵢ₊₁-φᵢ₋₁) — the central difference, bulk unchanged;
//   * immersed cut cell (solid neighbour): the closed face's φ is DECOUPLED (≈0, AC≈0 in the
//   operator), so
//     reading it corrupts the gradient — zeroing that face uses only the fluid-side (open)
//     gradient;
//   * domain-BC wall (Neumann, φ-ghost = interior): the closed wall face truly has ∂φ/∂n≈0, and
//   zeroing it
//     gives ½·g_open — identical to the previous central difference with the Neumann ghost (no
//     change there).
// Cut faces (0<o<1) keep their real gradient (the neighbour is fluid). Only the cell field is
// touched; the projection's face divergence-free guarantee is unaffected. phi ghosts + face
// openness must be filled first.
inline void projectCorrectCenter(CCField u, CCField v, CCField w, CCConst phi, CCConst ox,
                                 CCConst oy, CCConst oz, C3 e, int g) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::correct_center", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const double gm_x = (ox(i) > 1e-12) ? (phi(i) - phi(i - sx)) : 0.0;
        const double gp_x = (ox(i + sx) > 1e-12) ? (phi(i + sx) - phi(i)) : 0.0;
        const double gm_y = (oy(i) > 1e-12) ? (phi(i) - phi(i - sy)) : 0.0;
        const double gp_y = (oy(i + sy) > 1e-12) ? (phi(i + sy) - phi(i)) : 0.0;
        const double gm_z = (oz(i) > 1e-12) ? (phi(i) - phi(i - sz)) : 0.0;
        const double gp_z = (oz(i + sz) > 1e-12) ? (phi(i + sz) - phi(i)) : 0.0;
        u(i) -= 0.5 * (gm_x + gp_x);
        v(i) -= 0.5 * (gm_y + gp_y);
        w(i) -= 0.5 * (gm_z + gp_z);
      });
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_MAC_APPROX_PROJECTION_HPP