scalar_mg.hpp

Go to the documentation of this file.

// ScalarMG: the multigrid preconditioner of the cut-cell scalar path (doc/scalar_ibm_design.md
// §5.2, D10; WO-4 single phase). One instance per cut-cell scalar; it preconditions the SPD
// lumped-probe SURROGATE S of §4.3 — the BiCGStab of `scalar_krylov.hpp` runs on the true probe
// operator, and one V-cycle of this class is its z = M^-1 r.
//
// Level table (§5.2): VelocityMG::init / initMpi's rule verbatim — in-place coarsening per axis
// while the global dimension is even and >= 4 and (distributed) every rank's block origin and size
// are even on that axis, then the shared aspect rule `CutcellMG::mgChooseRatio` (theta = 2, read
// from setMetric's h'). No depth cap (the "full table"), no telescoping (§13 Q3). Level 0 is the
// scalar's own G = 2 block (its ghost exchange is the solver's, injected as `Fill`); coarse levels
// carry g = 1 and their own `GridHalo<double>`.
//
// Per level the surrogate is stored in FACE form, double: AC (diagonal) and AFX/AFY/AFZ, AFX(i) the
// off-diagonal coefficient of cell i's LOW x face (mac_pressure.hpp's face form). Level 0 (§4.3):
// AC = the WO-3 surrogate diagonal SAC (bands' diagonal + the lumped wall term), AF = -Lam w_a a,
// guarded toward a cell that is not an unknown exactly as `sco::buildBands` (so the level-0 sweeps
// are bitwise the WO-3 band sweeps). Coarse level L, coarse cell C (REDISCRETIZED, §5.2):
//   * face:  t_a = w_a(L) <Lam a>, w_a(L) = w_a / cfac_a^2, <Lam a> the plain average of the fine
//            sub-faces' products, coarsened recursively by `coarsenOpenAvg`'s cell body;
//   * mass:  m_C = restrictAvg of the children's m (kappa idt, + the lumped implicit outflow
//            under advection, WO-5);
//   * wall:  W_C = (1/N_L) sum over the level-0 facets under C of alpha G(s_phi) — the plain
//   average
//            of the level-0 wall terms cw, at the FINE probe distance (design Amendment A1; the
//            level's own s_L overshoots the coarse correction by ~2^L and diverged, log WO-4).
//            Robin and Dirichlet alike (Neumann: 0); in WO-7 the conjugate coupling the same way.
//            A direct gather, one thread per coarse cell over its level-0 descendants in fixed
//            order;
//   * pins:  C is pinned iff every child is (the unknown flags restricted with max); pinned rows
//   are
//            identity rows with zero correction;
//   * Dirichlet domain faces: 2 w_a(L) <Lam a_bf>, the boundary plane averaged level by level.
// Ghosts beyond a non-periodic global face are pinned on every level (the record clears the level-0
// ghost unknown flags; a coarse ghost's children are those ghosts), so a coarse correction is 0
// there before the trilinear prolongation reads it, and every face coefficient toward them is 0.
//
// Smoother: red-black Gauss-Seidel, colour = (gx + gy + gz) mod 2 from GLOBAL indices, an exchange
// before each colour — on a single rank on a device backend (WO-9b; flow's A3) the colour passes
// and the residual of a level with even inner dimensions read the periodic neighbour through the
// wrapped index instead (`wrapReads`); a coarse level's fill is ONE launch (`smg::fillShell`), no
// fence. V-cycle: pre 2 sweeps R -> B, residual (after a fresh exchange), restrictAvg, (singular:
// the coarse rhs mean removed), recurse, coarse ghosts filled, prolongAdd (trilinear), pinned cells
// re-zeroed, post 2 sweeps B -> R. Bottom: 16 sweeps (8 R -> B then 8 B -> R, so the cycle stays
// symmetric) plus the mean removal when singular. The level rule (§5.2): transient with kappa_A = 1
// + 4 dt' D' sum_a w_a < kFullTableKappa = 25 uses level 0 alone (2 + 2 sweeps, the WO-3
// preconditioner); otherwise, and always when steady, the full table. The cycle is a fixed linear
// operator.
//
// Steady mode WITH advection (design Amendment A2, §6.7; `Inputs::advective`, the same on every
// rank): the cycle runs on the ADVECTIVE surrogate S_adv, in band form (AC, AW..AT). Level 0 is
// the operator's 7 bands with SAC as the diagonal (no new storage). A coarse level keeps everything
// above (faces, wall, Dirichlet planes, pins; its mass carries the open-face outflow omega_open,
// not the interior outflow) and adds the summed positive parts of the sub-face fluxes, Qp/Qm per
// low face per unit level volume (the piecewise-constant Galerkin FOU): AC += the cell's outgoing
// Q, AW = AFX - Qp_x, AE = AFX(i+e_x) - Qm_x(i+e_x), ... The sweep (unknown cells of the colour),
// the residual, applySurrogate and the contraction instrument use the bands; the cycle structure,
// transfers, colours and exchanges are unchanged. Transient mode and steady mode at rest take the
// face-form path above, bitwise unchanged.
#ifndef PECLET_FLOW_SCALAR_MG_HPP
#define PECLET_FLOW_SCALAR_MG_HPP
 
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <functional>
#include <Kokkos_Core.hpp>
#include <limits>
#include <memory>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
#include "mac_cutcell_mg.hpp"  // coarsenOpenAvgCell, restrictAvg, prolongAdd, CutcellMG::mgChooseRatio
#include "mac_pressure.hpp"  // cutcellSmoothFaceCell, cutcellResidualFaceCell, cutcellApplyFaceCell
#include "policy.hpp"
#include "scalar_cutcell_geometry.hpp"  // ScalarFacetOverlay
 
namespace peclet::flow {
 
namespace smg {
 
struct GuardedProduct {
  CCConst ap, unk;
  double lam = 0.0;
  long st = 1;
  KOKKOS_INLINE_FUNCTION double operator()(long i) const {
    return (unk(i) > 0.5 && unk(i - st) > 0.5) ? lam * ap(i) : 0.0;
  }
};
 
struct MassField {
  CCConst kappa, unk, out;
  double idt = 0.0;
  bool hasOut = false;
  KOKKOS_INLINE_FUNCTION double operator()(long i) const {
    if (!(unk(i) > 0.5))
      return 0.0;
    return hasOut ? kappa(i) * idt + out(i) : kappa(i) * idt;
  }
};
 
template <class FV>
inline void restrictMax(CCField coarse, const FV& fine, C3 cext, C3 fext, int gc, int gf, C3 cinner,
                        C3 ratio) {
  ccFor3(
      "peclet::flow::smg_restrict_max", C3{0, 0, 0}, cinner,
      KOKKOS_LAMBDA(int icx, int icy, int icz) {
        const long fsy = fext.x, fsz = (long)fext.x * fext.y;
        double m = 0.0;
        for (int dz = 0; dz < ratio.z; ++dz)
          for (int dy = 0; dy < ratio.y; ++dy)
            for (int dx = 0; dx < ratio.x; ++dx) {
              const long fi = (long)(ratio.x * icx + dx + gf) +
                              (long)(ratio.y * icy + dy + gf) * fsy +
                              (long)(ratio.z * icz + dz + gf) * fsz;
              m = fine(fi) > m ? fine(fi) : m;
            }
        coarse((long)(icx + gc) + (long)(icy + gc) * cext.x +
               (long)(icz + gc) * (long)cext.x * cext.y) = m > 0.5 ? 1.0 : 0.0;
      });
}
 
inline void zeroPinned(CCField x, CCConst unk, C3 e, int g) {
  ccFor3(
      "peclet::flow::smg_zero_pinned", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
        if (!(unk(i) > 0.5))
          x(i) = 0.0;
      });
}
 
inline void zeroPlanes(CCField f, C3 e, int a, int lo, int hi) {
  if (hi <= lo)
    return;
  const int ext[3] = {e.x, e.y, e.z};
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long st[3] = {1, (long)e.x, (long)e.x * e.y};
  const long sa = st[a], sb = st[b], sc = st[c];
  const int nb = ext[b], nc = ext[c];
  CCExec space;
  auto row = KOKKOS_LAMBDA(long t) {  // one transverse row (pb, pc) through the planes
    const int pb = (int)(t % nb), pc = (int)(t / nb);
    for (int k = lo; k < hi; ++k)
      f((long)k * sa + (long)pb * sb + (long)pc * sc) = 0.0;
  };
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    if (hostRunSerial((long)nb * nc * (hi - lo))) {
      for (long t = 0; t < (long)nb * nc; ++t)
        row(t);
      return;
    }
  }
  Kokkos::parallel_for("peclet::flow::smg_zero_planes",
                       Kokkos::RangePolicy<CCExec>(space, 0, (long)nb * nc), row);
}
 
