PoPS/doxygen/face__flux_8hpp_source.html

#pragma once


#include <pops/core/model/physical_model.hpp>  // HasPrimitiveVars: optional primitive reconstruction

#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/fv/reconstruction.hpp>

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

#include <pops/numerics/spatial/primitives/state_access.hpp>


#include <stdexcept>  // require_reconstruction_ghosts: state without the stencil width -> clear error


namespace pops {


template <class Model>


POPS_HD inline typename Model::State rusanov_flux(const Model& m, const typename Model::State& UL,

                                                 const Aux& AL, const typename Model::State& UR,

                                                 const Aux& AR, int dir) {

  return RusanovFlux{}(m, UL, AL, UR, AR, dir);

}


template <class Model, class Limiter>


POPS_HD inline typename Model::State reconstruct(const Model& model, const ConstArray4& u, int i,

                                                int j, int dir, Real sgn, const Limiter& lim,

                                                bool prim) {

  if constexpr (HasPrimitiveVars<Model> && Limiter::n_ghost >= 2) {

    if (prim) {  // convert the stencil U->P, limit on P, convert back P->U

      using Prim = typename Model::Prim;

      const Prim P0 = model.to_primitive(load_state<Model>(u, i, j));

      Prim Pf{};

      if constexpr (Limiter::n_ghost == 2) {

        const Prim Pm =

            model.to_primitive(load_state<Model>(u, dir == 0 ? i - 1 : i, dir == 0 ? j : j - 1));

        const Prim Pp =

            model.to_primitive(load_state<Model>(u, dir == 0 ? i + 1 : i, dir == 0 ? j : j + 1));

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

          Pf[c] = P0[c] + sgn * Real(0.5) * lim(P0[c] - Pm[c], Pp[c] - P0[c]);

      } else {  // WENO5 on the 5-point stencil in primitive variables

        const int d = (sgn > Real(0)) ? 1 : -1;

        const Prim Pm2 = model.to_primitive(

            load_state<Model>(u, dir == 0 ? i - 2 * d : i, dir == 0 ? j : j - 2 * d));

        const Prim Pm1 =

            model.to_primitive(load_state<Model>(u, dir == 0 ? i - d : i, dir == 0 ? j : j - d));

        const Prim Pp1 =

            model.to_primitive(load_state<Model>(u, dir == 0 ? i + d : i, dir == 0 ? j : j + d));

        const Prim Pp2 = model.to_primitive(

            load_state<Model>(u, dir == 0 ? i + 2 * d : i, dir == 0 ? j : j + 2 * d));

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

          Pf[c] = weno5z(Pm2[c], Pm1[c], P0[c], Pp1[c], Pp2[c]);

      }

      return model.to_conservative(Pf);

    }

  }

  (void)model;

  (void)prim;

  typename Model::State s = load_state<Model>(u, i, j);

  if constexpr (Limiter::n_ghost == 2) {

    // MUSCL: per-component limited slope (order 2).

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

      const Real am = (dir == 0) ? u(i, j, c) - u(i - 1, j, c) : u(i, j, c) - u(i, j - 1, c);

      const Real ap = (dir == 0) ? u(i + 1, j, c) - u(i, j, c) : u(i, j + 1, c) - u(i, j, c);

      s[c] += sgn * Real(0.5) * lim(am, ap);

    }

  } else if constexpr (Limiter::n_ghost >= 3) {

    // WENO5 (order 5): face value from a 5-point stencil oriented by sgn

    // (sgn>0 -> +dir face; sgn<0 -> -dir face, reversed stencil). lim unused.

    (void)lim;

    const int d = (sgn > Real(0)) ? 1 : -1;

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

      if (dir == 0)

        s[c] = weno5z(u(i - 2 * d, j, c), u(i - d, j, c), u(i, j, c), u(i + d, j, c),

                      u(i + 2 * d, j, c));

      else

        s[c] = weno5z(u(i, j - 2 * d, c), u(i, j - d, c), u(i, j, c), u(i, j + d, c),

                      u(i, j + 2 * d, c));

    }

  }

  return s;

}


template <class Model, class Limiter>


POPS_HD inline typename Model::State reconstruct_pp(const Model& model, const ConstArray4& u, int i,

                                                   int j, int dir, Real sgn, const Limiter& lim,

                                                   bool prim, Real pos_floor, int pos_comp) {

  typename Model::State s = reconstruct<Model>(model, u, i, j, dir, sgn, lim, prim);

  zhang_shu_scale<Model>(s, u, i, j, pos_floor, pos_comp);

  return s;

}


