PoPS/doxygen/amr__reflux_8hpp_source.html

#pragma once


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

#include <pops/amr/hierarchy/refinement_ratio.hpp>

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

#include <pops/mesh/index/box2d.hpp>

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

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

#include <pops/numerics/spatial_operator.hpp>


#include <vector>


namespace pops {


static_assert(kAmrRefRatio == 2, "ratio-2-structural kernels below assume kAmrRefRatio == 2");


// xface_box / yface_box: provided by numerics/spatial_operator.hpp (included above),

// same face-box conventions. We do not redefine them here.


namespace detail {


// NAMED FUNCTORS (and not POPS_HD lambdas) for the AMR reflux device kernels. Same reasons as the

// rest of the mesh/time path (AmrSspRhsKernel, fill_boundary recipe): these kernels are first-

// instantiated from an external loader TU, where an extended lambda trips the device kernel emission

// under nvcc. Body strictly identical to the previous lambdas -> bit-identical CPU and device.


template <class Model>


struct RusanovFaceXKernel {

  Model m;

  ConstArray4 u, ax;

  Array4 F;


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

    typename Model::State UL{}, UR{};

    UL[0] = u(i - 1, j);

    UR[0] = u(i, j);

    F(i, j) = rusanov_flux(m, UL, load_aux<aux_comps<Model>()>(ax, i - 1, j), UR,

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

  }


};


template <class Model>


struct RusanovFaceYKernel {

  Model m;

  ConstArray4 u, ax;

  Array4 F;


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

    typename Model::State UL{}, UR{};

    UL[0] = u(i, j - 1);

    UR[0] = u(i, j);

    F(i, j) = rusanov_flux(m, UL, load_aux<aux_comps<Model>()>(ax, i, j - 1), UR,

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

  }


};


struct AdvanceFab1cKernel {

  Array4 uu;

  ConstArray4 FX, FY;

  double dx, dy, dt;


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

    uu(i, j) -= dt * ((FX(i + 1, j) - FX(i, j)) / dx + (FY(i, j + 1) - FY(i, j)) / dy);

  }


};


}  // namespace detail


// First-order Rusanov flux, 1 component, aux variable in space (Fab2D with

// 3 components [phi, gx, gy], ghosts filled), on a Fab2D.

// fx(i,j) = flux at the left face of cell i; fy(i,j) = bottom face of j.

template <class Model>


void compute_fluxes_1c(const Model& m, const Fab2D& U, const Fab2D& aux, Fab2D& fx, Fab2D& fy) {

  const ConstArray4 u = U.const_array();

  const ConstArray4 ax = aux.const_array();

  {

    Array4 F = fx.array();

    for_each_cell(fx.box(), detail::RusanovFaceXKernel<Model>{m, u, ax, F});

  }

  {

    Array4 F = fy.array();

    for_each_cell(fy.box(), detail::RusanovFaceYKernel<Model>{m, u, ax, F});

  }

}


// Explicit Euler: U -= dt div(F). The ghosts of U must be filled.

template <class Model>


void advance_fab_1c(const Model& m, Fab2D& U, const Fab2D& aux, double dx, double dy, double dt,

                    Fab2D& fx, Fab2D& fy) {

  compute_fluxes_1c(m, U, aux, fx, fy);

  Array4 uu = U.array();

  const ConstArray4 FX = fx.const_array();

  const ConstArray4 FY = fy.const_array();

  for_each_cell(U.box(), detail::AdvanceFab1cKernel{uu, FX, FY, dx, dy, dt});

}


// Periodic ghosts for a single Fab2D covering the domain.


inline void fill_periodic_fab(Fab2D& U, const Box2D& dom) {

  const int ng = U.n_ghost();

  const int nx = dom.nx(), ny = dom.ny();

  Array4 a = U.array();

  for (int j = dom.lo[1]; j <= dom.hi[1]; ++j)

    for (int g = 1; g <= ng; ++g) {

      a(dom.lo[0] - g, j) = a(dom.hi[0] - g + 1, j);

      a(dom.hi[0] + g, j) = a(dom.lo[0] + g - 1, j);

    }

  for (int i = dom.lo[0] - ng; i <= dom.hi[0] + ng; ++i)

    for (int g = 1; g <= ng; ++g) {

      a(i, dom.lo[1] - g) = a(i, dom.hi[1] - g + 1);

      a(i, dom.hi[1] + g) = a(i, dom.lo[1] + g - 1);

    }

}