inline void fillShell(CCField f, C3 e, C3 inner, int g, const bool per[3]) {
  CCExec space;
  const int nx = inner.x, ny = inner.y, nz = inner.z;
  const long exy = (long)e.x * e.y;
  const bool px = per[0], py = per[1], pz = per[2];
  CCField ff = f;
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    // host launch rule: rows (one index computation per row, a contiguous x loop inside)
    const long rZ = 2L * g * e.y, rY = (long)nz * 2 * g, rX = (long)nz * ny;
    auto row = [=](long t) {
      long dst, src;
      bool zero;
      if (t < rZ) {  // z-ghost plane row (y, zg), all x
        const int l = (int)(t / e.y), y = (int)(t - (long)l * e.y);
        const int z = l < g ? l : nz + l;
        const bool yg = y < g || y >= g + ny;
        const int sy = y < g ? y + ny : (y >= g + ny ? y - ny : y);
        const int sz = z < g ? z + nz : z - nz;
        zero = !pz || (yg && !py);
        dst = (long)y * e.x + (long)z * exy;
        src = (long)sy * e.x + (long)sz * exy;
      } else if (t < rZ + rY) {  // y-ghost row (yg, z) of the inner z range, all x
        const long t2 = t - rZ;
        const int zi = (int)(t2 / (2 * g)), l = (int)(t2 - (long)zi * 2 * g);
        const int y = l < g ? l : ny + l, z = g + zi;
        const int sy = y < g ? y + ny : y - ny;
        zero = !py;
        dst = (long)y * e.x + (long)z * exy;
        src = (long)sy * e.x + (long)z * exy;
      } else {  // the x-ghost cells of the inner row (y, z)
        const long t3 = t - rZ - rY;
        const int zi = (int)(t3 / ny), yi = (int)(t3 - (long)zi * ny);
        const long base = (long)(g + yi) * e.x + (long)(g + zi) * exy;
        for (int x = 0; x < g; ++x)
          ff(base + x) = px ? ff(base + x + nx) : 0.0;
        for (int x = g + nx; x < e.x; ++x)
          ff(base + x) = px ? ff(base + x - nx) : 0.0;
        return;
      }
      if (zero) {
        for (int x = 0; x < e.x; ++x)
          ff(dst + x) = 0.0;
        return;
      }
      for (int x = 0; x < g; ++x)
        ff(dst + x) = px ? ff(src + x + nx) : 0.0;
      for (int x = g; x < g + nx; ++x)
        ff(dst + x) = ff(src + x);
      for (int x = g + nx; x < e.x; ++x)
        ff(dst + x) = px ? ff(src + x - nx) : 0.0;
    };
    const long cells = 2L * g * exy + (long)nz * 2 * g * e.x + (long)nz * ny * 2 * g;
    if (hostRunSerial(cells)) {
      for (long t = 0; t < rZ + rY + rX; ++t)
        row(t);
      return;
    }
    Kokkos::parallel_for("peclet::flow::smg_fill_shell",
                         Kokkos::RangePolicy<CCExec>(space, 0, rZ + rY + rX), row);
    return;
  }
  const long nZ = 2L * g * exy, nY = (long)nz * 2 * g * e.x, nX = (long)nz * ny * 2 * g;
  Kokkos::parallel_for(
      "peclet::flow::smg_fill_shell", Kokkos::RangePolicy<CCExec>(space, 0, nZ + nY + nX),
      KOKKOS_LAMBDA(long t) {
        int x, y, z;
        if (t < nZ) {
          const int l = (int)(t / exy);
          const long r = t - (long)l * exy;
          x = (int)(r % e.x);
          y = (int)(r / e.x);
          z = l < g ? l : nz + l;
        } else if (t < nZ + nY) {
          const long t2 = t - nZ, rw = 2L * g * e.x;
          const long zi = t2 / rw, r = t2 - zi * rw;
          const int l = (int)(r / e.x);
          x = (int)(r % e.x);
          y = l < g ? l : ny + l;
          z = g + (int)zi;
        } else {
          const long t3 = t - nZ - nY, rw = (long)ny * 2 * g;
          const long zi = t3 / rw, r = t3 - zi * rw;
          const int l = (int)(r % (2 * g));
          x = l < g ? l : nx + l;
          y = g + (int)(r / (2 * g));
          z = g + (int)zi;
        }
        const bool gx = x < g || x >= g + nx, gy = y < g || y >= g + ny, gz = z < g || z >= g + nz;
        const long d = (long)x + (long)y * e.x + (long)z * exy;
        if ((gx && !px) || (gy && !py) || (gz && !pz)) {
          ff(d) = 0.0;
          return;
        }
        const int sx = x < g ? x + nx : (x >= g + nx ? x - nx : x);
        const int sy = y < g ? y + ny : (y >= g + ny ? y - ny : y);
        const int sz = z < g ? z + nz : (z >= g + nz ? z - nz : z);
        ff(d) = ff((long)sx + (long)sy * e.x + (long)sz * exy);
      });
}
 
struct Nb {
  long xp, xm, yp, ym, zp, zm;
};
KOKKOS_INLINE_FUNCTION Nb nbrs(int lx, int ly, int lz, C3 n, int g, long i, long sy, long sz,
                               int wrap) {
  Nb w;
  w.xm = ((wrap & 1) && lx == g) ? i + (long)(n.x - 1) : i - 1;
  w.xp = ((wrap & 1) && lx == g + n.x - 1) ? i - (long)(n.x - 1) : i + 1;
  w.ym = ((wrap & 2) && ly == g) ? i + (long)(n.y - 1) * sy : i - sy;
  w.yp = ((wrap & 2) && ly == g + n.y - 1) ? i - (long)(n.y - 1) * sy : i + sy;
  w.zm = ((wrap & 4) && lz == g) ? i + (long)(n.z - 1) * sz : i - sz;
  w.zp = ((wrap & 4) && lz == g + n.z - 1) ? i - (long)(n.z - 1) * sz : i + sz;
  return w;
}
 
KOKKOS_INLINE_FUNCTION Nb nbrsOn(int lx, int ly, int lz, C3 n, int g, long i, long sy, long sz,
                                 int wrap) {
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    (void)lx, (void)ly, (void)lz, (void)n, (void)g, (void)wrap;
    return Nb{i + 1, i - 1, i + sy, i - sy, i + sz, i - sz};
  } else {
    return nbrs(lx, ly, lz, n, g, i, sy, sz, wrap);
  }
}
 
template <class OpV>
inline void sweepColorFace(CCField x, CCConst rhs, OpV AC, OpV AFX, OpV AFY, OpV AFZ, C3 e, C3 og,
                           int g, int color, C3 n, int wrap) {
  CCExec space;
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    const int nyi = e.y - 2 * g, nzi = e.z - 2 * g;
    const long cells = (long)nyi * nzi * (e.x - 2 * g);
    auto pencil = KOKKOS_LAMBDA(long t) {
      const int ly = g + (int)(t % nyi), lz = g + (int)(t / nyi);
      const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
      const int P = (color + og.x + og.y + ly + og.z + lz) & 1;
      PECLET_FLOW_OMP_SIMD  // same-colour cells are independent (ccFor3's contract)
          for (int lx = g + ((P ^ (g & 1)) & 1); lx < e.x - g; lx += 2) {
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        cutcellSmoothFaceCell(x, rhs, AC, AFX, AFY, AFZ, i, sx, sy, sz, w.xp, w.xm, w.yp, w.ym,
                              w.zp, w.zm);
      }
    };
    if (hostRunSerial(cells)) {  // coarse MG level: the fork/join costs more than the sweep
      for (long t = 0; t < (long)nyi * nzi; ++t)
        pencil(t);
      return;
    }
    Kokkos::parallel_for("peclet::flow::smg_smooth",
                         Kokkos::RangePolicy<CCExec>(space, 0, (long)nyi * nzi), pencil);
    return;
  }
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::smg_smooth", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        if (((og.x + lx + og.y + ly + og.z + lz) & 1) != color)
          return;
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        cutcellSmoothFaceCell(x, rhs, AC, AFX, AFY, AFZ, i, sx, sy, sz, w.xp, w.xm, w.yp, w.ym,
                              w.zp, w.zm);
      });
}
 
inline void residualFace(CCField r, CCConst x, CCConst b, CCField AC, CCField AFX, CCField AFY,
                         CCField AFZ, C3 e, int g, C3 n, int wrap) {
  ccFor3(
      "peclet::flow::smg_residual", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        cutcellResidualFaceCell(r, x, b, AC, AFX, AFY, AFZ, i, sx, sy, sz, w.xp, w.xm, w.yp, w.ym,
                                w.zp, w.zm);
      });
}
 
// ---- the advective path (design Amendment A2, §6.7): band form ----------------------------------
 
struct GuardedPart {
  CCConst phi, unk;
  long st = 1;
  double sgn = 1.0;
  KOKKOS_INLINE_FUNCTION double operator()(long i) const {
    return (unk(i) > 0.5 && unk(i - st) > 0.5) ? Kokkos::fmax(sgn * phi(i), 0.0) : 0.0;
  }
};
 
template <class XV, class BV>
KOKKOS_INLINE_FUNCTION double bandOff(const XV& x, const BV& AW, const BV& AE, const BV& AS,
                                      const BV& AN, const BV& AB, const BV& AT, long i,
                                      const Nb& w) {
  return ((AW(i) * x(w.xm) + AE(i) * x(w.xp)) + (AS(i) * x(w.ym) + AN(i) * x(w.yp))) +
         (AB(i) * x(w.zm) + AT(i) * x(w.zp));
}
 
template <class BV, class UV>
inline void sweepColorBand(CCField x, CCConst rhs, BV AC, BV AW, BV AE, BV AS, BV AN, BV AB, BV AT,
                           UV unk, C3 e, C3 og, int g, int color, C3 n, int wrap) {
  CCExec space;
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    const int nyi = e.y - 2 * g, nzi = e.z - 2 * g;
    const long cells = (long)nyi * nzi * (e.x - 2 * g);
    auto pencil = KOKKOS_LAMBDA(long t) {
      const int ly = g + (int)(t % nyi), lz = g + (int)(t / nyi);
      const long sy = e.x, sz = (long)e.x * e.y;
      const int P = (color + og.x + og.y + ly + og.z + lz) & 1;
      for (int lx = g + ((P ^ (g & 1)) & 1); lx < e.x - g; lx += 2) {
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        if (!(unk(i) > 0.5))
          continue;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        x(i) = (rhs(i) - bandOff(x, AW, AE, AS, AN, AB, AT, i, w)) / AC(i);
      }
    };
    if (hostRunSerial(cells)) {  // coarse MG level: the fork/join costs more than the sweep
      for (long t = 0; t < (long)nyi * nzi; ++t)
        pencil(t);
      return;
    }
    Kokkos::parallel_for("peclet::flow::smg_smooth_band",
                         Kokkos::RangePolicy<CCExec>(space, 0, (long)nyi * nzi), pencil);
    return;
  }
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::smg_smooth_band", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        if (((og.x + lx + og.y + ly + og.z + lz) & 1) != color)
          return;
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        if (!(unk(i) > 0.5))
          return;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        x(i) = (rhs(i) - bandOff(x, AW, AE, AS, AN, AB, AT, i, w)) / AC(i);
      });
}
 
template <class XV, class OV>
KOKKOS_INLINE_FUNCTION double faceOff(const XV& x, const OV& AFX, const OV& AFY, const OV& AFZ,
                                      long i, long sx, long sy, long sz, const Nb& w) {
  return AFX(i + sx) * x(w.xp) + AFX(i) * x(w.xm) + AFY(i + sy) * x(w.yp) + AFY(i) * x(w.ym) +
         AFZ(i + sz) * x(w.zp) + AFZ(i) * x(w.zm);
}
 
