mac_ibm.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_MAC_IBM_HPP
#define PECLET_FLOW_MAC_IBM_HPP
 
#include <Kokkos_Core.hpp>
#include <type_traits>
 
#include "cut_cell_ibm.hpp"    // IbmOverlayT, ibmFillEntry
#include "mac_cutcell.hpp"     // peclet::flow::C3, peclet::flow::ccSampleExt, CCConst
#include "mac_cutcell_mg.hpp"  // MReal (G.6: the overlay follows the operator storage precision)
#include "policy.hpp"
 
namespace peclet::flow {
 
using mreal = MReal;  // matrix coefficient type; float unless -DPECLET_FLOW_OPERATOR_DOUBLE (G.6)
// MConst / ibmRbgsSweep stay FIXED float: a standalone kernel-level facility exercised by
// tests/kokkos/test_ibm.cpp with its own literal-float coefficients, independent of the solver's
// MReal -- unlike IbmOverlay below, it is not on IbmSolver's overlay-precision path (ibmRbgsSweep
// itself has no caller in flow_ibm.hpp), so tying it to `mreal` would only make the test's fixed
// float buffers fail to bind under -DPECLET_FLOW_OPERATOR_DOUBLE for no G.6 benefit.
using MConst = Kokkos::View<const float*, CCMem>;
// The canonical cut-cell overlay type: IbmOverlayT<Space, Real> instantiated at this build's
// operator precision. Bit-identical to the pre-G.6 hardcoded-float IbmOverlay when mreal = float.
using IbmOverlay = IbmOverlayT<IMem, mreal>;
// The ghost-projection overlay (ghost_projection.hpp) at this build's operator precision --
// defined here, not in ghost_projection.hpp, because mac_cutcell_mg.hpp #includes
// ghost_projection.hpp BEFORE its own `using MReal = ...`, so ghost_projection.hpp cannot name
// MReal/mreal itself. Bit-identical to the pre-G.6 hardcoded-float GpOverlay when mreal = float.
using GpOverlayMReal = GpOverlayReal<mreal>;
 
struct Off3 {
  float x, y, z;
};
 
// Cut cell = fluid centre with at least one solid axis neighbour.
template <class Real>
KOKKOS_INLINE_FUNCTION bool ibmIsCut(Real sc, const Real sn[6]) {
  if (sc <= Real(0))
    return false;
  for (int k = 0; k < 6; ++k)
    if (sn[k] < Real(0))
      return true;
  return false;
}
 
// Find cut cells over the inner block and build the Robust-Scaled overlay (port of ibm_count_ext_k
// + ibm_geometry_ext_k): per inner cell, gather the 7 staggered SDF samples; if cut, atomically
// claim a slot, set idMap[cell]=slot, and fill the overlay. counter/idMap are reset here. Returns
// the cut count (overlay arrays must be sized >= number of inner cells). bc_type 0=Dirichlet,
// 1=Neumann.
template <int SCHEME>
inline int buildIbmOverlay(CCConst sdf, C3 ext, int g, Off3 off, int bc_type, const IbmOverlay& ov,
                           Kokkos::View<int*, CCMem> idMap, Kokkos::View<int, CCMem> counter,
                           CCConst tx = CCConst(), CCConst ty = CCConst(), CCConst tz = CCConst(),
                           C3 nn = C3{0, 0, 0}, mreal slipLambda = mreal(0), int comp = -1,
                           Kokkos::View<int, CCMem> sandwichSkipped = Kokkos::View<int, CCMem>(),
                           mreal hpx = mreal(1), mreal hpy = mreal(1), mreal hpz = mreal(1)) {
  CCExec space;
  Kokkos::deep_copy(space, counter, 0);
  Kokkos::deep_copy(space, idMap, -1);
  const bool hasEx = tx.size() > 0;
  const bool hasSlip = (slipLambda > mreal(0)) && (comp >= 0) && (SCHEME == 0);
  const mreal lamBase = slipLambda;
  const int slipComp = comp;
  int* skipPtr = sandwichSkipped.data();
  if (sandwichSkipped.data() != nullptr)
    Kokkos::deep_copy(space, sandwichSkipped, 0);
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::ibm_build_overlay", MD(space, {g, g, g}, {ext.x - g, ext.y - g, ext.z - g}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long idx = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
        // G.6: (mreal) instead of a hard (float) -- a double build no longer narrows the SDF
        // sample before the overlay is built from it.
        const mreal sc = (mreal)ccSampleExt(sdf, ext, lx + off.x, ly + off.y, lz + off.z);
        const int d[6][3] = {{1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 0, -1}};
        mreal sn[6];
        for (int k = 0; k < 6; ++k)
          sn[k] = (mreal)ccSampleExt(sdf, ext, lx + d[k][0] + off.x, ly + d[k][1] + off.y,
                                     lz + d[k][2] + off.z);
        if (!ibmIsCut(sc, sn))
          return;
        const int slot = Kokkos::atomic_fetch_add(&counter(), 1);
        idMap(idx) = slot;
        // Optional analytic-SDF exact crossings (setExactCrossings): tx/ty/tz are INNER-grid
        // arrays, t_a(i) = exact crossing fraction from this component's staggered point i toward
        // i + e_a (NaN = no crossing). +a dir: theta = t_a(i); -a dir: theta = 1 - t_a(i - e_a),
        // periodic wrap on the inner grid.
        mreal thEx[6];
        if (hasEx) {
          const int ix = lx - g, iy = ly - g, iz = lz - g;
          auto wrap = [](int v, int n) {
            v %= n;
            return v < 0 ? v + n : v;
          };
          const CCConst* ta[3] = {&tx, &ty, &tz};
          for (int a = 0; a < 3; ++a) {
            const long ip = (long)ix + (long)iy * nn.x + (long)iz * (long)nn.x * nn.y;
            const int mx = a == 0 ? wrap(ix - 1, nn.x) : ix;
            const int my = a == 1 ? wrap(iy - 1, nn.y) : iy;
            const int mz = a == 2 ? wrap(iz - 1, nn.z) : iz;
            const long im = (long)mx + (long)my * nn.x + (long)mz * (long)nn.x * nn.y;
            thEx[2 * a] = (mreal)(*ta[a])(ip);
            thEx[2 * a + 1] = mreal(1) - (mreal)(*ta[a])(im);  // NaN propagates -> fallback
          }
        }
        const mreal hpv[3] = {hpx, hpy, hpz};  // h_a'; exactly 1 isotropic (trap 9)
        mreal lamAxis[3];
        if (hasSlip) {
          // grad(sdf) from the same 7 samples; k = 0/2/4 are the +a neighbours, 1/3/5 the -a ones.
          // Phase 2 (§6.1/§6.4): g_a = h_a'·n_a, so the PHYSICAL unit normal is g/h' renormalised.
          const mreal gx = mreal(0.5) * (sn[0] - sn[1]) / hpv[0],
                      gy = mreal(0.5) * (sn[2] - sn[3]) / hpv[1],
                      gz = mreal(0.5) * (sn[4] - sn[5]) / hpv[2];
          mreal gm = Kokkos::sqrt(gx * gx + gy * gy + gz * gz);
          if (gm < mreal(1e-6))
            gm = mreal(1e-6);
          const mreal nrm[3] = {gx / gm, gy / gm, gz / gm};
          mreal tang = mreal(1) - nrm[slipComp] * nrm[slipComp];
          if (tang < mreal(0))
            tang = mreal(0);
          const mreal lamEff = lamBase * tang;
          for (int a = 0; a < 3; ++a) {
            mreal na = Kokkos::fabs(nrm[a]);
            if (na < mreal(1e-3))
              na = mreal(1e-3);  // wall nearly parallel to this axis: cap at 1000x (still finite)
            lamAxis[a] = lamEff / (hpv[a] * na);
          }
        }
        ibmFillEntry<SCHEME>(ov, slot, (int)idx, sc, sn, bc_type, hasEx ? thEx : nullptr,
                             hasSlip ? lamAxis : nullptr, skipPtr);
      });
 
  int cnt = 0;
  Kokkos::deep_copy(cnt, counter);
  return cnt;
}
 
