scalar_cutcell_operator.hpp

Go to the documentation of this file.

// Cut-cell scalar operator: the stored true operator of doc/scalar_ibm_design.md §4.1 (single
// phase, WO-3) — the 7 bands, the facet overlay coefficients, the right-hand side, the lumped
// level-0 surrogate diagonal of §4.3 — and the per-facet flux evaluation the budget, the wall-flux
// getter and the facet diagnostics share. Block kernels on flow's G = 2 extended block; the Solver
// side (unit conversions, refusals, the solve, the getters) is `flow_ibm_scalars_cutcell.hpp`.
//
// The fluid row of unknown i, per unit FULL-cell volume V (§1.4, backward Euler; idt = 0 steady):
//
//   kappa_i idt (c_i - c_i^n) + sum_a sum_+- Lam w_a a (c_i - c_nb) + sum_phi cw_phi u_p,phi + B_i
//     = kappa_i S_i + sum_phi rw_phi + (Dirichlet domain faces: 2 Lam w_a a_bf g_bf)
//
// with cw = alpha G and rw = alpha G g (Dirichlet G = Lam/s, Robin G = 1/(s/Lam + 1/k)) or cw = 0
// and rw = alpha q (Neumann). The bands hold kappa idt, the face diffusion Lam w a and the
// Dirichlet-domain fold 2 Lam w a_bf; the overlay adds cw times the probe interpolation. Every row
// that is not a fluid unknown is an identity row (AC = 1, bands 0, rhs 0), so c stays 0 there.
//
// A face toward a cell that is not an unknown carries no band (the guard of ruling D-WO3-4). By the
// geometry's invariant (kappa = 0 => every aperture 0; a sealed cell has every aperture 0) it is a
// no-op on every interior face; at a non-periodic global face the ghost's unknown flag was cleared
// by the record, which is exactly the "band toward the ghost is 0" of §1.4's domain-face table.
#ifndef PECLET_FLOW_SCALAR_CUTCELL_OPERATOR_HPP
#define PECLET_FLOW_SCALAR_CUTCELL_OPERATOR_HPP
 
#include <cstdint>
#include <Kokkos_Core.hpp>
#include <limits>
#include <memory>
#include <string>
#include <utility>
#include <vector>
 
#include "mac_cutcell.hpp"  // CCField, CCConst, CCExec, CCMem, C3, ccFor3
#include "peclet/core/scheme/probe_flux.hpp"
#include "policy.hpp"
#include "scalar_cutcell_geometry.hpp"
#include "staggered_advection.hpp"  // sadv::tvd / sou / fou_flux, the legacy reconstructions
 
namespace peclet::flow {
 
class ScalarMG;  // scalar_mg.hpp (WO-4): the preconditioner, one per cut-cell scalar
 
struct ScalarWallSpec {
  int type = 0;              
  double value = 0.0;        
  double coefficient = 0.0;  
  // conjugate (§1.1, §8.1 `set_scalar_solid`): the solid's own diffusivity D_s (L^2/T), the
  // capacity ratio C_s, the partition coefficient K (= c_s / c at equilibrium) and the contact
  // resistance R_c (T/L, on psi). psi form: Lam_s = C_s D_s K, capacity C_s K.
  double solidD = 0.0, solidC = 1.0, solidK = 1.0, solidRc = 0.0;
};
 
struct ScalarCutState {
  // ---- settings (physical) ----
  ScalarWallSpec wallDefault;  
  std::vector<std::pair<int, ScalarWallSpec>> wallInstance;  
  double rtol = 1e-10;         
  int maxit = 200;             
  double sourceConst = 0.0;    
  CCField sourceField;         
  bool sourceIsField = false;  
  bool hasProfile[6] = {false, false, false, false, false, false};
  std::vector<double> profile[6];  
  int profileN[6][2] = {};         
  double meanGrad[3] = {0.0, 0.0, 0.0};
  // ---- the operator of the last advance / steady solve ----
  CCField SAC;                            
  Kokkos::View<double*, CCMem> cw, rw;    
  Kokkos::View<double*, CCMem> gFace[6];  
  bool dirFace[6] = {false, false, false, false, false, false};
  double lam = 0.0;                   
  double idt = 0.0;                   
  double dt = 0.0;                    
  double srcFactor = 0.0;             
  std::vector<int> wtType;            
  std::vector<double> wtK, wtG, wtQ;  
  long numUnknowns = 0;               
  bool steady = false, built = false, singular = false;
  bool meanGradOn = false;            
  double gPhys[3] = {0.0, 0.0, 0.0};  
  double gInt[3] = {0.0, 0.0, 0.0};   
  double ubar[3] = {0.0, 0.0, 0.0};   
  // ---- Krylov scratch and statistics ----
  CCField kr, krh, kp, kv, kt, kz, kz2;
  Kokkos::View<double*, CCMem> ks;              
  Kokkos::View<double*, Kokkos::HostSpace> pk;  
  CCField kb, kq;  
  int iterations = 0;
  double residual = 0.0;  
  bool converged = true;
  bool warnedNoConv = false;
  long warnedGeometryVersion = -1;
  bool warnedKoren = false;  
  // persistent device tables, refilled every advance (review finding 7: no per-advance allocation):
  // the wall table (sco::WallTable) and the conjugate material table (sco::MaterialTable)
  Kokkos::View<int*, CCMem> wtTypeV, mtConjV;
  Kokkos::View<double*, CCMem> wtKV, wtGV, wtQV, mtLamV, mtCKV, mtKV, mtRcV;
  double incompatibility = 0.0;  
  // ---- ScalarMG (§5.2): level table per geometry version / block, coefficients per build ----
  std::shared_ptr<ScalarMG> mg;
  long mgVersion = -1;  
  std::size_t mgN = 0;  
  int mgLevels = 1;     
  // ---- advection (WO-5, §6): the face flux of the last advance and its classification ----
  CCField
      phi[3];  
  CCField Phi[3];   
  CCField small;    
  CCField outflow;  
  CCField omegaOpen;
  double maxCellPeclet = 0.0;  
  bool openFace[6] = {false, false, false, false, false, false};  
  bool openInflow[6] = {false, false, false, false, false, false};  
  bool advecting = false;       
  long numSmall = 0;            
  long numImplicitFaces = 0;    
  long numFluxFaces = 0;        
  long numGuardedFaces = 0;     
  long numBackflowFaces = 0;    
  bool warnedBackflow = false;  
  bool warnedBottom = false;    
  double bulkCourant = 0.0;     
  // ---- conjugate (WO-7, §1.1, §1.3, §2.7, §3.4, §4): the solid field psi_s = c_s / K ----------
  // Allocated by the first set_scalar_solid; the two-field path runs only while some body is
  // conjugate (`conj` of the last build), so a scalar without a conjugate body takes the
  // single-phase path above, bitwise.
  CCField solid;      
  bool conj = false;  
  // the per-scalar solid record (§2.7, §3.4), rebuilt when the geometry or the conjugate set moves
  long solidVersion = -1;           
  std::vector<char> solidKey;       
  CCField mat;                      
  CCField sunk;                     
  Kokkos::View<double*, CCMem> sS;  
  Kokkos::View<int* [8], Kokkos::LayoutLeft, CCMem> offS;
  Kokkos::View<double* [8], Kokkos::LayoutLeft, CCMem> wS;
  Kokkos::View<std::uint8_t*, CCMem> rungS;
  Kokkos::View<std::uint8_t*, CCMem> conjF;  
  long numSolidUnknowns = 0;                 
  double solidVolume = 0.0;                  
  double sealedSolidVolume = 0.0;            
  long solidRungs[4] = {0, 0, 0, 0};         
  // the solid operator of the last build (§4.1: the solid phase's own 7 bands)
  CCField spx, spy, spz;  
  CCField ms;             
  CCField ASC, ASW, ASE, ASS, ASN, ASB, AST, bS, sOld;
  CCField SACs, W0;                 
  Kokkos::View<double*, CCMem> cc;  
  std::vector<double> mtLam, mtCK, mtK, mtRc;  
  double maxSolidD = 0.0;  
  CCField krS, krhS, kpS, kvS, ktS, kzS, kz2S, kbS, kqS;  
};
 
namespace sco {
 
KOKKOS_INLINE_FUNCTION double lowAp(const CCConst& sax, const CCConst& say, const CCConst& saz,
                                    int a, long i) {
  return a == 0 ? sax(i) : (a == 1 ? say(i) : saz(i));
}
 
inline void zeroNonUnknown(CCField c, CCConst unk, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_zero_nonunknown", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (!(unk(i) > 0.5))
          c(i) = 0.0;
      });
}
 
