mac_bc.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_MAC_BC_HPP
#define PECLET_FLOW_MAC_BC_HPP
 
#include <Kokkos_Core.hpp>
 
namespace peclet::flow {
 
using BExec = Kokkos::DefaultExecutionSpace;
using BMem = BExec::memory_space;
using BField = Kokkos::View<double*, BMem>;
 
struct B3 {
  int x, y, z;
};
 
namespace bcdetail {
KOKKOS_INLINE_FUNCTION void axisDims(B3 ext, int (&dims)[3], long (&strides)[3]) {
  dims[0] = ext.x;
  dims[1] = ext.y;
  dims[2] = ext.z;
  strides[0] = 1;
  strides[1] = ext.x;
  strides[2] = static_cast<long>(ext.x) * ext.y;
}
}  // namespace bcdetail
 
// Fill component comp (0=u,1=v,2=w) ghosts for one domain face (axis a, side s=0 low/1 high) with a
// scalar wall velocity. fold=1 drops the tangential wall face (ghost=0, implicit diffusion).
inline void bcVelocityComp(BField f, B3 ext, int g, int a, int s, int comp, double wall, int fold,
                           BField prof = BField(), int prof_nc = 0) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int na = dims[a];
  const int bf = (s == 0) ? g : (na - g);
  const bool hasProf =
      prof.extent(0) > 0;  // per-position inlet profile (resampled to the face grid)
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_vel", MD(space, {0, 0}, {dims[b], dims[c]}), KOKKOS_LAMBDA(int p0, int p1) {
        const long base = static_cast<long>(p0) * sb + static_cast<long>(p1) * sc;
        const double wc = hasProf ? prof((static_cast<long>(p0) * prof_nc + p1) * 3 + comp) : wall;
        auto at = [&](int ia) -> double& { return f(base + static_cast<long>(ia) * sa); };
        if (comp == a) {  // normal: Dirichlet face + odd reflection
          at(bf) = wc;
          if (s == 0)
            for (int ia = 0; ia < g; ++ia)
              at(ia) = 2.0 * wc - at(2 * bf - ia);
          else
            for (int ia = na - g + 1; ia < na; ++ia)
              at(ia) = 2.0 * wc - at(2 * bf - ia);
        } else if (fold) {  // tangential implicit: drop wall face
          if (s == 0)
            for (int ia = 0; ia < g; ++ia)
              at(ia) = 0.0;
          else
            for (int ia = na - g; ia < na; ++ia)
              at(ia) = 0.0;
        } else {  // tangential explicit: cell-centred reflection about bf-0.5
          if (s == 0)
            for (int ia = 0; ia < g; ++ia)
              at(ia) = 2.0 * wc - at(2 * bf - 1 - ia);
          else
            for (int ia = na - g; ia < na; ++ia)
              at(ia) = 2.0 * wc - at(2 * bf - 1 - ia);
        }
      });
}
 
// Collocated (cell-centered) velocity Dirichlet / no-slip ghost on one domain face. The wall sits
// at the boundary FACE (between the last inner cell and the first ghost), so EVERY component is
// reflected about it
// -- ghost = 2*wc - mirror(interior) makes the face-interpolated value equal `wc`. No
// normal/tangential split and no implicit fold: explicit reflection, re-imposed each smoother sweep
// (converges to no-slip). `wc` is the scalar `wall`, or a per-position value from `prof` (the
// resampled inlet profile, indexed (p0*prof_nc+p1)*3+comp over the face's perpendicular plane) when
// one is supplied (e.g. the BFS step inlet).
inline void bcVelocityColocated(BField f, B3 ext, int g, int a, int s, double wall, int comp = 0,
                                BField prof = BField(), int prof_nc = 0) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int na = dims[a];
  const bool hasProf = prof.extent(0) > 0;
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_vel_coloc", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        const long base = static_cast<long>(p0) * sb + static_cast<long>(p1) * sc;
        const double wc = hasProf ? prof((static_cast<long>(p0) * prof_nc + p1) * 3 + comp) : wall;
        auto at = [&](int ia) -> double& { return f(base + static_cast<long>(ia) * sa); };
        if (s == 0)
          for (int ia = 0; ia < g; ++ia)
            at(ia) = 2.0 * wc - at(2 * g - 1 - ia);  // mirror about g-1/2
        else
          for (int ia = na - g; ia < na; ++ia)
            at(ia) = 2.0 * wc - at(2 * (na - g) - 1 - ia);  // about na-g-1/2
      });
}
 
