PoPS/doxygen/cartesian__operator_8hpp_source.html

#pragma once


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

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

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

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

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

#include <pops/numerics/spatial/primitives/face_flux.hpp>     // reconstruct_pp, require_reconstruction_ghosts

#include <pops/numerics/spatial/primitives/positivity.hpp>    // detail::positivity_comp

#include <pops/numerics/spatial/primitives/state_access.hpp>  // load_state, load_aux, DiffusiveModel

#include <pops/numerics/spatial/primitives/wave_speed.hpp>    // fill_wave_speed_cache


namespace pops {


namespace detail {

//

// nvcc does not reliably emit the device kernel of a Model-template extended lambda first

// instantiated from an EXTERNAL TU through the std::function / host-lambda nesting of block_builder:

// the test passes on Serial and under compute-sanitizer but segfaults at runtime on Cuda (Heisenbug).

// A device-callable class does not have these instantiation-context restrictions. Body IDENTICAL to

// the former lambda -> residual BIT-IDENTICAL to add_block on CPU (and, targeted, on device).

template <class Limiter, class NumericalFlux, class Model>


struct AssembleRhsKernel {

  Model model;

  ConstArray4 u, ax;

  Array4 r;

  Real dx, dy;

  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 Aux Ac = load_aux<aux_comps<Model>()>(ax, i, j);

    const Aux Axm = load_aux<aux_comps<Model>()>(ax, i - 1, j);

    const Aux Axp = load_aux<aux_comps<Model>()>(ax, i + 1, j);

    const Aux Aym = load_aux<aux_comps<Model>()>(ax, i, j - 1);

    const Aux Ayp = load_aux<aux_comps<Model>()>(ax, i, j + 1);


    // x faces: reconstruction of the states on either side of each face

    const auto Lxm =

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

    const auto Rxm =

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

    const auto Lxp =

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

    const auto Rxp =

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

    const auto Fxm = nflux(model, Lxm, Axm, Rxm, Ac, 0);

    const auto Fxp = nflux(model, Lxp, Ac, Rxp, Axp, 0);


    // y faces

    const auto Lym =

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

    const auto Rym =

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

    const auto Lyp =

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

    const auto Ryp =

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

    const auto Fym = nflux(model, Lym, Aym, Rym, Ac, 1);

    const auto Fyp = nflux(model, Lyp, Ac, Ryp, Ayp, 1);


    const auto S = model.source(load_state<Model>(u, i, j), Ac);

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

      r(i, j, c) = S[c] - (Fxp[c] - Fxm[c]) / dx - (Fyp[c] - Fym[c]) / dy;


    // Parabolic (Fickian) term: +nu Lap(U), 5-point centered differences.

    // Guarded by DiffusiveModel: no effect (nor codegen) for a non-diffusive model.

    if constexpr (DiffusiveModel<Model>) {

      const Real nu = model.diffusivity();

      const Real idx2 = Real(1) / (dx * dx), idy2 = Real(1) / (dy * dy);

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

        r(i, j, c) += nu * ((u(i + 1, j, c) - 2 * u(i, j, c) + u(i - 1, j, c)) * idx2 +

                            (u(i, j + 1, c) - 2 * u(i, j, c) + u(i, j - 1, c)) * idy2);

    }

  }


};


}  // namespace detail


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


void assemble_rhs(const Model& model, const MultiFab& U, const MultiFab& aux, const Geometry& geom,

                  MultiFab& R, bool recon_prim = false, Real pos_floor = Real(0)) {

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

  const Real dx = geom.dx(), dy = geom.dy();

  const Limiter lim{};

  const NumericalFlux nflux{};

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

  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 r = R.fab(li).array();

    const Box2D v = R.box(li);

    for_each_cell(v, detail::AssembleRhsKernel<Limiter, NumericalFlux, Model>{

                         model, u, ax, r, dx, dy, lim, nflux, recon_prim, pos_floor, pos_comp});

  }

}


namespace detail {

template <class Limiter, class Model>


struct AssembleRhsHllCachedKernel {

  Model model;

  ConstArray4 u, ax, ws;

  Array4 r;

  Real dx, dy;

  Limiter lim;

  bool recon_prim;

  Real pos_floor = Real(0);

  int pos_comp = 0;


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

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

    const Aux Axm = load_aux<aux_comps<Model>()>(ax, i - 1, j);

    const Aux Axp = load_aux<aux_comps<Model>()>(ax, i + 1, j);

    const Aux Aym = load_aux<aux_comps<Model>()>(ax, i, j - 1);

    const Aux Ayp = load_aux<aux_comps<Model>()>(ax, i, j + 1);


    // x faces: reconstruction of the states on both sides of each face

    const auto Lxm =

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

    const auto Rxm =

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