inline void buildBands(CCField AC, CCField AW, CCField AE, CCField AS, CCField AN, CCField AB,
                       CCField AT, CCConst kappa, CCConst unk, CCConst sax, CCConst say,
                       CCConst saz, double lam, double idt, const double w[3], C3 e, int g) {
  const double lx = lam * w[0], ly = lam * w[1], lz = lam * w[2];
  ccFor3(
      "peclet::flow::sco_bands", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5)) {
          AC(i) = 1.0;
          AW(i) = AE(i) = AS(i) = AN(i) = AB(i) = AT(i) = 0.0;
          return;
        }
        const double tw = unk(i - sx) > 0.5 ? lx * sax(i) : 0.0;
        const double te = unk(i + sx) > 0.5 ? lx * sax(i + sx) : 0.0;
        const double ts = unk(i - sy) > 0.5 ? ly * say(i) : 0.0;
        const double tn = unk(i + sy) > 0.5 ? ly * say(i + sy) : 0.0;
        const double tb = unk(i - sz) > 0.5 ? lz * saz(i) : 0.0;
        const double tt = unk(i + sz) > 0.5 ? lz * saz(i + sz) : 0.0;
        AW(i) = -tw;
        AE(i) = -te;
        AS(i) = -ts;
        AN(i) = -tn;
        AB(i) = -tb;
        AT(i) = -tt;
        AC(i) = kappa(i) * idt + ((tw + te) + (ts + tn)) + (tb + tt);
      });
}
 
inline void buildRhs(CCField b, CCConst cOld, CCConst kappa, CCConst unk, double idt,
                     double srcConst, CCConst src, bool useField, double srcFactor, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_rhs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (!(unk(i) > 0.5)) {
          b(i) = 0.0;
          return;
        }
        const double S = (useField ? src(i) : srcConst) * srcFactor;
        b(i) = kappa(i) * (idt * cOld(i) + S);
      });
}
 
KOKKOS_INLINE_FUNCTION void faceTangents(int a, int& t1, int& t2) {
  t1 = (a == 0) ? 1 : 0;
  t2 = (a == 2) ? 1 : 2;
}
 
inline void dirichletFaceFold(CCField AC, CCField b, CCConst unk, CCConst saA,
                              Kokkos::View<const double*, CCMem> gv, double lamW, int a, int side,
                              C3 e, int g) {
  int t1, t2;
  faceTangents(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;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_dirichlet_face", 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;
        const double abf = side == 0 ? saA(i) : saA(i + sa);
        const double t = 2.0 * lamW * abf;
        AC(i) += t;
        b(i) += t * gv((long)j1 + (long)j2 * n1);
      });
  space.fence();
}
 
inline double dirichletFaceInflux(CCConst c, CCConst unk, CCConst saA,
                                  Kokkos::View<const double*, CCMem> gv, double lamW, double vol,
                                  int a, int side, C3 e, int g) {
  int t1, t2;
  faceTangents(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;
  Kokkos::View<double*, CCMem> q("peclet::flow::sco_face_influx", (long)n1 * n2);
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_dirichlet_face_influx", 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 long k = (long)j1 + (long)j2 * n1;
        if (!(unk(i) > 0.5)) {
          q(k) = 0.0;
          return;
        }
        const double abf = side == 0 ? saA(i) : saA(i + sa);
        q(k) = vol * (2.0 * lamW * abf) * (gv(k) - c(i));
      });
  space.fence();
  auto h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), q);
  double s = 0.0;
  for (long k = 0; k < (long)h.extent(0); ++k)
    s += h(k);
  return s;
}
 
inline void residualFluxForm(CCField r, CCConst c, CCConst b, CCConst kappa, CCConst unk,
                             CCConst sax, CCConst say, CCConst saz, double lam, double idt,
                             const double w[3], C3 e, int g) {
  const double lx = lam * w[0], ly = lam * w[1], lz = lam * w[2];
  ccFor3(
      "peclet::flow::sco_residual_flux", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5)) {
          r(i) = 0.0;
          return;
        }
        const double ci = c(i);
        const double tw = unk(i - sx) > 0.5 ? lx * sax(i) : 0.0;
        const double te = unk(i + sx) > 0.5 ? lx * sax(i + sx) : 0.0;
        const double ts = unk(i - sy) > 0.5 ? ly * say(i) : 0.0;
        const double tn = unk(i + sy) > 0.5 ? ly * say(i + sy) : 0.0;
        const double tb = unk(i - sz) > 0.5 ? lz * saz(i) : 0.0;
        const double tt = unk(i + sz) > 0.5 ? lz * saz(i + sz) : 0.0;
        const double flux = ((tw * (ci - c(i - sx)) + te * (ci - c(i + sx))) +
                             (ts * (ci - c(i - sy)) + tn * (ci - c(i + sy)))) +
                            (tb * (ci - c(i - sz)) + tt * (ci - c(i + sz)));
        r(i) = b(i) - kappa(i) * idt * ci - flux;
      });
}
 
inline void dirichletFaceResidual(CCField r, CCConst c, CCConst unk, CCConst saA, double lamW,
                                  int a, int side, C3 e, int g) {
  int t1, t2;
  faceTangents(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;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_dirichlet_face_residual", 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;
        const double abf = side == 0 ? saA(i) : saA(i + sa);
        r(i) -= 2.0 * lamW * abf * c(i);
      });
  space.fence();
}
 
struct WallTable {
  Kokkos::View<int*, CCMem> type;
  Kokkos::View<double*, CCMem> k, gval, q;
  int nb = 0;
};
 
inline void facetCoefficients(Kokkos::View<double*, CCMem> cw, Kokkos::View<double*, CCMem> rw,
                              const scg::ScalarFacetOverlay& fo, CCConst unk, const WallTable& wt,
                              double lam) {
  auto fcell = fo.cell;
  auto falpha = fo.alpha;
  auto fbody = fo.body;
  auto sF = fo.sF;
  auto type = wt.type;
  auto kk = wt.k;
  auto gg = wt.gval;
  auto qq = wt.q;
  const int nb = wt.nb;
  const double kInf = std::numeric_limits<double>::infinity();  // Dirichlet: G = Lam / s
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_facet_coeffs", Kokkos::RangePolicy<CCExec>(space, 0, fo.n),
      KOKKOS_LAMBDA(const long f) {
        if (!(unk(fcell(f)) > 0.5)) {
          cw(f) = 0.0;
          rw(f) = 0.0;
          return;
        }
        int bd = fbody(f);
        if (bd < 0 || bd >= nb)
          bd = 0;
        const int t = type(bd);
        const double al = falpha(f);
        if (t == 0) {
          cw(f) = 0.0;
          rw(f) = al * qq(bd);
          return;
        }
        if (t == 3) {  // conjugate (WO-7): the interface coupling is `cc`, not a wall term
          cw(f) = 0.0;
          rw(f) = 0.0;
          return;
        }
        const double G = peclet::core::scheme::wallConductance(sF(f), lam, t == 1 ? kInf : kk(bd));
        cw(f) = al * G;
        rw(f) = al * G * gg(bd);
      });
  space.fence();
}
 
inline void overlayApply(CCField y, CCConst x, const scg::ScalarFacetOverlay& fo,
                         Kokkos::View<const double*, CCMem> cw, bool subtract = false) {
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  auto offF = fo.offF;
  auto wF = fo.wF;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_overlay_apply", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        const long i = cutCell(c);
        double acc = 0.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f) {
          double up = 0.0;
          for (int k = 0; k < 8; ++k)
            up += wF(f, k) * x(i + (long)offF(f, k));
          acc += cw(f) * up;
        }
        if (subtract)
          y(i) -= acc;
        else
          y(i) += acc;
      });
}
 
