wetting_dynamic.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_WETTING_DYNAMIC_HPP
#define PECLET_FLOW_VOF_WETTING_DYNAMIC_HPP
 
#include <cmath>
#include <functional>
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
#include "policy.hpp"
#include "vof/advect_wy.hpp"
 
namespace peclet::flow::vof {
 
enum VofDynamicState : int {
  kVofDynNone = 0,       
  kVofDynStatic = 1,     
  kVofDynPinned = 2,     
  kVofDynAdvancing = 3,  
  kVofDynReceding = 4,   
  kVofDynStateCount = 5
};
 
KOKKOS_INLINE_FUNCTION double coxVoinovAngle(double thetaE, double ca, double logRatio,
                                             double thetaMin, double thetaMax) {
  const double t3 = thetaE * thetaE * thetaE + 9.0 * ca * logRatio;
  double th = t3 > 0.0 ? Kokkos::pow(t3, 1.0 / 3.0) : thetaMin;
  if (!(th == th))  // NaN guard: a caller-supplied NaN velocity must not reach the fill
    th = thetaE;
  return th < thetaMin ? thetaMin : (th > thetaMax ? thetaMax : th);
}
 
KOKKOS_INLINE_FUNCTION int vofHysteresisBase(double thetaApp, double thetaA, double thetaR,
                                             double& base) {
  if (thetaApp > thetaA) {
    base = thetaA;
    return kVofDynAdvancing;
  }
  if (thetaApp < thetaR) {
    base = thetaR;
    return kVofDynReceding;
  }
  base = thetaApp;
  return kVofDynPinned;
}
 
KOKKOS_INLINE_FUNCTION int vofDynamicContactAngle(double thetaApp, double uCl, double muL,
                                                  double sigma, double logRatio, bool hyst,
                                                  double thetaA, double thetaR, double thetaE,
                                                  double thetaMin, double thetaMax,
                                                  double& thetaOut, double& caOut) {
  caOut = sigma > 0.0 ? muL * uCl / sigma : 0.0;
  double base = thetaE;
  int state = kVofDynStatic;
  if (hyst) {
    state = vofHysteresisBase(thetaApp, thetaA, thetaR, base);
    if (state == kVofDynPinned) {
      // The fill is idempotent, so imposing the apparent angle reproduces the current interface
      // and the line does not move. No Cox-Voinov correction: the line is not moving.
      thetaOut = base < thetaMin ? thetaMin : (base > thetaMax ? thetaMax : base);
      return state;
    }
  }
  thetaOut = coxVoinovAngle(base, caOut, logRatio, thetaMin, thetaMax);
  return state;
}
 
KOKKOS_INLINE_FUNCTION bool vofWallTangent(const double mfIn[3], const double nwIn[3], double tEps,
                                           double that[3], double& cosApp, const VofMetric& g) {
  // Phase 3 (V4.3): the apparent angle and the in-wall direction are PHYSICAL, and both inputs
  // arrive as index-space gradients — `vofPhysNormal` is the pullback, exactly as in
  // `vofWettingPlane`. At the unit metric it divides by 1.0 and returns `m/|m|`, so the two
  // normalisations below are the ones this function always did, bit for bit. `that` comes back as
  // a PHYSICAL unit direction; the caller that steps on the index lattice pulls it back itself.
  double mh[3] = {mfIn[0], mfIn[1], mfIn[2]};
  if (!(vofPhysNormal(mfIn, g, mh) > 0.0)) {
    mh[0] = mfIn[0];
    mh[1] = mfIn[1];
    mh[2] = mfIn[2];
  }
  double nw[3] = {nwIn[0], nwIn[1], nwIn[2]};
  if (!(vofPhysNormal(nwIn, g, nw) > 0.0)) {
    nw[0] = nwIn[0];
    nw[1] = nwIn[1];
    nw[2] = nwIn[2];
  }
  cosApp = mh[0] * nw[0] + mh[1] * nw[1] + mh[2] * nw[2];
  double t[3] = {mh[0] - cosApp * nw[0], mh[1] - cosApp * nw[1], mh[2] - cosApp * nw[2]};
  const double tn = Kokkos::sqrt(t[0] * t[0] + t[1] * t[1] + t[2] * t[2]);
  if (tn < tEps) {
    that[0] = that[1] = that[2] = 0.0;
    return false;
  }
  for (int d = 0; d < 3; ++d)
    that[d] = t[d] / tn;
  return true;
}
 
KOKKOS_INLINE_FUNCTION bool vofWallTangent(const double mfIn[3], const double nwIn[3], double tEps,
                                           double that[3], double& cosApp) {
  return vofWallTangent(mfIn, nwIn, tEps, that, cosApp, VofMetric{});
}
 
KOKKOS_INLINE_FUNCTION double vofContactLineSpeed(const double u[3], const double that[3]) {
  return u[0] * that[0] + u[1] * that[1] + u[2] * that[2];
}
 
struct VofDynamicWetting {
  // ---- configuration (all angles in RADIANS) -------------------------------------------------
  bool dynamic = false;     
  bool hysteresis = false;  
  double thetaA = 0.0, thetaR = 0.0;
  VofMetric metric;
  double slip = 0.1;      
  double muLiquid = 0.0;  
  double sigma = 0.0;     
  double thetaMin = 1.0 * 3.14159265358979323846 / 180.0;
  double thetaMax = 179.0 * 3.14159265358979323846 / 180.0;
  bool smooth = true;     
  double pureEps = 1e-8;  
  double tanEps = 1e-6;   
 
