ghost_projection.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_GHOST_PROJECTION_HPP
#define PECLET_FLOW_GHOST_PROJECTION_HPP
 
#include <cstdint>
#include <Kokkos_Core.hpp>
 
#include "cut_cell_ibm.hpp"  // poly_D / poly_Nc / poly_N_nb (momentum IBM)
#include "mac_cutcell.hpp"   // CCField/CCConst, C3, CCExec
// The PURE per-face closure pieces (state cascade, closure weights, row fill) are shared with
// core's AMR octree band — lifted verbatim into peclet::core::scheme (float arithmetic identical
// to the definitions that lived here; Kokkos_Core.hpp is included above, so they compile as
// KOKKOS_INLINE_FUNCTION). The grid-specific parts (overlay SoA, periodic wrap, delta kernels,
// gpCenterGrad) stay here.
#include "gauge_exact_gradient.hpp"  // gpCenterGrad (moved out: the production scheme needs it)
#include "peclet/core/scheme/ghost_closure.hpp"
#include "policy.hpp"
 
namespace peclet::flow {
 
using peclet::core::scheme::GP_THETA_MIN;
using enum peclet::core::scheme::GpState;
using peclet::core::scheme::gpClassifyFace;
using peclet::core::scheme::GpFace;
using peclet::core::scheme::gpFillRow;
using peclet::core::scheme::gpOrderWeights;
using peclet::core::scheme::GpState;
 
template <class Space, class Real = float>
struct GpOverlayT {
  Kokkos::View<int*, Space> cell;               // packed INNER flat index x + y*nx + z*nx*ny
  Kokkos::View<Real*, Space> rescale;           // rho = min(1, min_f D_f) of the MATRIX weights
  Kokkos::View<int8_t*, Space> coupled;         // 1 if the row has any phi coupling at all
  Kokkos::View<int8_t*, Space> state;           // [slot*6+k]
  Kokkos::View<Real*, Space> th;                // [slot*6+k] (parity/diagnostics)
  Kokkos::View<Real*, Space> w_bc, w_n1, w_n2;  // [slot*6+k] RHS/diagnostic closure weights
  Kokkos::View<Real*, Space> wm_n1, wm_n2;      // [slot*6+k] matrix (implicit phi) weights
};
template <class Real>
using GpOverlayReal = GpOverlayT<CCMem, Real>;
using GpOverlay = GpOverlayReal<float>;
 
template <class Real = float>
inline GpOverlayReal<Real> gpMakeOverlay(long n) {
  GpOverlayReal<Real> ov;
  ov.cell = Kokkos::View<int*, CCMem>("gp_cell", n);
  ov.rescale = Kokkos::View<Real*, CCMem>("gp_rescale", n);
  ov.coupled = Kokkos::View<int8_t*, CCMem>("gp_coupled", n);
  ov.state = Kokkos::View<int8_t*, CCMem>("gp_state", 6 * n);
  ov.th = Kokkos::View<Real*, CCMem>("gp_th", 6 * n);
  ov.w_bc = Kokkos::View<Real*, CCMem>("gp_wbc", 6 * n);
  ov.w_n1 = Kokkos::View<Real*, CCMem>("gp_wn1", 6 * n);
  ov.w_n2 = Kokkos::View<Real*, CCMem>("gp_wn2", 6 * n);
  ov.wm_n1 = Kokkos::View<Real*, CCMem>("gp_wmn1", 6 * n);
  ov.wm_n2 = Kokkos::View<Real*, CCMem>("gp_wmn2", 6 * n);
  return ov;
}
 
KOKKOS_INLINE_FUNCTION int gpWrap(int v, int n) {
  v %= n;
  return v < 0 ? v + n : v;
}
 
template <class Real>
inline int buildGpOverlay(CCConst sdf, C3 ext, int g, C3 nn, const GpOverlayReal<Real>& ov,
                          Kokkos::View<int*, CCMem> idMap, Kokkos::View<int, CCMem> counter,
                          int matrixOrder = 2, int rhsOrder = 2, CCConst tx = CCConst(),
                          CCConst ty = CCConst(), CCConst tz = CCConst(), bool useGhost = false) {
  CCExec space;
  Kokkos::deep_copy(space, counter, 0);
  Kokkos::deep_copy(space, idMap, -1);
  const bool hasEx = tx.size() > 0;
  const bool ug = useGhost;
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::gp_build_overlay", MD(space, {0, 0, 0}, {nn.x, nn.y, nn.z}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        auto S = [&](int dx, int dy, int dz) {
          const long i = (long)((ug ? x + dx : gpWrap(x + dx, nn.x)) + g) +
                         (long)((ug ? y + dy : gpWrap(y + dy, nn.y)) + g) * ext.x +
                         (long)((ug ? z + dz : gpWrap(z + dz, nn.z)) + g) * (long)ext.x * ext.y;
          return (float)sdf(i);  // PRECISION-EXEMPT: feeds core::scheme::gpFillRow, float-hardcoded
        };
        const float sc = S(0, 0, 0);
        if (sc < 0.0f)
          return;  // solid-centered: decoupled row (phi = 0), not in the overlay
        float F[3][4], Cq[3][5];
        const int d[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
        for (int a = 0; a < 3; ++a) {
          for (int q = -2; q <= 2; ++q)
            Cq[a][q + 2] = S(q * d[a][0], q * d[a][1], q * d[a][2]);
          for (int m = -1; m <= 2; ++m)  // face i+m = mean of centers i+m-1, i+m
            F[a][m + 1] = 0.5f * (Cq[a][m + 1] + Cq[a][m + 2]);
        }
        // cheap pre-check: fully interior rows (all six 1st neighbors fluid on faces+centers)
        bool clean = true;
        for (int a = 0; a < 3; ++a)
          clean =
              clean && F[a][1] >= 0.0f && F[a][2] >= 0.0f && Cq[a][1] >= 0.0f && Cq[a][3] >= 0.0f;
        if (clean)
          return;
        const int inner = x + y * nn.x + z * nn.x * nn.y;
        // Optional analytic-SDF exact crossings (setExactCrossings): t_a(j) = exact crossing
        // fraction from same-component face point j toward j + e_a along axis a (NaN = none).
        // minus side: std theta = 1 - t_a(i);      sliver theta = 2 - t_a(i-1)
        // plus  side: std theta = t_a(i);          sliver theta = 1 + t_a(i+1)
        float exStd[6], exSliver[6];
        if (hasEx) {
          const CCConst* ta[3] = {&tx, &ty, &tz};
          for (int a = 0; a < 3; ++a) {
            auto T = [&](int m) {
              const int cx = a == 0 ? gpWrap(x + m, nn.x) : x;
              const int cy = a == 1 ? gpWrap(y + m, nn.y) : y;
              const int cz = a == 2 ? gpWrap(z + m, nn.z) : z;
              const long ii = (long)cx + (long)cy * nn.x + (long)cz * (long)nn.x * nn.y;
              return (float)(*ta[a])(ii);  // PRECISION-EXEMPT: ditto
            };
            exStd[2 * a + 1] = 1.0f - T(0);
            exSliver[2 * a + 1] = 2.0f - T(-1);
            exStd[2 * a] = T(0);
            exSliver[2 * a] = 1.0f + T(1);
          }
        }
        const int slot = Kokkos::atomic_fetch_add(&counter(), 1);
        if (!gpFillRow(ov, slot, inner, F, Cq, matrixOrder, rhsOrder, hasEx ? exStd : nullptr,
                       hasEx ? exSliver : nullptr)) {
          // all faces COUPLED after all (pre-check was conservative): release the slot lazily by
          // marking it inert (zero-weight row on its own cell).
          ov.cell(slot) = inner;
          ov.rescale(slot) = 1.0f;
          ov.coupled(slot) = 1;
          for (int k = 0; k < 6; ++k) {
            ov.state(slot * 6 + k) = GP_COUPLED;
            ov.th(slot * 6 + k) = 1.0f;
            ov.w_bc(slot * 6 + k) = 0.0f;
            ov.w_n1(slot * 6 + k) = 0.0f;
            ov.w_n2(slot * 6 + k) = 0.0f;
            ov.wm_n1(slot * 6 + k) = 0.0f;
            ov.wm_n2(slot * 6 + k) = 0.0f;
          }
        } else {
          idMap(inner) = slot;
        }
      });
  int cnt = 0;
  Kokkos::deep_copy(cnt, counter);
  return cnt;
}
 