inline void overlayRhs(CCField b, const scg::ScalarFacetOverlay& fo,
                       Kokkos::View<const double*, CCMem> rw) {
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_overlay_rhs", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        double acc = 0.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f)
          acc += rw(f);
        b(cutCell(c)) += acc;
      });
  space.fence();
}
 
inline void surrogateDiagonal(CCField SAC, CCConst AC, const scg::ScalarFacetOverlay& fo,
                              Kokkos::View<const double*, CCMem> cw) {
  Kokkos::deep_copy(CCExec(), SAC, AC);
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_surrogate_diag", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        double acc = 0.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f)
          acc += cw(f);
        SAC(cutCell(c)) += acc;
      });
  space.fence();
}
 
inline void facetFlux(Kokkos::View<double*, CCMem> up, Kokkos::View<double*, CCMem> qin, CCConst x,
                      const scg::ScalarFacetOverlay& fo, Kokkos::View<const double*, CCMem> cw,
                      Kokkos::View<const double*, CCMem> rw) {
  auto fcell = fo.cell;
  auto offF = fo.offF;
  auto wF = fo.wF;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_facet_flux", Kokkos::RangePolicy<CCExec>(space, 0, fo.n),
      KOKKOS_LAMBDA(const long f) {
        const long i = fcell(f);
        double u = 0.0;
        for (int k = 0; k < 8; ++k)
          u += wF(f, k) * x(i + (long)offF(f, k));
        up(f) = u;
        qin(f) = rw(f) - cw(f) * u;
      });
  space.fence();
}
 
// ---- inner-cell reductions (this rank; the Solver adds the MPI_Allreduce) ----------------------
 
// Each reduction is written once, with its result as a parameter: a host double (the Local
// form) or a device slot (the resident Krylov, scalar_krylov.hpp). Same policy, same functor, so
// the same value either way (WO-9b).
 
template <class R>
inline void dotReduce(CCConst a, CCConst b, C3 e, int g, R&& result) {
  ccReduce3(
      "peclet::flow::sco_dot", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        acc += a(i) * b(i);
      },
      std::forward<R>(result));
}
inline double dotLocal(CCConst a, CCConst b, C3 e, int g) {
  double s = 0.0;
  dotReduce(a, b, e, g, s);
  return s;
}
 
template <class R1, class R2>
inline void dot2Reduce(CCConst a, CCConst b, CCConst c, C3 e, int g, R1&& r1, R2&& r2) {
  CCExec space;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_dot2", MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        p1 += a(i) * b(i);
        p2 += c(i) * c(i);
      },
      std::forward<R1>(r1), std::forward<R2>(r2));
}
inline void dot2Local(CCConst a, CCConst b, CCConst c, C3 e, int g, double& ab, double& cc) {
  double s1 = 0.0, s2 = 0.0;
  dot2Reduce(a, b, c, e, g, s1, s2);
  ab = s1;
  cc = s2;
}
 
template <class Red>
inline void maxabsReduce(CCConst a, C3 e, int g, Red red) {
  ccReduce3(
      "peclet::flow::sco_maxabs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        const double v = Kokkos::fabs(a(i));
        if (v > acc)
          acc = v;
      },
      red);
}
inline double maxabsLocal(CCConst a, C3 e, int g) {
  double m = 0.0;
  maxabsReduce(a, e, g, Kokkos::Max<double>(m));
  return m;
}
 
inline void sumUnknownLocal(CCConst a, CCConst unk, C3 e, int g, double& sum, double& sumAbs) {
  CCExec space;
  double s1 = 0.0, s2 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_sum_unknown",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (unk(i) > 0.5) {
          p1 += a(i);
          p2 += Kokkos::fabs(a(i));
        }
      },
      s1, s2);
  sum = s1;
  sumAbs = s2;
}
 
inline void addOnUnknown(CCField a, CCConst unk, double d, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_add_unknown", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (unk(i) > 0.5)
          a(i) += d;
      });
}
 
inline void kappaMomentsLocal(CCConst c, CCConst kappa, CCConst unk, C3 e, int g, double& kc,
                              double& k) {
  CCExec space;
  double s1 = 0.0, s2 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_kappa_moments",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (unk(i) > 0.5) {
          p1 += kappa(i) * c(i);
          p2 += kappa(i);
        }
      },
      s1, s2);
  kc = s1;
  k = s2;
}
 
inline double sourceSumLocal(CCConst kappa, CCConst unk, double srcConst, CCConst src,
                             bool useField, double srcFactor, C3 e, int g) {
  double s = 0.0;
  ccReduce3(
      "peclet::flow::sco_source_sum", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (unk(i) > 0.5)
          acc += kappa(i) * ((useField ? src(i) : srcConst) * srcFactor);
      },
      s);
  return s;
}
 
// ---- advection and small cells (WO-5, §6) -------------------------------------------------------
//
// phi_a(i) is the volume flux through the LOW a-face of cell i per unit FULL-cell volume per
// internal time, positive along +a: open * vel with the projection's own openness and face velocity
// (Solver::scalarFaceFlux, §6.1). The velocity is an INDEX velocity (cells per tRef), so the A_a/V
// = 1/h'_a of F/V is already inside it (doc/anisotropic_metric.md §2). The flux out of cell i
// through its high a-face is phi_a(i + s_a), through its low a-face -phi_a(i).
 
inline void faceFluxes(CCField phi, CCConst vel, CCConst open, CCConst unk, int a, C3 e) {
  const long st[3] = {1, (long)e.x, (long)e.x * e.y};
  const long s = st[a];
  ccFor3(
      "peclet::flow::sco_face_flux", 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;
        const int xa = a == 0 ? x : (a == 1 ? y : z);
        if (xa < 1 || !(unk(i) > 0.5) || !(unk(i - s) > 0.5)) {
          phi(i) = 0.0;
          return;
        }
        phi(i) = open(i) * vel(i);
      });
}
 
inline void faceFluxCountsLocal(CCConst phi, CCConst vel, CCConst open, int a, C3 e, int g,
                                long& nFlux, long& nGuarded) {
  (void)a;
  CCExec space;
  long n1 = 0, n2 = 0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_flux_counts",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, long& p1, long& p2) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (phi(i) != 0.0)
          ++p1;
        else if (open(i) * vel(i) != 0.0)
          ++p2;
      },
      n1, n2);
  nFlux = n1;
  nGuarded = n2;
}
 
inline double maxAbsFluxLocal(CCConst phx, CCConst phy, CCConst phz, C3 e, int g) {
  double m = 0.0;
  ccReduce3(
      "peclet::flow::sco_max_flux", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        const double v = Kokkos::fmax(Kokkos::fabs(phx(i)),
                                      Kokkos::fmax(Kokkos::fabs(phy(i)), Kokkos::fabs(phz(i))));
        if (v > acc)
          acc = v;
      },
      Kokkos::Max<double>(m));
  return m;
}
 
inline double maxCellPecletLocal(CCConst phx, CCConst phy, CCConst phz, CCConst unk, double lam,
                                 const double w[3], C3 e, int g) {
  const double dx = lam * w[0], dy = lam * w[1], dz = lam * w[2];
  double m = 0.0;
  ccReduce3(
      "peclet::flow::sco_max_cell_peclet", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        const double P[3] = {unk(i - sx) > 0.5 ? Kokkos::fabs(phx(i)) / dx : 0.0,
                             unk(i - sy) > 0.5 ? Kokkos::fabs(phy(i)) / dy : 0.0,
                             unk(i - sz) > 0.5 ? Kokkos::fabs(phz(i)) / dz : 0.0};
        for (int k = 0; k < 3; ++k)
          if (P[k] > acc)  // (0 / 0 = NaN when lam = 0 and no flux: skipped)
            acc = P[k];
      },
      Kokkos::Max<double>(m));
  return m > 0.0 ? m : 0.0;
}
 
inline void clampNonNegative(CCField omega, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_clamp_omega", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (omega(i) < 0.0)
          omega(i) = 0.0;
      });
}
 
