fv_op.hpp

Go to the documentation of this file.

// core — device (Kokkos) Poisson operator + smoother on a BlockOctreeView.
//
// The device compute path for the AMR octree: the Laplacian matvec and a weighted-
// Jacobi smoother run as Kokkos::parallel_for over the leaf Views, using the
// device-callable face-neighbour walk of BlockOctreeView (block_octree_kokkos.hpp).
// Same arithmetic as the host BlockOctree::faceNeighbor, so the result is identical
// to the host operator (validated bit-for-bit in tests/test_amr_device_poisson_kokkos).
//
// Scope: the plain (uniform, openness-free) Laplacian over the leaves; a face with
// no neighbour (block edge) is skipped (homogeneous Neumann). Openness / graded /
// distributed-on-device build on this + the device halo.
//
// Requires a Kokkos build (MORTON_ENABLE_KOKKOS ⇒ MORTON_HD = KOKKOS_FUNCTION) and
#ifndef PECLET_AMR_FV_OP_HPP
#define PECLET_AMR_FV_OP_HPP
 
#include "peclet/amr/block_octree_view.hpp"
#include "peclet/amr/common.hpp"
#include "peclet/amr/face_csr.hpp"  // shared host+device FV (weight-CSR) row kernels
#include "peclet/core/common/view.hpp"
 