// Volume fraction theta = clamp(0.5 + sdf_sample, 0, 1) at the staggered point (lx+off, ...).
//
// PHASE 2 (anisotropic cells, doc/anisotropic_metric.md §6.4), RATE ONLY — this is the velocity
// multigrid's staircase classification, so it decides a coarse operator's coefficient and never the
// fixed point.  `sd` is a distance in hRef units and the staircase wants the distance in CELLS
// along the normal, i.e. `theta = clamp(0.5 + d' |m|)` with `|m| = sqrt(sum_a n_a^2 w_a)` the index
// distance per unit physical distance along the normal (n = the physical unit normal of §6.1, from
// the sampled gradient at the staggered point).  The isotropic `|m|` is a `sqrt` of a unit vector,
// which is NOT exactly 1, so the `!aniso` path runs the literal legacy expression (the same
// host-bool dispatch C1 used for the const-coefficient smoother).
// A degenerate gradient (|g/h'| below the guard) falls back to `|m| = 1`, i.e. to the legacy
// expression: with no direction there is no per-axis information to use, and a plateau of the
// sampled SDF must not be turned into theta = 0.5 on a deep cell.
inline void ibmVolfrac(CCField theta, CCConst sdf, C3 ext, Off3 off, bool aniso = false,
                       double hpx = 1.0, double hpy = 1.0, double hpz = 1.0, double wx = 1.0,
                       double wy = 1.0, double wz = 1.0) {
  CCExec space;
  using MD = MDRange3<CCExec>;
  if (!aniso) {
    Kokkos::parallel_for(
        "peclet::flow::ibm_volfrac", MD(space, {0, 0, 0}, {ext.x, ext.y, ext.z}),
        KOKKOS_LAMBDA(int lx, int ly, int lz) {
          const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
          const double sd = ccSampleExt(sdf, ext, lx + off.x, ly + off.y, lz + off.z);
          const double t = 0.5 + sd;
          theta(i) = t < 0.0 ? 0.0 : (t > 1.0 ? 1.0 : t);
        });
    return;
  }
  Kokkos::parallel_for(
      "peclet::flow::ibm_volfrac_aniso", MD(space, {0, 0, 0}, {ext.x, ext.y, ext.z}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
        const double px = lx + off.x, py = ly + off.y, pz = lz + off.z;
        const double sd = ccSampleExt(sdf, ext, px, py, pz);
        // g_a = h_a' n_a from the central difference at the SAME staggered point.
        const double gx =
            0.5 *
            (ccSampleExt(sdf, ext, px + 1.0, py, pz) - ccSampleExt(sdf, ext, px - 1.0, py, pz)) /
            hpx;
        const double gy =
            0.5 *
            (ccSampleExt(sdf, ext, px, py + 1.0, pz) - ccSampleExt(sdf, ext, px, py - 1.0, pz)) /
            hpy;
        const double gz =
            0.5 *
            (ccSampleExt(sdf, ext, px, py, pz + 1.0) - ccSampleExt(sdf, ext, px, py, pz - 1.0)) /
            hpz;
        const double gm = Kokkos::sqrt(gx * gx + gy * gy + gz * gz);
        double mm = 1.0;
        if (gm > 1e-12) {
          const double nx = gx / gm, ny = gy / gm, nz = gz / gm;
          mm = Kokkos::sqrt(nx * nx * wx + ny * ny * wy + nz * nz * wz);
        }
        const double t = 0.5 + sd * mm;
        theta(i) = t < 0.0 ? 0.0 : (t > 1.0 ? 1.0 : t);
      });
}
 