KOKKOS_INLINE_FUNCTION double cellOutflow(const CCConst& phx, const CCConst& phy,
                                          const CCConst& phz, long i, long sy, long sz) {
  return ((Kokkos::fmax(phx(i + 1), 0.0) + Kokkos::fmax(-phx(i), 0.0)) +
          (Kokkos::fmax(phy(i + sy), 0.0) + Kokkos::fmax(-phy(i), 0.0))) +
         (Kokkos::fmax(phz(i + sz), 0.0) + Kokkos::fmax(-phz(i), 0.0));
}
 
inline double bulkCourantLocal(CCConst phx, CCConst phy, CCConst phz, CCConst kappa, CCConst unk,
                               CCConst sax, CCConst say, CCConst saz, double dt, C3 e, int g) {
  double m = 0.0;
  ccReduce3(
      "peclet::flow::sco_bulk_courant", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5) || !(kappa(i) >= 1.0) || !(sax(i) >= 1.0) || !(sax(i + 1) >= 1.0) ||
            !(say(i) >= 1.0) || !(say(i + sy) >= 1.0) || !(saz(i) >= 1.0) || !(saz(i + sz) >= 1.0))
          return;
        const double v = dt * cellOutflow(phx, phy, phz, i, sy, sz);
        if (v > acc)
          acc = v;
      },
      Kokkos::Max<double>(m));
  return m;
}
 
inline void smallCells(CCField small, CCConst phx, CCConst phy, CCConst phz, CCConst kappa,
                       CCConst unk, double dt, double thr, C3 e, int g) {
  Kokkos::deep_copy(CCExec(), small, 0.0);
  ccFor3(
      "peclet::flow::sco_small_cells", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const double k = kappa(i);
        if (!(unk(i) > 0.5) || !(k < 1.0))
          return;
        small(i) = dt * cellOutflow(phx, phy, phz, i, sy, sz) > thr * k ? 1.0 : 0.0;
      });
}
 
inline long countSmallLocal(CCConst small, C3 e, int g) {
  long n = 0;
  ccReduce3(
      "peclet::flow::sco_count_small", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, long& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (small(i) > 0.5)
          ++acc;
      },
      n);
  return n;
}
 
inline long countImplicitLocal(CCConst phi, CCConst small, CCConst unk, int a, bool steady, C3 e,
                               int g) {
  const long st[3] = {1, (long)e.x, (long)e.x * e.y};
  const long s = st[a];
  long n = 0;
  ccReduce3(
      "peclet::flow::sco_count_implicit", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, long& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (phi(i) != 0.0 && unk(i) > 0.5 && unk(i - s) > 0.5 &&
            (steady || small(i) > 0.5 || small(i - s) > 0.5))
          ++acc;
      },
      n);
  return n;
}
 
inline void advectionBands(CCField AC, CCField AW, CCField AE, CCField AS, CCField AN, CCField AB,
                           CCField AT, CCField outflow, CCConst phx, CCConst phy, CCConst phz,
                           CCConst small, CCConst unk, bool steady, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_adv_bands", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5)) {
          outflow(i) = 0.0;
          return;
        }
        const bool si = small(i) > 0.5;
        // F_out per face: W, E, S, N, B, T; implicit flag per face
        const double F[6] = {-phx(i), phx(i + sx), -phy(i), phy(i + sy), -phz(i), phz(i + sz)};
        const long nb[6] = {i - sx, i + sx, i - sy, i + sy, i - sz, i + sz};
        double d[6], o[6];
        for (int f = 0; f < 6; ++f) {
          const bool imp = (steady || si || small(nb[f]) > 0.5) && unk(nb[f]) > 0.5;
          d[f] = imp ? Kokkos::fmax(F[f], 0.0) : 0.0;
          o[f] = imp ? Kokkos::fmin(F[f], 0.0) : 0.0;
        }
        const double om = ((d[0] + d[1]) + (d[2] + d[3])) + (d[4] + d[5]);
        outflow(i) = om;
        AC(i) += om;
        AW(i) += o[0];
        AE(i) += o[1];
        AS(i) += o[2];
        AN(i) += o[3];
        AB(i) += o[4];
        AT(i) += o[5];
      });
}
 
inline void explicitFaceFluxes(CCField Phi, CCConst phi, CCConst c, CCConst unk, CCConst small,
                               int a, int scheme, bool steady, C3 e, int g) {
  const long st[3] = {1, (long)e.x, (long)e.x * e.y};
  const long s = st[a];
  const C3 hi{a == 0 ? e.x - g + 1 : e.x - g, a == 1 ? e.y - g + 1 : e.y - g,
              a == 2 ? e.z - g + 1 : e.z - g};
  Kokkos::deep_copy(CCExec(), Phi, 0.0);
  ccFor3(
      "peclet::flow::sco_explicit_flux", C3{g, g, g}, hi, KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        const double v = phi(i);
        if (v == 0.0 || steady || small(i) > 0.5 || small(i - s) > 0.5 || !(unk(i) > 0.5) ||
            !(unk(i - s) > 0.5))
          return;
        const double cLL = c(i - 2 * s), cL = c(i - s), cR = c(i), cRR = c(i + s);
        bool fou = scheme == 0;
        if (!fou)
          fou = v > 0.0 ? !(unk(i - 2 * s) > 0.5 && unk(i - s) > 0.5 && unk(i) > 0.5)
                        : !(unk(i + s) > 0.5 && unk(i) > 0.5 && unk(i - s) > 0.5);
        Phi(i) =
            fou ? sadv::fou_flux(cL, cR, v)
                : (scheme == 2 ? sadv::sou(cLL, cL, cR, cRR, v) : sadv::tvd(cLL, cL, cR, cRR, v));
      });
}
 
inline void explicitAdvectionRhs(CCField b, CCConst Px, CCConst Py, CCConst Pz, CCConst unk, C3 e,
                                 int g) {
  ccFor3(
      "peclet::flow::sco_explicit_rhs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        b(i) -= ((Px(i + sx) - Px(i)) + (Py(i + sy) - Py(i))) + (Pz(i + sz) - Pz(i));
      });
}
 
inline void advectionResidual(CCField r, CCConst c, CCConst phx, CCConst phy, CCConst phz,
                              CCConst small, CCConst unk, bool steady, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_adv_residual", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        const bool si = small(i) > 0.5;
        const double F[6] = {-phx(i), phx(i + sx), -phy(i), phy(i + sy), -phz(i), phz(i + sz)};
        const long nb[6] = {i - sx, i + sx, i - sy, i + sy, i - sz, i + sz};
        const double ci = c(i);
        double q[6];
        for (int f = 0; f < 6; ++f) {
          const bool imp = (steady || si || small(nb[f]) > 0.5) && unk(nb[f]) > 0.5;
          q[f] = imp ? (F[f] > 0.0 ? F[f] * ci : F[f] * c(nb[f])) : 0.0;
        }
        r(i) -= ((q[0] + q[1]) + (q[2] + q[3])) + (q[4] + q[5]);
      });
}
 
// ---- open domain faces (WO-5b; §1.4 B_i, §6.5, ruling D-WO5-3) ---------------------------------
//
// A flow inflow / outflow global face carries the boundary-face flux the projection constrained,
// open * vel at the boundary face (the low face of the first inner cell, or the first GHOST index
// on the high side), captured right after the projection (Solver::scalarCaptureOpenFaceFlux) as a
// plane over the block's FULL transverse extent, index j1 + j2 * ext[t1] (faceTangents order).
// openFaceFluxes writes it into phi_a at the boundary face, so Out_i, C_bulk and the small flags
// of §6.3 see it; the interior kernels above skip every face toward a non-unknown, so its row is
// carried here alone, per §1.4's domain-face table:
//   inflow   b_i += F_in g_bf (g_bf the scalar's Dirichlet value; always on the rhs);
//   outflow  F_out c_i (upwind on the inner cell: the zero-gradient exit), explicit at c^n on the
//            rhs, or implicit (AC += F_out) when the inner cell is small, and always when steady.
// F_out is the volume flux per unit V out of the inner cell through that face: -phi on the low
// side, +phi on the high side.
 