template <bool Band, class XV, class BV>
KOKKOS_INLINE_FUNCTION void twoPhaseCell(const XV& x, const XV& rhs, const BV& AC, const BV& AFX,
                                         const BV& AFY, const BV& AFZ, const BV& AW, const BV& AE,
                                         const BV& AS, const BV& AN, const BV& AB, const BV& AT,
                                         const XV& xs, const XV& rhss, const BV& ACs,
                                         const BV& AFXs, const BV& AFYs, const BV& AFZs,
                                         const BV& W, const CCConst& unkF, long i, long sx, long sy,
                                         long sz, const Nb& w) {
  const double wc = W(i);
  const double af = AC(i), as = ACs(i);
  if (wc > 0.0) {
    const double nf = Band ? bandOff(x, AW, AE, AS, AN, AB, AT, i, w)
                           : faceOff(x, AFX, AFY, AFZ, i, sx, sy, sz, w);
    const double rf = rhs(i) - nf;
    const double rs = rhss(i) - faceOff(xs, AFXs, AFYs, AFZs, i, sx, sy, sz, w);
    const double det = af * as - wc * wc;
    if (det > 1e-300 * af * as) {
      x(i) = (as * rf + wc * rs) / det;
      xs(i) = (wc * rf + af * rs) / det;
    } else {
      x(i) = (rf + wc * xs(i)) / af;
      xs(i) = (rs + wc * x(i)) / as;
    }
    return;
  }
  if constexpr (Band) {
    if (unkF(i) > 0.5)
      x(i) = (rhs(i) - bandOff(x, AW, AE, AS, AN, AB, AT, i, w)) / af;
  } else {
    if (!(af < 1e-30))
      x(i) = (rhs(i) - faceOff(x, AFX, AFY, AFZ, i, sx, sy, sz, w)) / af;
  }
  if (!(as < 1e-30))
    xs(i) = (rhss(i) - faceOff(xs, AFXs, AFYs, AFZs, i, sx, sy, sz, w)) / as;
}
 
template <bool Band>
inline void sweepColorTwo(CCField x, CCField rhs, CCField AC, CCField AFX, CCField AFY, CCField AFZ,
                          CCField AW, CCField AE, CCField AS, CCField AN, CCField AB, CCField AT,
                          CCField xs, CCField rhss, CCField ACs, CCField AFXs, CCField AFYs,
                          CCField AFZs, CCField W, CCConst unkF, C3 e, C3 og, int g, int color,
                          C3 n, int wrap) {
  CCExec space;
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    const int nyi = e.y - 2 * g, nzi = e.z - 2 * g;
    const long cells = (long)nyi * nzi * (e.x - 2 * g);
    auto pencil = KOKKOS_LAMBDA(long t) {
      const int ly = g + (int)(t % nyi), lz = g + (int)(t / nyi);
      const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
      const int P = (color + og.x + og.y + ly + og.z + lz) & 1;
      for (int lx = g + ((P ^ (g & 1)) & 1); lx < e.x - g; lx += 2) {
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        twoPhaseCell<Band>(x, rhs, AC, AFX, AFY, AFZ, AW, AE, AS, AN, AB, AT, xs, rhss, ACs, AFXs,
                           AFYs, AFZs, W, unkF, i, sx, sy, sz, w);
      }
    };
    if (hostRunSerial(cells)) {
      for (long t = 0; t < (long)nyi * nzi; ++t)
        pencil(t);
      return;
    }
    Kokkos::parallel_for("peclet::flow::smg_smooth_two",
                         Kokkos::RangePolicy<CCExec>(space, 0, (long)nyi * nzi), pencil);
    return;
  }
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::smg_smooth_two", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        if (((og.x + lx + og.y + ly + og.z + lz) & 1) != color)
          return;
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        twoPhaseCell<Band>(x, rhs, AC, AFX, AFY, AFZ, AW, AE, AS, AN, AB, AT, xs, rhss, ACs, AFXs,
                           AFYs, AFZs, W, unkF, i, sx, sy, sz, w);
      });
}
 
template <bool Band>
inline void residualTwo(CCField r, CCConst x, CCConst b, CCField AC, CCField AFX, CCField AFY,
                        CCField AFZ, CCField AW, CCField AE, CCField AS, CCField AN, CCField AB,
                        CCField AT, CCField rs, CCConst xs, CCConst bs, CCField ACs, CCField AFXs,
                        CCField AFYs, CCField AFZs, CCField W, C3 e, int g, C3 n, int wrap) {
  ccFor3(
      "peclet::flow::smg_residual_two", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        const double ax = Band ? AC(i) * x(i) + bandOff(x, AW, AE, AS, AN, AB, AT, i, w)
                               : AC(i) * x(i) + faceOff(x, AFX, AFY, AFZ, i, sx, sy, sz, w);
        const double axs = ACs(i) * xs(i) + faceOff(xs, AFXs, AFYs, AFZs, i, sx, sy, sz, w);
        const double wc = W(i);
        r(i) = (b(i) - ax) + wc * xs(i);
        rs(i) = (bs(i) - axs) + wc * x(i);
      });
}
 
inline void residualBand(CCField r, CCConst x, CCConst b, CCField AC, CCField AW, CCField AE,
                         CCField AS, CCField AN, CCField AB, CCField AT, C3 e, int g, C3 n,
                         int wrap) {
  ccFor3(
      "peclet::flow::smg_residual_band", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)lx + (long)ly * sy + (long)lz * sz;
        const Nb w = nbrsOn(lx, ly, lz, n, g, i, sy, sz, wrap);
        r(i) = b(i) - (AC(i) * x(i) + bandOff(x, AW, AE, AS, AN, AB, AT, i, w));
      });
}
 
}  // namespace smg
 
class ScalarMG {
 public:
  static constexpr int G0 = 2;  
  static constexpr int GC = 1;  
  static constexpr double kFullTableKappa = 25.0;
  using Fill = std::function<void(CCField)>;
 
  struct Level {
    C3 ext{0, 0, 0}, inner{0, 0, 0}, ratio{1, 1, 1}, cfac{1, 1, 1};
    C3 og{0, 0, 0};    
    C3 gdim{0, 0, 0};  
    int g = GC;
    std::size_t n = 0;
    CCField x, rhs, res;        
    CCField AC, AFX, AFY, AFZ;  
    CCField unk;                
    CCField mass, px, py, pz;   
    Kokkos::View<double*, CCMem> plane[6];  
    double nUnk = 0.0;                      
    CCField AW, AE, AS, AN, AB, AT;
    CCField qp[3], qm[3];
    CCField xs, rhss, ress, ACs, AFXs, AFYs, AFZs, unks, masss, pxs, pys, pzs, W;
    CCField wallG, wallCG;
    double nUnkS = 0.0;  
#ifdef PECLET_FLOW_MPI
    std::shared_ptr<GridHaloTopology<3>> halo;
    std::shared_ptr<GridHalo<double>> dev;
#endif
  };
 
  void setMetric(const double w[3], const double hp[3]) {
    for (int a = 0; a < 3; ++a) {
      w_[a] = w[a];
      hp_[a] = hp[a];
    }
    aniso_ = !(hp_[0] == 1.0 && hp_[1] == 1.0 && hp_[2] == 1.0);
  }
  void setPeriodic(const bool per[3]) {
    for (int a = 0; a < 3; ++a)
      per_[a] = per[a];
  }
 
  void setLevel0Wrap(bool on) { fill0Wraps_ = on; }
 
  void init(int nx, int ny, int nz, Fill fill0) {
    lv_.clear();
    distributed_ = false;
    fill0_ = std::move(fill0);
    C3 inner{nx, ny, nz}, cf{1, 1, 1};
    for (int L = 0; L < kMaxLevels; ++L) {
      Level v;
      v.g = L == 0 ? G0 : GC;
      v.inner = inner;
      v.gdim = inner;
      v.ext = C3{inner.x + 2 * v.g, inner.y + 2 * v.g, inner.z + 2 * v.g};
      v.cfac = cf;
      v.n = (std::size_t)v.ext.x * v.ext.y * v.ext.z;
      const bool canA[3] = {can(inner.x), can(inner.y), can(inner.z)};
      const C3 ratio = chooseRatio(canA, cf);
      const C3 next{ratio.x == 2 ? inner.x / 2 : inner.x, ratio.y == 2 ? inner.y / 2 : inner.y,
                    ratio.z == 2 ? inner.z / 2 : inner.z};
      v.ratio = ratio;
      allocate(v, L);
      lv_.push_back(v);
      if (next.x == inner.x && next.y == inner.y && next.z == inner.z)
        break;
      inner = next;
      cf = C3{cf.x * ratio.x, cf.y * ratio.y, cf.z * ratio.z};
    }
  }
#ifdef PECLET_FLOW_MPI
  void initMpi(const peclet::core::decomp::BlockDecomposer<3>& dec0, MPI_Comm comm, C3 og0,
               Fill fill0) {
    lv_.clear();
    distributed_ = true;
    comm_ = comm;
    fill0_ = std::move(fill0);
    int rank = 0;
    MPI_Comm_rank(comm, &rank);
    std::array<bool, 3> per{true, true, true};  // periodic everywhere; non-periodic ghosts zeroed
    const auto& g0 = dec0.globalSize();
    C3 gs{(int)g0[0], (int)g0[1], (int)g0[2]}, cf{1, 1, 1};
    auto evenOn = [](const peclet::core::decomp::BlockDecomposer<3>& d, int ax) {
      for (std::size_t b = 0; b < d.numBlocks(); ++b) {
        const auto blk = d.block(b);
        if (blk.origin[ax] % 2 != 0 || blk.size[ax] % 2 != 0)
          return false;
      }
      return true;
    };
    peclet::core::decomp::BlockDecomposer<3> dec = dec0;
    for (int L = 0; L < kMaxLevels; ++L) {
      Level v;
      v.g = L == 0 ? G0 : GC;
      if (L == 0) {
        const auto blk = dec.block((std::size_t)rank);
        v.inner = C3{(int)blk.size[0], (int)blk.size[1], (int)blk.size[2]};
        v.ext = C3{v.inner.x + 2 * G0, v.inner.y + 2 * G0, v.inner.z + 2 * G0};
        v.og = og0;
      } else {
        v.halo = std::make_shared<GridHaloTopology<3>>();
        v.halo->buildTopology(dec, rank, v.g, per, comm);
        v.dev = std::make_shared<GridHalo<double>>();
        v.dev->init(*v.halo);
        v.dev->setLabel("smg L" + std::to_string(L));
        const auto& idx = v.halo->indexer();
        const auto eg = idx.sizeInclGhost(), ino = idx.sizeInner(), oig = idx.originInclGhost();
        v.ext = {(int)eg[0], (int)eg[1], (int)eg[2]};
        v.inner = {(int)ino[0], (int)ino[1], (int)ino[2]};
        v.og = {(int)oig[0] + v.g, (int)oig[1] + v.g, (int)oig[2] + v.g};
      }
      v.gdim = gs;
      v.cfac = cf;
      v.n = (std::size_t)v.ext.x * v.ext.y * v.ext.z;
      const bool canA[3] = {can(gs.x) && evenOn(dec, 0), can(gs.y) && evenOn(dec, 1),
                            can(gs.z) && evenOn(dec, 2)};
      const C3 ratio = chooseRatio(canA, cf);
      const C3 next{ratio.x == 2 ? gs.x / 2 : gs.x, ratio.y == 2 ? gs.y / 2 : gs.y,
                    ratio.z == 2 ? gs.z / 2 : gs.z};
      v.ratio = ratio;
      allocate(v, L);
      lv_.push_back(v);
      if (next.x == gs.x && next.y == gs.y && next.z == gs.z)
        break;
      gs = next;
      cf = C3{cf.x * ratio.x, cf.y * ratio.y, cf.z * ratio.z};
      dec = dec.coarsened(peclet::core::IVec<3>{ratio.x, ratio.y, ratio.z});
    }
  }
#endif
 