    const auto Lxp =

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

    const auto Rxp =

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

    // face signal speeds = union (min of the lo, max of the hi) of the two adjacent cells: the LEFT

    // cell is always the lower-index neighbor (same operands/order as hll_speeds).

    const Real sLxm = ws(i - 1, j, 0) < ws(i, j, 0) ? ws(i - 1, j, 0) : ws(i, j, 0);

    const Real sRxm = ws(i - 1, j, 1) > ws(i, j, 1) ? ws(i - 1, j, 1) : ws(i, j, 1);

    const Real sLxp = ws(i, j, 0) < ws(i + 1, j, 0) ? ws(i, j, 0) : ws(i + 1, j, 0);

    const Real sRxp = ws(i, j, 1) > ws(i + 1, j, 1) ? ws(i, j, 1) : ws(i + 1, j, 1);

    const auto Fxm = hll_flux_with_speeds(model, Lxm, Axm, Rxm, Ac, 0, sLxm, sRxm);

    const auto Fxp = hll_flux_with_speeds(model, Lxp, Ac, Rxp, Axp, 0, sLxp, sRxp);


    // y faces (components 2 = lo_y, 3 = hi_y of the scratch)

    const auto Lym =

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

    const auto Rym =

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

    const auto Lyp =

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

    const auto Ryp =

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

    const Real sLym = ws(i, j - 1, 2) < ws(i, j, 2) ? ws(i, j - 1, 2) : ws(i, j, 2);

    const Real sRym = ws(i, j - 1, 3) > ws(i, j, 3) ? ws(i, j - 1, 3) : ws(i, j, 3);

    const Real sLyp = ws(i, j, 2) < ws(i, j + 1, 2) ? ws(i, j, 2) : ws(i, j + 1, 2);

    const Real sRyp = ws(i, j, 3) > ws(i, j + 1, 3) ? ws(i, j, 3) : ws(i, j + 1, 3);

    const auto Fym = hll_flux_with_speeds(model, Lym, Aym, Rym, Ac, 1, sLym, sRym);

    const auto Fyp = hll_flux_with_speeds(model, Lyp, Ac, Ryp, Ayp, 1, sLyp, sRyp);


    const auto S = model.source(load_state<Model>(u, i, j), Ac);

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

      r(i, j, c) = S[c] - (Fxp[c] - Fxm[c]) / dx - (Fyp[c] - Fym[c]) / dy;


    // Parabolic (Fickian) term: identical to AssembleRhsKernel, guarded by DiffusiveModel.

    if constexpr (DiffusiveModel<Model>) {

      const Real nu = model.diffusivity();

      const Real idx2 = Real(1) / (dx * dx), idy2 = Real(1) / (dy * dy);

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

        r(i, j, c) += nu * ((u(i + 1, j, c) - 2 * u(i, j, c) + u(i - 1, j, c)) * idx2 +

                            (u(i, j + 1, c) - 2 * u(i, j, c) + u(i, j - 1, c)) * idy2);

    }

  }


};


}  // namespace detail


template <class Limiter = NoSlope, class Model>


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

                             const Geometry& geom, MultiFab& R, MultiFab& cache,

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

  const Real dx = geom.dx(), dy = geom.dy();

  const Limiter lim{};

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

  fill_wave_speed_cache(model, U, aux, cache);

  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 ws = cache.fab(li).const_array();

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

    const Box2D v = R.box(li);

    for_each_cell(v, detail::AssembleRhsHllCachedKernel<Limiter, Model>{

                         model, u, ax, ws, r, dx, dy, lim, recon_prim, pos_floor, pos_comp});

  }

}


}  // namespace pops

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::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::DiffusiveModel
DiffusiveModel: optional concept for models with isotropic scalar diffusion.
Definition state_access.hpp:38

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

face_flux.hpp
Face-state reconstruction and face fluxes of the Cartesian spatial operator.

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
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::assemble_rhs
void assemble_rhs(const Model &model, const MultiFab &U, const MultiFab &aux, const Geometry &geom, MultiFab &R, bool recon_prim=false, Real pos_floor=Real(0))
assemble_rhs<Limiter,NumericalFlux>: residual R = -div Fhat + S over all boxes.
Definition cartesian_operator.hpp:107

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

pops::assemble_rhs_hll_cached
void assemble_rhs_hll_cached(const Model &model, const MultiFab &U, const MultiFab &aux, const Geometry &geom, MultiFab &R, MultiFab &cache, bool recon_prim=false, Real pos_floor=Real(0))
assemble_rhs_hll_cached<Limiter>: residual R = -div Fhat + S at the HLL flux, wave speeds PRE-COMPUTE...
Definition cartesian_operator.hpp:205

