PoPS/doxygen/polar__operator_8hpp_source.html

#pragma once


#include <pops/core/state/state.hpp>

#include <pops/core/foundation/types.hpp>

#include <pops/mesh/layout/box_array.hpp>

#include <pops/mesh/storage/fab2d.hpp>

#include <pops/mesh/execution/for_each.hpp>

#include <pops/mesh/geometry/geometry.hpp>  // PolarGeometry

#include <pops/mesh/storage/multifab.hpp>

#include <pops/numerics/fv/numerical_flux.hpp>

#include <pops/numerics/fv/reconstruction.hpp>

#include <pops/numerics/spatial_operator.hpp>  // reconstruct<>, load_state/load_aux (REUSED verbatim)


#include <concepts>

#include <utility>

#include <vector>


namespace pops {


namespace detail {


// The source term is OPTIONAL in polar transport (Phase 1): a pure transport brick

// (ExBVelocityPolar) does not expose source(), only a composite model (CompositeModel) does. We

// detect it and fall back to 0 otherwise -> the polar operator works just as well on the standalone

// brick as on a composite model, without requiring source() (the Cartesian path, for its part,

// assumes a composite model).

template <class M>


concept PolarHasSource = requires(const M m, const typename M::State u, const Aux a) {

  { m.source(u, a) } -> std::convertible_to<typename M::State>;

};


template <class Model>


POPS_HD inline typename Model::State polar_source(const Model& m, const typename Model::State& u,

                                                 const Aux& a) {

  if constexpr (PolarHasSource<Model>)

    return m.source(u, a);

  else

    return typename Model::State{};

}


template <class M>


concept PolarHasGeomSource = requires(const M m, const typename M::State u, Real r) {

  { m.polar_geom_source(u, r) } -> std::convertible_to<typename M::State>;

};


template <class Model>


POPS_HD inline typename Model::State polar_geom_source(const Model& m,

                                                      const typename Model::State& u, Real r) {

  if constexpr (PolarHasGeomSource<Model>)

    return m.polar_geom_source(u, r);

  else

    return typename Model::State{};

}


// RADIAL FACE FLUX kernel (dir 0): numerical flux at the radial face i (between i-1 and i), ALREADY

// weighted by the face radius r_face(i). We store r_{i-1/2} Fr_{i-1/2} so the divergence is a

// simple difference. NAMED FUNCTOR (device-clean cross-TU).

template <class Limiter, class NumericalFlux, class Model>


struct PolarFaceFluxRKernel {

  Model model;

  ConstArray4 u, ax;

  Array4 fr;       // output: r_face(i) * Fr at the radial face i (ncomp components)

  Real r_min, dr;  // radial geometry (r_face(i) = r_min + i*dr)

  Limiter lim;

  NumericalFlux nflux;

  bool recon_prim;

  // Optional RADIAL WALL (no-penetration). wall_radial == false (default): no effect, boundary flux

  // computed like the interior (BIT-IDENTICAL to the history: MMS, azimuthal conservation). true:

  // the radial flux at BOTH physical boundary faces (i = i_lo_face = lo, i = i_hi_face = hi+1) is

  // forced to ZERO -> the radial term telescopes EXACTLY (each interior face is shared, the

  // boundaries no longer count) -> mass Sum n r dr dtheta conserved to machine precision, whatever

  // v_r (solid wall).

  bool wall_radial;

  int i_lo_face,

      i_hi_face;  // FACE indices of physical boundaries (lo and hi+1); ignored if !wall_radial

  Real pos_floor = Real(0);

  int pos_comp = 0;