// Zero-gradient (Neumann) ghost on one domain face: every ghost layer = the boundary-adjacent inner
// cell (cell-centered). Used for the collocated pressure increment phi at walls so the
// cell-centered correction carries no spurious normal acceleration (the same role pressureBcGhost
// plays for P in the predictor).
inline void bcNeumannGhost(BField f, B3 ext, int g, int a, int s) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int na = dims[a];
  const int bic = (s == 0) ? g : (na - g - 1);
  const int lo = (s == 0) ? 0 : (na - g), hi = (s == 0) ? (g - 1) : (na - 1);
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_neumann_ghost", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        const long base = static_cast<long>(p0) * sb + static_cast<long>(p1) * sc;
        const double pin = f(base + static_cast<long>(bic) * sa);
        for (int ia = lo; ia <= hi; ++ia)
          f(base + static_cast<long>(ia) * sa) = pin;
      });
}
 
// Zero-gradient (Neumann) outflow velocity ghost for component comp on one face.
inline void bcOutflowComp(BField f, B3 ext, int g, int a, int s, int comp, int fold) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int na = dims[a];
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_outflow", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        const long base = static_cast<long>(p0) * sb + static_cast<long>(p1) * sc;
        auto at = [&](int ia) -> double& { return f(base + static_cast<long>(ia) * sa); };
        if (s == 0) {
          const int src = (comp == a) ? g + 1 : g;
          const int last = (comp == a) ? g : g - 1;
          const double v = fold ? 0.0 : at(src);
          for (int ia = 0; ia <= last; ++ia)
            at(ia) = v;
        } else {
          const int src = na - g - 1;
          const double v = fold ? 0.0 : at(src);
          for (int ia = na - g; ia < na; ++ia)
            at(ia) = v;
        }
      });
}
 
// Implicit-diffusion face-fold accumulation at the boundary-adjacent inner cell.
inline void bcDiffusionFold(BField dcorr, BField brhs, B3 ext, int g, int a, int s, double dval,
                            double bval) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int bic = (s == 0) ? g : (dims[a] - g - 1);
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_fold", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        const long i =
            static_cast<long>(p0) * sb + static_cast<long>(p1) * sc + static_cast<long>(bic) * sa;
        dcorr(i) += dval;
        brhs(i) += bval;
      });
}
 
// Hold the pressure ghost at 0 on an outflow face (Dirichlet p=0; the open face couples to it).
inline void bcZeroPressureGhost(BField phi, B3 ext, int g, int a, int s) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int na = dims[a];
  const int lo = (s == 0) ? 0 : (na - g), hi = (s == 0) ? (g - 1) : (na - 1);
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_zero_p_ghost", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        const long base = (long)p0 * sb + (long)p1 * sc;
        for (int ia = lo; ia <= hi; ++ia)
          phi(base + (long)ia * sa) = 0.0;
      });
}
 
// Projection correction of the high-side outflow normal face (index na-g) that correct_k misses:
// f -= phi[bf] - phi[bf-sa] (with the Dirichlet ghost phi[bf]=0 -> += phi_inner).
// (correct_outflow_k.)
inline void bcCorrectOutflow(BField f, BField phi, B3 ext, int g, int a) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_correct_outflow", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        const long bf = (long)p0 * sb + (long)p1 * sc + (long)(dims[a] - g) * sa;
        f(bf) -= phi(bf) - phi(bf - sa);
      });
}
 
// Set the a-component face openness on a domain face to `val` (Neumann wall/inflow -> 0; the
// periodic fill would otherwise wrap the wrong value into an outflow face from the opposite
// boundary -> set it open = 1).
inline void bcSetOpenness(BField oa, B3 ext, int g, int a, int s, double val) {
  BExec space;
  int dims[3];
  long strides[3];
  bcdetail::axisDims(ext, dims, strides);
  const int b = (a + 1) % 3, c = (a + 2) % 3;
  const long sa = strides[a], sb = strides[b], sc = strides[c];
  const int bf = (s == 0) ? g : (dims[a] - g);
  using MD = Kokkos::MDRangePolicy<BExec, Kokkos::Rank<2>>;
  Kokkos::parallel_for(
      "peclet::flow::bc_setopen", MD(space, {0, 0}, {dims[b], dims[c]}),
      KOKKOS_LAMBDA(int p0, int p1) {
        oa(static_cast<long>(p0) * sb + static_cast<long>(p1) * sc + static_cast<long>(bf) * sa) =
            val;
      });
}
inline void bcZeroOpenness(BField oa, B3 ext, int g, int a, int s) {
  bcSetOpenness(oa, ext, g, a, s, 0.0);
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_MAC_BC_HPP