pops::fill_wave_speed_cache
void fill_wave_speed_cache(const Model &model, const MultiFab &U, const MultiFab &aux, MultiFab &cache)
fill_wave_speed_cache: fills the per-cell wave speed scratch (lo_x, hi_x, lo_y, hi_y).
Definition wave_speed.hpp:277

pops::hll_flux_with_speeds
POPS_HD Model::State hll_flux_with_speeds(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir, Real sL, Real sR)
hll_flux_with_speeds: HLL flux from ALREADY estimated signal speeds (sL, sR).
Definition numerical_flux.hpp:96

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

positivity.hpp
Zhang-Shu positivity limiter and Density-role resolution.

state_access.hpp
Model/state/aux access layer of the Cartesian spatial operator.

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::ConstArray4
READ-only handle (const counterpart of Array4): same layout and same contract (POD device-copyable,...
Definition fab2d.hpp:44

pops::Geometry
Cartesian geometry of a level: index domain + physical bounds [xlo, xhi] x [ylo, yhi].
Definition geometry.hpp:20

pops::Geometry::dy
POPS_HD Real dy() const
Grid spacing in y (= (yhi - ylo) / domain.ny()). POPS_HD.
Definition geometry.hpp:33

pops::Geometry::dx
POPS_HD Real dx() const
Grid spacing in x (= (xhi - xlo) / domain.nx()). POPS_HD.
Definition geometry.hpp:31

pops::detail::AssembleRhsHllCachedKernel
AssembleRhsHllCachedKernel: kernel of the residual R = -div Fhat + S for the HLL flux with wave speed...
Definition cartesian_operator.hpp:132

pops::detail::AssembleRhsHllCachedKernel::ax
ConstArray4 ax
Definition cartesian_operator.hpp:134

pops::detail::AssembleRhsHllCachedKernel::lim
Limiter lim
Definition cartesian_operator.hpp:137

pops::detail::AssembleRhsHllCachedKernel::recon_prim
bool recon_prim
Definition cartesian_operator.hpp:138

pops::detail::AssembleRhsHllCachedKernel::u
ConstArray4 u
Definition cartesian_operator.hpp:134

pops::detail::AssembleRhsHllCachedKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition cartesian_operator.hpp:141

pops::detail::AssembleRhsHllCachedKernel::ws
ConstArray4 ws
Definition cartesian_operator.hpp:134

pops::detail::AssembleRhsHllCachedKernel::model
Model model
Definition cartesian_operator.hpp:133

pops::detail::AssembleRhsHllCachedKernel::dx
Real dx
Definition cartesian_operator.hpp:136

pops::detail::AssembleRhsHllCachedKernel::dy
Real dy
Definition cartesian_operator.hpp:136

pops::detail::AssembleRhsHllCachedKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition cartesian_operator.hpp:139

pops::detail::AssembleRhsHllCachedKernel::r
Array4 r
Definition cartesian_operator.hpp:135

pops::detail::AssembleRhsHllCachedKernel::pos_comp
int pos_comp
Density role component (resolved by the host caller)
Definition cartesian_operator.hpp:140

pops::detail::AssembleRhsKernel
AssembleRhsKernel<Limiter,NumericalFlux,Model>: device kernel of the central residual of assemble_rhs...
Definition cartesian_operator.hpp:41

pops::detail::AssembleRhsKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition cartesian_operator.hpp:49

pops::detail::AssembleRhsKernel::dx
Real dx
Definition cartesian_operator.hpp:45

pops::detail::AssembleRhsKernel::ax
ConstArray4 ax
Definition cartesian_operator.hpp:43

pops::detail::AssembleRhsKernel::nflux
NumericalFlux nflux
Definition cartesian_operator.hpp:47

pops::detail::AssembleRhsKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition cartesian_operator.hpp:51

pops::detail::AssembleRhsKernel::model
Model model
Definition cartesian_operator.hpp:42

pops::detail::AssembleRhsKernel::recon_prim
bool recon_prim
Definition cartesian_operator.hpp:48

pops::detail::AssembleRhsKernel::pos_comp
int pos_comp
component of the Density role (resolved by the host caller)
Definition cartesian_operator.hpp:50

pops::detail::AssembleRhsKernel::u
ConstArray4 u
Definition cartesian_operator.hpp:43

pops::detail::AssembleRhsKernel::r
Array4 r
Definition cartesian_operator.hpp:44

pops::detail::AssembleRhsKernel::lim
Limiter lim
Definition cartesian_operator.hpp:46

pops::detail::AssembleRhsKernel::dy
Real dy
Definition cartesian_operator.hpp:45

POPS_HD
#define POPS_HD
Definition types.hpp:25

wave_speed.hpp
Global wave-speed / step-bound reductions and the HLL wave-speed cache.