star_elimination.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_STAR_ELIMINATION_HPP
#define PECLET_FLOW_STAR_ELIMINATION_HPP
 
#include <Kokkos_Core.hpp>
 
#include "mac_cutcell.hpp"  // CCField/CCConst, C3, CCExec, CCMem
#include "policy.hpp"
 
namespace peclet::flow {
 
struct StarOverlay {
  Kokkos::View<int*, CCMem> cell;
  // [slot*6+k]. DOUBLE since P2: this is the masked face openness, and the openness source
  // (Level::ox/oy/oz) is double — storing it float made the star row a float-rounded copy of a
  // coefficient the solver already holds exactly. Gate-off reads it back through (float) so the
  // default path is bitwise what it was; see starAval().
  Kokkos::View<double*, CCMem> a;
};
 
inline StarOverlay starMakeOverlay(long n) {
  StarOverlay ov;
  ov.cell = Kokkos::View<int*, CCMem>("star_cell", n);
  ov.a = Kokkos::View<double*, CCMem>("star_a", 6 * n);
  return ov;
}
 
KOKKOS_INLINE_FUNCTION double starAval(const StarOverlay& ov, long k, bool exact) {
  const double v = ov.a(k);
  // Deliberate bit-compat gate (not an operator-precision leak): OFF reproduces the float-rounded
  // value the overlay used to store (bitwise unchanged default path), ON returns the exact double.
  return exact ? v : (double)(float)v;  // PRECISION-EXEMPT: bit-compat gate, see above
}
 
KOKKOS_INLINE_FUNCTION int starWrap(int v, int n) {
  v %= n;
  return v < 0 ? v + n : v;
}
 
inline int buildStarOverlay(CCConst sdf, CCConst ox, CCConst oy, CCConst oz, C3 ext, int g, C3 nn,
                            const StarOverlay& ov, Kokkos::View<int, CCMem> counter) {
  CCExec space;
  Kokkos::deep_copy(counter, 0);
  const bool fill = ov.cell.extent(0) > 0;
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::star_build", MD(space, {0, 0, 0}, {nn.x, nn.y, nn.z}),
      KOKKOS_LAMBDA(int x, int y, int z) {
        const long sx = 1, sy = ext.x, sz = (long)ext.x * ext.y;
        const long st[3] = {sx, sy, sz};
        const long i = (long)(x + g) + (long)(y + g) * sy + (long)(z + g) * sz;
        if (sdf(i) >= 0.0)
          return;  // fluid-centered: not eliminated
        CCConst oa[3] = {ox, oy, oz};
        double av[6];
        double D = 0.0;
        for (int a2 = 0; a2 < 3; ++a2) {
          const double op = oa[a2](i + st[a2]), om = oa[a2](i);  // +side face / -side face
          const bool fp = sdf(i + st[a2]) >= 0.0, fm = sdf(i - st[a2]) >= 0.0;
          av[2 * a2] = (fp && op > 0.0) ? op : 0.0;
          av[2 * a2 + 1] = (fm && om > 0.0) ? om : 0.0;
          D += av[2 * a2] + av[2 * a2 + 1];
        }
        if (D <= 0.0)
          return;  // no open fluid face: fully decoupled solid cell
        const int slot = Kokkos::atomic_fetch_add(&counter(), 1);
        if (fill) {
          ov.cell(slot) = x + y * nn.x + z * nn.x * nn.y;
          for (int k = 0; k < 6; ++k)
            ov.a(slot * 6 + k) = av[k];
        }
      });
  space.fence();
  auto h = Kokkos::create_mirror_view(counter);
  Kokkos::deep_copy(h, counter);
  return h();
}
 