inline void gpBinaryOpenness(CCField ox, CCField oy, CCField oz, CCConst sdf, C3 ext) {
  CCExec space;
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::gp_binary_openness", MD(space, {0, 0, 0}, {ext.x, ext.y, ext.z}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sy = ext.x, sz = (long)ext.x * ext.y;
        const long i = (long)x + (long)y * sy + (long)z * sz;
        auto cl = [](int v, int n) { return v < 0 ? 0 : (v >= n ? n - 1 : v); };
        auto S = [&](int xx, int yy, int zz) {
          return sdf((long)cl(xx, ext.x) + (long)cl(yy, ext.y) * sy + (long)cl(zz, ext.z) * sz);
        };
        const double sc = S(x, y, z);
        const double sw = S(x - 1, y, z), ss = S(x, y - 1, z), sb = S(x, y, z - 1);
        ox(i) = (0.5 * (sw + sc) >= 0.0 && sw >= 0.0 && sc >= 0.0) ? 1.0 : 0.0;
        oy(i) = (0.5 * (ss + sc) >= 0.0 && ss >= 0.0 && sc >= 0.0) ? 1.0 : 0.0;
        oz(i) = (0.5 * (sb + sc) >= 0.0 && sb >= 0.0 && sc >= 0.0) ? 1.0 : 0.0;
      });
}
 