  POPS_HD void operator()(int i, int j) const {

    const Real rf = r_min + i * dr;  // r_face(i) (positive on the annulus: r_min >= 0, i >= 0)

    if (wall_radial && (i == i_lo_face || i == i_hi_face)) {

      for (int c = 0; c < Model::n_vars; ++c)

        fr(i, j, c) = Real(0);  // wall: zero radial flux

      return;

    }

    // Reconstructed states on either side of the radial face i (REUSES Cartesian reconstruct_pp<>,

    // dir == 0). L = extrapolation from cell i-1 toward its + face; R = from cell i toward its

    // - face.

    const auto L =

        reconstruct_pp<Model>(model, u, i - 1, j, 0, +1, lim, recon_prim, pos_floor, pos_comp);

    const auto Rr =

        reconstruct_pp<Model>(model, u, i, j, 0, -1, lim, recon_prim, pos_floor, pos_comp);

    const auto F = nflux(model, L, load_aux<aux_comps<Model>()>(ax, i - 1, j), Rr,

                         load_aux<aux_comps<Model>()>(ax, i, j), 0);

    for (int c = 0; c < Model::n_vars; ++c)

      fr(i, j, c) = rf * F[c];

  }


};


// AZIMUTHAL FACE FLUX kernel (dir 1): numerical flux at the theta face j (between j-1 and j). NO

// weighting by r (the azimuthal metric 1/r is applied in cell, cf. PolarAssembleRhsKernel).

template <class Limiter, class NumericalFlux, class Model>


struct PolarFaceFluxThetaKernel {

  Model model;

  ConstArray4 u, ax;

  Array4 ft;  // output: Ftheta at the azimuthal face j (ncomp components)

  Limiter lim;

  NumericalFlux nflux;

  bool recon_prim;

  Real pos_floor = Real(0);

  int pos_comp = 0;


  POPS_HD void operator()(int i, int j) const {

    const auto L =

        reconstruct_pp<Model>(model, u, i, j - 1, 1, +1, lim, recon_prim, pos_floor, pos_comp);

    const auto Rr =

        reconstruct_pp<Model>(model, u, i, j, 1, -1, lim, recon_prim, pos_floor, pos_comp);

    const auto F = nflux(model, L, load_aux<aux_comps<Model>()>(ax, i, j - 1), Rr,

                         load_aux<aux_comps<Model>()>(ax, i, j), 1);

    for (int c = 0; c < Model::n_vars; ++c)

      ft(i, j, c) = F[c];

  }


};


// Polar residual assembly kernel in cell (i, j):

//   R = S - (1/r_i) (Fr_weighted_{i+1} - Fr_weighted_i) / dr - (1/r_i) (Ftheta_{j+1} - Ftheta_j) / dtheta

// Fr_weighted = r_face * Fr (already produced by PolarFaceFluxRKernel). NAMED FUNCTOR.

template <class Model>


struct PolarAssembleRhsKernel {

  Model model;

  ConstArray4 u, ax, fr, ft;  // state, aux, r-weighted radial flux, azimuthal flux

  Array4 r;                   // output: residual

  Real r_min, dr, dtheta;


