probe_flux.hpp

Go to the documentation of this file.

 
#ifndef PECLET_CORE_SCHEME_PROBE_FLUX_HPP
#define PECLET_CORE_SCHEME_PROBE_FLUX_HPP
 
#include <cmath>
 
#include "peclet/core/scheme/cut_cell_closure.hpp"  // PECLET_CORE_CC_HD
 
namespace peclet::core::scheme {
 
inline constexpr int kProbeR0 = 0;
inline constexpr int kProbeR1a = 1;
inline constexpr int kProbeR1b = 2;
inline constexpr int kProbeR2 = 3;
 
inline constexpr double kProbeSigma = 1.1;
inline constexpr double kProbeSigmaLong = 1.5;
inline constexpr double kProbeSigmaLast = 0.5;
inline constexpr double kProbeWeightEps = 1e-12;
inline constexpr double kProbeClearance = 0.5;
inline constexpr int kProbeReach = 2;
inline constexpr double kProbeGapMin = 0.2;
inline constexpr double kProbeMinKeptWeight = 0.5;
 
struct ProbeResult {
  int rung;
  double s;
  int n;
  int off[8][3];
  double w[8];
};
 
PECLET_CORE_CC_HD double probeSupport(const double n[3], const double h[3]) {
  const double ax = n[0] < 0.0 ? -n[0] : n[0];
  const double ay = n[1] < 0.0 ? -n[1] : n[1];
  const double az = n[2] < 0.0 ? -n[2] : n[2];
  return 0.5 * (ax * h[0] + ay * h[1] + az * h[2]);
}
 
namespace detail {
 
PECLET_CORE_CC_HD double probeFloor(double x) {
#ifdef KOKKOS_INLINE_FUNCTION
  return Kokkos::floor(x);
#else
  return std::floor(x);
#endif
}
 
}  // namespace detail
 
PECLET_CORE_CC_HD void trilinearStencil(const double xi[3], int base[3], double w[8]) {
  double f[3];
  for (int a = 0; a < 3; ++a) {
    const double b = detail::probeFloor(xi[a]);
    base[a] = static_cast<int>(b);
    f[a] = xi[a] - b;
  }
  for (int k = 0; k < 8; ++k) {
    const double wx = (k & 1) ? f[0] : 1.0 - f[0];
    const double wy = ((k >> 1) & 1) ? f[1] : 1.0 - f[1];
    const double wz = ((k >> 2) & 1) ? f[2] : 1.0 - f[2];
    w[k] = wx * wy * wz;
  }
}
 
namespace detail {
 
// The candidate probe at distance s: index-space point xi, its trilinear stencil, and V3.
PECLET_CORE_CC_HD bool probeCandidate(const double xw[3], const double n[3], const double h[3],
                                      double s, double xi[3], int base[3], double w[8]) {
  for (int a = 0; a < 3; ++a)
    xi[a] = (xw[a] + s * n[a]) / h[a];
  trilinearStencil(xi, base, w);
  for (int a = 0; a < 3; ++a)
    if (base[a] < -kProbeReach || base[a] + 1 > kProbeReach)
      return false;  // V3
  return true;
}
 
// V1: every cell with weight > kProbeWeightEps is an unknown of the phase.
template <class Lookup>
PECLET_CORE_CC_HD bool probeAllUnknown(const int base[3], const double w[8], const Lookup& L) {
  for (int k = 0; k < 8; ++k)
    if (w[k] > kProbeWeightEps &&
        !L.unknown(base[0] + (k & 1), base[1] + ((k >> 1) & 1), base[2] + ((k >> 2) & 1)))
      return false;
  return true;
}
 
PECLET_CORE_CC_HD void probeFillAll(ProbeResult& r, int rung, double s, const int base[3],
                                    const double w[8]) {
  r.rung = rung;
  r.s = s;
  r.n = 8;
  for (int k = 0; k < 8; ++k) {
    r.off[k][0] = base[0] + (k & 1);
    r.off[k][1] = base[1] + ((k >> 1) & 1);
    r.off[k][2] = base[2] + ((k >> 2) & 1);
    r.w[k] = w[k];
  }
}
 
}  // namespace detail
 
template <class Lookup>
PECLET_CORE_CC_HD ProbeResult buildProbe(const double xw[3], const double n[3], const double h[3],
                                         const Lookup& L) {
  ProbeResult r;
  const double S = probeSupport(n, h);
  double xi[3], w[8];
  int base[3];
 
  // R0: s0 = 1.1 S, trilinear, V1 ∧ V2 ∧ V3.
  const double s0 = kProbeSigma * S;
  const bool reach0 = detail::probeCandidate(xw, n, h, s0, xi, base, w);
  double phi0 = 0.0;
  if (reach0) {
    phi0 = L.phi(xi);
    const bool v1 = detail::probeAllUnknown(base, w, L);
    const bool v2 = phi0 >= kProbeClearance * s0;
    if (v1 && v2) {
      detail::probeFillAll(r, kProbeR0, s0, base, w);
      return r;
    }
    // R1a: R0 failed V1 only -> a longer probe, s = 1.5 S.
    if (v2) {
      const double s1 = kProbeSigmaLong * S;
      if (detail::probeCandidate(xw, n, h, s1, xi, base, w) &&
          detail::probeAllUnknown(base, w, L) && L.phi(xi) >= kProbeClearance * s1) {
        detail::probeFillAll(r, kProbeR1a, s1, base, w);
        return r;
      }
    }
    // R1b: probe the estimated gap midpoint, s = ½ w with w = s0 + φ(p0), dropping the cells
    // that are not unknowns and renormalizing.
    const double gap = s0 + phi0;
    if (gap > kProbeGapMin * S) {
      const double s1 = 0.5 * gap;
      if (detail::probeCandidate(xw, n, h, s1, xi, base, w)) {
        int m = 0;
        double sum = 0.0;
        for (int k = 0; k < 8; ++k) {
          const int o0 = base[0] + (k & 1), o1 = base[1] + ((k >> 1) & 1), o2 = base[2] + (k >> 2);
          if (!L.unknown(o0, o1, o2))
            continue;
          r.off[m][0] = o0;
          r.off[m][1] = o1;
          r.off[m][2] = o2;
          r.w[m] = w[k];
          sum += w[k];
          ++m;
        }
        if (sum >= kProbeMinKeptWeight) {
          for (int k = 0; k < m; ++k)
            r.w[k] /= sum;
          r.rung = kProbeR1b;
          r.s = s1;
          r.n = m;
          return r;
        }
      }
    }
  }
 
  // R2 (last resort, first order): the cell itself at s = 0.5 S.
  r.rung = kProbeR2;
  r.s = kProbeSigmaLast * S;
  r.n = 1;
  r.off[0][0] = r.off[0][1] = r.off[0][2] = 0;
  r.w[0] = 1.0;
  return r;
}
 
PECLET_CORE_CC_HD double wallConductance(double s, double L, double k) {
  constexpr double kDoubleMax = 1.7976931348623157e308;
  if (!(k > 0.0))
    return 0.0;
  if (k > kDoubleMax)
    return L / s;
  return 1.0 / (s / L + 1.0 / k);
}
 
PECLET_CORE_CC_HD double interfaceConductance(double sf, double Lf, double Rc, double ss,
                                              double Ls) {
  return 1.0 / (sf / Lf + Rc + ss / Ls);
}
 
}  // namespace peclet::core::scheme
 
#endif  // PECLET_CORE_SCHEME_PROBE_FLUX_HPP