energy_advect.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_ENERGY_ADVECT_HPP
#define PECLET_FLOW_VOF_ENERGY_ADVECT_HPP
 
#include <functional>
#include <Kokkos_Core.hpp>
#include <stdexcept>
 
#include "policy.hpp"
#include "vof/advect_wy.hpp"
#include "vof/momentum_advect.hpp"  // vofMinmod
#include "vof/phase_change.hpp"
 
namespace peclet::flow::vof {
 
class VofEnergyAdvector {
 public:
  void init(const WyAdvector& w, double rcpG, double rcpL) {
    if (!(rcpG > 0.0) || !(rcpL > 0.0))
      throw std::invalid_argument(
          "peclet::flow::vof::VofEnergyAdvector: both phase rho*c_p must be > 0");
    n_ = w.inner();
    e_ = w.extent();
    g_ = w.ghost();
    len_ = static_cast<long>(e_.x) * e_.y * e_.z;
    rcpG_ = rcpG;
    rcpL_ = rcpL;
    T_ = SField("vof::energyT", len_);
    Tw_ = SField("vof::energyTw", len_);
  }
  bool initialized() const { return len_ > 0; }
  void setPhaseRcp(double g, double l) {
    rcpG_ = g;
    rcpL_ = l;
  }
  double phaseRcpG() const { return rcpG_; }
  double phaseRcpL() const { return rcpL_; }
  SField temperature() const { return T_; }
 
  bool energyMuscl = false;
 
  std::function<void(SField)> exchange;
 
  void advect(WyAdvector& w, double dt, long step) {
    if (!initialized())
      throw std::runtime_error("peclet::flow::vof::VofEnergyAdvector: init() was never called");
    if (!exchange)
      throw std::runtime_error("peclet::flow::vof::VofEnergyAdvector: no `exchange` hook");
    if (w.hasGeometry())
      throw std::runtime_error(
          "peclet::flow::vof::VofEnergyAdvector: cut-cell geometry is not composed with the "
          "consistent energy transport at this rung (enable_phase_change already refuses a solid)");
    w.requireExchange();
    w.checkCourant(dt);
    const double dth = dt / w.h();
    w.freezeDilationFlag();
    w.timedExchange([&] { exchange(T_); });
    const int* perm = kWySweepPerm[static_cast<int>(step % 6)];
    for (int s = 0; s < 3; ++s) {
      const int d = perm[s];
      w.reconstruct();
      w.computeFluxes(d, dth);
      w.applySweep(d, dth);     // C^{(s+1)}: the capacity the recovery divides by
      energyUpdate(w, d, dth);  // T from the pre-sweep T, the fluxes, and that new C (swaps T_)
      w.exchangeTimed(w.colour());
      w.timedExchange([&] { exchange(T_); });  // regenerates every ghost of the swapped-in buffer
    }
  }
 
  void extrema(double& lo, double& hi) const {
    const I3 e = e_, n = n_;
    const int g = g_;
    SField T = T_;
    double mn = 0.0, mx = 0.0;
    using MD = MDRange3<SExec>;
    MD pol(SExec(), {g, g, g}, {g + n.x, g + n.y, g + n.z});
    Kokkos::parallel_reduce(
        "vof::energy::min", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& a) { a = Kokkos::fmin(a, T(L3(x, y, z, e))); },
        Kokkos::Min<double>(mn));
    Kokkos::parallel_reduce(
        "vof::energy::max", pol,
        KOKKOS_LAMBDA(int x, int y, int z, double& a) { a = Kokkos::fmax(a, T(L3(x, y, z, e))); },
        Kokkos::Max<double>(mx));
    Kokkos::fence();
    lo = mn;
    hi = mx;
  }
 
  // nvcc requires members holding extended device lambdas to be public.
 
  void energyUpdate(WyAdvector& w, int d, double dth) {
    const I3 e = e_, n = n_;
    const int g = g_;
    const long sd =
        d == 0 ? 1 : (d == 1 ? static_cast<long>(e_.x) : static_cast<long>(e_.x) * e_.y);
    // The update is OUT OF PLACE: every cell reads its two upwind neighbours' temperatures, so
    // writing `T` in place is a read-write race between threads (measured: it moved the P2
    // interface position by 0.6 % and turned an order-2.8 ladder into order 0.5). `Tw_` takes the
    // new values and the two handles are swapped; `exchange` then regenerates every ghost.
    SField T = T_, Tw = Tw_, fl = w.faceFlux(), u = w.faceVel(d), c = w.colour();
    const double rg = rcpG_, rl = rcpL_;
    const double floor = 1e-12 * (rg < rl ? rg : rl);
    const bool muscl = energyMuscl;
    Kokkos::parallel_for(
        "vof::energy::update", 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);
          const double aP = u(i) * dth, aM = u(i - sd) * dth;
          const double FP = fl(i), FM = fl(i - sd);
          const double phiP = rg * (aP - FP) + rl * FP;
          const double phiM = rg * (aM - FM) + rl * FM;
          const double ti = T(i);
          double tP, tM;
          if (!muscl) {
            tP = (aP > 0.0) ? ti : T(i + sd);
            tM = (aM > 0.0) ? T(i - sd) : ti;
          } else {
            // the flux slab of a face occupies the fraction |a| of the DONOR next to that face, so
            // its centroid sits (1/2 - |a|/2) from the donor centre TOWARDS the face
            const double wP = 0.5 - 0.5 * Kokkos::fabs(aP), wM = 0.5 - 0.5 * Kokkos::fabs(aM);
            if (aP > 0.0)
              tP = ti + wP * vofMinmod(ti - T(i - sd), T(i + sd) - ti);
            else
              tP = T(i + sd) - wP * vofMinmod(T(i + sd) - ti, T(i + 2 * sd) - T(i + sd));
            if (aM > 0.0)
              tM = T(i - sd) + wM * vofMinmod(T(i - sd) - T(i - 2 * sd), ti - T(i - sd));
            else
              tM = ti - wM * vofMinmod(ti - T(i - sd), T(i + sd) - ti);
          }
          const double cl = Kokkos::fmin(Kokkos::fmax(c(i), 0.0), 1.0);  // C AFTER this sweep
          const double r = pcPhaseMix(rg, rl, cl);
          Tw(i) = ti + (phiM * (tM - ti) - phiP * (tP - ti)) / (r > floor ? r : floor);
        });
    Kokkos::fence();
    SField tmp = T_;
    T_ = Tw_;
    Tw_ = tmp;
  }
 
 private:
  I3 n_{0, 0, 0}, e_{0, 0, 0};
  int g_ = 0;
  long len_ = 0;
  double rcpG_ = 1.0, rcpL_ = 1.0;
  SField T_, Tw_;
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_ENERGY_ADVECT_HPP