namespace detail {

template <class Limiter>


inline void require_reconstruction_ghosts(const MultiFab& U) {

  if (U.n_grow() < Limiter::n_ghost)

    throw std::runtime_error(

        "spatial operator: the state must carry at least Limiter::n_ghost ghost layers "

        "(the reconstruction stencil reads i+-Limiter::n_ghost at the edge of the valid box); "

        "allocate the state MultiFab with this number of ghosts.");

}


}  // namespace detail


inline Box2D xface_box(const Box2D& v) {

  return Box2D{{v.lo[0], v.lo[1]}, {v.hi[0] + 1, v.hi[1]}};

}


inline Box2D yface_box(const Box2D& v) {

  return Box2D{{v.lo[0], v.lo[1]}, {v.hi[0], v.hi[1] + 1}};

}


namespace detail {

template <class Limiter, class NumericalFlux, class Model>


struct FaceFluxXKernel {

  Model model;

  ConstArray4 u, ax;

  Array4 fx;

  Real dx;

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

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

    if constexpr (DiffusiveModel<Model>) {

      const Real nu = model.diffusivity();

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

        fx(i, j, c) += -nu * (u(i, j, c) - u(i - 1, j, c)) / dx;

    }

  }


};


template <class Limiter, class NumericalFlux, class Model>


struct FaceFluxYKernel {

  Model model;

  ConstArray4 u, ax;

  Array4 fy;

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

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

    if constexpr (DiffusiveModel<Model>) {

      const Real nu = model.diffusivity();

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

        fy(i, j, c) += -nu * (u(i, j, c) - u(i, j - 1, c)) / dy;

    }

  }


};


}  // namespace detail


//

// compute_face_fluxes: writes the numerical fluxes at the FACES (Fx at faces normal to x,

// Fy at y), BEFORE divergence. This is the brick the AMR reflux needs (it accumulates the

// fine fluxes and subtracts the coarse flux at the coarse-fine interfaces; assemble_rhs

// itself computes -div F directly and discards the face fluxes).

//

// Conventions: Fx(i,j) = flux at the face between cells (i-1,j) and (i,j), i in [lo..hi+1].

// Fy(i,j) = flux between (i,j-1) and (i,j), j in [lo..hi+1]. Same reconstruction (Limiter)

// and numerical flux (NumericalFlux) as assemble_rhs, so

//   r(i,j) = S - (Fx(i+1,j)-Fx(i,j))/dx - (Fy(i,j+1)-Fy(i,j))/dy

// gives back EXACTLY the assemble_rhs residual. Fx, Fy sized by the caller (xface_box/yface_box

// boxes, ncomp = Model::n_vars, 0 ghost). Device-callable.

//

// DIFFUSION on AMR (milestone 4): for a DiffusiveModel, we add the FACE Fickian flux

// F_diff = -nu (u_R - u_L)/h (centered gradient at the face, cell values). Its divergence

// -(Fx(i+1)-Fx(i))/dx gives back EXACTLY +nu Lap(u) of assemble_rhs, but treated as a FLUX:

// the AMR reflux therefore sees it, and the diffusion stays conservative at the coarse-fine

// interfaces (otherwise a direct Laplacian would be ignored by the reflux). dx/dy = step of

// the LEVEL (passed by the caller; 0 by default, not read for a non-diffusive model -> the

// hyperbolic path is strictly bit-identical).

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