template <class Real>
inline void gpApplyDelta(CCField y, CCConst x, const GpOverlayReal<Real>& ov, int nOv, C3 nn,
                         C3 extY, int gbY, C3 extX, int gbX, bool useGhost = false) {
  if (nOv <= 0)
    return;
  CCExec space;
  const bool ug = useGhost;
  Kokkos::parallel_for(
      "peclet::flow::gp_apply_delta", Kokkos::RangePolicy<CCExec>(space, 0, nOv),
      KOKKOS_LAMBDA(int s) {
        const int inner = ov.cell(s);
        const int ix = inner % nn.x, iy = (inner / nn.x) % nn.y, iz = inner / (nn.x * nn.y);
        auto X = [&](int a, int q) {  // phi at cell offset q along axis a
          const int cx = a == 0 ? (ug ? ix + q : gpWrap(ix + q, nn.x)) : ix;
          const int cy = a == 1 ? (ug ? iy + q : gpWrap(iy + q, nn.y)) : iy;
          const int cz = a == 2 ? (ug ? iz + q : gpWrap(iz + q, nn.z)) : iz;
          return x((long)(cx + gbX) + (long)(cy + gbX) * extX.x +
                   (long)(cz + gbX) * (long)extX.x * extX.y);
        };
        const long r =
            (long)(ix + gbY) + (long)(iy + gbY) * extY.x + (long)(iz + gbY) * (long)extY.x * extY.y;
        double delta = 0.0;
        for (int k = 0; k < 6; ++k) {
          const int8_t st = ov.state(s * 6 + k);
          if (st != GP_QUAD && st != GP_LIN)
            continue;
          const int a = k / 2;
          const int sgn = (k & 1) ? -1 : 1;  // odd k = minus side
          const int mn = (k & 1) ? 1 : 0;    // near-face relative index
          const int mf = (k & 1) ? 2 : -1;   // far-face relative index
          const double w1 = ov.wm_n1(s * 6 + k), w2 = ov.wm_n2(s * 6 + k);
          delta += sgn * w1 * (X(a, mn - 1) - X(a, mn));
          if (st == GP_QUAD && w2 != 0.0)
            delta += sgn * w2 * (X(a, mf - 1) - X(a, mf));
        }
        const double rho = ov.rescale(s);
        y(r) = ov.coupled(s) ? rho * (y(r) + delta) : 0.0;
      });
}
 
template <class Real>
inline void gpDivergDelta(CCField d, CCConst u, CCConst v, CCConst w, const GpOverlayReal<Real>& ov,
                          int nOv, C3 nn, C3 extb, int gb, bool useGhost = false) {
  if (nOv <= 0)
    return;
  CCExec space;
  const bool ug = useGhost;
  Kokkos::parallel_for(
      "peclet::flow::gp_diverg_delta", Kokkos::RangePolicy<CCExec>(space, 0, nOv),
      KOKKOS_LAMBDA(int s) {
        const int inner = ov.cell(s);
        const int ix = inner % nn.x, iy = (inner / nn.x) % nn.y, iz = inner / (nn.x * nn.y);
        auto U = [&](int a, int m) {  // face-field value at face index i+m along axis a
          const int cx = a == 0 ? (ug ? ix + m : gpWrap(ix + m, nn.x)) : ix;
          const int cy = a == 1 ? (ug ? iy + m : gpWrap(iy + m, nn.y)) : iy;
          const int cz = a == 2 ? (ug ? iz + m : gpWrap(iz + m, nn.z)) : iz;
          const long i =
              (long)(cx + gb) + (long)(cy + gb) * extb.x + (long)(cz + gb) * (long)extb.x * extb.y;
          return a == 0 ? u(i) : (a == 1 ? v(i) : w(i));
        };
        const long r =
            (long)(ix + gb) + (long)(iy + gb) * extb.x + (long)(iz + gb) * (long)extb.x * extb.y;
        double dd = d(r);
        for (int k = 0; k < 6; ++k) {
          const int8_t st = ov.state(s * 6 + k);
          if (st == GP_COUPLED)
            continue;
          const int a = k / 2;
          const int sgn = (k & 1) ? -1 : 1;
          const int mg = (k & 1) ? 0 : 1;  // the closed face's own index
          const int mn = (k & 1) ? 1 : 0;
          const int mf = (k & 1) ? 2 : -1;
          if (st == GP_EXPLICIT) {
            dd += sgn * U(a, mg);  // sliver without crossing: explicit u* flux
            continue;
          }
          if (st == GP_BC_ONLY)
            continue;  // u_bc = 0 (v1); w_bc kept in the overlay for moving walls
          double val = ov.w_n1(s * 6 + k) * U(a, mn);  // + w_bc*u_bc (= 0)
          if (st == GP_QUAD)
            val += ov.w_n2(s * 6 + k) * U(a, mf);
          dd += sgn * val;
        }
        d(r) = ov.coupled(s) ? ov.rescale(s) * dd : 0.0;
      });
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_GHOST_PROJECTION_HPP