gauge_exact_gradient.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_GAUGE_EXACT_GRADIENT_HPP
#define PECLET_FLOW_GAUGE_EXACT_GRADIENT_HPP
 
#include <Kokkos_Core.hpp>
 
#include "mac_cutcell.hpp"  // CCField/CCConst, C3, CCExec
 
namespace peclet::flow {
 
inline void gpCenterGrad(CCField out, CCConst p, CCConst sdf, int axis, C3 e, int g,
                         bool grad2a = false) {
  CCExec space;
  using MD = Kokkos::MDRangePolicy<CCExec, Kokkos::Rank<3>>;
  Kokkos::parallel_for(
      "peclet::flow::gp_center_grad", MD(space, {g, g, g}, {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 long sa = (axis == 0) ? 1 : (axis == 1) ? sy : sz;
        if (sdf(i) < 0.0) {
          out(i) = 0.0;
          return;
        }
        const bool am = sdf(i - sa) >= 0.0, ap = sdf(i + sa) >= 0.0;
        // the +/-3 cell of the grad2a stencil can sit one past the g-wide ghost pad; the axis
        // coordinate guard falls back to the 2-point form there (never reads out of the block)
        const int ca = (axis == 0) ? x : (axis == 1) ? y : z;
        const int ea = (axis == 0) ? e.x : (axis == 1) ? e.y : e.z;
        double gr;
        if (am && ap)
          gr = 0.5 * (p(i + sa) - p(i - sa));
        else if (ap) {
          if (grad2a)
            gr = (ca + 3 < ea && sdf(i + 2 * sa) >= 0.0 && sdf(i + 3 * sa) >= 0.0)
                     ? 0.5 * (-2.0 * p(i) + p(i + sa) + 2.0 * p(i + 2 * sa) - p(i + 3 * sa))
                     : (p(i + sa) - p(i));
          else
            gr = (sdf(i + 2 * sa) >= 0.0) ? 0.5 * (-3.0 * p(i) + 4.0 * p(i + sa) - p(i + 2 * sa))
                                          : (p(i + sa) - p(i));
        } else if (am) {
          if (grad2a)
            gr = (ca - 3 >= 0 && sdf(i - 2 * sa) >= 0.0 && sdf(i - 3 * sa) >= 0.0)
                     ? 0.5 * (2.0 * p(i) - p(i - sa) - 2.0 * p(i - 2 * sa) + p(i - 3 * sa))
                     : (p(i) - p(i - sa));
          else
            gr = (sdf(i - 2 * sa) >= 0.0) ? 0.5 * (3.0 * p(i) - 4.0 * p(i - sa) + p(i - 2 * sa))
                                          : (p(i) - p(i - sa));
        } else
          gr = 0.0;
        out(i) = gr;
      });
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_GAUGE_EXACT_GRADIENT_HPP