// Fine ghosts by injection from the coarse (ratio 2), interpolated in time

// between the old coarse state (frac=0) and the new one (frac=1): space-time

// FillPatch for the Berger-Oliger subcycling.


inline void fill_fine_ghosts_t(Fab2D& Uf, const Fab2D& Uco, const Fab2D& Ucn, double frac) {

  const ConstArray4 co = Uco.const_array();

  const ConstArray4 cn = Ucn.const_array();

  Array4 f = Uf.array();

  const Box2D g = Uf.grown_box();

  const Box2D v = Uf.box();

  auto coarsen = [](int x) { return (x >= 0) ? x / 2 : -((-x + 1) / 2); };

  for (int j = g.lo[1]; j <= g.hi[1]; ++j)

    for (int i = g.lo[0]; i <= g.hi[0]; ++i)

      if (!v.contains(i, j)) {

        const int ci = coarsen(i), cj = coarsen(j);

        f(i, j) = (1 - frac) * co(ci, cj) + frac * cn(ci, cj);

      }

}


// Fine -> coarse average over the covered region (ratio 2).


inline void average_down_fab(const Fab2D& Uf, Fab2D& Uc, int CI0, int CI1, int CJ0, int CJ1) {

  const ConstArray4 f = Uf.const_array();

  Array4 c = Uc.array();

  for (int J = CJ0; J <= CJ1; ++J)

    for (int I = CI0; I <= CI1; ++I)

      c(I, J) = 0.25 * (f(2 * I, 2 * J) + f(2 * I + 1, 2 * J) + f(2 * I, 2 * J + 1) +

                        f(2 * I + 1, 2 * J + 1));

}


// One conservative 2-level step with Berger-Oliger subcycling (the fine does

// r=2 substeps of dt/2) and reflux. Fine region = coarse cells

// [CI0..CI1] x [CJ0..CJ1] (strictly interior), refined into

// [2CI0..2CI1+1] x [...]. dom = periodic coarse domain.

//

// Flux register: we accumulate the coarse flux (x dt) and the sum of the fine

// fluxes (x dt/2 per substep, spatially averaged) at the 4 faces of the fine

// region, then we correct the adjacent coarse cells by their difference.

template <class Model>