  POPS_HD void operator()(int i, int j) const {

    const Real ri = r_min + (i + Real(0.5)) * dr;  // r_cell(i)

    const Real inv_r = Real(1) / ri;

    const Aux Ac = load_aux<aux_comps<Model>()>(ax, i, j);

    const auto Us = load_state<Model>(u, i, j);

    const auto S = polar_source<Model>(model, Us, Ac);

    // GEOMETRIC curvature term (centrifugal + cross curvature): not captured by the conservative

    // divergence in the rotating local basis (e_r, e_theta). Zero for scalar bricks (polar ExB =

    // bit-identical); non-zero for a polar fluid (IsothermalFluxPolar). See PolarHasGeomSource /

    // IsothermalFluxPolar for the derivation. r_cell(i) > 0 (annulus).

    const auto Sg = polar_geom_source<Model>(model, Us, ri);

    for (int c = 0; c < Model::n_vars; ++c) {

      const Real div_r = (fr(i + 1, j, c) - fr(i, j, c)) / dr;      // d_r(r Fr) discrete

      const Real div_t = (ft(i, j + 1, c) - ft(i, j, c)) / dtheta;  // d_theta(Ftheta) discrete

      r(i, j, c) = S[c] + Sg[c] - inv_r * (div_r + div_t);

    }

  }


};


// Radial / azimuthal FACE boxes (cf. Cartesian xface_box/yface_box): radial faces = nr+1 x ntheta,

// azimuthal faces = nr x ntheta+1. Reuses the Cartesian helpers (same face index conventions) to

// avoid duplication.

}  // namespace detail


template <class Limiter = NoSlope, class NumericalFlux = RusanovFlux, class Model>


void assemble_rhs_polar(const Model& model, const MultiFab& U, const MultiFab& aux,

                        const PolarGeometry& geom, MultiFab& R, bool recon_prim = false,

                        bool wall_radial = false, Real pos_floor = Real(0)) {

  // STATE-GHOST WIDTH: exactly Limiter::n_ghost, like the Cartesian operator. The polar face kernels

  // (PolarFaceFluxRKernel / PolarFaceFluxThetaKernel) reuse reconstruct_pp<> VERBATIM at the SAME

  // i-1/i (radial) and j-1/j (azimuthal) offsets over the SAME face boxes (xface_box/yface_box, up

  // to hi+1) as compute_face_fluxes -> identical reconstruction stencil (i+-Limiter::n_ghost at the

  // edge of the valid box). The polar metric (r_face = r_min + i*dr, r_cell, 1/r) is computed from

  // INDICES, never read from U, so it adds NO state-ghost width; aux is read at i+-1 only (1 ghost,

  // narrower). HOST-only guard, BEFORE the pass-1/pass-2 loops -- never inside a kernel.

  detail::require_reconstruction_ghosts<Limiter>(U);  // state ghosts >= stencil (otherwise OOB)

  const int pos_comp = detail::positivity_comp<Model>(pos_floor);

  const Real r_min = geom.r_min, dr = geom.dr(), dtheta = geom.dtheta();

  // Physical radial boundary faces (wall): r_min at the lo face of the index domain, r_max at the

  // hi+1 face. geom.domain is the GLOBAL domain (System's PolarGeometry covers the whole annulus,

  // mono-box).

  const int i_lo_face = geom.domain.lo[0];

  const int i_hi_face = geom.domain.hi[0] + 1;

  const Limiter lim{};

  const NumericalFlux nflux{};

  // FACE BoxArrays (cf. Cartesian compute_face_fluxes): radial faces = surroundingNodes in x

  // (xface_box), azimuthal faces in y (yface_box). Same face index conventions -> fr(i, .) is the

  // radial face between i-1 and i, ft(., j) the azimuthal face between j-1 and j.

  std::vector<Box2D> rfaces, tfaces;

  rfaces.reserve(U.box_array().size());

  tfaces.reserve(U.box_array().size());

  for (const Box2D& b : U.box_array().boxes()) {

    rfaces.push_back(xface_box(b));

    tfaces.push_back(yface_box(b));

  }

  MultiFab Fr(BoxArray(std::move(rfaces)), U.dmap(), Model::n_vars, 0);

  MultiFab Ft(BoxArray(std::move(tfaces)), U.dmap(), Model::n_vars, 0);

  for (int li = 0; li < U.local_size(); ++li) {

    const ConstArray4 u = U.fab(li).const_array();

    const ConstArray4 ax = aux.fab(li).const_array();

    Array4 fr = Fr.fab(li).array();

    Array4 ft = Ft.fab(li).array();

    const Box2D v = R.box(li);

    // Radial faces: i in [lo..hi+1], j in [lo..hi] (cf. xface_box).

    for_each_cell(xface_box(v), detail::PolarFaceFluxRKernel<Limiter, NumericalFlux, Model>{

                                    model, u, ax, fr, r_min, dr, lim, nflux, recon_prim,

                                    wall_radial, i_lo_face, i_hi_face, pos_floor, pos_comp});

    // Azimuthal faces: i in [lo..hi], j in [lo..hi+1] (cf. yface_box).

    for_each_cell(yface_box(v), detail::PolarFaceFluxThetaKernel<Limiter, NumericalFlux, Model>{

                                    model, u, ax, ft, lim, nflux, recon_prim, pos_floor, pos_comp});

  }

  for (int li = 0; li < U.local_size(); ++li) {

    const ConstArray4 u = U.fab(li).const_array();

    const ConstArray4 ax = aux.fab(li).const_array();

    const ConstArray4 fr = Fr.fab(li).const_array();

    const ConstArray4 ft = Ft.fab(li).const_array();

    Array4 r = R.fab(li).array();

    const Box2D v = R.box(li);

    for_each_cell(

        v, detail::PolarAssembleRhsKernel<Model>{model, u, ax, fr, ft, r, r_min, dr, dtheta});

  }

}


}  // namespace pops