inline void openFaceFluxes(CCField phi, Kokkos::View<const double*, CCMem> plane, CCConst unk,
                           int a, int side, C3 e, int g) {
  int t1, t2;
  faceTangents(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 long sa = st[a], s1 = st[t1], s2 = st[t2];
  const int m1 = ext[t1], m2 = ext[t2];
  const long bfa = side == 0 ? g : ext[a] - g;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_open_face_flux", MDRange2<CCExec>(space, {0, 0}, {m1, m2}),
      KOKKOS_LAMBDA(int j1, int j2) {
        const long bf = bfa * sa + (long)j1 * s1 + (long)j2 * s2;
        const long in = side == 0 ? bf : bf - sa;
        phi(bf) = unk(in) > 0.5 ? plane((long)j1 + (long)j2 * m1) : 0.0;
      });
  space.fence();
}
 
inline void openFaceAdvection(CCField AC, CCField b, CCField outflow, CCConst cOld, CCConst phi,
                              CCConst small, CCConst unk, Kokkos::View<const double*, CCMem> gv,
                              bool inflow, int a, int side, bool steady, C3 e, int g,
                              CCField omega = CCField(), bool wOmega = false) {
  int t1, t2;
  faceTangents(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;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_open_face_adv", 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;
        const double Fout = side == 0 ? -phi(i) : phi(i + sa);
        if (inflow) {
          b(i) -= Fout * gv((long)j1 + (long)j2 * n1);
        } else if (steady || small(i) > 0.5) {
          AC(i) += Fout;
          outflow(i) += Fout;
          if (wOmega)
            omega(i) += Fout;
        } else {
          b(i) -= Fout * cOld(i);
        }
      });
  space.fence();
}
 
inline void openFaceResidual(CCField r, CCConst c, CCConst phi, CCConst small, CCConst unk,
                             bool inflow, int a, int side, bool steady, C3 e, int g) {
  if (inflow)
    return;
  int t1, t2;
  faceTangents(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;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_open_face_residual", 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) || !(steady || small(i) > 0.5))
          return;
        const double Fout = side == 0 ? -phi(i) : phi(i + sa);
        r(i) -= Fout * c(i);
      });
  space.fence();
}
 
inline double openFaceInflux(CCConst c, CCConst cOld, CCConst phi, CCConst small, CCConst unk,
                             Kokkos::View<const double*, CCMem> gv, bool inflow, double vol, int a,
                             int side, bool steady, C3 e, int g) {
  int t1, t2;
  faceTangents(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;
  Kokkos::View<double*, CCMem> q("peclet::flow::sco_open_face_influx", (long)n1 * n2);
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_open_face_influx", 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 long k = (long)j1 + (long)j2 * n1;
        if (!(unk(i) > 0.5)) {
          q(k) = 0.0;
          return;
        }
        const double Fout = side == 0 ? -phi(i) : phi(i + sa);
        const double cv = inflow ? gv(k) : ((steady || small(i) > 0.5) ? c(i) : cOld(i));
        q(k) = -vol * (Fout * cv);
      });
  space.fence();
  auto h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), q);
  double s = 0.0;
  for (long k = 0; k < (long)h.extent(0); ++k)
    s += h(k);
  return s;
}
 
inline void openFaceCountsLocal(CCConst phi, CCConst small, CCConst unk, bool inflow, int a,
                                int side, bool steady, C3 e, int g, long& nFlux, long& nImpl,
                                long& nBack) {
  int t1, t2;
  faceTangents(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;
  CCExec space;
  long p = 0, q = 0, r = 0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_open_face_counts", MDRange2<CCExec>(space, {0, 0}, {n1, n2}),
      KOKKOS_LAMBDA(int j1, int j2, long& c1, long& c2, long& c3) {
        const long i = (long)aInner * sa + (long)(j1 + g) * s1 + (long)(j2 + g) * s2;
        const double v = side == 0 ? phi(i) : phi(i + sa);
        if (!(unk(i) > 0.5) || v == 0.0)
          return;
        if (side == 1)
          ++c1;
        if (!inflow && (steady || small(i) > 0.5))
          ++c2;
        if (!inflow && (side == 0 ? v > 0.0 : v < 0.0))  // F_out < 0
          ++c3;
      },
      p, q, r);
  nFlux = p;
  nImpl = q;
  nBack = r;
}
 
inline double maxFacetLocal(Kokkos::View<const double*, CCMem> cw) {
  double m = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_max_cw", Kokkos::RangePolicy<CCExec>(0, (long)cw.extent(0)),
      KOKKOS_LAMBDA(const long f, double& acc) {
        if (cw(f) > acc)
          acc = cw(f);
      },
      Kokkos::Max<double>(m));
  return m;
}
 
// ---- mean-gradient (closure) mode (WO-6, §1.5) --------------------------------------------------
//
// The unknown is theta, c = G.x + theta. Wherever the discretization evaluates c at a point x it
// uses theta + G.x and moves the G.x part to the right-hand side. Face differences use the TRUE
// displacement +-h'_a e_a (never a wrapped index difference); a probe stencil cell or a domain-face
// centre uses its global position, unwrapped relative to the cut cell (an extended-block index
// past a periodic face is simply its unwrapped position).
 
struct MeanGradPosition {
  double G[3] = {0.0, 0.0, 0.0};    
  double org[3] = {0.0, 0.0, 0.0};  
  double h[3] = {1.0, 1.0, 1.0};    
  long base[3] = {0, 0, 0};         
  long sy = 1, sz = 1;              
  KOKKOS_INLINE_FUNCTION double operator()(long j) const {
    const long z = j / sz, y = (j - z * sz) / sy, x = j - z * sz - y * sy;
    const double px = org[0] + ((double)(base[0] + x) + 0.5) * h[0];
    const double py = org[1] + ((double)(base[1] + y) + 0.5) * h[1];
    const double pz = org[2] + ((double)(base[2] + z) + 0.5) * h[2];
    return (G[0] * px + G[1] * py) + G[2] * pz;
  }
};
 
inline void meanGradientFaceRhs(CCField b, CCConst unk, CCConst sax, CCConst say, CCConst saz,
                                const double gf[3], C3 e, int g) {
  const double gx = gf[0], gy = gf[1], gz = gf[2];
  ccFor3(
      "peclet::flow::sco_mg_face_rhs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        const double aw = unk(i - sx) > 0.5 ? sax(i) : 0.0;
        const double ae = unk(i + sx) > 0.5 ? sax(i + sx) : 0.0;
        const double as = unk(i - sy) > 0.5 ? say(i) : 0.0;
        const double an = unk(i + sy) > 0.5 ? say(i + sy) : 0.0;
        const double ab = unk(i - sz) > 0.5 ? saz(i) : 0.0;
        const double at = unk(i + sz) > 0.5 ? saz(i + sz) : 0.0;
        b(i) += (gx * (ae - aw) + gy * (an - as)) + gz * (at - ab);
      });
}
 
inline void meanGradientProbeShift(Kokkos::View<double*, CCMem> rw,
                                   const scg::ScalarFacetOverlay& fo,
                                   Kokkos::View<const double*, CCMem> cw,
                                   const MeanGradPosition& pos) {
  auto fcell = fo.cell;
  auto offF = fo.offF;
  auto wF = fo.wF;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_mg_probe_shift", Kokkos::RangePolicy<CCExec>(space, 0, fo.n),
      KOKKOS_LAMBDA(const long f) {
        if (cw(f) == 0.0)
          return;
        const long i = fcell(f);
        double gxp = 0.0;
        for (int k = 0; k < 8; ++k)
          gxp += wF(f, k) * pos(i + (long)offF(f, k));
        rw(f) -= cw(f) * gxp;
      });
  space.fence();
}
 
KOKKOS_INLINE_FUNCTION void cellVelocity(const CCConst& phx, const CCConst& phy, const CCConst& phz,
                                         long i, long sy, long sz, const double hh[3],
                                         double u[3]) {
  u[0] = hh[0] * (phx(i + 1) + phx(i));
  u[1] = hh[1] * (phy(i + sy) + phy(i));
  u[2] = hh[2] * (phz(i + sz) + phz(i));
}
 