// Solid mask: 1 where the staggered SDF point is inside the solid OR exactly ON the wall
// (sd<=0), else 0.
//
// WHY `<=` AND NOT `<` (WO-V6b part A). `sdf == 0` at a velocity DOF happens exactly when a flat
// wall sits on a cell FACE (an integer coordinate): the staggered sample is the mean of the two
// cell-centre values +-1/2. Every other consumer of the same field calls that point NON-fluid --
// `ibmIsCut` (`sc <= 0` => not a cut cell, so no Dirichlet closure is built), `ibmCleanFluidMask`
// (`sc <= 0` => solid), `ccFractionCore`/`ccFaceOpen` (`sd <= 0` => the face is CLOSED). With `<`
// this mask alone called it fluid, so the DOF was left unpinned AND unclosed: an unconstrained
// normal velocity sitting on the wall, invisible to the projection (its face is closed, so it
// never enters the divergence) but read by the advection operator of its fluid neighbours. A
// driven two-phase slit then diverges geometrically while every diagnostic stays healthy
// (WO-V7 finding 1). With `<=` the DOF is pinned to the wall datum, which for a DOF lying exactly
// ON the wall is the exact Dirichlet value, and the neighbouring fluid DOF's plain interior
// stencil becomes a wall-resolved discretization.
inline void ibmSolidMask(CCField mask, CCConst sdf, C3 ext, Off3 off) {
  CCExec space;
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::ibm_solid", MD(space, {0, 0, 0}, {ext.x, ext.y, ext.z}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
        const double sd = ccSampleExt(sdf, ext, lx + off.x, ly + off.y, lz + off.z);
        mask(i) = (sd <= 0.0) ? 1.0 : 0.0;
      });
}
 