  bool active() const { return dynamic || hysteresis; }
  double logRatio() const { return dynamic ? -std::log(slip) : 0.0; }
 
  // ---- state ---------------------------------------------------------------------------------
  bool allocated() const { return base_.extent(0) > 0; }
  SField base() const { return base_; }
  SField imposed() const { return imposed_; }     
  SField apparent() const { return app_; }        
  SField speed() const { return uclS_; }          
  SField capillary() const { return ca_; }        
  SField stateField() const { return stateD_; }   
  SField cellVel(int d) const { return uc_[d]; }  
 
  void allocate(long len) {
    if (base_.extent(0) == (std::size_t)len)
      return;
    base_ = SField("vof::dyn::base", len);
    imposed_ = SField("vof::dyn::imposed", len);
    app_ = SField("vof::dyn::apparent", len);
    ucl_ = SField("vof::dyn::uclRaw", len);
    uclS_ = SField("vof::dyn::ucl", len);
    wgt_ = SField("vof::dyn::valid", len);
    ca_ = SField("vof::dyn::ca", len);
    stateD_ = SField("vof::dyn::state", len);
    for (int d = 0; d < 3; ++d) {
      uc_[d] = SField("vof::dyn::cellVel", len);
      th_[d] = SField("vof::dyn::tangent", len);
    }
  }
 