inline void cellVelocitySumLocal(CCConst phx, CCConst phy, CCConst phz, CCConst unk,
                                 const double hh[3], C3 e, int g, double s[3]) {
  CCExec space;
  const double h0 = hh[0], h1 = hh[1], h2 = hh[2];
  double a = 0.0, b = 0.0, c = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_mg_cell_velocity_sum",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2, double& p3) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        const double hq[3] = {h0, h1, h2};
        double u[3];
        cellVelocity(phx, phy, phz, i, sy, sz, hq, u);
        p1 += u[0];
        p2 += u[1];
        p3 += u[2];
      },
      a, b, c);
  s[0] = a;
  s[1] = b;
  s[2] = c;
}
 
inline void meanGradientAdvectionRhs(CCField b, CCConst phx, CCConst phy, CCConst phz,
                                     CCConst kappa, CCConst unk, const double gi[3],
                                     const double hh[3], const double ub[3], C3 e, int g) {
  const double g0 = gi[0], g1 = gi[1], g2 = gi[2], h0 = hh[0], h1 = hh[1], h2 = hh[2], u0 = ub[0],
               u1 = ub[1], u2 = ub[2];
  ccFor3(
      "peclet::flow::sco_mg_adv_rhs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        const double hq[3] = {h0, h1, h2};
        double u[3];
        cellVelocity(phx, phy, phz, i, sy, sz, hq, u);
        const double k = kappa(i);
        b(i) -= (g0 * (u[0] - k * u0) + g1 * (u[1] - k * u1)) + g2 * (u[2] - k * u2);
      });
}
 
inline void meanFluxLocal(CCConst th, CCConst kappa, CCConst unk, CCConst sax, CCConst say,
                          CCConst saz, CCConst phx, CCConst phy, CCConst phz, bool advecting,
                          double lam, const double gi[3], const double hp[3], const double ub[3],
                          C3 e, int g, double diff[3], double adv[3]) {
  CCExec space;
  const double g0 = gi[0], g1 = gi[1], g2 = gi[2], hp0 = hp[0], hp1 = hp[1], hp2 = hp[2];
  const double h0 = 0.5 * hp[0], h1 = 0.5 * hp[1], h2 = 0.5 * hp[2], u0 = ub[0], u1 = ub[1],
               u2 = ub[2];
  double d0 = 0.0, d1 = 0.0, d2 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_mean_flux_diff",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2, double& p3) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        if (unk(i - 1) > 0.5)
          p1 += -lam * sax(i) * ((th(i) - th(i - 1)) + g0 * hp0) / hp0;
        if (unk(i - sy) > 0.5)
          p2 += -lam * say(i) * ((th(i) - th(i - sy)) + g1 * hp1) / hp1;
        if (unk(i - sz) > 0.5)
          p3 += -lam * saz(i) * ((th(i) - th(i - sz)) + g2 * hp2) / hp2;
      },
      d0, d1, d2);
  diff[0] = d0;
  diff[1] = d1;
  diff[2] = d2;
  adv[0] = adv[1] = adv[2] = 0.0;
  if (!advecting)
    return;
  double a0 = 0.0, a1 = 0.0, a2 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_mean_flux_adv",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2, double& p3) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(unk(i) > 0.5))
          return;
        const double hq[3] = {h0, h1, h2};
        double u[3];
        cellVelocity(phx, phy, phz, i, sy, sz, hq, u);
        const double k = kappa(i), t = th(i);
        p1 += (u[0] - k * u0) * t;
        p2 += (u[1] - k * u1) * t;
        p3 += (u[2] - k * u2) * t;
      },
      a0, a1, a2);
  adv[0] = a0;
  adv[1] = a1;
  adv[2] = a2;
}
 
// ---- conjugate transport (WO-7; §1.1, §1.3, §1.4, §2.7, §3.4, §4) -------------------------------
//
// Two fields on one grid: the fluid c and the solid psi_s = c_s / K (the psi form of §1.1, in which
// the interface conditions are flux continuity and q = (psi_s - psi_f)/R_c). Per phase the 7 bands
// of §4.1; the solid's hold C_s K kappa_s idt and the face diffusion Lam_face w_a a_s, a_s = 1 -
// a^snap, on faces whose two cells are solid unknowns (Lam_face the harmonic mean of the two
// cells' Lam_s when their materials differ). The overlay adds, per conjugate facet, cc = alpha G_c
// times (u_pf - u_ps) to the fluid row and its negative to the solid row of the same cell: one
// number entering two rows with opposite signs (§1.6.1).
 
struct MaterialTable {
  Kokkos::View<int*, CCMem> conj;
  Kokkos::View<double*, CCMem> lam, ck, kp, rc;
  int nb = 0;
};
 
KOKKOS_INLINE_FUNCTION int bodyIndex(double m, int nb) {
  const int b = (int)m;
  return (b < 0 || b >= nb) ? 0 : b;
}
 
KOKKOS_INLINE_FUNCTION double solidFaceLam(const Kokkos::View<double*, CCMem>& lam, int ma,
                                           int mb) {
  const double la = lam(ma);
  if (ma == mb)
    return la;
  const double lb = lam(mb);
  return (la + lb) > 0.0 ? 2.0 * la * lb / (la + lb) : 0.0;
}
 
inline void buildSolidMaterial(CCField mat, CCField sunk, CCConst kappa, CCConst sax, CCConst say,
                               CCConst saz, const scg::ScalarFacetOverlay& fo,
                               Kokkos::View<const int*, CCMem> owner, bool hasOwner,
                               const MaterialTable& mt, C3 e, int g, double vol[2]) {
  const int nx = e.x - 2 * g, ny = e.y - 2 * g;
  Kokkos::deep_copy(CCExec(), mat, -1.0);
  Kokkos::deep_copy(CCExec(), sunk, 0.0);  // first: the conjugate-facet marker
  ccFor3(
      "peclet::flow::sco_solid_owner", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        const long ii = (long)(x - g) + (long)(y - g) * nx + (long)(z - g) * (long)nx * ny;
        mat(i) = hasOwner ? (double)owner(ii) : 0.0;
      });
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  auto alpha = fo.alpha;
  auto body = fo.body;
  auto conj = mt.conj;
  const int nb = mt.nb;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_solid_facet_material", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        int best = -1, bestAny = -1;
        double aBest = -1.0, aAny = -1.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f) {
          const int b = bodyIndex((double)body(f), nb);
          if (alpha(f) > aAny) {
            aAny = alpha(f);
            bestAny = b;
          }
          if (conj(b) != 0 && alpha(f) > aBest) {
            aBest = alpha(f);
            best = b;
          }
        }
        mat(cutCell(c)) = (double)(best >= 0 ? best : bestAny);
        sunk(cutCell(c)) = best >= 0 ? 1.0 : 0.0;  // the cell has a conjugate facet
      });
  space.fence();
  double v0 = 0.0, v1 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_solid_flags",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& pv, double& ps) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        const double m = mat(i);
        const bool cand = kappa(i) < 1.0 && m >= 0.0 && conj(bodyIndex(m, nb)) != 0;
        const double asSum =
            (((1.0 - sax(i)) + (1.0 - sax(i + 1))) + ((1.0 - say(i)) + (1.0 - say(i + sy)))) +
            ((1.0 - saz(i)) + (1.0 - saz(i + sz)));
        const bool u = cand && (asSum > 0.0 || sunk(i) > 0.5);
        if (cand) {
          pv += 1.0 - kappa(i);
          if (!u)
            ps += 1.0 - kappa(i);
        }
        sunk(i) = u ? 1.0 : 0.0;
      },
      v0, v1);
  vol[0] = v0;
  vol[1] = v1;
}
 
