advect_recon.hpp

Go to the documentation of this file.

// core — shared, backend-agnostic high-order advection face reconstruction.
//
// The SOU / Koren-TVD reconstruction of the advected face value from the two upwind cells (upup,
// up) and the downwind cell (down) is the numerically delicate part of the advection scheme —
// exactly the kind of formula that must not drift between the serial host solver (AmrFlow::hoFace,
// flow.hpp) and the Kokkos device solver (deferredSou / advectExplicit, flow.hpp). This is that
// formula, written ONCE as a MORTON_HD function (host- and device-callable; empty MORTON_HD in the
// pure-C++ build). min/max/abs are expressed with branches rather than std::/Kokkos:: math so the
// body is identical on both backends (and bit-identical to the previous std::fmin/fmax form for the
// finite velocity/scalar data the solver produces).
#ifndef PECLET_AMR_ADVECT_RECON_HPP
#define PECLET_AMR_ADVECT_RECON_HPP
 
#include <morton/morton.hpp>
 
#include "peclet/amr/common.hpp"
#ifndef MORTON_HD
#define MORTON_HD
#endif
 
namespace peclet::amr {
 
MORTON_HD inline double korenPsi(double r) {
  const double t = (1.0 + 2.0 * r) / 3.0;
  double m = (2.0 * r < t) ? (2.0 * r) : t;
  if (m > 2.0)
    m = 2.0;
  return (m > 0.0) ? m : 0.0;
}
 
MORTON_HD inline double hoFaceValue(double upup, double up, double down, int scheme) {
  if (scheme == 0)
    return 1.5 * up - 0.5 * upup;  // SOU
  const double den = down - up;
  const double aden = (den < 0.0) ? -den : den;
  const double r = (aden < 1e-10) ? 0.0 : (up - upup) / den;
  return up + 0.5 * korenPsi(r) * den;  // Koren TVD
}
 
MORTON_HD inline double hoFaceValueSeam(double upup, double up, double upStar, double downStar,
                                        double d1, double dUU, double dD, int scheme) {
  const double gUp = (up - upup) / dUU;
  if (scheme == 0)
    return upStar + d1 * gUp;  // SOU
  const double den = downStar - upStar;
  const double aden = (den < 0.0) ? -den : den;
  const double r = (aden < 1e-10) ? 0.0 : gUp * dD / den;
  return upStar + (d1 / dD) * korenPsi(r) * den;  // Koren TVD
}
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_ADVECT_RECON_HPP