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
 
namespace peclet::flow {
 
struct StarOverlay {
  Kokkos::View<int*, CCMem> cell;
  Kokkos::View<float*, CCMem> a;  // [slot*6+k]
};
 
inline StarOverlay starMakeOverlay(long n) {
  StarOverlay ov;
  ov.cell = Kokkos::View<int*, CCMem>("star_cell", n);
  ov.a = Kokkos::View<float*, CCMem>("star_a", 6 * n);
  return ov;
}
 
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 = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  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};
        float 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) ? (float)op : 0.0f;
          av[2 * a2 + 1] = (fm && om > 0.0) ? (float)om : 0.0f;
          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) {
  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 xv[6];
        for (int k = 0; k < 6; ++k) {
          const double a = ov.a(s * 6 + k);
          if (a <= 0.0)
            continue;
          xv[k] = x(idx(nb[k][0], nb[k][1], nb[k][2], extX, gX));
          D += a;
          num += a * xv[k];
        }
        const double phibar = num / D;
        for (int k = 0; k < 6; ++k) {
          const double a = ov.a(s * 6 + k);
          if (a <= 0.0)
            continue;
          Kokkos::atomic_add(&y(idx(nb[k][0], nb[k][1], nb[k][2], extY, gY)),
                             a * (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) {
  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 = ov.a(s * 6 + k);
          if (a <= 0.0)
            continue;
          D += a;
          num += a * phi(idxP(nb[k][0], nb[k][1], nb[k][2]));
        }
        const double phibar = num / D;
        CCField fa[3] = {uf, vf, wf};
        for (int k = 0; k < 6; ++k) {
          if (ov.a(s * 6 + k) <= 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])) += phibar;
          } else {
            // - side: face is s's own LOW face; s is the HIGH cell -> -= phibar
            fa[a2](idxF(ix, iy, iz)) -= phibar;
          }
        }
      });
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_STAR_ELIMINATION_HPP