// Clean-fluid-interior mask: 1 only at fluid cells with no solid neighbour (not cut, not solid).
inline void ibmCleanFluidMask(CCField m, CCConst sdf, C3 ext, Off3 off) {
  CCExec space;
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::ibm_clean", MD(space, {0, 0, 0}, {ext.x, ext.y, ext.z}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long i = (long)lx + (long)ly * ext.x + (long)lz * (long)ext.x * ext.y;
        const mreal sc = (mreal)ccSampleExt(sdf, ext, lx + off.x, ly + off.y, lz + off.z);
        const int d[6][3] = {{1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 0, -1}};
        mreal sn[6];
        for (int k = 0; k < 6; ++k)
          sn[k] = (mreal)ccSampleExt(sdf, ext, lx + d[k][0] + off.x, ly + d[k][1] + off.y,
                                     lz + d[k][2] + off.z);
        const bool solid = (sc <= mreal(0));
        m(i) = (solid || ibmIsCut(sc, sn)) ? 0.0 : 1.0;
      });
}
 
// One Red-Black sweep of the variable-coefficient stencil: x[i] = (b[i] - sum(A_off*x_nbr)) /
// A_C[i]. float matrix coeffs promote to double; solid cells pinned to 0. Call colour 0 then 1.
// Host backends take a line-sweep form: one task per (y,z) pencil, contiguous stride-2 x-loop
// starting at the colour's parity — no per-cell parity branch, unit-stride reads. Bit-identical to
// the MDRange form (same-colour cells are independent, so update order within a colour cannot
// change the result); the device keeps the MDRange form untouched.
template <class MC>
inline void ibmRbgsStencilColor(CCField x, CCConst b, MC AC, MC AW, MC AE, MC AS, MC AN, MC AB,
                                MC AT, CCConst solidmask, C3 ext, C3 og, int g, int color) {
  CCExec space;
  const bool hasMask = (solidmask.extent(0) != 0);
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    const int nyi = ext.y - 2 * g, nzi = ext.z - 2 * g;
    Kokkos::parallel_for(
        "peclet::flow::ibm_rbgs", Kokkos::RangePolicy<CCExec>(space, 0, (long)nyi * nzi),
        KOKKOS_LAMBDA(long t) {
          const int ly = g + (int)(t % nyi), lz = g + (int)(t / nyi);
          const long sx = 1, sy = ext.x, sz = (long)ext.x * ext.y;
          const int P = (color + og.x + og.y + ly + og.z + lz) & 1;
          for (int lx = g + ((P ^ (g & 1)) & 1); lx < ext.x - g; lx += 2) {
            const long i = (long)lx + (long)ly * sy + (long)lz * sz;
            if (hasMask && solidmask(i) > 0.5) {
              x(i) = 0.0;
              continue;
            }
            const double ac = AC(i);
            if (Kokkos::fabs(ac) < 1e-30)
              continue;
            const double s = (double)AE(i) * x(i + sx) + (double)AW(i) * x(i - sx) +
                             (double)AN(i) * x(i + sy) + (double)AS(i) * x(i - sy) +
                             (double)AT(i) * x(i + sz) + (double)AB(i) * x(i - sz);
            x(i) = (b(i) - s) / ac;
          }
        });
    return;
  }
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::ibm_rbgs", MD(space, {g, g, g}, {ext.x - g, ext.y - g, ext.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 = ext.x, sz = (long)ext.x * ext.y;
        const long i = (long)lx + (long)ly * ext.x + (long)lz * sz;
        if (hasMask && solidmask(i) > 0.5) {
          x(i) = 0.0;
          return;
        }
        const double ac = AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double s = (double)AE(i) * x(i + sx) + (double)AW(i) * x(i - sx) +
                         (double)AN(i) * x(i + sy) + (double)AS(i) * x(i - sy) +
                         (double)AT(i) * x(i + sz) + (double)AB(i) * x(i - sz);
        x(i) = (b(i) - s) / ac;
      });
}
 