inline void buildSolidProbes(Kokkos::View<double*, CCMem> sS,
                             Kokkos::View<int* [8], Kokkos::LayoutLeft, CCMem> offS,
                             Kokkos::View<double* [8], Kokkos::LayoutLeft, CCMem> wS,
                             Kokkos::View<std::uint8_t*, CCMem> rungS,
                             Kokkos::View<const std::uint8_t*, CCMem> conjF,
                             const scg::ScalarFacetOverlay& fo, CCConst sdf, CCConst sunk, C3 e,
                             const double hp[3]) {
  CCExec space;
  const long sy = e.x, sz = (long)e.x * e.y;
  const double h0 = hp[0], h1 = hp[1], h2 = hp[2];
  auto fcell = fo.cell;
  auto fnormal = fo.normal;
  auto fcen = fo.centroid;
  Kokkos::parallel_for(
      "peclet::flow::sco_solid_probes", Kokkos::RangePolicy<CCExec>(space, 0, fo.n),
      KOKKOS_LAMBDA(const long f) {
        if (conjF(f) == 0) {
          sS(f) = 0.0;
          rungS(f) = 0;
          for (int k = 0; k < 8; ++k) {
            offS(f, k) = 0;
            wS(f, k) = 0.0;
          }
          return;
        }
        const double h[3] = {h0, h1, h2};
        const double xw[3] = {fcen(f, 0), fcen(f, 1), fcen(f, 2)};
        const double n[3] = {-fnormal(f, 0), -fnormal(f, 1), -fnormal(f, 2)};
        const scg::BlockLookup L{sdf, sunk, fcell(f), sy, sz, -1.0};
        const peclet::core::scheme::ProbeResult r = peclet::core::scheme::buildProbe(xw, n, h, L);
        sS(f) = r.s;
        rungS(f) = (std::uint8_t)r.rung;
        for (int k = 0; k < 8; ++k) {
          const bool used = k < r.n;
          offS(f, k) =
              used ? (int)((long)r.off[k][0] + (long)r.off[k][1] * sy + (long)r.off[k][2] * sz) : 0;
          wS(f, k) = used ? r.w[k] : 0.0;
        }
      });
  space.fence();
}
 
inline void solidFaceProducts(CCField spx, CCField spy, CCField spz, CCConst sunk, CCConst mat,
                              CCConst sax, CCConst say, CCConst saz, const MaterialTable& mt,
                              C3 e) {
  auto lam = mt.lam;
  const int nb = mt.nb;
  ccFor3(
      "peclet::flow::sco_solid_face_products", C3{0, 0, 0}, e, KOKKOS_LAMBDA(int x, int y, int z) {
        const long st[3] = {1, (long)e.x, (long)e.x * e.y};
        const long i = (long)x + (long)y * st[1] + (long)z * st[2];
        const int xa[3] = {x, y, z};
        const bool ui = sunk(i) > 0.5;
        double P[3];
        for (int a = 0; a < 3; ++a) {
          const long j = i - st[a];
          if (xa[a] < 1 || !ui || !(sunk(j) > 0.5)) {
            P[a] = 0.0;
            continue;
          }
          const double ap = a == 0 ? sax(i) : (a == 1 ? say(i) : saz(i));
          P[a] = solidFaceLam(lam, bodyIndex(mat(i), nb), bodyIndex(mat(j), nb)) * (1.0 - ap);
        }
        spx(i) = P[0];
        spy(i) = P[1];
        spz(i) = P[2];
      });
}
 
inline void solidMass(CCField ms, CCConst kappa, CCConst sunk, CCConst mat, const MaterialTable& mt,
                      double idt, C3 e, int g) {
  auto ck = mt.ck;
  const int nb = mt.nb;
  Kokkos::deep_copy(CCExec(), ms, 0.0);
  ccFor3(
      "peclet::flow::sco_solid_mass", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (sunk(i) > 0.5)
          ms(i) = ck(bodyIndex(mat(i), nb)) * (1.0 - kappa(i)) * idt;
      });
}
 
inline void solidBands(CCField AC, CCField AW, CCField AE, CCField AS, CCField AN, CCField AB,
                       CCField AT, CCConst ms, CCConst sunk, CCConst spx, CCConst spy, CCConst spz,
                       const double w[3], C3 e, int g) {
  const double wx = w[0], wy = w[1], wz = w[2];
  ccFor3(
      "peclet::flow::sco_solid_bands", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(sunk(i) > 0.5)) {
          AC(i) = 1.0;
          AW(i) = AE(i) = AS(i) = AN(i) = AB(i) = AT(i) = 0.0;
          return;
        }
        const double tw = wx * spx(i), te = wx * spx(i + sx);
        const double ts = wy * spy(i), tn = wy * spy(i + sy);
        const double tb = wz * spz(i), tt = wz * spz(i + sz);
        AW(i) = -tw;
        AE(i) = -te;
        AS(i) = -ts;
        AN(i) = -tn;
        AB(i) = -tb;
        AT(i) = -tt;
        AC(i) = ms(i) + (((tw + te) + (ts + tn)) + (tb + tt));
      });
}
 
inline void solidRhs(CCField bs, CCConst ms, CCConst sOld, C3 e, int g) {
  ccFor3(
      "peclet::flow::sco_solid_rhs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        bs(i) = ms(i) * sOld(i);
      });
}
 
inline void conjugateCoefficients(Kokkos::View<double*, CCMem> cc,
                                  const scg::ScalarFacetOverlay& fo,
                                  Kokkos::View<const std::uint8_t*, CCMem> conjF,
                                  Kokkos::View<const double*, CCMem> sS, CCConst unk, CCConst sunk,
                                  const MaterialTable& mt, double lamF) {
  auto fcell = fo.cell;
  auto falpha = fo.alpha;
  auto fbody = fo.body;
  auto sF = fo.sF;
  auto lam = mt.lam;
  auto rc = mt.rc;
  const int nb = mt.nb;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_conj_coeffs", Kokkos::RangePolicy<CCExec>(space, 0, fo.n),
      KOKKOS_LAMBDA(const long f) {
        const long i = fcell(f);
        if (conjF(f) == 0 || !(unk(i) > 0.5) || !(sunk(i) > 0.5)) {
          cc(f) = 0.0;
          return;
        }
        const int b = bodyIndex((double)fbody(f), nb);
        cc(f) = falpha(f) *
                peclet::core::scheme::interfaceConductance(sF(f), lamF, rc(b), sS(f), lam(b));
      });
  space.fence();
}
 
inline void conjOverlayApply(CCField yf, CCField ys, CCConst xf, CCConst xs,
                             const scg::ScalarFacetOverlay& fo,
                             Kokkos::View<const double*, CCMem> cw,
                             Kokkos::View<const double*, CCMem> cc,
                             Kokkos::View<int* [8], Kokkos::LayoutLeft, CCMem> offS,
                             Kokkos::View<double* [8], Kokkos::LayoutLeft, CCMem> wS,
                             bool subtract = false) {
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  auto offF = fo.offF;
  auto wF = fo.wF;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_conj_overlay_apply", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        const long i = cutCell(c);
        double af = 0.0, as = 0.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f) {
          double upf = 0.0;
          for (int k = 0; k < 8; ++k)
            upf += wF(f, k) * xf(i + (long)offF(f, k));
          af += cw(f) * upf;
          if (cc(f) != 0.0) {
            double ups = 0.0;
            for (int k = 0; k < 8; ++k)
              ups += wS(f, k) * xs(i + (long)offS(f, k));
            const double q = cc(f) * (upf - ups);
            af += q;
            as -= q;
          }
        }
        if (subtract) {
          yf(i) -= af;
          ys(i) -= as;
        } else {
          yf(i) += af;
          ys(i) += as;
        }
      });
}
 
inline void conjSurrogate(CCField SACf, CCField SACs, CCField W0, CCConst ACf, CCConst ACs,
                          const scg::ScalarFacetOverlay& fo, Kokkos::View<const double*, CCMem> cw,
                          Kokkos::View<const double*, CCMem> cc) {
  Kokkos::deep_copy(CCExec(), SACf, ACf);
  Kokkos::deep_copy(CCExec(), SACs, ACs);
  Kokkos::deep_copy(CCExec(), W0, 0.0);
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_conj_surrogate", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        double aw = 0.0, ac = 0.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f) {
          aw += cw(f) + cc(f);
          ac += cc(f);
        }
        const long i = cutCell(c);
        SACf(i) += aw;
        SACs(i) += ac;
        W0(i) = ac;
      });
  space.fence();
}
 
