PoPS/doxygen/masked__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


namespace pops {


// ============================================================================

// DOMAIN MASK (T2 effort, conservative, OPT-IN -- default path untouched)

// ============================================================================

// The mask makes the FV transport aware of an ACTIVE sub-domain (e.g. a bounded disk-shaped region).

// Convention: mask(i, j) >= 0.5 -> ACTIVE cell, otherwise INACTIVE. A face is OPEN (normal flux

// computed) if BOTH adjacent cells are active; it is CLOSED (normal flux set to ZERO) if at least

// one is inactive. Zeroing the normal flux at active/inactive faces makes the step CONSERVATIVE

// over the active sub-domain: no mass crosses the boundary, so the total mass over the active cells

// is conserved to machine precision (telescoping internal fluxes, zero boundary fluxes). This is the

// FV counterpart of the conducting wall (which only acts on the elliptic part).

//

// The residual is written ONLY on the active cells; an inactive cell keeps its residual at 0

// (the caller does not advance it). This header does NOT wire this path into System::step: it

// provides the mask-aware brick, exercised directly by the tests and, eventually, behind the

// active-sub-domain opt-in.


namespace detail {


POPS_HD inline bool mask_active(const ConstArray4& mask, int i, int j) {

  return mask(i, j, 0) >= Real(0.5);

}


template <class Limiter, class NumericalFlux, class Model>


struct AssembleRhsMaskedKernel {

  Model model;

  ConstArray4 u, ax, mask;

  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 {

    if (!mask_active(mask, i,

                     j)) {  // cell outside the active sub-domain: zero residual, not advanced

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

        r(i, j, c) = Real(0);

      return;

    }

    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 on either side, numerical flux, THEN mask gate (closed face

    // -> zero normal flux) -- an inactive neighbor cell closes the face between it and (i, j).

    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);

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

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

    if (!mask_active(mask, i - 1, j))

      Fxm = typename Model::State{};

    if (!mask_active(mask, i + 1, j))

      Fxp = typename Model::State{};


    // 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);

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

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

    if (!mask_active(mask, i, j - 1))

      Fym = typename Model::State{};

    if (!mask_active(mask, i, j + 1))

      Fyp = typename Model::State{};


    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;

  }


};


}  // namespace detail


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


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

                         const MultiFab& mask, 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();

    const ConstArray4 mk = mask.fab(li).const_array();

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

    const Box2D v = R.box(li);

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

                         model, u, ax, mk, r, dx, dy, lim, nflux, 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

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::detail::mask_active
POPS_HD bool mask_active(const ConstArray4 &mask, int i, int j)
Activity indicator of a cell from a 0/1 cell-centered mask (>= 0.5 -> active).
Definition masked_operator.hpp:44

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::assemble_rhs_masked
void assemble_rhs_masked(const Model &model, const MultiFab &U, const MultiFab &aux, const MultiFab &mask, const Geometry &geom, MultiFab &R, bool recon_prim=false, Real pos_floor=Real(0))
assemble_rhs_masked<Limiter,NumericalFlux>: residual R = -div Fhat + S RESTRICTED to a 0/1 cell-cente...
Definition masked_operator.hpp:133

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::AssembleRhsMaskedKernel
AssembleRhsMaskedKernel: variant of AssembleRhsKernel AWARE of a domain mask.
Definition masked_operator.hpp:59

pops::detail::AssembleRhsMaskedKernel::nflux
NumericalFlux nflux
Definition masked_operator.hpp:65

pops::detail::AssembleRhsMaskedKernel::recon_prim
bool recon_prim
Definition masked_operator.hpp:66

pops::detail::AssembleRhsMaskedKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition masked_operator.hpp:69

pops::detail::AssembleRhsMaskedKernel::u
ConstArray4 u
Definition masked_operator.hpp:61

pops::detail::AssembleRhsMaskedKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition masked_operator.hpp:67

pops::detail::AssembleRhsMaskedKernel::r
Array4 r
Definition masked_operator.hpp:62

pops::detail::AssembleRhsMaskedKernel::pos_comp
int pos_comp
component of the Density role (resolved by the host caller)
Definition masked_operator.hpp:68

pops::detail::AssembleRhsMaskedKernel::ax
ConstArray4 ax
Definition masked_operator.hpp:61

pops::detail::AssembleRhsMaskedKernel::model
Model model
Definition masked_operator.hpp:60

pops::detail::AssembleRhsMaskedKernel::dx
Real dx
Definition masked_operator.hpp:63

pops::detail::AssembleRhsMaskedKernel::mask
ConstArray4 mask
Definition masked_operator.hpp:61

pops::detail::AssembleRhsMaskedKernel::lim
Limiter lim
Definition masked_operator.hpp:64

pops::detail::AssembleRhsMaskedKernel::dy
Real dy
Definition masked_operator.hpp:63

POPS_HD
#define POPS_HD
Definition types.hpp:25