void amr_step_2level(const Model& m, Fab2D& Uc, const Box2D& dom, double dxc, double dyc, Fab2D& Uf,

                     int CI0, int CI1, int CJ0, int CJ1, const Fab2D& auxc, const Fab2D& auxf,

                     double dt) {

  const int r = kAmrRefRatio;

  const double dxf = dxc / kAmrRefRatio, dyf = dyc / kAmrRefRatio, dtf = dt / r;

  const int nJ = CJ1 - CJ0 + 1, nI = CI1 - CI0 + 1;


  const Fab2D Uc_old = Uc;  // coarse state at time t (for temporal interp)


  // --- coarse fluxes (before update) at the 4 faces of the fine region ---

  fill_periodic_fab(Uc, dom);

  Fab2D fxc(xface_box(Uc.box()), 1, 0), fyc(yface_box(Uc.box()), 1, 0);

  compute_fluxes_1c(m, Uc, auxc, fxc, fyc);

  const ConstArray4 FXc = fxc.const_array();

  const ConstArray4 FYc = fyc.const_array();

  std::vector<double> cL(nJ), cR(nJ), cB(nI), cT(nI);

  for (int J = CJ0; J <= CJ1; ++J) {

    cL[J - CJ0] = FXc(CI0, J);

    cR[J - CJ0] = FXc(CI1 + 1, J);

  }

  for (int I = CI0; I <= CI1; ++I) {

    cB[I - CI0] = FYc(I, CJ0);

    cT[I - CI0] = FYc(I, CJ1 + 1);

  }


  advance_fab_1c(m, Uc, auxc, dxc, dyc, dt, fxc, fyc);  // Uc becomes the "t+dt" state


  // --- fine subcycling: r substeps, accumulation of fine fluxes (x dtf) ---

  std::vector<double> fL(nJ, 0), fR(nJ, 0), fB(nI, 0), fT(nI, 0);

  Fab2D fxf(xface_box(Uf.box()), 1, 0), fyf(yface_box(Uf.box()), 1, 0);

  for (int s = 0; s < r; ++s) {

    fill_fine_ghosts_t(Uf, Uc_old, Uc, double(s) / r);  // time-interpolated BC

    compute_fluxes_1c(m, Uf, auxf, fxf, fyf);

    const ConstArray4 FXf = fxf.const_array();

    const ConstArray4 FYf = fyf.const_array();

    for (int J = CJ0; J <= CJ1; ++J) {

      fL[J - CJ0] += 0.5 * (FXf(2 * CI0, 2 * J) + FXf(2 * CI0, 2 * J + 1)) * dtf;

      fR[J - CJ0] += 0.5 * (FXf(2 * CI1 + 2, 2 * J) + FXf(2 * CI1 + 2, 2 * J + 1)) * dtf;

    }

    for (int I = CI0; I <= CI1; ++I) {

      fB[I - CI0] += 0.5 * (FYf(2 * I, 2 * CJ0) + FYf(2 * I + 1, 2 * CJ0)) * dtf;

      fT[I - CI0] += 0.5 * (FYf(2 * I, 2 * CJ1 + 2) + FYf(2 * I + 1, 2 * CJ1 + 2)) * dtf;

    }

    advance_fab_1c(m, Uf, auxf, dxf, dyf, dtf, fxf, fyf);

  }


  average_down_fab(Uf, Uc, CI0, CI1, CJ0, CJ1);  // sync of the covered cells


  // --- reflux: coarse flux (x dt) replaced by sum of fine fluxes (x dtf) ---

  Array4 c = Uc.array();

  for (int J = CJ0; J <= CJ1; ++J) {

    c(CI0 - 1, J) -= (fL[J - CJ0] - cL[J - CJ0] * dt) / dxc;

    c(CI1 + 1, J) += (fR[J - CJ0] - cR[J - CJ0] * dt) / dxc;

  }

  for (int I = CI0; I <= CI1; ++I) {

    c(I, CJ0 - 1) -= (fB[I - CI0] - cB[I - CI0] * dt) / dyc;

    c(I, CJ1 + 1) += (fT[I - CI0] - cT[I - CI0] * dt) / dyc;

  }

}


}  // namespace pops

box2d.hpp
Box2D: the integer index space of a 2D cell-centered Cartesian grid.

pops::Fab2D
Single-grid data on a Box2D: VALID box + ng ghost layers, ncomp components, component-slow layout.
Definition fab2d.hpp:59

pops::Fab2D::grown_box
const Box2D & grown_box() const
Grown box (valid + ng ghosts) = actual memory footprint.
Definition fab2d.hpp:76

pops::Fab2D::box
const Box2D & box() const
VALID box (without ghosts).
Definition fab2d.hpp:74

pops::Fab2D::n_ghost
int n_ghost() const
Number of ghost layers.
Definition fab2d.hpp:80

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

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

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::fill_periodic_fab
void fill_periodic_fab(Fab2D &U, const Box2D &dom)
Definition amr_reflux.hpp:119

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::compute_fluxes_1c
void compute_fluxes_1c(const Model &m, const Fab2D &U, const Fab2D &aux, Fab2D &fx, Fab2D &fy)
Definition amr_reflux.hpp:94

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::average_down_fab
void average_down_fab(const Fab2D &Uf, Fab2D &Uc, int CI0, int CI1, int CJ0, int CJ1)
Definition amr_reflux.hpp:154