inline void starApplyDelta(CCField y, CCConst x, const StarOverlay& ov, int nOv, C3 nn, C3 extY,
                           int gY, C3 extX, int gX, bool exact) {
  if (nOv <= 0)
    return;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::star_apply", 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 idx = [&](int cx, int cy, int cz, C3 ext, int gb) {
          return (long)(starWrap(cx, nn.x) + gb) + (long)(starWrap(cy, nn.y) + gb) * ext.x +
                 (long)(starWrap(cz, nn.z) + gb) * (long)ext.x * ext.y;
        };
        const int nb[6][3] = {{ix + 1, iy, iz}, {ix - 1, iy, iz}, {ix, iy + 1, iz},
                              {ix, iy - 1, iz}, {ix, iy, iz + 1}, {ix, iy, iz - 1}};
        double D = 0.0, num = 0.0;
        double av[6], xv[6];
        for (int k = 0; k < 6; ++k) {
          av[k] = starAval(ov, (long)s * 6 + k, exact);
          if (av[k] <= 0.0)
            continue;
          xv[k] = x(idx(nb[k][0], nb[k][1], nb[k][2], extX, gX));
          D += av[k];
          num += av[k] * xv[k];
        }
        if (exact) {
          // FLUX FORM (P2, the P1 trick applied to the star row). Algebraically identical to the
          // phibar form below --
          //   a_k (x_k - phibar) = (a_k/D)(D x_k - sum_j a_j x_j) = (a_k/D) sum_j a_j (x_k - x_j)
          // -- but every term is a coefficient times a DIFFERENCE, so a constant x is annihilated
          // BITWISE. The phibar form is not: phibar = sum(a x)/sum(a) does not return x exactly
          // for constant x even in double (the two sums round independently, then divide), so the
          // star delta leaked a nonzero row contribution into A*1 that P1's exact 7-point apply
          // could not cancel. Composing the two is the whole point of the rung: the level-0
          // operator is the 7-point part PLUS this delta, and A*1 = 0 needs both halves exact.
          for (int k = 0; k < 6; ++k) {
            if (av[k] <= 0.0)
              continue;
            double acc = 0.0;
            for (int j = 0; j < 6; ++j)
              if (av[j] > 0.0)
                acc += av[j] * (xv[k] - xv[j]);
            Kokkos::atomic_add(&y(idx(nb[k][0], nb[k][1], nb[k][2], extY, gY)), (av[k] / D) * acc);
          }
          return;
        }
        const double phibar = num / D;
        for (int k = 0; k < 6; ++k) {
          if (av[k] <= 0.0)
            continue;
          Kokkos::atomic_add(&y(idx(nb[k][0], nb[k][1], nb[k][2], extY, gY)),
                             av[k] * (xv[k] - phibar));
        }
      });
}
 
inline void starCorrectFaces(CCField uf, CCField vf, CCField wf, CCConst phi, const StarOverlay& ov,
                             int nOv, C3 nn, C3 ext, int g, C3 extP, int gP, bool exact,
                             double wx = 1.0, double wy = 1.0, double wz = 1.0) {
  if (nOv <= 0)
    return;
  CCExec space;
  Kokkos::parallel_for(
      "peclet::flow::star_correct_faces", 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 idxP = [&](int cx, int cy, int cz) {
          return (long)(starWrap(cx, nn.x) + gP) + (long)(starWrap(cy, nn.y) + gP) * extP.x +
                 (long)(starWrap(cz, nn.z) + gP) * (long)extP.x * extP.y;
        };
        auto idxF = [&](int cx, int cy, int cz) {
          return (long)(starWrap(cx, nn.x) + g) + (long)(starWrap(cy, nn.y) + g) * ext.x +
                 (long)(starWrap(cz, nn.z) + g) * (long)ext.x * ext.y;
        };
        const int nb[6][3] = {{ix + 1, iy, iz}, {ix - 1, iy, iz}, {ix, iy + 1, iz},
                              {ix, iy - 1, iz}, {ix, iy, iz + 1}, {ix, iy, iz - 1}};
        double D = 0.0, num = 0.0;
        for (int k = 0; k < 6; ++k) {
          const double a = starAval(ov, (long)s * 6 + k, exact);
          if (a <= 0.0)
            continue;
          D += a;
          num += a * phi(idxP(nb[k][0], nb[k][1], nb[k][2]));
        }
        // phibar itself IS the wanted quantity here (the eliminated cell's reconstructed value),
        // not an operator row, so there is no difference form to take -- only the aperture
        // precision matters, and it must match the operator that was actually solved.
        const double phibar = num / D;
        CCField fa[3] = {uf, vf, wf};
        const double wv[3] = {wx, wy, wz};
        for (int k = 0; k < 6; ++k) {
          if (starAval(ov, (long)s * 6 + k, exact) <= 0.0)
            continue;
          const int a2 = k / 2;
          if ((k & 1) == 0) {
            // + side: face is the LOW face of the + neighbour; s is the LOW cell -> += phibar
            fa[a2](idxF(nb[k][0], nb[k][1], nb[k][2])) += wv[a2] * phibar;
          } else {
            // - side: face is s's own LOW face; s is the HIGH cell -> -= phibar
            fa[a2](idxF(ix, iy, iz)) -= wv[a2] * phibar;
          }
        }
      });
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_STAR_ELIMINATION_HPP