box_array.hpp
BoxArray: the set of boxes tiling a level (disjoint, covering).

pops::BoxArray
Ordered list of boxes tiling a level.
Definition box_array.hpp:22

pops::BoxArray::size
int size() const
Number of boxes in the tiling.
Definition box_array.hpp:42

pops::BoxArray::boxes
const std::vector< Box2D > & boxes() const
View on the underlying vector (element-by-element equality = same boxes AND same order).
Definition box_array.hpp:46

pops::Fab2D::const_array
ConstArray4 const_array() const
READ handle (POD device-copyable) over this Fab. Valid as long as the Fab lives.
Definition fab2d.hpp:96

pops::Fab2D::array
Array4 array()
WRITE handle (POD device-copyable) over this Fab. Valid as long as the Fab lives.
Definition fab2d.hpp:91

pops::MultiFab
Field distributed over a level: decomposition (BoxArray) + distribution (DistributionMapping) + ncomp...
Definition multifab.hpp:33

pops::MultiFab::fab
Fab2D & fab(int li)
Local fab at index li (0 <= li < local_size()), for writing.
Definition multifab.hpp:67

pops::MultiFab::dmap
const DistributionMapping & dmap() const
GLOBAL distribution (owner rank per box).
Definition multifab.hpp:58

pops::MultiFab::box_array
const BoxArray & box_array() const
GLOBAL decomposition of the level (all boxes, all ranks).
Definition multifab.hpp:56

pops::MultiFab::box
const Box2D & box(int li) const
VALID box of local fab li.
Definition multifab.hpp:71

pops::MultiFab::local_size
int local_size() const
Number of fabs OWNED by this rank (bound on local indices).
Definition multifab.hpp:65

pops::detail::PolarHasGeomSource
PolarHasGeomSource<M>: internal concept – true if M exposes polar_geom_source(u, r) -> State.
Definition polar_operator.hpp:85

pops::detail::PolarHasSource
PolarHasSource<M>: internal concept – true if M exposes source(U, aux) -> State.
Definition polar_operator.hpp:61