namespace peclet::amr {
 
template <int Dim, unsigned Bits>
void laplacian(BlockOctreeView<Dim, Bits> dev, View<double> x, View<double> y, double inv) {
  const Index n = dev.numLeaves();
  Kokkos::parallel_for(
      "amr::device_laplacian", n, KOKKOS_LAMBDA(const Index i) {
        double s = 0.0;
        for (int axis = 0; axis < Dim; ++axis)
          for (int dir = -1; dir <= 1; dir += 2) {
            Index j = dev.faceNeighbor(i, axis, dir);
            if (j >= 0)
              s += x(j) - x(i);
          }
        y(i) = inv * s;
      });
}
 
template <int Dim, unsigned Bits>
void jacobiSweep(BlockOctreeView<Dim, Bits> dev, View<double> x, View<const double> b,
                 View<double> lx, double inv, double omega) {
  const Index n = dev.numLeaves();
  const double diag = 2.0 * Dim * inv;
  // lx = L x = ∇² x
  Kokkos::parallel_for(
      "amr::device_jacobi_apply", n, KOKKOS_LAMBDA(const Index i) {
        double s = 0.0;
        for (int axis = 0; axis < Dim; ++axis)
          for (int dir = -1; dir <= 1; dir += 2) {
            Index j = dev.faceNeighbor(i, axis, dir);
            if (j >= 0)
              s += x(j) - x(i);
          }
        lx(i) = inv * s;
      });
  Kokkos::parallel_for(
      "amr::device_jacobi_update", n,
      KOKKOS_LAMBDA(const Index i) { x(i) += omega * (lx(i) - b(i)) / diag; });
}
 
// ===========================================================================
// Consistent FV (graded) operator as a precomputed face CSR.
//
// The plain operator above ignores 2:1-interface geometry (it treats every face
// as same-level, spacing h0). The *consistent* conservative FV Laplacian is
//   (L u)_i = invVol_i · Σ_{f ∈ [start_i, start_{i+1})} w_f · (u[nbr_f] − u_i),
// with w_f = openness · A_f/d_f and the coarse side of a 2:1 interface carrying
// one face per fine sub-neighbour — exactly the host AmrPoisson enumeration
// (forEachFaceNeighbor). The CSR is built once on the host in that same face
// order, so the device operator is bit-identical to the host operator (and the
// graded V-cycle converges instead of stalling on the inconsistent plain op).
// ===========================================================================
struct FvOp {
  View<double> invVol;    
  View<Index> faceStart;  
  View<Index> faceNbr;    
  View<double> faceW;     
  View<double> bcDiag;    
  Index n = 0;
  // Helmholtz generalisation: the applied operator is H = c0·I + cD·L (c0=0, cD=1 ⇒ the
  // pure Laplacian L, the default — bit-exact unchanged). Setting c0=idiag, cD=−μ turns the
  // existing hierarchy into the momentum operator idiag·I − μ∇², so Multigrid can be
  // used as a (non-singular when c0≠0) preconditioner for the momentum BiCGStab. The L path
  // is bit-exact preserved: c0·u+cD·(L) with c0=0,cD=1 is L exactly in IEEE arithmetic, and
  // the Jacobi point-solve branches on (c0==0 && cD==1) to keep the original expression.
  double c0 = 0.0;
  double cD = 1.0;
};
 
inline FvCsrOpT<View<const double>, View<const Index>> fvView(const FvOp& op) {
  FvCsrOpT<View<const double>, View<const Index>> v;
  v.n = op.n;
  v.invVol = op.invVol;
  v.coef = op.faceW;
  v.start = op.faceStart;
  v.nbr = op.faceNbr;
  v.bcDiag = op.bcDiag;
  v.c0 = op.c0;
  v.cD = op.cD;
  return v;
}
 
inline void applyFv(const FvOp& op, View<const double> u, View<double> Lu) {
  const auto A = fvView(op);
  Kokkos::parallel_for(
      "amr::fv_apply", op.n, KOKKOS_LAMBDA(const Index i) { Lu(i) = fvApplyRow(A, i, u); });
}
 
inline void residualFv(const FvOp& op, View<const double> u, View<const double> rhs,
                       View<double> res) {
  const auto A = fvView(op);
  Kokkos::parallel_for(
      "amr::fv_residual", op.n,
      KOKKOS_LAMBDA(const Index i) { res(i) = rhs(i) - fvApplyRow(A, i, u); });
}
 
inline void jacobiFv(const FvOp& op, View<double> u, View<const double> rhs, View<double> tmp,
                     double omega) {
  const auto A = fvView(op);
  Kokkos::parallel_for(
      "amr::fv_jacobi_compute", op.n,
      KOKKOS_LAMBDA(const Index i) { tmp(i) = fvPointSolve(A, i, u, rhs(i), u(i)); });
  Kokkos::parallel_for(
      "amr::fv_jacobi_update", op.n,
      KOKKOS_LAMBDA(const Index i) { u(i) = (1.0 - omega) * u(i) + omega * tmp(i); });
}
 
inline void removeMeanFv(const FvOp& op, View<double> u) {
  auto invVol = op.invVol;
  auto fs = op.faceStart;
  auto fw = op.faceW;
  auto bc = op.bcDiag;
  double sx = 0.0, sv = 0.0;
  // Numerator and denominator share the same active-cell diagonal d = bc + Σfw; fold them into a
  // single multi-output reduction so that diagonal is walked once instead of twice. Each output
  // keeps its original term verbatim (u/invVol, 1/invVol) and its own Sum reducer over the same
  // range, so the result is bit-identical to the two separate reduces. (The subtract pass below
  // recomputes d a third time but is a parallel_for, where reloading a cached mask would cost the
  // same as the walk.)
  Kokkos::parallel_reduce(
      "amr::fv_rmean", op.n,
      KOKKOS_LAMBDA(const Index i, double& an, double& ad) {
        double d = bc(i);
        for (Index k = fs(i); k < fs(i + 1); ++k)
          d += fw(k);
        if (d > 1e-30) {
          an += u(i) / invVol(i);  // V_i = 1/invVol_i
          ad += 1.0 / invVol(i);
        }
      },
      sx, sv);
  if (sv <= 0.0)
    return;
  const double m = sx / sv;
  Kokkos::parallel_for(
      "amr::fv_rmean_sub", op.n, KOKKOS_LAMBDA(const Index i) {
        double d = bc(i);
        for (Index k = fs(i); k < fs(i + 1); ++k)
          d += fw(k);
        if (d > 1e-30)
          u(i) -= m;
      });
}
 
inline void quadDelta(View<const Index> qStart, View<const Index> qSlot, View<const double> qCoef,
                      View<const double> u, View<double> dq, Index n) {
  Kokkos::parallel_for(
      "amr::fv_quad_delta", n, KOKKOS_LAMBDA(const Index i) {
        double acc = 0.0;
        for (Index k = qStart(i); k < qStart(i + 1); ++k)
          acc += qCoef(k) * u(qSlot(k));
        dq(i) = acc;
      });
}
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_FV_OP_HPP