  int levels() const { return (int)lv_.size(); }
  int levelsUsed() const { return nUse_; }
  const Level& level(int L) const { return lv_[L]; }
  std::vector<C3> levelRatios() const {
    std::vector<C3> r;
    for (const auto& v : lv_)
      r.push_back(v.ratio);
    return r;
  }
 
  struct Inputs {
    double lam = 0.0;       
    double idt = 0.0;       
    bool fullTable = true;  
    bool singular = false;  
    CCField SAC;            
    CCConst kappa, unknown, sax, say, saz;
    const scg::ScalarFacetOverlay* fac = nullptr;
    Kokkos::View<const double*, CCMem> facetW;  
    bool dirFace[6] = {false, false, false, false, false, false};
    CCConst outflow;
    bool hasOutflow = false;
    bool advective = false;
    CCField AW, AE, AS, AN, AB, AT;  
    CCConst phi[3];                  
    CCConst omegaOpen;               
    bool twoPhase = false;
    CCField SACs, W0;
    CCConst sunk, spx, spy, spz, ms0;
    Kokkos::View<const double*, CCMem> facetC;
  };
 
  void build(const Inputs& in) {
    singular_ = in.singular;
    adv_ = in.advective;
    nUse_ = in.fullTable ? (int)lv_.size() : 1;
    Level& l0 = lv_[0];
    l0.AC = in.SAC;
    if (adv_) {  // A2: the operator's bands, no new storage
      l0.AW = in.AW;
      l0.AE = in.AE;
      l0.AS = in.AS;
      l0.AN = in.AN;
      l0.AB = in.AB;
      l0.AT = in.AT;
    }
    l0.unk = CCField();  // level 0's flag is the geometry's (const) view
    unk0_ = in.unknown;
    buildLevel0Faces(in);
    l0.nUnk = countUnknown(0);
    two_ = in.twoPhase;
    if (two_) {  // WO-7: the solid phase and the coupling on level 0
      l0.ACs = in.SACs;
      l0.W = in.W0;
      unkS0_ = in.sunk;
      buildLevel0SolidFaces(in);
      l0.nUnkS = countUnknownS(0);
    }
    for (int f = 0; f < 6; ++f)
      dirFace_[f] = in.dirFace[f];
    if (nUse_ > 1) {
      const bool newLists = buildCutLists(in);
      wallReuse_ = !newLists && wallKeep(in);
      for (int f = 0; f < 6; ++f)
        if (dirFace_[f] && touches(l0, f))
          buildPlane0(in, f);
      for (int L = 1; L < nUse_; ++L)
        buildCoarse(in, L);
    }
  }
 
  void apply(CCField z, CCField r) {
    Level& l0 = lv_[0];
    l0.x = z;
    l0.rhs = r;
    Kokkos::deep_copy(CCExec(), z, 0.0);
    if (nUse_ == 1) {
      for (int k = 0; k < 2; ++k) {
        sweepColor(0, 0);
        sweepColor(0, 1);
      }
      for (int k = 0; k < 2; ++k) {
        sweepColor(0, 1);
        sweepColor(0, 0);
      }
      return;
    }
    vcycle(0);
  }
 
  void applyTwo(CCField zf, CCField zs, CCField rf, CCField rs) {
    Level& l0 = lv_[0];
    l0.x = zf;
    l0.rhs = rf;
    l0.xs = zs;
    l0.rhss = rs;
    Kokkos::deep_copy(CCExec(), zf, 0.0);
    Kokkos::deep_copy(CCExec(), zs, 0.0);
    if (nUse_ == 1) {
      for (int k = 0; k < 2; ++k) {
        sweepColor(0, 0);
        sweepColor(0, 1);
      }
      for (int k = 0; k < 2; ++k) {
        sweepColor(0, 1);
        sweepColor(0, 0);
      }
      return;
    }
    vcycleTwo(0);
  }
 