fab2d.hpp
Fab2D: single-grid data on a Box2D (in-house equivalent of AMReX's FArrayBox); Array4 / ConstArray4: ...

for_each.hpp
for_each_cell and reductions: the parallelism SEAM over the cells of a Box2D; sync_host / sync_device...

geometry.hpp
Geometry: index-space (Box2D) <-> Cartesian physical-space mapping; PolarGeometry: SIBLING for a glob...

multifab.hpp
MultiFab: a field DISTRIBUTED over a level (equivalent of AMReX's MultiFab).

pops::detail::polar_geom_source
POPS_HD Model::State polar_geom_source(const Model &m, const typename Model::State &u, Real r)
polar_geom_source<Model>: returns m.polar_geom_source(u, r) if PolarHasGeomSource<Model>,...
Definition polar_operator.hpp:93

pops::detail::polar_source
POPS_HD Model::State polar_source(const Model &m, const typename Model::State &u, const Aux &a)
polar_source<Model>: returns m.source(u, a) if PolarHasSource<Model>, otherwise the zero state.
Definition polar_operator.hpp:69

pops
Definition amr_hierarchy.hpp:29

pops::for_each_cell
void for_each_cell(const Box2D &b, F f)
Applies f to EACH cell (i, j) of box b (bounds inclusive), via Kokkos::parallel_for (Serial / OpenMP ...
Definition for_each.hpp:138

pops::Real
double Real
Definition types.hpp:30

pops::load_aux
POPS_HD Aux load_aux(const ConstArray4 &a, int i, int j)
load_aux<NComp>: reads NComp components of the auxiliary from an Array4 at (i,j).
Definition state_access.hpp:141

pops::assemble_rhs_polar
void assemble_rhs_polar(const Model &model, const MultiFab &U, const MultiFab &aux, const PolarGeometry &geom, MultiFab &R, bool recon_prim=false, bool wall_radial=false, Real pos_floor=Real(0))
assemble_rhs_polar<Limiter, NumericalFlux>: R = -div_polar F* + S on a PolarGeometry.
Definition polar_operator.hpp:230

pops::xface_box
Box2D xface_box(const Box2D &v)
xface_box / yface_box: face boxes normal to x (resp.
Definition face_flux.hpp:147

pops::yface_box
Box2D yface_box(const Box2D &v)
Definition face_flux.hpp:150

numerical_flux.hpp
Single-interface numerical flux policies: Rusanov, HLL, HLLC, Roe.

reconstruction.hpp
Interface reconstruction policies: MUSCL limiters and WENO5-Z.

spatial_operator.hpp
Cartesian spatial operator: assembles R(U, aux) = -div F + S over the cells of a level.

state.hpp
Pointwise types of the physics layer: StateVec<N> (conserved state) and Aux (auxiliary fields from th...

pops::Array4
WRITE POD handle (raw pointer + strides) over a Fab2D buffer, indexed by (i, j, c) IN GLOBAL INDICES ...
Definition fab2d.hpp:29

pops::Aux
POINTWISE auxiliary fields shared with the physics: single coupling channel.
Definition state.hpp:123

pops::Box2D
2D integer index space, cell-centered.
Definition box2d.hpp:37

pops::Box2D::hi
int hi[2]
Definition box2d.hpp:39

pops::Box2D::lo
int lo[2]
Definition box2d.hpp:38

pops::ConstArray4
READ-only handle (const counterpart of Array4): same layout and same contract (POD device-copyable,...
Definition fab2d.hpp:44

pops::PolarGeometry
Definition geometry.hpp:59

pops::PolarGeometry::dr
POPS_HD Real dr() const
Definition geometry.hpp:70

pops::PolarGeometry::r_min
Real r_min
Definition geometry.hpp:61

pops::PolarGeometry::domain
Box2D domain
nx() = nr (radial cells), ny() = ntheta (azimuthal cells)
Definition geometry.hpp:60

pops::PolarGeometry::dtheta
POPS_HD Real dtheta() const
Definition geometry.hpp:71

pops::detail::PolarAssembleRhsKernel
PolarAssembleRhsKernel: device kernel of the polar residual in cell (i,j).
Definition polar_operator.hpp:187

pops::detail::PolarAssembleRhsKernel::ax
ConstArray4 ax
Definition polar_operator.hpp:189

pops::detail::PolarAssembleRhsKernel::model
Model model
Definition polar_operator.hpp:188

pops::detail::PolarAssembleRhsKernel::r
Array4 r
Definition polar_operator.hpp:190

pops::detail::PolarAssembleRhsKernel::u
ConstArray4 u
Definition polar_operator.hpp:189

pops::detail::PolarAssembleRhsKernel::fr
ConstArray4 fr
Definition polar_operator.hpp:189

pops::detail::PolarAssembleRhsKernel::dtheta
Real dtheta
Definition polar_operator.hpp:191

pops::detail::PolarAssembleRhsKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_operator.hpp:192

pops::detail::PolarAssembleRhsKernel::dr
Real dr
Definition polar_operator.hpp:191

pops::detail::PolarAssembleRhsKernel::ft
ConstArray4 ft
Definition polar_operator.hpp:189

pops::detail::PolarAssembleRhsKernel::r_min
Real r_min
Definition polar_operator.hpp:191

pops::detail::PolarFaceFluxRKernel
Definition polar_operator.hpp:111

pops::detail::PolarFaceFluxRKernel::nflux
NumericalFlux nflux
Definition polar_operator.hpp:117

pops::detail::PolarFaceFluxRKernel::ax
ConstArray4 ax
Definition polar_operator.hpp:113

pops::detail::PolarFaceFluxRKernel::r_min
Real r_min
Definition polar_operator.hpp:115

pops::detail::PolarFaceFluxRKernel::u
ConstArray4 u
Definition polar_operator.hpp:113

pops::detail::PolarFaceFluxRKernel::fr
Array4 fr
Definition polar_operator.hpp:114

pops::detail::PolarFaceFluxRKernel::wall_radial
bool wall_radial
Definition polar_operator.hpp:125

pops::detail::PolarFaceFluxRKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_operator.hpp:130

pops::detail::PolarFaceFluxRKernel::model
Model model
Definition polar_operator.hpp:112

pops::detail::PolarFaceFluxRKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition polar_operator.hpp:128

pops::detail::PolarFaceFluxRKernel::lim
Limiter lim
Definition polar_operator.hpp:116

pops::detail::PolarFaceFluxRKernel::i_lo_face
int i_lo_face
Definition polar_operator.hpp:126

pops::detail::PolarFaceFluxRKernel::dr
Real dr
Definition polar_operator.hpp:115

pops::detail::PolarFaceFluxRKernel::i_hi_face
int i_hi_face
Definition polar_operator.hpp:127

pops::detail::PolarFaceFluxRKernel::recon_prim
bool recon_prim
Definition polar_operator.hpp:118

pops::detail::PolarFaceFluxRKernel::pos_comp
int pos_comp
component of the Density role (resolved by the host caller)
Definition polar_operator.hpp:129

pops::detail::PolarFaceFluxThetaKernel
PolarFaceFluxThetaKernel: device kernel of the flux at the azimuthal face j.
Definition polar_operator.hpp:158

pops::detail::PolarFaceFluxThetaKernel::pos_comp
int pos_comp
component of the Density role (resolved by the host caller)
Definition polar_operator.hpp:166

pops::detail::PolarFaceFluxThetaKernel::lim
Limiter lim
Definition polar_operator.hpp:162

pops::detail::PolarFaceFluxThetaKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_operator.hpp:167

pops::detail::PolarFaceFluxThetaKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition polar_operator.hpp:165

pops::detail::PolarFaceFluxThetaKernel::recon_prim
bool recon_prim
Definition polar_operator.hpp:164

pops::detail::PolarFaceFluxThetaKernel::model
Model model
Definition polar_operator.hpp:159

pops::detail::PolarFaceFluxThetaKernel::ax
ConstArray4 ax
Definition polar_operator.hpp:160

pops::detail::PolarFaceFluxThetaKernel::ft
Array4 ft
Definition polar_operator.hpp:161

pops::detail::PolarFaceFluxThetaKernel::nflux
NumericalFlux nflux
Definition polar_operator.hpp:163

pops::detail::PolarFaceFluxThetaKernel::u
ConstArray4 u
Definition polar_operator.hpp:160

types.hpp
Base scalar types and the POPS_HD macro (host+device portability).

POPS_HD
#define POPS_HD
Definition types.hpp:25