  void measure(const WyAdvector& adv) {
    const I3 e = adv.extent(), n = adv.inner();
    const int g = adv.ghost();
    const long sx = 1, sy = e.x, sz = static_cast<long>(e.x) * e.y;
    SField c = adv.colour(), sdf = adv.wallSdf();
    SField mx = adv.wettingNormal(0), my = adv.wettingNormal(1), mz = adv.wettingNormal(2);
    SField ux = uc_[0], uy = uc_[1], uz = uc_[2];
    SField ucl = ucl_, wg = wgt_, ap = app_, tx = th_[0], ty = th_[1], tz = th_[2];
    UCField kk = adv.cellKind();
    const double pe = pureEps, te = tanEps;
    const VofMetric gme = metric;  // `g` is the ghost width in this scope
    Kokkos::parallel_for(
        "vof::dyn::measure", MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          ucl(i) = 0.0;
          wg(i) = 0.0;
          ap(i) = 0.0;
          tx(i) = ty(i) = tz(i) = 0.0;
          if (vofGhostDepth(x, y, z, g, n.x, n.y, n.z) > 2)
            return;
          if (kk(i) != kVofSolid)
            return;
          double nw[3] = {0.5 * (sdf(i + sx) - sdf(i - sx)), 0.5 * (sdf(i + sy) - sdf(i - sy)),
                          0.5 * (sdf(i + sz) - sdf(i - sz))};
          // Phase 3 (V4.2): the WALK follows the physical wall normal expressed in index steps,
          // `H^-1 n_w ~ H^-2 grad_xi(sdf)`, with ONE normalisation — identical to `nw/|nw|` at the
          // unit metric. `nw` itself stays the raw index gradient: `vofWallTangent` pulls it back.
          double wk[3] = {nw[0] / (gme.h[0] * gme.h[0]), nw[1] / (gme.h[1] * gme.h[1]),
                          nw[2] / (gme.h[2] * gme.h[2])};
          const double nn = Kokkos::sqrt(wk[0] * wk[0] + wk[1] * wk[1] + wk[2] * wk[2]);
          if (!(nn > 1e-12))
            return;
          for (int d = 0; d < 3; ++d)
            wk[d] /= nn;
          for (int step = 1; step <= 4; ++step) {
            const int fx = x + static_cast<int>(Kokkos::round(step * wk[0]));
            const int fy = y + static_cast<int>(Kokkos::round(step * wk[1]));
            const int fz = z + static_cast<int>(Kokkos::round(step * wk[2]));
            if (fx < 0 || fy < 0 || fz < 0 || fx >= e.x || fy >= e.y || fz >= e.z)
              return;
            const long fi = L3(fx, fy, fz, e);
            if (kk(fi) != kVofFluid)
              continue;
            double mAcc[3] = {0.0, 0.0, 0.0}, uAcc[3] = {0.0, 0.0, 0.0};
            int cnt = 0;
            const double cf = c(fi);
            const double ml = Kokkos::fabs(mx(fi)) + Kokkos::fabs(my(fi)) + Kokkos::fabs(mz(fi));
            if (cf > pe && cf < 1.0 - pe && ml > 0.0) {
              mAcc[0] = mx(fi);
              mAcc[1] = my(fi);
              mAcc[2] = mz(fi);
              uAcc[0] = ux(fi);
              uAcc[1] = uy(fi);
              uAcc[2] = uz(fi);
              cnt = 1;
            } else {
              // The anchor column carries no interface: average over its MIXED fluid neighbours,
              // the same set `solidBandFillPassWetting`'s neighbour branch averages over.
              const bool inner =
                  fx >= 1 && fy >= 1 && fz >= 1 && fx + 1 < e.x && fy + 1 < e.y && fz + 1 < e.z;
              if (inner)
                for (int kz = -1; kz <= 1; ++kz)
                  for (int ky = -1; ky <= 1; ++ky)
                    for (int kx = -1; kx <= 1; ++kx) {
                      const long gi = L3(fx + kx, fy + ky, fz + kz, e);
                      if (kk(gi) != kVofFluid)
                        continue;
                      const double cg = c(gi);
                      if (cg <= pe || cg >= 1.0 - pe)
                        continue;
                      const double gm[3] = {mx(gi), my(gi), mz(gi)};
                      const double gl =
                          Kokkos::fabs(gm[0]) + Kokkos::fabs(gm[1]) + Kokkos::fabs(gm[2]);
                      if (!(gl > 0.0))
                        continue;
                      const double g2 = Kokkos::sqrt(gm[0] * gm[0] + gm[1] * gm[1] + gm[2] * gm[2]);
                      for (int d = 0; d < 3; ++d)
                        mAcc[d] += gm[d] / g2;
                      uAcc[0] += ux(gi);
                      uAcc[1] += uy(gi);
                      uAcc[2] += uz(gi);
                      ++cnt;
                    }
            }
            if (cnt == 0)
              return;
            const double inv = 1.0 / static_cast<double>(cnt);
            // `u_a = h_a v_a` (the index -> physical velocity map, Phase 1 `velToPhys(a)` without
            // the tRef that cancels in `U_cl/|u|`): the contact-line speed is a PHYSICAL projection
            // on a PHYSICAL direction. Identity at the unit metric.
            const double uu[3] = {uAcc[0] * inv * gme.h[0], uAcc[1] * inv * gme.h[1],
                                  uAcc[2] * inv * gme.h[2]};
            double that[3], cosApp;
            const bool ok = vofWallTangent(mAcc, nw, te, that, cosApp, gme);
            ap(i) = Kokkos::acos(cosApp < -1.0 ? -1.0 : (cosApp > 1.0 ? 1.0 : cosApp));
            wg(i) = 1.0;
            if (ok) {
              // The stored tangent is what `impose()` STEPS along on the index lattice, so it is
              // the index pullback `H^-1 t_hat`, renormalized once.
              //
              // RULE B(3) GUARD (`flow/doc/anisotropic_vof.md` §1) — the ONLY one on the VoF half.
              // `that` is already an L2-unit vector, but its norm is 1.0 only to round-off, so a
              // second normalisation would perturb the stored bits on the isotropic path. The
              // pullback is the identity there, so it is skipped outright.
              double ti[3] = {that[0], that[1], that[2]};
              if (!gme.isotropic()) {
                for (int d = 0; d < 3; ++d)
                  ti[d] = that[d] / gme.h[d];
                const double tin = Kokkos::sqrt(ti[0] * ti[0] + ti[1] * ti[1] + ti[2] * ti[2]);
                if (tin > 0.0)
                  for (int d = 0; d < 3; ++d)
                    ti[d] /= tin;
              }
              tx(i) = ti[0];
              ty(i) = ti[1];
              tz(i) = ti[2];
              ucl(i) = vofContactLineSpeed(uu, that);
            }
            return;
          }
        });
    Kokkos::fence();
  }
 
