wetting.hpp

Go to the documentation of this file.

 
#ifndef PECLET_CORE_VOF_WETTING_HPP
#define PECLET_CORE_VOF_WETTING_HPP
 
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
#include "peclet/core/vof/curvature.hpp"  // plicPolygon + polygonAreaCentroid (the PLIC centroid)
#include "peclet/core/vof/plic.hpp"
 
namespace peclet::core::vof {
 
enum VofWettingPivot : int {
  kVofPivotVolume = 0,      
  kVofPivotInterface = 1,   
  kVofPivotWallNormal = 2,  
  kVofPivotContactLine = 3  
};
 
enum VofWettingBranch : int {
  kVofWetNone = 0,       
  kVofWetTheta = 1,      
  kVofWetNeighbour = 2,  
  kVofWetPure = 3,       
  kVofWetParallel = 4,   
  kVofWetNeutral = 5,    
  kVofWetCount = 6
};
 
KOKKOS_INLINE_FUNCTION bool youngsNormalFluidOnly(const double c[27], const unsigned char fluid[27],
                                                  double m[3]) {
  for (int d = 0; d < 3; ++d) {
    double s[3] = {0.0, 0.0, 0.0}, w[3] = {0.0, 0.0, 0.0};  // planes at offset -1, 0, +1 in d
    for (int q = 0; q < 27; ++q) {
      if (!fluid[q])
        continue;
      const int o[3] = {q % 3 - 1, (q / 3) % 3 - 1, q / 9 - 1};
      const int t1 = o[(d + 1) % 3], t2 = o[(d + 2) % 3];
      const double ww = (2.0 - (t1 < 0 ? -t1 : t1)) * (2.0 - (t2 < 0 ? -t2 : t2));
      const int p = o[d] + 1;
      s[p] += ww * c[q];
      w[p] += ww;
    }
    if (w[0] > 0.0 && w[2] > 0.0)
      m[d] = s[0] / w[0] - s[2] / w[2];
    else if (w[0] > 0.0 && w[1] > 0.0)
      m[d] = 2.0 * (s[0] / w[0] - s[1] / w[1]);
    else if (w[2] > 0.0 && w[1] > 0.0)
      m[d] = 2.0 * (s[1] / w[1] - s[2] / w[2]);
    else
      m[d] = 0.0;
  }
  return plicNormalizeL1(m) >= 1e-30;
}
 
KOKKOS_INLINE_FUNCTION void vofPlicCentroid(const double m[3], double alpha, double p[3]) {
  double v[8][3];
  const int nv = plicPolygon(m[0], m[1], m[2], alpha, v);
  if (nv >= 3) {
    double area = 0.0;
    polygonAreaCentroid(v, nv, p, area);
    if (area > 0.0)
      return;
  }
  const double nn = m[0] * m[0] + m[1] * m[1] + m[2] * m[2];
  const double t = nn > 0.0 ? (alpha - 0.5 * (m[0] + m[1] + m[2])) / nn : 0.0;
  p[0] = 0.5 + t * m[0];
  p[1] = 0.5 + t * m[1];
  p[2] = 0.5 + t * m[2];
}
 
KOKKOS_INLINE_FUNCTION int vofWettingPlane(const double mfIn[3], double cf, const double nwIn[3],
                                           double cosT, double sinT, double sdfF, int pivot,
                                           double tEps, double mth[3], double& alphaTh,
                                           double& cosApp, const VofMetric& g) {
  // (a) L2-unit frames — PHYSICAL (Phase 3, `flow/doc/anisotropic_vof.md` §6 decision V4.1).
  //
  // theta is an angle between two directions, so the rotation below has to happen in physical
  // space: both inputs arrive as INDEX-space gradients (a colour Youngs gradient and a central
  // difference of the SDF on the index lattice), and an index normal `m` is the physical direction
  // `H^-1 m` renormalized. `vofPhysNormal` is that pullback, and at `g.h = {1,1,1}` it divides by
  // 1.0 and returns `m/|m|` — the same two expressions this function had, bit for bit.
  //
  // The plane is handed back to `plicAlpha`/`plicVolume` in INDEX space (`m_idx = H m_theta`),
  // because those routines live on the unit cube; they renormalize internally, so no L1 pass is
  // needed and the isotropic path multiplies by 1.0.
  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];
  }
 
  // (b) rotate to the prescribed angle about the contact line's direction: keep the AZIMUTH of the
  //     interface (its in-wall component, which a fluid-only stencil measures well) and replace the
  //     angle to the wall by the prescribed one.
  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]);
  int branch = kVofWetTheta;
  if (tn < tEps) {
    mth[0] = mh[0];
    mth[1] = mh[1];
    mth[2] = mh[2];
    branch = kVofWetParallel;
  } else {
    for (int d = 0; d < 3; ++d)
      mth[d] = cosT * nw[d] + sinT * t[d] / tn;
  }
 
  // (c) back to INDEX space (m_idx = H m_theta), then anchor it on the unit cube.
  vofIndexNormal(mth, g, mth);
  if (pivot == kVofPivotVolume) {
    // The plane of normal m_theta whose liquid volume in the anchor cell is exactly C_f. No pivot
    // point is needed and the poorly-measured wall-normal part of m_f is never read.
    alphaTh = plicAlpha(mth[0], mth[1], mth[2], cf);
    return branch;
  }
  // The three ablations all reconstruct the FLUID cell with m_f's own normal first. `mfIn` is
  // already the index normal the unit-cube reconstruction wants.
  double mfl[3] = {mfIn[0], mfIn[1], mfIn[2]};
  plicNormalizeL1(mfl);
  double pf[3];
  vofPlicCentroid(mfl, plicAlpha(mfl[0], mfl[1], mfl[2], cf), pf);
  double c0[3] = {pf[0], pf[1], pf[2]};
  if (pivot == kVofPivotWallNormal || pivot == kVofPivotContactLine) {
    // |grad sdf| = 1, so the signed distance at p_f is the cell-centre value plus the normal
    // displacement. (Cell centre = (0.5,0.5,0.5) in the local frame; lengths are cells.)
    // `pf` is an index position on the unit cube; the wall normal is physical, so the
    // displacement is mapped before the dot product (identity at g.h = {1,1,1}).
    const double dpi[3] = {pf[0] - 0.5, pf[1] - 0.5, pf[2] - 0.5};
    double dpp[3];
    g.toPhys(dpi, dpp);
    const double sp = sdfF + nw[0] * dpp[0] + nw[1] * dpp[1] + nw[2] * dpp[2];
    // A physical displacement `s * n` is the index displacement `H^-1 (s n)` (identity at
    // g.h = {1,1,1}: `x / 1.0 == x`).
    if (pivot == kVofPivotWallNormal) {
      for (int d = 0; d < 3; ++d)
        c0[d] = pf[d] - sp * nw[d] / g.h[d];
    } else if (tn >= tEps) {
      // Move within the FLUID plane along the in-plane part of n_w until the wall is reached.
      double q[3] = {nw[0] - cosApp * mh[0], nw[1] - cosApp * mh[1], nw[2] - cosApp * mh[2]};
      const double qn = Kokkos::sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2]);
      if (qn > tEps)
        for (int d = 0; d < 3; ++d)
          c0[d] = pf[d] - (sp / qn) * (q[d] / qn) / g.h[d];
    }
  }
  alphaTh = mth[0] * c0[0] + mth[1] * c0[1] + mth[2] * c0[2];
  return branch;
}
 
KOKKOS_INLINE_FUNCTION int vofWettingPlane(const double mfIn[3], double cf, const double nwIn[3],
                                           double cosT, double sinT, double sdfF, int pivot,
                                           double tEps, double mth[3], double& alphaTh,
                                           double& cosApp) {
  return vofWettingPlane(mfIn, cf, nwIn, cosT, sinT, sdfF, pivot, tEps, mth, alphaTh, cosApp,
                         VofMetric{});
}
 
KOKKOS_INLINE_FUNCTION double vofWettingFraction(const double mth[3], double alphaTh,
                                                 const int ds[3]) {
  const double a = alphaTh - (mth[0] * ds[0] + mth[1] * ds[1] + mth[2] * ds[2]);
  const double v = plicVolume(mth[0], mth[1], mth[2], a);
  return v < 0.0 ? 0.0 : (v > 1.0 ? 1.0 : v);
}
 
}  // namespace peclet::core::vof
 
#endif  // PECLET_CORE_VOF_WETTING_HPP