// ibmRbgsStencilColor + a fused max|Δx| reduction over the swept colour (values already in
// registers — no extra memory pass). The momentum tolerance stop runs this as the SECOND colour of
// a sweep: one colour's max increment is the convergence proxy. Same math as the plain sweep.
template <class MC>
inline double ibmRbgsStencilColorDu(CCField x, CCConst b, MC AC, MC AW, MC AE, MC AS, MC AN, MC AB,
                                    MC AT, CCConst solidmask, C3 ext, C3 og, int g, int color) {
  CCExec space;
  const bool hasMask = (solidmask.extent(0) != 0);
  double du = 0.0;
  if constexpr (std::is_same_v<typename CCExec::memory_space, Kokkos::HostSpace>) {
    const int nyi = ext.y - 2 * g, nzi = ext.z - 2 * g;
    Kokkos::parallel_reduce(
        "peclet::flow::ibm_rbgs_du", Kokkos::RangePolicy<CCExec>(space, 0, (long)nyi * nzi),
        KOKKOS_LAMBDA(long t, double& m) {
          const int ly = g + (int)(t % nyi), lz = g + (int)(t / nyi);
          const long sx = 1, sy = ext.x, sz = (long)ext.x * ext.y;
          const int P = (color + og.x + og.y + ly + og.z + lz) & 1;
          for (int lx = g + ((P ^ (g & 1)) & 1); lx < ext.x - g; lx += 2) {
            const long i = (long)lx + (long)ly * sy + (long)lz * sz;
            if (hasMask && solidmask(i) > 0.5) {
              x(i) = 0.0;
              continue;
            }
            const double ac = AC(i);
            if (Kokkos::fabs(ac) < 1e-30)
              continue;
            const double s = (double)AE(i) * x(i + sx) + (double)AW(i) * x(i - sx) +
                             (double)AN(i) * x(i + sy) + (double)AS(i) * x(i - sy) +
                             (double)AT(i) * x(i + sz) + (double)AB(i) * x(i - sz);
            const double xn = (b(i) - s) / ac;
            const double d = Kokkos::fabs(xn - x(i));
            if (d > m)
              m = d;
            x(i) = xn;
          }
        },
        Kokkos::Max<double>(du));
    return du;
  }
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_reduce(
      "peclet::flow::ibm_rbgs_du", MD(space, {g, g, g}, {ext.x - g, ext.y - g, ext.z - g}),
      KOKKOS_LAMBDA(int lx, int ly, int lz, double& m) {
        if (((og.x + lx + og.y + ly + og.z + lz) & 1) != color)
          return;
        const long sx = 1, sy = ext.x, sz = (long)ext.x * ext.y;
        const long i = (long)lx + (long)ly * ext.x + (long)lz * sz;
        if (hasMask && solidmask(i) > 0.5) {
          x(i) = 0.0;
          return;
        }
        const double ac = AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double s = (double)AE(i) * x(i + sx) + (double)AW(i) * x(i - sx) +
                         (double)AN(i) * x(i + sy) + (double)AS(i) * x(i - sy) +
                         (double)AT(i) * x(i + sz) + (double)AB(i) * x(i - sz);
        const double xn = (b(i) - s) / ac;
        const double d = Kokkos::fabs(xn - x(i));
        if (d > m)
          m = d;
        x(i) = xn;
      },
      Kokkos::Max<double>(du));
  return du;
}
 