  void impose(const WyAdvector& adv) {
    const I3 e = adv.extent(), n = adv.inner();
    const int g = adv.ghost();
    SField th = adv.contactAngle();
    SField bs = base_, ap = app_, ucl = ucl_, wg = wgt_, us = uclS_, ca = ca_, st = stateD_;
    SField tx = th_[0], ty = th_[1], tz = th_[2], imp = imposed_;
    UCField kk = adv.cellKind();
    const bool hyst = hysteresis, sm = smooth;
    const double ta = thetaA, tr = thetaR, mu = muLiquid, sg = sigma, lr = logRatio();
    const double tmin = thetaMin, tmax = thetaMax;
    Kokkos::parallel_for(
        "vof::dyn::impose", MDRange3<SExec>(SExec(), {0, 0, 0}, {e.x, e.y, e.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(x, y, z, e);
          us(i) = 0.0;
          ca(i) = 0.0;
          st(i) = static_cast<double>(kVofDynNone);
          imp(i) = bs(i);
          if (vofGhostDepth(x, y, z, g, n.x, n.y, n.z) > 2)
            return;
          if (kk(i) != kVofSolid) {
            th(i) = bs(i);
            return;
          }
          if (!(wg(i) > 0.0)) {
            th(i) = bs(i);  // no contact-line data: the static base stands
            return;
          }
          double su = ucl(i), sw = 1.0;
          const double td[3] = {tx(i), ty(i), tz(i)};
          if (sm && (td[0] != 0.0 || td[1] != 0.0 || td[2] != 0.0)) {
            for (int s = -1; s <= 1; s += 2) {
              const int dx = static_cast<int>(Kokkos::round(s * td[0]));
              const int dy = static_cast<int>(Kokkos::round(s * td[1]));
              const int dz = static_cast<int>(Kokkos::round(s * td[2]));
              if (dx == 0 && dy == 0 && dz == 0)
                continue;
              const int px = x + dx, py = y + dy, pz = z + dz;
              if (px < 0 || py < 0 || pz < 0 || px >= e.x || py >= e.y || pz >= e.z)
                continue;
              const long j = L3(px, py, pz, e);
              if (!(wg(j) > 0.0))
                continue;
              su += ucl(j);
              sw += 1.0;
            }
          }
          const double u = su / sw;
          double thOut, caOut;
          const int state = vofDynamicContactAngle(ap(i), u, mu, sg, lr, hyst, ta, tr, bs(i), tmin,
                                                   tmax, thOut, caOut);
          us(i) = u;
          ca(i) = caOut;
          st(i) = static_cast<double>(state);
          imp(i) = thOut;
          th(i) = thOut;
        });
    Kokkos::fence();
  }
 
  SField uclRaw() const { return ucl_; }
  SField valid() const { return wgt_; }
 
  struct Census {
    long cells[kVofDynStateCount] = {0, 0, 0, 0, 0};
    double meanImposedDeg = 0.0, meanApparentDeg = 0.0, maxCa = 0.0, maxSpeed = 0.0;
    long contactCells = 0;
  };
  Census census(const WyAdvector& adv) const {
    Census r;
    if (!allocated())
      return r;
    const I3 e = adv.extent(), n = adv.inner();
    const int g = adv.ghost();
    SField wg = wgt_, imp = imposed_, ap = app_, ca = ca_, st = stateD_, us = uclS_;
    UCField kk = adv.cellKind();
    using MD = MDRange3<SExec>;
    MD pol(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z});
    for (int s = 0; s < kVofDynStateCount; ++s) {
      long acc = 0;
      const double sv = static_cast<double>(s);
      Kokkos::parallel_reduce(
          "vof::dyn::census_state", pol,
          KOKKOS_LAMBDA(int x, int y, int z, long& a) {
            const long i = L3(x, y, z, e);
            if (kk(i) == kVofSolid && wg(i) > 0.0 && st(i) == sv)
              ++a;
          },
          acc);
      r.cells[s] = acc;
    }
    double si = 0.0, sa = 0.0, mc = 0.0, ms = 0.0;
    long cnt = 0;
    Kokkos::parallel_reduce(
        "vof::dyn::census_mean", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& ai, double& aa, long& ac) {
          const long i = L3(x, y, z, e);
          if (kk(i) == kVofSolid && wg(i) > 0.0) {
            ai += imp(i);
            aa += ap(i);
            ++ac;
          }
        },
        si, sa, cnt);
    Kokkos::parallel_reduce(
        "vof::dyn::census_maxca", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& am) {
          const long i = L3(x, y, z, e);
          if (kk(i) == kVofSolid && wg(i) > 0.0) {
            const double v = Kokkos::fabs(ca(i));
            if (v > am)
              am = v;
          }
        },
        Kokkos::Max<double>(mc));
    Kokkos::parallel_reduce(
        "vof::dyn::census_maxu", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& as) {
          const long i = L3(x, y, z, e);
          if (kk(i) == kVofSolid && wg(i) > 0.0) {
            const double w = Kokkos::fabs(us(i));
            if (w > as)
              as = w;
          }
        },
        Kokkos::Max<double>(ms));
    Kokkos::fence();
    const double toDeg = 180.0 / 3.14159265358979323846;
    r.contactCells = cnt;
    r.meanImposedDeg = cnt ? si / static_cast<double>(cnt) * toDeg : 0.0;
    r.meanApparentDeg = cnt ? sa / static_cast<double>(cnt) * toDeg : 0.0;
    r.maxCa = mc;
    r.maxSpeed = ms;
    return r;
  }
 
 private:
  SField base_, imposed_, app_, ucl_, uclS_, wgt_, ca_, stateD_;
  SField uc_[3], th_[3];
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_WETTING_DYNAMIC_HPP