interface_area_field.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_INTERFACE_AREA_FIELD_HPP
#define PECLET_FLOW_VOF_INTERFACE_AREA_FIELD_HPP
 
#include <Kokkos_Core.hpp>
#include <stdexcept>
 
#include "mac_stencils.hpp"
#include "policy.hpp"
#include "vof/advect_wy.hpp"
#include "vof/curvature.hpp"
#include "vof/curvature_field.hpp"  // vofIsInterface
#include "vof/interface_area.hpp"
 
namespace peclet::flow::vof {
 
class VofInterfaceArea {
 public:
  struct Stats {
    long interfacial = 0;
    long hf = 0;          
    long hfMixed = 0;     
    long pv = 0;          
    long noEstimate = 0;  
    double area = 0.0;    
  };
 
  void init(int nx, int ny, int nz, int ghost) {
    if (ghost < kHfColumn / 2)
      throw std::invalid_argument("peclet::flow::vof::VofInterfaceArea: ghost width must be >= 3");
    if (nx < 1 || ny < 1 || nz < 1)
      throw std::invalid_argument("peclet::flow::vof::VofInterfaceArea: empty block");
    n_ = I3{nx, ny, nz};
    g_ = ghost;
    e_ = I3{nx + 2 * ghost, ny + 2 * ghost, nz + 2 * ghost};
    len_ = static_cast<long>(e_.x) * e_.y * e_.z;
    mx_ = SField("vof::area::mx", len_);
    my_ = SField("vof::area::my", len_);
    mz_ = SField("vof::area::mz", len_);
    alpha_ = SField("vof::area::alpha", len_);
    area_ = SField("vof::area::A", len_);
    branch_ = SField("vof::area::branch", len_);
  }
 
  bool ready() const { return area_.extent(0) != 0; }
  SField area() const { return area_; }
  SField branch() const { return branch_; }
 
  // ---- tunables: the SAME knobs and defaults the curvature cascade runs with ----------------
  VofMetric metric;
  double weightWidth = kPvWeightWidth;
  double monoTol = 1e-6;
  double pureEps = 0.0;
  double cosMin = 0.2;
  double interfaceEps = 0.0;
 
  Stats compute(SField c, int mode) {
    if (!ready())
      throw std::runtime_error("peclet::flow::vof::VofInterfaceArea::compute: init() not called");
    reconstructPlanes(c);
    heightPass(c, mode);
    fallbackPass(c, mode);
    return census();
  }
 
  // public for the same nvcc reason `VofCurvature`'s passes are.
  void reconstructPlanes(SField c) {
    const I3 e = e_, n = n_;
    const int g = g_, gr = kPvHalf;
    const long sy = e_.x, sz = static_cast<long>(e_.x) * e_.y;
    SField mx = mx_, my = my_, mz = mz_, al = alpha_;
    const double ieps = interfaceEps;
    Kokkos::parallel_for(
        "vof::area::planes",
        MDRange3<SExec>(SExec(), {g - gr, g - gr, g - gr},
                        {g + n.x + gr, g + n.y + gr, g + n.z + gr}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          if (!vofIsInterface(c(i), ieps)) {
            mx(i) = 0.0;
            my(i) = 0.0;
            mz(i) = 0.0;
            al(i) = 0.0;
            return;
          }
          wyReconstructCell(c, i, sy, sz, mx, my, mz, al);
        });
    Kokkos::fence();
  }
 
  void heightPass(SField c, int mode) {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long st[3] = {1, e_.x, static_cast<long>(e_.x) * e_.y};
    SField mx = mx_, my = my_, mz = mz_, al = alpha_, ar = area_, br = branch_;
    const double mtol = monoTol, ieps = interfaceEps;
    const double peps = pureEps;
    const VofMetric gm = metric;  // `g` is the ghost width in this scope
    const int md = mode;
    Kokkos::parallel_for(
        "vof::area::hf", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          ar(i) = 0.0;
          if (!vofIsInterface(c(i), ieps)) {
            br(i) = static_cast<double>(kCurvNone);
            return;
          }
          const double am[3] = {Kokkos::fabs(mx(i)), Kokkos::fabs(my(i)), Kokkos::fabs(mz(i))};
          int ord[3] = {0, 1, 2};
          for (int p = 1; p < 3; ++p)
            for (int q = p; q > 0 && am[ord[q]] > am[ord[q - 1]]; --q) {
              const int t = ord[q];
              ord[q] = ord[q - 1];
              ord[q - 1] = t;
            }
          for (int t = 0; t < 3; ++t) {
            const int d = ord[t];
            const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
            const long sd = st[d], s1 = st[d1], s2 = st[d2];
            double hh[9];
            int orient0 = 0;
            bool ok = true;
            for (int q = 0; q < 3 && ok; ++q)
              for (int p = 0; p < 3 && ok; ++p) {
                const long base = i + (p - 1) * s1 + (q - 1) * s2;
                double col[kHfColumn];
                for (int k = 0; k < kHfColumn; ++k)
                  col[k] = c(base + (k - kHfColumn / 2) * sd);
                double h;
                int orient;
                if (!hfColumnHeight(col, kHfColumn, h, orient, mtol, peps)) {
                  ok = false;
                  break;
                }
                if (orient0 == 0)
                  orient0 = orient;
                else if (orient != orient0) {
                  ok = false;
                  break;
                }
                hh[p + 3 * q] = h;
              }
            if (!ok)
              continue;
            if (md == kAreaFootprint) {
              ar(i) = hfFootprintArea(hh, gm, d);
            } else {
              double ns[3];
              hfSurfaceNormal(hh, d, gm, ns);
              ar(i) = interfaceAreaFromNormal(md, mx(i), my(i), mz(i), al(i), c(i), ns, gm);
            }
            br(i) = static_cast<double>(t == 0 ? kCurvHf : kCurvHfMixed);
            return;
          }
          br(i) = -1.0;
        });
    Kokkos::fence();
  }
 