pops::coarsen
BoxArray coarsen(const BoxArray &ba, int r)
Coarsens each box of the BoxArray by a ratio r (coarsen box by box, order preserved).
Definition refinement.hpp:39

pops::amr_step_2level
void amr_step_2level(const Model &m, Fab2D &Uc, const Box2D &dom, double dxc, double dyc, Fab2D &Uf, int CI0, int CI1, int CJ0, int CJ1, const Fab2D &auxc, const Fab2D &auxf, double dt)
Definition amr_reflux.hpp:172

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

pops::kAmrRefRatio
constexpr int kAmrRefRatio
The native AMR refinement ratio between two consecutive levels.
Definition refinement_ratio.hpp:30

pops::advance_fab_1c
void advance_fab_1c(const Model &m, Fab2D &U, const Fab2D &aux, double dx, double dy, double dt, Fab2D &fx, Fab2D &fy)
Definition amr_reflux.hpp:109

pops::fill_fine_ghosts_t
void fill_fine_ghosts_t(Fab2D &Uf, const Fab2D &Uco, const Fab2D &Ucn, double frac)
Definition amr_reflux.hpp:138

refinement_ratio.hpp
Single source of truth for the native AMR refinement ratio.

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::Box2D
2D integer index space, cell-centered.
Definition box2d.hpp:37

pops::Box2D::nx
POPS_HD int nx() const
Width (direction 0). POPS_HD (called from Geometry::dx() in a device kernel).
Definition box2d.hpp:50

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

pops::Box2D::contains
bool contains(int i, int j) const
true if cell (i, j) is inside the box (lo/hi bounds inclusive).
Definition box2d.hpp:61

pops::Box2D::ny
POPS_HD int ny() const
Height (direction 1). POPS_HD (called from Geometry::dy() in a device kernel).
Definition box2d.hpp:52

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::detail::AdvanceFab1cKernel
Explicit Euler, 1 component: U -= dt div(F) on cell (i,j).
Definition amr_reflux.hpp:79

pops::detail::AdvanceFab1cKernel::FY
ConstArray4 FY
Definition amr_reflux.hpp:81

pops::detail::AdvanceFab1cKernel::dx
double dx
Definition amr_reflux.hpp:82

pops::detail::AdvanceFab1cKernel::uu
Array4 uu
Definition amr_reflux.hpp:80

pops::detail::AdvanceFab1cKernel::dy
double dy
Definition amr_reflux.hpp:82

pops::detail::AdvanceFab1cKernel::dt
double dt
Definition amr_reflux.hpp:82

pops::detail::AdvanceFab1cKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition amr_reflux.hpp:83

pops::detail::AdvanceFab1cKernel::FX
ConstArray4 FX
Definition amr_reflux.hpp:81

pops::detail::RusanovFaceXKernel
Rusanov flux at the left face (x axis) of cell (i,j).
Definition amr_reflux.hpp:50

pops::detail::RusanovFaceXKernel::m
Model m
Definition amr_reflux.hpp:51

pops::detail::RusanovFaceXKernel::F
Array4 F
Definition amr_reflux.hpp:53

pops::detail::RusanovFaceXKernel::ax
ConstArray4 ax
Definition amr_reflux.hpp:52

pops::detail::RusanovFaceXKernel::u
ConstArray4 u
Definition amr_reflux.hpp:52

pops::detail::RusanovFaceXKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition amr_reflux.hpp:54

pops::detail::RusanovFaceYKernel
Rusanov flux at the bottom face (y axis) of cell (i,j).
Definition amr_reflux.hpp:65

pops::detail::RusanovFaceYKernel::u
ConstArray4 u
Definition amr_reflux.hpp:67

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

pops::detail::RusanovFaceYKernel::ax
ConstArray4 ax
Definition amr_reflux.hpp:67

pops::detail::RusanovFaceYKernel::m
Model m
Definition amr_reflux.hpp:66

pops::detail::RusanovFaceYKernel::F
Array4 F
Definition amr_reflux.hpp:68

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

POPS_HD
#define POPS_HD
Definition types.hpp:25