face_csr.hpp

Go to the documentation of this file.

// core — shared, backend-agnostic kernels for an assembled face-CSR operator.
//
// Lifted verbatim from the AMR tree (peclet/core/amr/face_csr.hpp, 2026-09-10; QUALITY_PLAN G.2) —
// the AMR package and voro's mesh optimiser both consume it. The "host" serial reference solver
// (peclet::amr::AmrCutCell / oracle::AmrFlow, pure C++20) and the "device" Kokkos solver
// (MomentumOp + csr_operator.hpp) used to carry two independent encodings of the SAME per-cell
// arithmetic — a real drift risk. This header is the single source of that arithmetic: the
// assembled operator
//
//     (A u)_i = diag_i · u_i + Σ_{k ∈ [start_i, start_{i+1})} coef_k · u[nbr_k]
//                              ( + optional implicit-FOU advection over a second CSR )
//
// expressed once as PECLET_HD row kernels over a *templated accessor*, so the SAME body serves
//   - the host: a plain serial loop, arrays = raw pointers (no Kokkos dependency — PECLET_HD is
//     plain `inline` in a non-Kokkos, non-CUDA build, so these compile as ordinary functions); and
//   - the device: inside Kokkos::parallel_for, arrays = Kokkos::View.
//
// Mesh-agnostic: it names no octree types and assumes nothing beyond the CSR, so it is equally the
// kernel layer for a Voronoi-cell (unstructured polyhedral) operator.
#ifndef PECLET_CORE_SOLVER_FACE_CSR_HPP
#define PECLET_CORE_SOLVER_FACE_CSR_HPP
 
#include "peclet/core/common/types.hpp"
 
// PECLET_HD (common/portable.hpp): KOKKOS_INLINE_FUNCTION under a Kokkos build, __host__ __device__
// inline under nvcc/hipcc, plain inline otherwise — so the row kernels are device-callable when
// compiled for the device and ordinary host functions in the pure-C++ build.
#include "peclet/core/common/portable.hpp"
 
namespace peclet::core::solver {
 
template <class T>
struct HostArr {
  const T* p = nullptr;
  HostArr() = default;
  explicit HostArr(const T* ptr) : p(ptr) {}
  PECLET_HD T operator()(Index i) const { return p[i]; }
};
 
template <class D, class I>
struct FaceCsrOpT {
  Index n = 0;
  D diag, coef;
  I start, nbr;
  bool hasAdv = false;
  D advDiag, advCoef;
  I advStart, advNbr;
};
 
template <class Op, class U>
PECLET_HD double faceCsrApplyRow(const Op& op, Index i, const U& u) {
  double acc = op.diag(i) * u(i);
  for (Index k = op.start(i); k < op.start(i + 1); ++k)
    acc += op.coef(k) * u(op.nbr(k));
  if (op.hasAdv) {
    acc += op.advDiag(i) * u(i);
    for (Index k = op.advStart(i); k < op.advStart(i + 1); ++k)
      acc += op.advCoef(k) * u(op.advNbr(k));
  }
  return acc;
}
 
template <class Op, class U>
PECLET_HD void faceCsrOffDiag(const Op& op, Index i, const U& u, double& off, double& d) {
  off = 0.0;
  d = op.diag(i);
  for (Index k = op.start(i); k < op.start(i + 1); ++k)
    off += op.coef(k) * u(op.nbr(k));
  if (op.hasAdv) {
    for (Index k = op.advStart(i); k < op.advStart(i + 1); ++k)
      off += op.advCoef(k) * u(op.advNbr(k));
    d += op.advDiag(i);
  }
}
 
PECLET_HD double faceCsrPointUpdate(double b_i, double off, double d, double uOld, double omega) {
  const double nu = (d != 0.0) ? (b_i - off) / d : uOld;
  return (1.0 - omega) * uOld + omega * nu;
}
 
// ---------------------------------------------------------------------------
// The conservative FV (graded-Laplacian) operator — the *pressure* form. Unlike the general
// momentum operator above it stores per-face conductances `w` and a per-cell `invVol`, with the
// diagonal *derived* (the symmetric Laplacian Σ w·(u_nbr − u_i)); a Helmholtz generalisation
//   H u = c0·u + cD·( invVol·( Σ w·(u_nbr − u_i) − bcDiag·u ) )
// (c0=0, cD=1 ⇒ the pure FV Laplacian L) keeps the L path bit-exact. Same body for host
// (AmrPoisson) and device (poisson.hpp FvOp).
// ---------------------------------------------------------------------------
 
template <class D, class I>
struct FvCsrOpT {
  Index n = 0;
  D invVol, coef;  // coef = per-face conductance w = openness·A_f/d_f
  I start, nbr;
  D bcDiag;  // homogeneous-Dirichlet boundary diagonal (0 when periodic)
  double c0 = 0.0, cD = 1.0;
};
 
template <class Op, class U>
PECLET_HD double fvApplyRow(const Op& op, Index i, const U& u) {
  const double ui = u(i);
  double acc = 0.0;
  for (Index k = op.start(i); k < op.start(i + 1); ++k)
    acc += op.coef(k) * (u(op.nbr(k)) - ui);
  return op.c0 * ui + op.cD * (op.invVol(i) * (acc - op.bcDiag(i) * ui));
}
 
template <class Op, class U>
PECLET_HD double fvPointSolve(const Op& op, Index i, const U& u, double rhs_i, double uOld) {
  double sumOff = 0.0, sw = 0.0;
  for (Index k = op.start(i); k < op.start(i + 1); ++k) {
    sumOff += op.coef(k) * u(op.nbr(k));
    sw += op.coef(k);
  }
  const double swbc = sw + op.bcDiag(i);
  if (op.c0 == 0.0 && op.cD == 1.0)
    return (swbc != 0.0) ? (sumOff - rhs_i / op.invVol(i)) / swbc : uOld;
  const double Hii = op.c0 - op.cD * op.invVol(i) * swbc;
  const double Hoff = op.cD * op.invVol(i) * sumOff;
  return (Hii != 0.0) ? (rhs_i - Hoff) / Hii : uOld;
}
 
}  // namespace peclet::core::solver
 
#endif  // PECLET_CORE_SOLVER_FACE_CSR_HPP