  void fallbackPass(SField c, int mode) {
    const I3 e = e_, n = n_;
    const int g = g_, gr = kPvHalf;
    SField mx = mx_, my = my_, mz = mz_, al = alpha_, ar = area_, br = branch_;
    const double dW = weightWidth * metric.maxH(), cmin = cosMin, ieps = interfaceEps;
    const VofMetric gm = metric;  // `g` is the ghost width in this scope
    const int md = mode;
    Kokkos::parallel_for(
        "vof::area::pv", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          if (br(i) >= 0.0)
            return;
          const double m0 = mx(i), m1 = my(i), m2 = mz(i);
          const double mi[3] = {m0, m1, m2};
          double nn[3] = {0.0, 0.0, 0.0};
          // PHYSICAL frame (Phase 3, V2.4); identity at the unit metric.
          if (!(vofPhysNormalInv(mi, gm, nn) > 0.0)) {
            ar(i) = 0.0;
            br(i) = static_cast<double>(kCurvNoEstimate);
            return;
          }
          double t1[3], t2[3];
          curvFrame(nn, t1, t2);
          double v[8][3], ctr[3], area;
          const int nv = plicPolygon(m0, m1, m2, al(i), v);
          polygonAreaCentroid(v, nv, ctr, area);
          const double org[3] = {ctr[0] - 0.5, ctr[1] - 0.5, ctr[2] - 0.5};
          PvFit fit;
          pvFitInit(fit);
          for (int oz = -gr; oz <= gr; ++oz)
            for (int oy = -gr; oy <= gr; ++oy)
              for (int ox = -gr; ox <= gr; ++ox) {
                const long j = L3(x + ox, y + oy, z + oz, e);
                if (!vofIsInterface(c(j), ieps))
                  continue;
                const double off[3] = {static_cast<double>(ox), static_cast<double>(oy),
                                       static_cast<double>(oz)};
                pvFitAdd(fit, mx(j), my(j), mz(j), al(j), off, org, t1, t2, nn, dW, cmin, gm);
              }
          double a[6];
          bool red = false;
          if (!pvFitSolve(fit, a, red)) {
            // no cascade geometry at all: keep the MYC PLIC area rather than emit a zero. Loud in
            // the census, exactly as `kCurvNoEstimate` is in the curvature.
            ar(i) = plicAreaMetric(m0, m1, m2, al(i), gm);
            br(i) = static_cast<double>(kCurvNoEstimate);
            return;
          }
          double ns[3];
          pvSurfaceNormal(a, t1, t2, nn, ns);
          // The PV branch has no column and therefore no tiling footprint of its own; under
          // `kAreaFootprint` it takes variant B (the plane rebuilt on the paraboloid's normal),
          // which is the closest thing to it that a single cell can supply. The branch census
          // says how much of the total that is.
          ar(i) = interfaceAreaFromNormal(md == kAreaFootprint ? kAreaNormal : md, m0, m1, m2,
                                          al(i), c(i), ns, gm);
          br(i) = static_cast<double>(red ? kCurvPvReduced : kCurvPv);
        });
    Kokkos::fence();
  }
 
  Stats census() const {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField br = branch_, ar = area_;
    Stats s;
    Kokkos::parallel_reduce(
        "vof::area::census", MDRange3<SExec>(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z}),
        KOKKOS_LAMBDA(int x, int y, int z, long& ni, long& n1, long& n2, long& n3, long& n4,
                      double& acc) {
          const long i = L3(x, y, z, e);
          const int b = static_cast<int>(br(i));
          if (b == kCurvNone)
            return;
          ++ni;
          acc += ar(i);
          if (b == kCurvHf)
            ++n1;
          else if (b == kCurvHfMixed)
            ++n2;
          else if (b == kCurvPv || b == kCurvPvReduced)
            ++n3;
          else
            ++n4;
        },
        s.interfacial, s.hf, s.hfMixed, s.pv, s.noEstimate, s.area);
    Kokkos::fence();
    return s;
  }
 
 private:
  I3 n_{0, 0, 0}, e_{0, 0, 0};
  int g_ = 0;
  long len_ = 0;
  SField mx_, my_, mz_, alpha_, area_, branch_;
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_INTERFACE_AREA_FIELD_HPP