  void applySurrogateTwo(CCField yf, CCField ys, CCField xf, CCField xs) {
    Level& l0 = lv_[0];
    fill(0, xf);
    fill(0, xs);
    const C3 e = l0.ext;
    const int g = l0.g;
    CCField AC = l0.AC, AFX = l0.AFX, AFY = l0.AFY, AFZ = l0.AFZ, ACs = l0.ACs, AFXs = l0.AFXs,
            AFYs = l0.AFYs, AFZs = l0.AFZs, W = l0.W, AW = l0.AW, AE = l0.AE, AS = l0.AS,
            AN = l0.AN, AB = l0.AB, AT = l0.AT;
    const bool band = adv_;
    CCConst x = xf, s = xs;
    ccFor3(
        "peclet::flow::smg_apply_two", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
          const long i = (long)lx + (long)ly * sy + (long)lz * sz;
          const smg::Nb w{i + sx, i - sx, i + sy, i - sy, i + sz, i - sz};
          const double ax = band ? AC(i) * x(i) + smg::bandOff(x, AW, AE, AS, AN, AB, AT, i, w)
                                 : AC(i) * x(i) + smg::faceOff(x, AFX, AFY, AFZ, i, sx, sy, sz, w);
          const double axs = ACs(i) * s(i) + smg::faceOff(s, AFXs, AFYs, AFZs, i, sx, sy, sz, w);
          yf(i) = ax - W(i) * s(i);
          ys(i) = axs - W(i) * x(i);
        });
  }
 
  void applySurrogate(CCField y, CCField x) {
    Level& l0 = lv_[0];
    fill(0, x);
    const C3 e = l0.ext;
    const int g = l0.g;
    if (adv_) {  // A2: S_adv at level 0 = the operator's band matvec with SAC as the diagonal
      applyCutcellOp(y, CCConst(x), l0.AC, l0.AW, l0.AE, l0.AS, l0.AN, l0.AB, l0.AT, e, g);
      return;
    }
    CCField AC = l0.AC, AFX = l0.AFX, AFY = l0.AFY, AFZ = l0.AFZ;
    CCConst xx = x;
    ccFor3(
        "peclet::flow::smg_apply", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
          const long i = (long)lx + (long)ly * sy + (long)lz * sz;
          y(i) = cutcellApplyFaceCell(xx, AC, AFX, AFY, AFZ, i, sx, sy, sz, i + sx, i - sx, i + sy,
                                      i - sy, i + sz, i - sz);
        });
  }
 
  std::vector<double> contraction(int iters, CCField e, CCField t, CCField z) {
    Level& l0 = lv_[0];
    const C3 ext = l0.ext, og = l0.og;
    const int g = l0.g;
    CCConst unk = unk0_;
    ccFor3(
        "peclet::flow::smg_seed", C3{g, g, g}, C3{ext.x - g, ext.y - g, ext.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
          std::uint64_t h = (std::uint64_t)(og.x + lx - g) * 0x9E3779B97F4A7C15ull ^
                            (std::uint64_t)(og.y + ly - g) * 0xC2B2AE3D27D4EB4Full ^
                            (std::uint64_t)(og.z + lz - g) * 0x165667B19E3779F9ull;
          h ^= h >> 29;
          h *= 0xBF58476D1CE4E5B9ull;
          h ^= h >> 32;
          e(i) = unk(i) > 0.5 ? (double)(h >> 11) * (1.0 / 9007199254740992.0) - 0.5 : 0.0;
        });
    if (singular_)
      removeMean(0, e);
    std::vector<double> ratios;
    double prev = norm2(e);
    for (int k = 0; k < iters; ++k) {
      applySurrogate(t, e);
      apply(z, t);
      if (singular_)
        removeMean(0, z);
      ccFor3(
          "peclet::flow::smg_err_update", C3{g, g, g}, C3{ext.x - g, ext.y - g, ext.z - g},
          KOKKOS_LAMBDA(int lx, int ly, int lz) {
            const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
            e(i) -= z(i);
          });
      if (singular_)  // keep the iterate off the null space (round-off would otherwise grow there)
        removeMean(0, e);
      const double cur = norm2(e);
      ratios.push_back(prev > 0.0 ? cur / prev : 0.0);
      if (!(cur > 0.0))
        break;
      // renormalize (keeps the iterate O(1) over many iterations)
      const double s = 1.0 / cur;
      ccFor3(
          "peclet::flow::smg_err_scale", C3{g, g, g}, C3{ext.x - g, ext.y - g, ext.z - g},
          KOKKOS_LAMBDA(int lx, int ly, int lz) {
            const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
            e(i) *= s;
          });
      prev = 1.0;
    }
    return ratios;
  }
 
  // ---- (public for nvcc extended lambdas) ----
  void vcycle(int L) {
    Level& lv = lv_[L];
    if (L + 1 == nUse_) {  // bottom: 16 sweeps, symmetric order
      for (int k = 0; k < kBottom / 2; ++k)
        sweep(L, true);
      for (int k = 0; k < kBottom / 2; ++k)
        sweep(L, false);
      if (singular_)
        removeMean(L, lv.x);
      return;
    }
    for (int k = 0; k < kPre; ++k)
      sweep(L, true);
    // the smoother leaves the ghosts one colour stale: a fresh fill, or the fused reads
    const int wr = wrapReads(L);
    if (wr < 0)
      fill(L, lv.x);
    if (adv_)
      smg::residualBand(lv.res, CCConst(lv.x), CCConst(lv.rhs), lv.AC, lv.AW, lv.AE, lv.AS, lv.AN,
                        lv.AB, lv.AT, lv.ext, lv.g, lv.inner, wr < 0 ? 0 : wr);
    else
      smg::residualFace(lv.res, CCConst(lv.x), CCConst(lv.rhs), lv.AC, lv.AFX, lv.AFY, lv.AFZ,
                        lv.ext, lv.g, lv.inner, wr < 0 ? 0 : wr);
    Level& cs = lv_[L + 1];
    restrictZeroX(cs, lv, cs.rhs, lv.res, cs.x);
    if (singular_)
      removeMean(L + 1, cs.rhs);
    vcycle(L + 1);
    fill(L + 1, cs.x);  // coarse ghosts: exchanged, 0 beyond a non-periodic global face
    prolongAdd(lv.x, CCConst(cs.x), lv.ext, cs.ext, lv.g, cs.g, lv.inner, lv.ratio);
    smg::zeroPinned(lv.x, unkOf(L), lv.ext, lv.g);
    for (int k = 0; k < kPost; ++k)
      sweep(L, false);
  }
 
  void vcycleTwo(int L) {
    Level& lv = lv_[L];
    if (L + 1 == nUse_) {
      for (int k = 0; k < kBottom / 2; ++k)
        sweep(L, true);
      for (int k = 0; k < kBottom / 2; ++k)
        sweep(L, false);
      if (singular_)
        removeMeanTwo(L, lv.x, lv.xs);
      return;
    }
    for (int k = 0; k < kPre; ++k)
      sweep(L, true);
    const int wr = wrapReads(L);
    if (wr < 0) {
      fill(L, lv.x);
      fill(L, lv.xs);
    }
    const int w = wr < 0 ? 0 : wr;
    if (adv_)
      smg::residualTwo<true>(lv.res, CCConst(lv.x), CCConst(lv.rhs), lv.AC, lv.AFX, lv.AFY, lv.AFZ,
                             lv.AW, lv.AE, lv.AS, lv.AN, lv.AB, lv.AT, lv.ress, CCConst(lv.xs),
                             CCConst(lv.rhss), lv.ACs, lv.AFXs, lv.AFYs, lv.AFZs, lv.W, lv.ext,
                             lv.g, lv.inner, w);
    else
      smg::residualTwo<false>(lv.res, CCConst(lv.x), CCConst(lv.rhs), lv.AC, lv.AFX, lv.AFY, lv.AFZ,
                              lv.AW, lv.AE, lv.AS, lv.AN, lv.AB, lv.AT, lv.ress, CCConst(lv.xs),
                              CCConst(lv.rhss), lv.ACs, lv.AFXs, lv.AFYs, lv.AFZs, lv.W, lv.ext,
                              lv.g, lv.inner, w);
    Level& cs = lv_[L + 1];
    restrictZeroX(cs, lv, cs.rhs, lv.res, cs.x);
    restrictZeroX(cs, lv, cs.rhss, lv.ress, cs.xs);
    if (singular_)
      removeMeanTwo(L + 1, cs.rhs, cs.rhss);
    vcycleTwo(L + 1);
    fill(L + 1, cs.x);
    fill(L + 1, cs.xs);
    prolongAdd(lv.x, CCConst(cs.x), lv.ext, cs.ext, lv.g, cs.g, lv.inner, lv.ratio);
    prolongAdd(lv.xs, CCConst(cs.xs), lv.ext, cs.ext, lv.g, cs.g, lv.inner, lv.ratio);
    smg::zeroPinned(lv.x, unkOf(L), lv.ext, lv.g);
    smg::zeroPinned(lv.xs, unkSOf(L), lv.ext, lv.g);
    for (int k = 0; k < kPost; ++k)
      sweep(L, false);
  }
 
  void restrictZeroX(Level& cs, Level& lv, CCField crhs, CCField fres, CCField cx) {
    if (!distributed_) {
      restrictAvgZeroX(crhs, cx, CCConst(fres), cs.ext, lv.ext, cs.g, lv.g, cs.inner, lv.ratio);
      return;
    }
    restrictAvg(crhs, CCConst(fres), cs.ext, lv.ext, cs.g, lv.g, cs.inner, lv.ratio);
    Kokkos::deep_copy(CCExec(), cx, 0.0);
  }
 
  void sweepColor(int L, int rb) {
    Level& lv = lv_[L];
    // WO-9b: on a single-rank level with every inner dimension even the pass reads the periodic
    // neighbours directly (wrapSmooth), so no fill; otherwise the fill before each colour
    const int ws = wrapSmooth(L);
    const int w = ws < 0 ? 0 : ws;
    if (ws < 0)
      fill(L, lv.x);
    const int kc = (rb + 3 * lv.g) & 1;
    if (two_) {  // WO-7: the coupled two-field update (one more exchange: the solid's)
      if (ws < 0)
        fill(L, lv.xs);
      if (adv_)
        smg::sweepColorTwo<true>(lv.x, lv.rhs, lv.AC, lv.AFX, lv.AFY, lv.AFZ, lv.AW, lv.AE, lv.AS,
                                 lv.AN, lv.AB, lv.AT, lv.xs, lv.rhss, lv.ACs, lv.AFXs, lv.AFYs,
                                 lv.AFZs, lv.W, unkOf(L), lv.ext, lv.og, lv.g, kc, lv.inner, w);
      else
        smg::sweepColorTwo<false>(lv.x, lv.rhs, lv.AC, lv.AFX, lv.AFY, lv.AFZ, lv.AW, lv.AE, lv.AS,
                                  lv.AN, lv.AB, lv.AT, lv.xs, lv.rhss, lv.ACs, lv.AFXs, lv.AFYs,
                                  lv.AFZs, lv.W, unkOf(L), lv.ext, lv.og, lv.g, kc, lv.inner, w);
      return;
    }
    if (adv_) {
      smg::sweepColorBand(lv.x, CCConst(lv.rhs), lv.AC, lv.AW, lv.AE, lv.AS, lv.AN, lv.AB, lv.AT,
                          unkOf(L), lv.ext, lv.og, lv.g, kc, lv.inner, w);
      return;
    }
    smg::sweepColorFace(lv.x, CCConst(lv.rhs), lv.AC, lv.AFX, lv.AFY, lv.AFZ, lv.ext, lv.og, lv.g,
                        kc, lv.inner, w);
  }
 
  void sweep(int L, bool fwd) {
    sweepColor(L, fwd ? 0 : 1);
    sweepColor(L, fwd ? 1 : 0);
  }
 
  void fill(int L, CCField f) {
    if (L == 0) {
      fill0_(f);
      return;
    }
    Level& lv = lv_[L];
#ifdef PECLET_FLOW_MPI
    if (distributed_) {
      lv.dev->exchange(f);
      zeroNonPeriodicGhosts(lv, f);
      return;
    }
#endif
    smg::fillShell(f, lv.ext, lv.inner, lv.g, per_);  // WO-9b: one launch, no fence
  }
 
  int wrapReads(int L) const {
    // Device backends only: on a host backend the per-cell wrapped-index arithmetic costs more than
    // the fill it saves — it keeps the x-row loop from vectorizing (128^3 bed advance, 1 thread:
    // the smoother 700 -> 1131 ms and the residual 80 -> 124 ms per advance, against 44 ms of
    // level-0 fills), so host passes fill first (`fill` is one launch per coarse level anyway).
    if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>)
      return -1;
    if (distributed_)
      return -1;
    if (L == 0)
      return fill0Wraps_ ? 7 : -1;
    return (per_[0] ? 1 : 0) | (per_[1] ? 2 : 0) | (per_[2] ? 4 : 0);
  }
  int wrapSmooth(int L) const {
    const Level& lv = lv_[L];
    if (lv.inner.x % 2 != 0 || lv.inner.y % 2 != 0 || lv.inner.z % 2 != 0)
      return -1;
    return wrapReads(L);
  }
 
  bool touches(const Level& lv, int f) const {
    const int a = f / 2;
    const int o = (a == 0) ? lv.og.x : (a == 1) ? lv.og.y : lv.og.z;
    const int n = (a == 0) ? lv.inner.x : (a == 1) ? lv.inner.y : lv.inner.z;
    const int gn = (a == 0) ? lv.gdim.x : (a == 1) ? lv.gdim.y : lv.gdim.z;
    return (f % 2 == 0) ? (o == 0) : (o + n == gn);
  }
 
 public:  // nvcc: a member enclosing an extended device lambda must be public (data stays private)
  static constexpr int kMaxLevels = 64;  
  static constexpr int kPre = 2, kPost = 2, kBottom = 16;
 
  static bool can(int d) { return (d % 2 == 0) && (d / 2 >= 2); }
  C3 chooseRatio(const bool canA[3], C3 cf) const {
    const double H[3] = {hp_[0] * (double)cf.x, hp_[1] * (double)cf.y, hp_[2] * (double)cf.z};
    return CutcellMG::mgChooseRatio(H, canA, aniso_, kAspectTheta);
  }
  static constexpr double kAspectTheta = 2.0;
 
  void allocate(Level& v, int L) {
    if (L == 0) {
      v.res = CCField("smg_res0", v.n);
      v.AFX = CCField("smg_afx0", v.n);
      v.AFY = CCField("smg_afy0", v.n);
      v.AFZ = CCField("smg_afz0", v.n);
      return;
    }
    v.x = CCField("smg_x", v.n);
    v.rhs = CCField("smg_rhs", v.n);
    v.res = CCField("smg_res", v.n);
    v.AC = CCField("smg_ac", v.n);
    v.AFX = CCField("smg_afx", v.n);
    v.AFY = CCField("smg_afy", v.n);
    v.AFZ = CCField("smg_afz", v.n);
    v.unk = CCField("smg_unk", v.n);
    v.mass = CCField("smg_mass", v.n);
    v.px = CCField("smg_px", v.n);
    v.py = CCField("smg_py", v.n);
    v.pz = CCField("smg_pz", v.n);
  }
 
  CCConst unkOf(int L) const { return L == 0 ? unk0_ : CCConst(lv_[L].unk); }
  CCConst unkSOf(int L) const { return L == 0 ? unkS0_ : CCConst(lv_[L].unks); }
 
  void allocateSolid(Level& v, int L) {
    if (v.ress.extent(0) == v.n)
      return;
    v.ress = CCField("smg_ress", v.n);
    v.AFXs = CCField("smg_afxs", v.n);
    v.AFYs = CCField("smg_afys", v.n);
    v.AFZs = CCField("smg_afzs", v.n);
    if (L == 0)
      return;
    v.xs = CCField("smg_xs", v.n);
    v.rhss = CCField("smg_rhss", v.n);
    v.ACs = CCField("smg_acs", v.n);
    v.unks = CCField("smg_unks", v.n);
    v.masss = CCField("smg_masss", v.n);
    v.pxs = CCField("smg_pxs", v.n);
    v.pys = CCField("smg_pys", v.n);
    v.pzs = CCField("smg_pzs", v.n);
    v.W = CCField("smg_w", v.n);
  }
 
  double countUnknownS(int L) const {
    const Level& lv = lv_[L];
    const C3 e = lv.ext;
    const int g = lv.g;
    CCConst unk = unkSOf(L);
    double s = 0.0;
    ccReduce3(
        "peclet::flow::smg_count_s", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz, double& acc) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          acc += unk(i) > 0.5 ? 1.0 : 0.0;
        },
        Kokkos::Sum<double>(s));
    allSum(&s, 1);
    return s;
  }
 
  void removeMeanTwo(int L, CCField f, CCField s) {
    const Level& lv = lv_[L];
    const double nTot = lv.nUnk + lv.nUnkS;
    if (!(nTot > 0.0))
      return;
    const C3 e = lv.ext;
    const int g = lv.g;
    CCConst uf = unkOf(L), us = unkSOf(L);
    double sum = 0.0;
    ccReduce3(
        "peclet::flow::smg_mean_two", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz, double& acc) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          acc += (uf(i) > 0.5 ? f(i) : 0.0) + (us(i) > 0.5 ? s(i) : 0.0);
        },
        Kokkos::Sum<double>(sum));
    allSum(&sum, 1);
    const double mean = sum / nTot;
    ccFor3(
        "peclet::flow::smg_mean_two_sub", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          if (uf(i) > 0.5)
            f(i) -= mean;
          if (us(i) > 0.5)
            s(i) -= mean;
        });
  }
 
  void buildLevel0SolidFaces(const Inputs& in) {
    Level& l0 = lv_[0];
    allocateSolid(l0, 0);
    const C3 e = l0.ext;
    const double wx = w_[0], wy = w_[1], wz = w_[2];
    CCConst px = in.spx, py = in.spy, pz = in.spz;
    CCField AFX = l0.AFXs, AFY = l0.AFYs, AFZ = l0.AFZs;
    ccFor3(
        "peclet::flow::smg_faces0_solid", C3{0, 0, 0}, e, KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
          AFX(i) = -(wx * px(i));
          AFY(i) = -(wy * py(i));
          AFZ(i) = -(wz * pz(i));
        });
  }
 
  void buildCoarseSolid(const Inputs& in, int L) {
    Level& c = lv_[L];
    Level& fin = lv_[L - 1];
    const C3 ratio = fin.ratio;
    if (L == 1)
      smg::restrictMax(c.unks, in.sunk, c.ext, fin.ext, c.g, fin.g, c.inner, ratio);
    else
      smg::restrictMax(c.unks, CCConst(fin.unks), c.ext, fin.ext, c.g, fin.g, c.inner, ratio);
    c.nUnkS = countUnknownS(L);
    {
      CCField px = c.pxs, py = c.pys, pz = c.pzs;
      const C3 ce = c.ext, fe = fin.ext;
      const int gc = c.g, gf = fin.g;
      const CCConst fx = L == 1 ? in.spx : CCConst(fin.pxs),
                    fy = L == 1 ? in.spy : CCConst(fin.pys),
                    fz = L == 1 ? in.spz : CCConst(fin.pzs);
      ccFor3(
          "peclet::flow::smg_coarsen_faces_solid", C3{0, 0, 0}, c.inner,
          KOKKOS_LAMBDA(int icx, int icy, int icz) {
            coarsenOpenAvgCell(px, py, pz, fx, fy, fz, ce, fe, gc, gf, ratio, icx, icy, icz);
          });
      fill(L, c.pxs);
      fill(L, c.pys);
      fill(L, c.pzs);
      const CCField P[3] = {c.pxs, c.pys, c.pzs};
      for (int a = 0; a < 3; ++a)
        if (!per_[a] && touches(c, 2 * a))
          smg::zeroPlanes(P[a], c.ext, a, 0, c.g + 1);
    }
    restrictAvg(c.masss, L == 1 ? in.ms0 : CCConst(fin.masss), c.ext, fin.ext, c.g, fin.g, c.inner,
                ratio);
    const C3 e = c.ext, ci = c.inner, cf = c.cfac;
    const int g = c.g;
    const double wx = w_[0] / ((double)cf.x * cf.x), wy = w_[1] / ((double)cf.y * cf.y),
                 wz = w_[2] / ((double)cf.z * cf.z);
    CCField AC = c.ACs, AFX = c.AFXs, AFY = c.AFYs, AFZ = c.AFZs, W = c.W;
    CCConst px = c.pxs, py = c.pys, pz = c.pzs, unk = c.unks, mass = c.masss;
    ccFor3(
        "peclet::flow::smg_coarse_faces_solid", C3{0, 0, 0}, e, KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
          AFX(i) = -(wx * px(i));
          AFY(i) = -(wy * py(i));
          AFZ(i) = -(wz * pz(i));
        });
    ccFor3(
        "peclet::flow::smg_coarse_diag_solid", C3{0, 0, 0}, ci,
        KOKKOS_LAMBDA(int icx, int icy, int icz) {
          const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
          const long i = (long)(icx + g) + (long)(icy + g) * sy + (long)(icz + g) * sz;
          if (!(unk(i) > 0.5)) {
            AC(i) = 1.0;
            AFX(i) = AFY(i) = AFZ(i) = 0.0;
            W(i) = 0.0;
            return;
          }
          const double tw = wx * px(i), te = wx * px(i + sx);
          const double ts = wy * py(i), tn = wy * py(i + sy);
          const double tb = wz * pz(i), tt = wz * pz(i + sz);
          AC(i) = (mass(i) + (((tw + te) + (ts + tn)) + (tb + tt))) + W(i);
        });
  }
 
  void zeroNonPeriodicGhosts(Level& lv, CCField f) {
    const int ext[3] = {lv.ext.x, lv.ext.y, lv.ext.z};
    for (int a = 0; a < 3; ++a) {
      if (per_[a])
        continue;
      if (touches(lv, 2 * a))
        smg::zeroPlanes(f, lv.ext, a, 0, lv.g);
      if (touches(lv, 2 * a + 1))
        smg::zeroPlanes(f, lv.ext, a, ext[a] - lv.g, ext[a]);
    }
  }
 
  // ---- reductions ----
  void allSum(double* v, int k) const {
#ifdef PECLET_FLOW_MPI
    if (distributed_) {
      std::vector<double> o((std::size_t)k);
      MPI_Allreduce(v, o.data(), k, MPI_DOUBLE, MPI_SUM, comm_);
      for (int j = 0; j < k; ++j)
        v[j] = o[(std::size_t)j];
    }
#endif
    (void)v;
    (void)k;
  }
  double countUnknown(int L) const {
    const Level& lv = lv_[L];
    const C3 e = lv.ext;
    const int g = lv.g;
    CCConst unk = unkOf(L);
    double s = 0.0;
    ccReduce3(
        "peclet::flow::smg_count", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz, double& acc) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          acc += unk(i) > 0.5 ? 1.0 : 0.0;
        },
        Kokkos::Sum<double>(s));
    allSum(&s, 1);
    return s;
  }
  void removeMean(int L, CCField f) {
    const Level& lv = lv_[L];
    if (!(lv.nUnk > 0.0))
      return;
    const C3 e = lv.ext;
    const int g = lv.g;
    CCConst unk = unkOf(L);
    double s = 0.0;
    ccReduce3(
        "peclet::flow::smg_mean", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz, double& acc) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          if (unk(i) > 0.5)
            acc += f(i);
        },
        Kokkos::Sum<double>(s));
    allSum(&s, 1);
    const double mean = s / lv.nUnk;
    ccFor3(
        "peclet::flow::smg_mean_sub", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          if (unk(i) > 0.5)
            f(i) -= mean;
        });
  }
  double norm2(CCField f) const {
    const Level& lv = lv_[0];
    const C3 e = lv.ext;
    const int g = lv.g;
    double s = 0.0;
    ccReduce3(
        "peclet::flow::smg_norm", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz, double& acc) {
          const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
          acc += f(i) * f(i);
        },
        Kokkos::Sum<double>(s));
    allSum(&s, 1);
    return std::sqrt(s);
  }
 
  // ---- build ----
  void buildLevel0Faces(const Inputs& in) {
    Level& l0 = lv_[0];
    const C3 e = l0.ext;
    const double lx = in.lam * w_[0], ly = in.lam * w_[1], lz = in.lam * w_[2];
    CCConst unk = in.unknown, sax = in.sax, say = in.say, saz = in.saz;
    CCField AFX = l0.AFX, AFY = l0.AFY, AFZ = l0.AFZ;
    ccFor3(
        "peclet::flow::smg_faces0", C3{0, 0, 0}, e, 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 bool u = unk(i) > 0.5;
          const double tw = (x > 0 && u && unk(i - sx) > 0.5) ? lx * sax(i) : 0.0;
          const double ts = (y > 0 && u && unk(i - sy) > 0.5) ? ly * say(i) : 0.0;
          const double tb = (z > 0 && u && unk(i - sz) > 0.5) ? lz * saz(i) : 0.0;
          AFX(i) = -tw;
          AFY(i) = -ts;
          AFZ(i) = -tb;
        });
  }
 
  bool buildCutLists(const Inputs& in) {
    const auto& fo = *in.fac;
    if (cutKey_ == fo.cutCell.data() && cutKeyN_ == fo.nCut && (int)cutStart_.size() == nUse_)
      return false;
    const Level& l0 = lv_[0];
    const C3 e0 = l0.ext, n0 = l0.inner;
    const int g0 = l0.g;
    auto hc = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), fo.cutCell);
    // the inner cut rows by level-0 linear index (the rows of one cell are unique)
    std::vector<std::pair<long, int>> rows;
    rows.reserve((std::size_t)fo.nCut);
    for (long r = 0; r < fo.nCut; ++r) {
      const long i0 = hc(r);
      const int x = (int)(i0 % e0.x) - g0, y = (int)((i0 / e0.x) % e0.y) - g0,
                z = (int)(i0 / ((long)e0.x * e0.y)) - g0;
      if (x >= 0 && x < n0.x && y >= 0 && y < n0.y && z >= 0 && z < n0.z)
        rows.emplace_back(i0, (int)r);
    }
    std::sort(rows.begin(), rows.end());
    cutStart_.assign((std::size_t)nUse_, Kokkos::View<int*, CCMem>());
    cutRows_.assign((std::size_t)nUse_, Kokkos::View<int*, CCMem>());
    for (int L = 1; L < nUse_; ++L) {
      const Level& c = lv_[L];
      const C3 ci = c.inner, cf = c.cfac;
      const long nc = (long)ci.x * ci.y * ci.z;
      std::vector<int> start((std::size_t)nc + 1, 0), list(rows.size());
      auto coarseOf = [&](long i0) {
        const int x = (int)(i0 % e0.x) - g0, y = (int)((i0 / e0.x) % e0.y) - g0,
                  z = (int)(i0 / ((long)e0.x * e0.y)) - g0;
        return (long)(x / cf.x) + (long)(y / cf.y) * ci.x + (long)(z / cf.z) * ci.x * ci.y;
      };
      for (const auto& q : rows)
        ++start[(std::size_t)coarseOf(q.first) + 1];
      for (long k = 0; k < nc; ++k)
        start[(std::size_t)k + 1] += start[(std::size_t)k];
      std::vector<int> pos(start.begin(), start.end() - 1);
      for (const auto& q : rows)  // stable: increasing level-0 index within each coarse cell
        list[(std::size_t)pos[(std::size_t)coarseOf(q.first)]++] = q.second;
      Kokkos::View<int*, CCMem> ds("smg_cut_start", (std::size_t)nc + 1),
          dl("smg_cut_rows", std::max<std::size_t>(list.size(), 1));
      Kokkos::deep_copy(
          ds, Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
                  start.data(), start.size()));
      if (!list.empty())
        Kokkos::deep_copy(
            Kokkos::subview(dl, std::make_pair((std::size_t)0, list.size())),
            Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
                list.data(), list.size()));
      cutStart_[(std::size_t)L] = ds;
      cutRows_[(std::size_t)L] = dl;
    }
    cutKey_ = fo.cutCell.data();
    cutKeyN_ = fo.nCut;
    return true;
  }
 
  bool wallKeep(const Inputs& in) {
    const auto& fo = *in.fac;
    const bool two = in.twoPhase;
    bool keep = wallValid_ && wallTwo_ == two && wallUse_ == nUse_ &&
                wallUnk_ == in.unknown.data() && cwSnap_.extent(0) == (std::size_t)fo.n &&
                (!two || ccSnap_.extent(0) == (std::size_t)fo.n);
    if (keep && fo.n > 0) {
      auto a = in.facetW, as = Kokkos::View<const double*, CCMem>(cwSnap_);
      auto c = two ? in.facetC : in.facetW,
           cs = two ? Kokkos::View<const double*, CCMem>(ccSnap_) : as;
      long diff = 0;
      Kokkos::parallel_reduce(
          "peclet::flow::smg_wall_same", Kokkos::RangePolicy<CCExec>(CCExec(), 0, fo.n),
          KOKKOS_LAMBDA(const long f, long& d) {
            auto same = [](double u, double w) {
              return u == w && Kokkos::signbit(u) == Kokkos::signbit(w);
            };
            if (!same(a(f), as(f)) || (two && !same(c(f), cs(f))))
              ++d;
          },
          diff);
      keep = diff == 0;
    }
    if (!keep) {  // the gather below runs: snapshot its inputs
      if (cwSnap_.extent(0) != (std::size_t)fo.n)
        cwSnap_ = Kokkos::View<double*, CCMem>("smg_cw_snap", (std::size_t)fo.n);
      Kokkos::deep_copy(CCExec(), cwSnap_, in.facetW);
      if (two) {
        if (ccSnap_.extent(0) != (std::size_t)fo.n)
          ccSnap_ = Kokkos::View<double*, CCMem>("smg_cc_snap", (std::size_t)fo.n);
        Kokkos::deep_copy(CCExec(), ccSnap_, in.facetC);
      }
      wallValid_ = true;
      wallTwo_ = two;
      wallUse_ = nUse_;
      wallUnk_ = in.unknown.data();
    }
    return keep;
  }
 
  void buildPlane0(const Inputs& in, int f) {
    Level& l0 = lv_[0];
    const int a = f / 2, side = f % 2;
    int t1, t2;
    tangents(a, t1, t2);
    const int ext[3] = {l0.ext.x, l0.ext.y, l0.ext.z};
    const long st[3] = {1, (long)l0.ext.x, (long)l0.ext.x * l0.ext.y};
    const int g = l0.g;
    const int n1 = ext[t1] - 2 * g, n2 = ext[t2] - 2 * g;
    const long sa = st[a], s1 = st[t1], s2 = st[t2];
    const int aInner = side == 0 ? g : ext[a] - g - 1;
    if (l0.plane[f].extent(0) != (std::size_t)n1 * n2)
      l0.plane[f] = Kokkos::View<double*, CCMem>("smg_plane0", (long)n1 * n2);
    auto pl = l0.plane[f];
    CCConst unk = in.unknown;
    CCConst saA = a == 0 ? in.sax : (a == 1 ? in.say : in.saz);
    const double lam = in.lam;
    CCExec space;
    Kokkos::parallel_for(
        "peclet::flow::smg_plane0", MDRange2<CCExec>(space, {0, 0}, {n1, n2}),
        KOKKOS_LAMBDA(int j1, int j2) {
          const long i = (long)aInner * sa + (long)(j1 + g) * s1 + (long)(j2 + g) * s2;
          const double abf = side == 0 ? saA(i) : saA(i + sa);
          pl((long)j1 + (long)j2 * n1) = unk(i) > 0.5 ? lam * abf : 0.0;
        });
  }
 
  static void tangents(int a, int& t1, int& t2) {
    t1 = (a == 0) ? 1 : 0;
    t2 = (a == 2) ? 1 : 2;
  }
 
  void buildCoarse(const Inputs& in, int L) {
    Level& c = lv_[L];
    Level& fin = lv_[L - 1];
    const C3 ratio = fin.ratio;
    // pins: max of the children's unknown flags
    if (L == 1)
      smg::restrictMax(c.unk, in.unknown, c.ext, fin.ext, c.g, fin.g, c.inner, ratio);
    else
      smg::restrictMax(c.unk, CCConst(fin.unk), c.ext, fin.ext, c.g, fin.g, c.inner, ratio);
    c.nUnk = countUnknown(L);
    // face products <Lam a>: coarsenOpenAvg's cell body on the fine products
    {
      CCField px = c.px, py = c.py, pz = c.pz;
      const C3 ce = c.ext, fe = fin.ext;
      const int gc = c.g, gf = fin.g;
      if (L == 1) {
        const long sy0 = fe.x, sz0 = (long)fe.x * fe.y;
        const smg::GuardedProduct fx{in.sax, in.unknown, in.lam, 1};
        const smg::GuardedProduct fy{in.say, in.unknown, in.lam, sy0};
        const smg::GuardedProduct fz{in.saz, in.unknown, in.lam, sz0};
        ccFor3(
            "peclet::flow::smg_coarsen_faces", C3{0, 0, 0}, c.inner,
            KOKKOS_LAMBDA(int icx, int icy, int icz) {
              coarsenOpenAvgCell(px, py, pz, fx, fy, fz, ce, fe, gc, gf, ratio, icx, icy, icz);
            });
      } else {
        CCConst fx = fin.px, fy = fin.py, fz = fin.pz;
        ccFor3(
            "peclet::flow::smg_coarsen_faces", C3{0, 0, 0}, c.inner,
            KOKKOS_LAMBDA(int icx, int icy, int icz) {
              coarsenOpenAvgCell(px, py, pz, fx, fy, fz, ce, fe, gc, gf, ratio, icx, icy, icz);
            });
      }
      fill(L, c.px);
      fill(L, c.py);
      fill(L, c.pz);
      // the low domain face of a non-periodic axis: no coupling toward the pinned ghost
      const CCField P[3] = {c.px, c.py, c.pz};
      for (int a = 0; a < 3; ++a)
        if (!per_[a] && touches(c, 2 * a))
          smg::zeroPlanes(P[a], c.ext, a, 0, c.g + 1);
    }
    // mass
    if (L == 1) {
      // A2: on the advective path the mass carries the open-face outflow omega_open (the interior
      // outflow is Q's); otherwise WO-5's lumped outflow, or nothing (WO-4) — those two verbatim
      const bool hasOut = in.advective || in.hasOutflow;
      const CCConst out = in.advective ? in.omegaOpen : (in.hasOutflow ? in.outflow : in.kappa);
      const smg::MassField m0{in.kappa, in.unknown, out, in.idt, hasOut};
      CCField cm = c.mass;
      const C3 ce = c.ext, fe = fin.ext, ci = c.inner;
      const int gc = c.g, gf = fin.g;
      ccFor3(
          "peclet::flow::smg_restrict_mass", C3{0, 0, 0}, ci,
          KOKKOS_LAMBDA(int icx, int icy, int icz) {
            restrictAvgCell(cm, m0, ce, fe, gc, gf, ratio, icx, icy, icz);
          });
    } else {
      restrictAvg(c.mass, CCConst(fin.mass), c.ext, fin.ext, c.g, fin.g, c.inner, ratio);
    }
    // Dirichlet boundary planes, averaged from the finer level's
    for (int f = 0; f < 6; ++f)
      if (dirFace_[f] && touches(c, f))
        coarsenPlane(L, f);
    // the operator
    if (two_) {
      allocateSolid(c, L);
      Kokkos::deep_copy(CCExec(), c.W, 0.0);
    }
    assembleCoarse(in, L);
    if (two_)  // WO-7: the solid phase, W gathered by assembleCoarse
      buildCoarseSolid(in, L);
    if (adv_) {  // A2: the coarse advection, then the 7 bands
      coarsenAdvection(in, L);
      assembleBands(L);
    }
  }
 
  void coarsenAdvection(const Inputs& in, int L) {
    Level& c = lv_[L];
    Level& fin = lv_[L - 1];
    if (c.qp[0].extent(0) != c.n)
      for (int a = 0; a < 3; ++a) {
        c.qp[a] = CCField("smg_qp", c.n);
        c.qm[a] = CCField("smg_qm", c.n);
      }
    const C3 ratio = fin.ratio, ce = c.ext, fe = fin.ext;
    const int gc = c.g, gf = fin.g;
    const double ir[3] = {1.0 / (double)ratio.x, 1.0 / (double)ratio.y, 1.0 / (double)ratio.z};
    for (int sgn = 0; sgn < 2; ++sgn) {
      CCField qx = sgn == 0 ? c.qp[0] : c.qm[0], qy = sgn == 0 ? c.qp[1] : c.qm[1],
              qz = sgn == 0 ? c.qp[2] : c.qm[2];
      const double irx = ir[0], iry = ir[1], irz = ir[2];
      if (L == 1) {
        const long sy0 = fe.x, sz0 = (long)fe.x * fe.y;
        const double sg = sgn == 0 ? 1.0 : -1.0;
        const smg::GuardedPart fx{in.phi[0], in.unknown, 1, sg};
        const smg::GuardedPart fy{in.phi[1], in.unknown, sy0, sg};
        const smg::GuardedPart fz{in.phi[2], in.unknown, sz0, sg};
        ccFor3(
            "peclet::flow::smg_coarsen_adv", C3{0, 0, 0}, c.inner,
            KOKKOS_LAMBDA(int icx, int icy, int icz) {
              coarsenOpenAvgCell(qx, qy, qz, fx, fy, fz, ce, fe, gc, gf, ratio, icx, icy, icz);
              const long ci =
                  (long)(icx + gc) + (long)(icy + gc) * ce.x + (long)(icz + gc) * (long)ce.x * ce.y;
              qx(ci) *= irx;
              qy(ci) *= iry;
              qz(ci) *= irz;
            });
      } else {
        const CCConst fx = sgn == 0 ? fin.qp[0] : fin.qm[0], fy = sgn == 0 ? fin.qp[1] : fin.qm[1],
                      fz = sgn == 0 ? fin.qp[2] : fin.qm[2];
        ccFor3(
            "peclet::flow::smg_coarsen_adv", C3{0, 0, 0}, c.inner,
            KOKKOS_LAMBDA(int icx, int icy, int icz) {
              coarsenOpenAvgCell(qx, qy, qz, fx, fy, fz, ce, fe, gc, gf, ratio, icx, icy, icz);
              const long ci =
                  (long)(icx + gc) + (long)(icy + gc) * ce.x + (long)(icz + gc) * (long)ce.x * ce.y;
              qx(ci) *= irx;
              qy(ci) *= iry;
              qz(ci) *= irz;
            });
      }
    }
    for (int a = 0; a < 3; ++a) {
      fill(L, c.qp[a]);
      fill(L, c.qm[a]);
    }
    for (int a = 0; a < 3; ++a)
      if (!per_[a] && touches(c, 2 * a)) {
        smg::zeroPlanes(c.qp[a], c.ext, a, 0, c.g + 1);
        smg::zeroPlanes(c.qm[a], c.ext, a, 0, c.g + 1);
      }
  }
 
  void assembleBands(int L) {
    Level& c = lv_[L];
    if (c.AW.extent(0) != c.n) {
      c.AW = CCField("smg_aw", c.n);
      c.AE = CCField("smg_ae", c.n);
      c.AS = CCField("smg_as", c.n);
      c.AN = CCField("smg_an", c.n);
      c.AB = CCField("smg_ab", c.n);
      c.AT = CCField("smg_at", c.n);
    }
    const C3 e = c.ext;
    const int g = c.g;
    CCField AC = c.AC, AW = c.AW, AE = c.AE, AS = c.AS, AN = c.AN, AB = c.AB, AT = c.AT;
    CCConst AFX = c.AFX, AFY = c.AFY, AFZ = c.AFZ, unk = c.unk;
    CCConst qpx = c.qp[0], qpy = c.qp[1], qpz = c.qp[2], qmx = c.qm[0], qmy = c.qm[1],
            qmz = c.qm[2];
    Kokkos::deep_copy(CCExec(), AW, 0.0);
    Kokkos::deep_copy(CCExec(), AE, 0.0);
    Kokkos::deep_copy(CCExec(), AS, 0.0);
    Kokkos::deep_copy(CCExec(), AN, 0.0);
    Kokkos::deep_copy(CCExec(), AB, 0.0);
    Kokkos::deep_copy(CCExec(), AT, 0.0);
    ccFor3(
        "peclet::flow::smg_coarse_bands", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
          const long i = (long)lx + (long)ly * sy + (long)lz * sz;
          if (!(unk(i) > 0.5))
            return;
          AC(i) =
              AC(i) + ((qmx(i) + qpx(i + sx)) + (qmy(i) + qpy(i + sy))) + (qmz(i) + qpz(i + sz));
          AW(i) = AFX(i) - qpx(i);
          AE(i) = AFX(i + sx) - qmx(i + sx);
          AS(i) = AFY(i) - qpy(i);
          AN(i) = AFY(i + sy) - qmy(i + sy);
          AB(i) = AFZ(i) - qpz(i);
          AT(i) = AFZ(i + sz) - qmz(i + sz);
        });
  }
 
  void coarsenPlane(int L, int f) {
    Level& c = lv_[L];
    Level& fin = lv_[L - 1];
    const int a = f / 2;
    int t1, t2;
    tangents(a, t1, t2);
    const int cin[3] = {c.inner.x, c.inner.y, c.inner.z};
    const int fin3[3] = {fin.inner.x, fin.inner.y, fin.inner.z};
    const int rr[3] = {fin.ratio.x, fin.ratio.y, fin.ratio.z};
    const int n1 = cin[t1], n2 = cin[t2], f1 = fin3[t1];
    const int r1 = rr[t1], r2 = rr[t2];
    if (c.plane[f].extent(0) != (std::size_t)n1 * n2)
      c.plane[f] = Kokkos::View<double*, CCMem>("smg_plane", (long)n1 * n2);
    auto pc = c.plane[f];
    auto pf = fin.plane[f];
    CCExec space;
    Kokkos::parallel_for(
        "peclet::flow::smg_coarsen_plane", MDRange2<CCExec>(space, {0, 0}, {n1, n2}),
        KOKKOS_LAMBDA(int j1, int j2) {
          double s = 0.0;
          for (int p = 0; p < r1; ++p)
            for (int q = 0; q < r2; ++q)
              s += pf((long)(r1 * j1 + p) + (long)(r2 * j2 + q) * f1);
          pc((long)j1 + (long)j2 * n1) = s / (double)(r1 * r2);
        });
  }
 
  void assembleCoarse(const Inputs& in, int L) {
    Level& c = lv_[L];
    const C3 e = c.ext, ci = c.inner, cf = c.cfac;
    const int g = c.g;
    const double wx = w_[0] / ((double)cf.x * cf.x), wy = w_[1] / ((double)cf.y * cf.y),
                 wz = w_[2] / ((double)cf.z * cf.z);
    CCField AC = c.AC, AFX = c.AFX, AFY = c.AFY, AFZ = c.AFZ;
    CCConst px = c.px, py = c.py, pz = c.pz;
    ccFor3(
        "peclet::flow::smg_coarse_faces", C3{0, 0, 0}, e, KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
          AFX(i) = -(wx * px(i));
          AFY(i) = -(wy * py(i));
          AFZ(i) = -(wz * pz(i));
        });
    // the wall term: the plain average of the level-0 wall terms (fine probe distance, A1)
    const double invN = 1.0 / ((double)cf.x * cf.y * cf.z);
    auto cstart = cutStart_[(std::size_t)L];
    auto crows = cutRows_[(std::size_t)L];
    auto cutCell = in.fac->cutCell;
    CCConst unk0 = in.unknown;
    CCConst unk = c.unk, mass = c.mass;
    auto cfs = in.fac->cellFacetStart;
    auto cw = in.facetW;
    // WO-7 (A1 for the interface): the coupling W_c is the plain average of the level-0 cc at the
    // fine probe distances, gathered in the same order; it enters the fluid diagonal and is stored
    // as the level's W (the solid diagonal and the off-diagonal -W follow in buildCoarseSolid).
    const bool two = two_;
    auto ccv = in.facetC;
    CCField Wl = c.W;
    if (c.wallG.extent(0) != c.n)  // the gathered W (and W_c) of this level, kept across builds
      c.wallG = CCField("smg_wall_gathered", c.n);
    if (two && c.wallCG.extent(0) != c.n)
      c.wallCG = CCField("smg_wallc_gathered", c.n);
    CCField wg = c.wallG, wcg = c.wallCG;
    const bool reuse = wallReuse_;
    ccFor3(
        "peclet::flow::smg_coarse_diag", C3{0, 0, 0}, ci, KOKKOS_LAMBDA(int icx, int icy, int icz) {
          const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
          const long i = (long)(icx + g) + (long)(icy + g) * sy + (long)(icz + g) * sz;
          if (!(unk(i) > 0.5)) {
            AC(i) = 1.0;
            AFX(i) = AFY(i) = AFZ(i) = 0.0;
            return;
          }
          double W = 0.0, Wc = 0.0;
          if (reuse) {  // wallKeep: the inputs of the last gather are unchanged
            W = wg(i);
            if (two)
              Wc = wcg(i);
          } else {
            // the cut descendants in increasing level-0 index = the (dz, dy, dx) order
            const long ic = (long)icx + (long)icy * ci.x + (long)icz * ci.x * ci.y;
            for (int k = cstart(ic); k < cstart(ic + 1); ++k) {
              const int row = crows(k);
              if (!(unk0(cutCell(row)) > 0.5))
                continue;
              for (int f = cfs(row); f < cfs(row + 1); ++f) {
                W += cw(f);
                if (two)
                  Wc += ccv(f);
              }
            }
            W *= invN;
            if (two)
              Wc *= invN;
            wg(i) = W;
            if (two)
              wcg(i) = Wc;
          }
          const double tw = wx * px(i), te = wx * px(i + sx);
          const double ts = wy * py(i), tn = wy * py(i + sy);
          const double tb = wz * pz(i), tt = wz * pz(i + sz);
          if (two) {
            AC(i) = ((mass(i) + (((tw + te) + (ts + tn)) + (tb + tt))) + W) + Wc;
            Wl(i) = Wc;
          } else {
            AC(i) = (mass(i) + (((tw + te) + (ts + tn)) + (tb + tt))) + W;
          }
        });
    // Dirichlet domain faces: AC += 2 w_a(L) <Lam a_bf> on the unknown boundary cells
    const double wa[3] = {wx, wy, wz};
    for (int f = 0; f < 6; ++f) {
      if (!(dirFace_[f] && touches(c, f)))
        continue;
      const int a = f / 2, side = f % 2;
      int t1, t2;
      tangents(a, t1, t2);
      const int ext[3] = {e.x, e.y, e.z};
      const long st[3] = {1, (long)e.x, (long)e.x * e.y};
      const int n1 = ext[t1] - 2 * g, n2 = ext[t2] - 2 * g;
      const long sa = st[a], s1 = st[t1], s2 = st[t2];
      const int aInner = side == 0 ? g : ext[a] - g - 1;
      const double tw2 = 2.0 * wa[a];
      auto pl = c.plane[f];
      CCExec space;
      Kokkos::parallel_for(
          "peclet::flow::smg_coarse_dirichlet", MDRange2<CCExec>(space, {0, 0}, {n1, n2}),
          KOKKOS_LAMBDA(int j1, int j2) {
            const long i = (long)aInner * sa + (long)(j1 + g) * s1 + (long)(j2 + g) * s2;
            if (!(unk(i) > 0.5))
              return;
            AC(i) += tw2 * pl((long)j1 + (long)j2 * n1);
          });
    }
  }
 
 private:
  std::vector<Level> lv_;
  Fill fill0_;
  bool fill0Wraps_ = false;  
  CCConst unk0_;
  CCConst unkS0_;     
  bool two_ = false;  
  std::vector<Kokkos::View<int*, CCMem>> cutStart_, cutRows_;
  const void* cutKey_ = nullptr;
  long cutKeyN_ = -1;
  Kokkos::View<double*, CCMem> cwSnap_, ccSnap_;
  bool wallValid_ = false, wallTwo_ = false, wallReuse_ = false;
  int wallUse_ = 0;
  const void* wallUnk_ = nullptr;
  double w_[3] = {1.0, 1.0, 1.0};
  double hp_[3] = {1.0, 1.0, 1.0};
  bool aniso_ = false;
  bool per_[3] = {true, true, true};
  bool dirFace_[6] = {false, false, false, false, false, false};
  bool singular_ = false;
  bool adv_ = false;  
  int nUse_ = 1;
  bool distributed_ = false;
#ifdef PECLET_FLOW_MPI
  MPI_Comm comm_ = MPI_COMM_NULL;
#endif
};
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_SCALAR_MG_HPP