// Box-restricted variant for the distributed overlap smoother: sweeps [rlo,rhi) only, skipping
// [slo,shi) (pass slo==shi to skip nothing). Interior-then-shell split of one colour's sweep is
// bit-identical to the full sweep — a colour's cells never read same-colour cells (see
// cutcellSmoothColorBox in mac_pressure.hpp).
template <class MC>
inline void ibmRbgsStencilColorBox(CCField x, CCConst b, MC AC, MC AW, MC AE, MC AS, MC AN, MC AB,
                                   MC AT, CCConst solidmask, C3 ext, C3 og, int color, C3 rlo,
                                   C3 rhi, C3 slo, C3 shi) {
  if (rhi.x <= rlo.x || rhi.y <= rlo.y || rhi.z <= rlo.z)
    return;
  CCExec space;
  const bool hasMask = (solidmask.extent(0) != 0);
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::ibm_rbgs_box", MD(space, {rlo.x, rlo.y, rlo.z}, {rhi.x, rhi.y, rhi.z}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        if (lx >= slo.x && lx < shi.x && ly >= slo.y && ly < shi.y && lz >= slo.z && lz < shi.z)
          return;
        if (((og.x + lx + og.y + ly + og.z + lz) & 1) != color)
          return;
        const long sx = 1, sy = ext.x, sz = (long)ext.x * ext.y;
        const long i = (long)lx + (long)ly * ext.x + (long)lz * sz;
        if (hasMask && solidmask(i) > 0.5) {
          x(i) = 0.0;
          return;
        }
        const double ac = AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double s = (double)AE(i) * x(i + sx) + (double)AW(i) * x(i - sx) +
                         (double)AN(i) * x(i + sy) + (double)AS(i) * x(i - sy) +
                         (double)AT(i) * x(i + sz) + (double)AB(i) * x(i - sz);
        x(i) = (b(i) - s) / ac;
      });
}
 
// ibmRbgsStencilColorBox + the fused max|Δx| reduction (the box form of ibmRbgsStencilColorDu),
// for the distributed overlap of the momentum tolerance stop: interior-then-shell partial maxima
// combine by max, which is order-independent — max(interior, shell) is bit-identical to the
// one-pass reduction.
template <class MC>
inline double ibmRbgsStencilColorDuBox(CCField x, CCConst b, MC AC, MC AW, MC AE, MC AS, MC AN,
                                       MC AB, MC AT, CCConst solidmask, C3 ext, C3 og, int color,
                                       C3 rlo, C3 rhi, C3 slo, C3 shi) {
  if (rhi.x <= rlo.x || rhi.y <= rlo.y || rhi.z <= rlo.z)
    return 0.0;
  CCExec space;
  const bool hasMask = (solidmask.extent(0) != 0);
  double du = 0.0;
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_reduce(
      "peclet::flow::ibm_rbgs_du_box", MD(space, {rlo.x, rlo.y, rlo.z}, {rhi.x, rhi.y, rhi.z}),
      KOKKOS_LAMBDA(int lx, int ly, int lz, double& m) {
        if (lx >= slo.x && lx < shi.x && ly >= slo.y && ly < shi.y && lz >= slo.z && lz < shi.z)
          return;
        if (((og.x + lx + og.y + ly + og.z + lz) & 1) != color)
          return;
        const long sx = 1, sy = ext.x, sz = (long)ext.x * ext.y;
        const long i = (long)lx + (long)ly * ext.x + (long)lz * sz;
        if (hasMask && solidmask(i) > 0.5) {
          x(i) = 0.0;
          return;
        }
        const double ac = AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double s = (double)AE(i) * x(i + sx) + (double)AW(i) * x(i - sx) +
                         (double)AN(i) * x(i + sy) + (double)AS(i) * x(i - sy) +
                         (double)AT(i) * x(i + sz) + (double)AB(i) * x(i - sz);
        const double xn = (b(i) - s) / ac;
        const double d = Kokkos::fabs(xn - x(i));
        if (d > m)
          m = d;
        x(i) = xn;
      },
      Kokkos::Max<double>(du));
  return du;
}
 
inline void ibmRbgsSweep(CCField x, CCConst b, MConst AC, MConst AW, MConst AE, MConst AS,
                         MConst AN, MConst AB, MConst AT, CCConst solidmask, C3 ext, C3 og, int g) {
  ibmRbgsStencilColor(x, b, AC, AW, AE, AS, AN, AB, AT, solidmask, ext, og, g, 0);
  ibmRbgsStencilColor(x, b, AC, AW, AE, AS, AN, AB, AT, solidmask, ext, og, g, 1);
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_MAC_IBM_HPP