phase_change.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_PHASE_CHANGE_HPP
#define PECLET_FLOW_VOF_PHASE_CHANGE_HPP
 
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
#include "vof/plic.hpp"
 
namespace peclet::flow::vof {
 
// ---------------------------------------------------------------------------------------------
// (1) the PLIC polygon area
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double plicArea(double mx, double my, double mz, double alpha) {
  const double ax = Kokkos::fabs(mx), ay = Kokkos::fabs(my), az = Kokkos::fabs(mz);
  const double s = ax + ay + az;
  if (!(s > 0.0))
    return 0.0;  // degenerate normal: no polygon
 
  double beta = alpha;  // mirror into the positive octant (plic.hpp `plicVolume` convention)
  if (mx < 0.0)
    beta += ax;
  if (my < 0.0)
    beta += ay;
  if (mz < 0.0)
    beta += az;
  const double a1 = beta / s;
  if (!(a1 > 0.0) || !(a1 < 1.0))
    return 0.0;  // the plane misses the cell
 
  const double n1 = Kokkos::fmin(Kokkos::fmin(ax, ay), az) / s;
  const double n3 = Kokkos::fmax(Kokkos::fmax(ax, ay), az) / s;
  const double n2 = Kokkos::fmax(1.0 - n1 - n3, 0.0);
  const double nrm = Kokkos::sqrt(n1 * n1 + n2 * n2 + n3 * n3);  // |m|_2 in the normalized frame
  const double w = Kokkos::fmin(a1, 1.0 - a1);
 
  if (!(n2 > 0.0))
    return nrm;  // axis-aligned plane (n3 = 1): the polygon is a whole cell face, V'(w) = 1
 
  if (!(n1 > 0.0)) {  // 2-D: a rectangle. V'(w) = (w > n2 ? n2 : w)/(n2 n3); <w-n3> == 0 for w<=1/2
    return nrm * ((w > n2) ? n2 : w) / (n2 * n3);
  }
 
  // full 3-D, cancellation-free (see the note above)
  const double A1 = (w > n1) ? n1 * (2.0 * w - n1) : w * w;
  const double u2 = w - n2, u3 = w - n3;
  const double A2 = (u2 > n1) ? n1 * (n1 - 2.0 * u2) : ((u2 > 0.0) ? -u2 * u2 : 0.0);
  const double A3 = (u3 > n1) ? n1 * (n1 - 2.0 * u3) : ((u3 > 0.0) ? -u3 * u3 : 0.0);
  return nrm * (A1 + A2 + A3) / (2.0 * n1 * n2 * n3);
}
 
KOKKOS_INLINE_FUNCTION double plicAreaMetric(double mx, double my, double mz, double alpha,
                                             const VofMetric& g) {
  const double a = plicArea(mx, my, mz, alpha);
  if (!(a > 0.0))
    return a;
  const double m[3] = {mx, my, mz};
  double n[3];
  const double sh = vofPhysNormal(m, g, n);
  return a * g.det() * sh;
}
 
// ---------------------------------------------------------------------------------------------
// (2) plane geometry in the cell frame
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double pcUnitNormal(double mx, double my, double mz, double n[3]) {
  const double q = Kokkos::sqrt(mx * mx + my * my + mz * mz);
  if (q > 0.0) {
    n[0] = mx / q;
    n[1] = my / q;
    n[2] = mz / q;
  }
  return q;
}
 
KOKKOS_INLINE_FUNCTION double pcCentreDistance(double mx, double my, double mz, double alpha) {
  const double q = Kokkos::sqrt(mx * mx + my * my + mz * mz);
  if (!(q > 0.0))
    return 0.0;
  return (0.5 * (mx + my + mz) - alpha) / q;
}
 
KOKKOS_INLINE_FUNCTION double pcOffsetDistance(double phic, const double n[3], double dx, double dy,
                                               double dz) {
  return phic + dx * n[0] + dy * n[1] + dz * n[2];
}
 
// ---------------------------------------------------------------------------------------------
// (3) the one-sided, pure-cell weighted least-squares normal gradient (Malan et al. 2021)
// ---------------------------------------------------------------------------------------------
 
struct PcGradFit {
  double num = 0.0;   
  double den = 0.0;   
  double num2 = 0.0;  
  double s3 = 0.0;    
  double s4 = 0.0;    
  int n = 0;          
};
 
KOKKOS_INLINE_FUNCTION double pcGradWeight(double dx, double dy, double dz, const double n[3]) {
  const double d2 = dx * dx + dy * dy + dz * dz;
  if (!(d2 > 0.0))
    return 0.0;
  const double dn = dx * n[0] + dy * n[1] + dz * n[2];
  return (dn * dn) / (d2 * Kokkos::sqrt(d2));
}
 
KOKKOS_INLINE_FUNCTION void pcGradAdd(PcGradFit& f, double w, double phi, double T, double Tg) {
  const double dT = T - Tg, p2 = phi * phi;
  f.num += w * phi * dT;
  f.den += w * p2;
  f.num2 += w * p2 * dT;
  f.s3 += w * p2 * phi;
  f.s4 += w * p2 * p2;
  ++f.n;
}
 
KOKKOS_INLINE_FUNCTION double pcGradSolve(const PcGradFit& f) {
  return (f.den > 0.0) ? (f.num / f.den) : 0.0;
}
 
KOKKOS_INLINE_FUNCTION double pcCurvedDistance(double phi, double dx, double dy, double dz,
                                               const double n[3], double kappa) {
  const double dn = dx * n[0] + dy * n[1] + dz * n[2];
  const double rho2 = Kokkos::fmax(dx * dx + dy * dy + dz * dz - dn * dn, 0.0);
  return phi + 0.25 * kappa * rho2;
}
 
KOKKOS_INLINE_FUNCTION double pcGradSolve2(const PcGradFit& f) {
  const double det = f.den * f.s4 - f.s3 * f.s3;
  if (f.n < 2 || !(Kokkos::fabs(det) > 1e-12 * f.den * f.s4))
    return pcGradSolve(f);
  return (f.num * f.s4 - f.num2 * f.s3) / det;
}
 
// ---------------------------------------------------------------------------------------------
// (4) the physics rules
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double pcMassFlux(double kg, double gradG, double kl, double gradL,
                                         double hlv) {
  return (kg * gradG - kl * gradL) / hlv;
}
 
KOKKOS_INLINE_FUNCTION double pcRegressVolume(double mdot, double area, double dt, double rhoL) {
  return mdot * area * dt / rhoL;
}
 
KOKKOS_INLINE_FUNCTION double pcDivSource(double mdot, double area, double rhoG, double rhoL) {
  return mdot * area * (1.0 / rhoG - 1.0 / rhoL);
}
 
KOKKOS_INLINE_FUNCTION double pcInterfaceTemperature(double Tsat, double mdot, double Rint) {
  return Tsat + mdot * Rint;
}
 
KOKKOS_INLINE_FUNCTION double pcCarriedValue(double Tg, double grad, double phic) {
  return Tg + grad * phic;
}
 
KOKKOS_INLINE_FUNCTION double pcGfmTheta(double phiC, double nd, double s, double thetaMin,
                                         double thetaMax) {
  const double den = Kokkos::fabs(nd);
  const double num = Kokkos::fabs(phiC - s * nd);
  if (!(num < thetaMax * den))  // includes den == 0: the limit, reached continuously
    return thetaMax;
  return Kokkos::fmax(num / den, thetaMin);
}
 
KOKKOS_INLINE_FUNCTION double pcFaceConductivity(double ki, double kj, bool mi, bool mj) {
  if (mi && !mj)
    return kj;
  if (mj && !mi)
    return ki;
  return 0.5 * (ki + kj);
}
 
KOKKOS_INLINE_FUNCTION double pcPhaseMix(double vg, double vl, double c) {
  return vg * (1.0 - c) + vl * c;
}
 
KOKKOS_INLINE_FUNCTION bool pcPushWeights(const double n[3], double sgn, const bool avail[3],
                                          int step[3], double w[3]) {
  double tot = 0.0, tota = 0.0;
  for (int d = 0; d < 3; ++d) {
    const double p = sgn * n[d];
    w[d] = p * p;
    step[d] = (p > 0.0) ? 1 : ((p < 0.0) ? -1 : 0);
    tot += w[d];
    if (avail[d])
      tota += w[d];
  }
  const bool ok = tota > 0.0;
  for (int d = 0; d < 3; ++d) {
    if (ok)
      w[d] = avail[d] ? w[d] / tota : 0.0;
    else
      w[d] = (tot > 0.0) ? w[d] / tot : 0.0;
  }
  return ok;
}
 
// ---------------------------------------------------------------------------------------------
// (5) WO-P3g — the SECOND-ORDER interfacial energy operator
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double pcGfmThetaK(double phiC, double nd, double s, double kappa,
                                          double thetaMin, double thetaMax) {
  if (kappa == 0.0)
    return pcGfmTheta(phiC, nd, s, thetaMin, thetaMax);
  const double den = Kokkos::fabs(nd);
  const double num = Kokkos::fabs(phiC - s * nd + 0.25 * kappa * (1.0 - nd * nd));
  if (!(num < thetaMax * den))
    return thetaMax;
  return Kokkos::fmax(num / den, thetaMin);
}
 
KOKKOS_INLINE_FUNCTION double pcGfmThetaKAniso(double phiC, const double n[3], int d, double s,
                                               double kappa, double thetaMin, double thetaMax,
                                               const VofMetric& g) {
  const double hd = g.h[d];
  const double nd = n[d];
  const double den = Kokkos::fabs(nd) * hd;
  double num = Kokkos::fabs(phiC - s * hd * nd);
  if (kappa != 0.0)
    num = Kokkos::fabs(phiC - s * hd * nd + 0.25 * kappa * hd * hd * (1.0 - nd * nd));
  if (!(num < thetaMax * den))  // includes den == 0: the limit, reached continuously
    return thetaMax;
  return Kokkos::fmax(num / den, thetaMin);
}
 
struct PcGfmRow {
  double aGamma;   
  double aBehind;  
};
KOKKOS_INLINE_FUNCTION PcGfmRow pcGfmRow(double theta, bool behind, int order) {
  if (order < 2 || !behind)
    return PcGfmRow{1.0 / theta, 1.0};
  return PcGfmRow{2.0 / ((1.0 + theta) * theta), 2.0 / (1.0 + theta)};
}
 
KOKKOS_INLINE_FUNCTION double pcOperatorMassFlux(double q, double coefSum, double area, double hlv,
                                                 double Rint) {
  const double den = hlv * area + Rint * coefSum;
  return (den > 0.0) ? q / den : 0.0;
}
 
KOKKOS_INLINE_FUNCTION bool pcIsInterfacial(double c, double eps) {
  return c > eps && c < 1.0 - eps;
}
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_PHASE_CHANGE_HPP