inline void conjMeanGradientShift(CCField bf, CCField bs, const scg::ScalarFacetOverlay& fo,
                                  Kokkos::View<const double*, CCMem> cc,
                                  Kokkos::View<int* [8], Kokkos::LayoutLeft, CCMem> offS,
                                  Kokkos::View<double* [8], Kokkos::LayoutLeft, CCMem> wS,
                                  const MeanGradPosition& pos) {
  auto cutCell = fo.cutCell;
  auto cfs = fo.cellFacetStart;
  auto offF = fo.offF;
  auto wF = fo.wF;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_conj_mg_shift", Kokkos::RangePolicy<CCExec>(space, 0, fo.nCut),
      KOKKOS_LAMBDA(const long c) {
        const long i = cutCell(c);
        double d = 0.0;
        for (int f = cfs(c); f < cfs(c + 1); ++f) {
          if (cc(f) == 0.0)
            continue;
          double gf = 0.0, gs = 0.0;
          for (int k = 0; k < 8; ++k) {
            gf += wF(f, k) * pos(i + (long)offF(f, k));
            gs += wS(f, k) * pos(i + (long)offS(f, k));
          }
          d += cc(f) * (gs - gf);
        }
        bf(i) += d;
        bs(i) -= d;
      });
  space.fence();
}
 
inline void solidMeanGradientFaceRhs(CCField bs, CCConst sunk, CCConst spx, CCConst spy,
                                     CCConst spz, const double gs[3], C3 e, int g) {
  const double gx = gs[0], gy = gs[1], gz = gs[2];
  ccFor3(
      "peclet::flow::sco_solid_mg_face_rhs", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(sunk(i) > 0.5))
          return;
        bs(i) += (gx * (spx(i + sx) - spx(i)) + gy * (spy(i + sy) - spy(i))) +
                 gz * (spz(i + sz) - spz(i));
      });
}
 
inline void solidResidualFluxForm(CCField r, CCConst x, CCConst b, CCConst ms, CCConst sunk,
                                  CCConst spx, CCConst spy, CCConst spz, const double w[3], C3 e,
                                  int g) {
  const double wx = w[0], wy = w[1], wz = w[2];
  ccFor3(
      "peclet::flow::sco_solid_residual_flux", 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 (!(sunk(i) > 0.5)) {
          r(i) = 0.0;
          return;
        }
        const double ci = x(i);
        const double flux =
            ((wx * spx(i) * (ci - x(i - sx)) + wx * spx(i + sx) * (ci - x(i + sx))) +
             (wy * spy(i) * (ci - x(i - sy)) + wy * spy(i + sy) * (ci - x(i + sy)))) +
            (wz * spz(i) * (ci - x(i - sz)) + wz * spz(i + sz) * (ci - x(i + sz)));
        r(i) = b(i) - ms(i) * ci - flux;
      });
}
 
inline void solidMomentsLocal(CCConst xs, CCConst kappa, CCConst sunk, CCConst mat,
                              const MaterialTable& mt, C3 e, int g, double& m, double& w) {
  auto ck = mt.ck;
  const int nb = mt.nb;
  CCExec space;
  double s1 = 0.0, s2 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_solid_moments",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        if (sunk(i) > 0.5) {
          const double wk = ck(bodyIndex(mat(i), nb)) * (1.0 - kappa(i));
          p1 += wk * xs(i);
          p2 += wk;
        }
      },
      s1, s2);
  m = s1;
  w = s2;
}
 
inline void conjFacetFlux(Kokkos::View<double*, CCMem> up, Kokkos::View<double*, CCMem> us,
                          Kokkos::View<double*, CCMem> qin, CCConst xf, CCConst xs,
                          const scg::ScalarFacetOverlay& fo, Kokkos::View<const double*, CCMem> cw,
                          Kokkos::View<const double*, CCMem> rw,
                          Kokkos::View<const double*, CCMem> cc,
                          Kokkos::View<int* [8], Kokkos::LayoutLeft, CCMem> offS,
                          Kokkos::View<double* [8], Kokkos::LayoutLeft, CCMem> wS,
                          const MeanGradPosition& pos, bool mgOn) {
  auto fcell = fo.cell;
  auto offF = fo.offF;
  auto wF = fo.wF;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::sco_conj_facet_flux", Kokkos::RangePolicy<CCExec>(space, 0, fo.n),
      KOKKOS_LAMBDA(const long f) {
        const long i = fcell(f);
        double u = 0.0, v = 0.0;
        for (int k = 0; k < 8; ++k)
          u += wF(f, k) * xf(i + (long)offF(f, k));
        for (int k = 0; k < 8; ++k)
          v += wS(f, k) * xs(i + (long)offS(f, k));
        up(f) = u;
        us(f) = v;
        if (cc(f) != 0.0) {
          // in mean-gradient mode the probes hold theta: the flux is that of theta + G.x (§1.5)
          double d = v - u;
          if (mgOn) {
            double gf = 0.0, gs = 0.0;
            for (int k = 0; k < 8; ++k) {
              gf += wF(f, k) * pos(i + (long)offF(f, k));
              gs += wS(f, k) * pos(i + (long)offS(f, k));
            }
            d += gs - gf;
          }
          qin(f) = cc(f) * d;
        } else {
          qin(f) = rw(f) - cw(f) * u;
        }
      });
  space.fence();
}
 
template <class R>
inline void dotTwoReduce(CCConst af, CCConst bf, CCConst as, CCConst bs, C3 e, int g, R&& result) {
  ccReduce3(
      "peclet::flow::sco_dot_two", C3{g, g, g}, C3{e.x - g, e.y - g, e.z - g},
      KOKKOS_LAMBDA(int x, int y, int z, double& acc) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        acc += af(i) * bf(i) + as(i) * bs(i);
      },
      std::forward<R>(result));
}
inline double dotTwoLocal(CCConst af, CCConst bf, CCConst as, CCConst bs, C3 e, int g) {
  double s = 0.0;
  dotTwoReduce(af, bf, as, bs, e, g, s);
  return s;
}
 
template <class R1, class R2>
inline void dot2TwoReduce(CCConst af, CCConst bf, CCConst cf, CCConst as, CCConst bs, CCConst cs,
                          C3 e, int g, R1&& r1, R2&& r2) {
  CCExec space;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_dot2_two", MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2) {
        const long i = (long)x + (long)y * e.x + (long)z * (long)e.x * e.y;
        p1 += af(i) * bf(i) + as(i) * bs(i);
        p2 += cf(i) * cf(i) + cs(i) * cs(i);
      },
      std::forward<R1>(r1), std::forward<R2>(r2));
}
inline void dot2TwoLocal(CCConst af, CCConst bf, CCConst cf, CCConst as, CCConst bs, CCConst cs,
                         C3 e, int g, double& ab, double& cc) {
  double s1 = 0.0, s2 = 0.0;
  dot2TwoReduce(af, bf, cf, as, bs, cs, e, g, s1, s2);
  ab = s1;
  cc = s2;
}
 
inline void solidMeanFluxLocal(CCConst th, CCConst sunk, CCConst spx, CCConst spy, CCConst spz,
                               const double gi[3], const double hp[3], C3 e, int g,
                               double diff[3]) {
  CCExec space;
  const double g0 = gi[0], g1 = gi[1], g2 = gi[2], hp0 = hp[0], hp1 = hp[1], hp2 = hp[2];
  double d0 = 0.0, d1 = 0.0, d2 = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::sco_solid_mean_flux",
      MDRange3<CCExec>(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int x, int y, int z, double& p1, double& p2, double& p3) {
        const long sy = e.x, sz = (long)e.x * e.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        if (!(sunk(i) > 0.5))
          return;
        if (spx(i) != 0.0)
          p1 += -spx(i) * ((th(i) - th(i - 1)) + g0 * hp0) / hp0;
        if (spy(i) != 0.0)
          p2 += -spy(i) * ((th(i) - th(i - sy)) + g1 * hp1) / hp1;
        if (spz(i) != 0.0)
          p3 += -spz(i) * ((th(i) - th(i - sz)) + g2 * hp2) / hp2;
      },
      d0, d1, d2);
  diff[0] = d0;
  diff[1] = d1;
  diff[2] = d2;
}
 
}  // namespace sco
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_SCALAR_CUTCELL_OPERATOR_HPP