void compute_face_fluxes(const Model& model, const MultiFab& U, const MultiFab& aux, MultiFab& Fx,

                         MultiFab& Fy, Real dx = 0, Real dy = 0, bool recon_prim = false,

                         Real pos_floor = Real(0)) {

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

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

    Array4 fy = Fy.fab(li).array();

    const Box2D v = U.box(li);

    for_each_cell(xface_box(v),

                  detail::FaceFluxXKernel<Limiter, NumericalFlux, Model>{

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

    for_each_cell(yface_box(v),

                  detail::FaceFluxYKernel<Limiter, NumericalFlux, Model>{

                      model, u, ax, fy, 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::n_grow
int n_grow() const
Number of ghost layers.
Definition multifab.hpp:62

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

pops::HasPrimitiveVars
OPTIONAL extension of a PhysicalModel: primitive variables + cons<->prim conversions.
Definition physical_model.hpp:167

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::require_reconstruction_ghosts
void require_reconstruction_ghosts(const MultiFab &U)
require_reconstruction_ghosts<Limiter>: STRUCTURAL ENTRY GUARD of the FV spatial operators.
Definition face_flux.hpp:133

pops
Definition amr_hierarchy.hpp:29

pops::rusanov_flux
POPS_HD Model::State rusanov_flux(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir)
rusanov_flux: free compat, delegates to RusanovFlux{} (policy of numerical_flux.hpp).
Definition face_flux.hpp:37

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::compute_face_fluxes
void compute_face_fluxes(const Model &model, const MultiFab &U, const MultiFab &aux, MultiFab &Fx, MultiFab &Fy, Real dx=0, Real dy=0, bool recon_prim=false, Real pos_floor=Real(0))
compute_face_fluxes<Limiter,NumericalFlux>: writes the face fluxes BEFORE divergence.
Definition face_flux.hpp:249

pops::reconstruct_pp
POPS_HD Model::State reconstruct_pp(const Model &model, const ConstArray4 &u, int i, int j, int dir, Real sgn, const Limiter &lim, bool prim, Real pos_floor, int pos_comp)
reconstruct_pp: reconstruct + zhang_shu_scale positivity limiter on the returned state.
Definition face_flux.hpp:114

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::weno5z
POPS_HD Real weno5z(Real vm2, Real vm1, Real v0, Real vp1, Real vp2)
weno5z: WENO5-Z reconstruction (Borges 2008) at one interface, on a 5-point stencil.
Definition reconstruction.hpp:70

pops::reconstruct
POPS_HD Model::State reconstruct(const Model &model, const ConstArray4 &u, int i, int j, int dir, Real sgn, const Limiter &lim, bool prim)
reconstruct<Model,Limiter>: face value at (i,j) extrapolated in direction dir.
Definition face_flux.hpp:51

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.

physical_model.hpp
C++20 concepts defining the contract of the physics layer.

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

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

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::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::RusanovFlux
RusanovFlux (local Lax-Friedrichs): robust flux, compatible with any minimal PhysicalModel.
Definition numerical_flux.hpp:55

pops::detail::FaceFluxXKernel
FaceFluxXKernel: device kernel for the flux at the radial x face (between i-1 and i).
Definition face_flux.hpp:161

pops::detail::FaceFluxXKernel::recon_prim
bool recon_prim
Definition face_flux.hpp:168

pops::detail::FaceFluxXKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition face_flux.hpp:169

pops::detail::FaceFluxXKernel::u
ConstArray4 u
Definition face_flux.hpp:163

pops::detail::FaceFluxXKernel::ax
ConstArray4 ax
Definition face_flux.hpp:163

pops::detail::FaceFluxXKernel::model
Model model
Definition face_flux.hpp:162

pops::detail::FaceFluxXKernel::lim
Limiter lim
Definition face_flux.hpp:166

pops::detail::FaceFluxXKernel::fx
Array4 fx
Definition face_flux.hpp:164

pops::detail::FaceFluxXKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition face_flux.hpp:171

pops::detail::FaceFluxXKernel::dx
Real dx
Definition face_flux.hpp:165

pops::detail::FaceFluxXKernel::pos_comp
int pos_comp
component of the Density role (resolved by the host caller)
Definition face_flux.hpp:170

pops::detail::FaceFluxXKernel::nflux
NumericalFlux nflux
Definition face_flux.hpp:167

pops::detail::FaceFluxYKernel
FaceFluxYKernel: device kernel for the flux at the y face (between j-1 and j).
Definition face_flux.hpp:191

pops::detail::FaceFluxYKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition face_flux.hpp:201

pops::detail::FaceFluxYKernel::lim
Limiter lim
Definition face_flux.hpp:196

pops::detail::FaceFluxYKernel::pos_comp
int pos_comp
component of the Density role (resolved by the host caller)
Definition face_flux.hpp:200

pops::detail::FaceFluxYKernel::recon_prim
bool recon_prim
Definition face_flux.hpp:198

pops::detail::FaceFluxYKernel::model
Model model
Definition face_flux.hpp:192

pops::detail::FaceFluxYKernel::ax
ConstArray4 ax
Definition face_flux.hpp:193

pops::detail::FaceFluxYKernel::u
ConstArray4 u
Definition face_flux.hpp:193

pops::detail::FaceFluxYKernel::fy
Array4 fy
Definition face_flux.hpp:194

pops::detail::FaceFluxYKernel::pos_floor
Real pos_floor
Zhang-Shu positivity limiter (<= 0: inactive, bit-identical)
Definition face_flux.hpp:199

pops::detail::FaceFluxYKernel::nflux
NumericalFlux nflux
Definition face_flux.hpp:197

pops::detail::FaceFluxYKernel::dy
Real dy
Definition face_flux.hpp:195

POPS_HD
#define POPS_HD
Definition types.hpp:25