mac_cheb_momentum.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_MAC_CHEB_MOMENTUM_HPP
#define PECLET_FLOW_MAC_CHEB_MOMENTUM_HPP
 
#include <Kokkos_Core.hpp>
 
#include "mac_cutcell_mg.hpp"  // FPC
#include "mac_ibm.hpp"         // CCField/CCConst/CCExec/C3
#include "policy.hpp"
 
namespace peclet::flow {
 
template <class MC>
inline void ibmStencilJacobiBounds(MC AC, MC AW, MC AE, MC AS, MC AN, MC AB, MC AT,
                                   CCConst solidmask, C3 e, int g, double& lo, double& hi) {
  CCExec space;
  const bool hasMask = (solidmask.extent(0) != 0);
  using MD = MDRange3<CCExec>;
  const MD pol(space, {g, g, g}, {e.x - g, e.y - g, e.z - g});
  double rlo = 0.0, rhi = 0.0;
  Kokkos::parallel_reduce(
      "peclet::flow::cheb_bounds_lo", pol,
      KOKKOS_LAMBDA(int lx, int ly, int lz, double& m) {
        const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
        if (hasMask && solidmask(i) > 0.5)
          return;
        const double ac = (double)AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double off = (double)AW(i) + (double)AE(i) + (double)AS(i) + (double)AN(i) +
                           (double)AB(i) + (double)AT(i);
        const double v = (ac + off) / ac;  // row sum / diagonal
        if (v < m)
          m = v;
      },
      Kokkos::Min<double>(rlo));
  Kokkos::parallel_reduce(
      "peclet::flow::cheb_bounds_hi", pol,
      KOKKOS_LAMBDA(int lx, int ly, int lz, double& m) {
        const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
        if (hasMask && solidmask(i) > 0.5)
          return;
        const double ac = (double)AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double off = (double)AW(i) + (double)AE(i) + (double)AS(i) + (double)AN(i) +
                           (double)AB(i) + (double)AT(i);
        const double v = (ac - off) / ac;
        if (v > m)
          m = v;
      },
      Kokkos::Max<double>(rhi));
  lo = rlo;
  hi = rhi;
}
 
inline void ibmChebUpdate(CCField x, CCField d, CCConst r, FPC AC, CCConst solidmask, C3 e, int g,
                          double alpha, double beta) {
  CCExec space;
  const bool hasMask = (solidmask.extent(0) != 0);
  using MD = MDRange3<CCExec>;
  Kokkos::parallel_for(
      "peclet::flow::cheb_update", MD(space, {g, g, g}, {e.x - g, e.y - g, e.z - g}),
      KOKKOS_LAMBDA(int lx, int ly, int lz) {
        const long i = (long)lx + (long)ly * e.x + (long)lz * (long)e.x * e.y;
        if (hasMask && solidmask(i) > 0.5) {
          d(i) = 0.0;
          x(i) = 0.0;
          return;
        }
        const double ac = (double)AC(i);
        if (Kokkos::fabs(ac) < 1e-30)
          return;
        const double dn = beta * d(i) + alpha * r(i) / ac;
        d(i) = dn;
        x(i) += dn;
      });
}
 
struct ChebCoeffs {
  double theta, delta, sigma, rho;
  ChebCoeffs(double lo, double hi)
      : theta(0.5 * (hi + lo)), delta(0.5 * (hi - lo)), sigma(theta / delta), rho(1.0 / sigma) {}
  void next(int it, double& alpha, double& beta) {
    if (it == 0) {
      alpha = 1.0 / theta;
      beta = 0.0;
      return;
    }
    const double rPrev = rho;
    rho = 1.0 / (2.0 * sigma - rPrev);
    alpha = 2.0 * rho / delta;
    beta = rho * rPrev;
  }
};
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_MAC_CHEB_MOMENTUM_HPP