PoPS/doxygen/amr__regrid__coupler_8hpp_source.html

#pragma once


#include <pops/amr/tagging/cluster.hpp>  // berger_rigoutsos, ClusterParams

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

#include <pops/amr/regridding/regrid.hpp>   // tag_cells, grow_tags

#include <pops/amr/tagging/tag_box.hpp>  // TagBox

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

#include <pops/numerics/time/amr/reflux/amr_reflux_mf.hpp>  // AmrLevelMP, mf_find_box

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

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

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

#include <pops/mesh/execution/for_each.hpp>  // device_fence (barrier after async parallel_copy under Cuda)

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

#include <pops/mesh/layout/refinement.hpp>  // coarsen_index

#include <pops/parallel/comm.hpp>    // n_ranks (explicit include, no longer an indirect path)


#include <algorithm>

#include <utility>

#include <vector>


namespace pops {


// ------------------------------------------------------------------------------------------------

// REFACTOR (docs/AMR_REGRID_UNION_TAGS_DESIGN.md section 6): amr_regrid_finest split into TWO

// responsibilities so the multi-block UNION regrid (AmrRuntime::regrid) can re-grid

// SEVERAL fields on a SINGLE layout imposed from outside (the same one for all blocks):

//   (1) regrid_compute_fine_layout: tags of a parent -> grow -> all_reduce_or (if distributed) ->

//       berger_rigoutsos -> nesting clamp -> (fine BoxArray, DistributionMapping). This is the

//       COMPUTATION of the layout, done ONCE on the union of the tags by the multi-block path.

//   (2) regrid_field_on_layout: takes an IMPOSED (fb, dmap) and rebuilds ONE fine MultiFab (parent

//       interp + fine carry-over), exactly the BODY of the old amr_regrid_finest, but without the

//       layout computation. Called PER BLOCK on the same union layout.

// amr_regrid_finest stays the CHAINING of (1) then (2) on a SINGLE block -> the single-block path

// (AmrCouplerMP::regrid) stays BIT-IDENTICAL (same operations, same order).

// ------------------------------------------------------------------------------------------------


inline std::pair<BoxArray, DistributionMapping> regrid_compute_fine_layout(

    TagBox grown, const Box2D& pdom, int pk, int margin, bool coarse_replicated = true) {

  const int PNX = pdom.nx(), PNY = pdom.ny();

  // DISTRIBUTED COARSE (pk == 0 && !coarse_replicated): each rank only tagged its LOCAL boxes.

  // Global OR before the clustering -> all ranks start from the SAME tag grid (otherwise the

  // fine BoxArray would differ per rank -> incompatible dmaps, MPI desynchronized). Replicated: the

  // grid is already complete on each rank (no-op, all_reduce_or would be the identity).

  if (pk == 0 && !coarse_replicated)

    all_reduce_or_inplace(grown.t.data(), static_cast<int>(grown.t.size()));

  std::vector<Box2D> cl = berger_rigoutsos(grown, ClusterParams{});

  std::vector<Box2D> fb;  // fine patches (fine level coords = parent x2)

  for (Box2D b : cl) {

    b.lo[0] = std::max(b.lo[0], margin);

    b.lo[1] = std::max(b.lo[1], margin);

    b.hi[0] = std::min(b.hi[0], PNX - 1 - margin);

    b.hi[1] = std::min(b.hi[1], PNY - 1 - margin);

    if (b.hi[0] < b.lo[0] || b.hi[1] < b.lo[1])

      continue;

    fb.push_back(Box2D{{2 * b.lo[0], 2 * b.lo[1]}, {2 * b.hi[0] + 1, 2 * b.hi[1] + 1}});

  }

  if (fb.empty())

    return {BoxArray{}, DistributionMapping{}};  // nothing to refine

  BoxArray ba(fb);

  return {ba, DistributionMapping(static_cast<int>(ba.size()), n_ranks())};

}


inline MultiFab regrid_field_on_layout(const BoxArray& fb, const DistributionMapping& dmap,

                                       const MultiFab& par, const MultiFab& old, int pk, int ngf,

                                       bool coarse_replicated = true) {

  MultiFab nU(fb, dmap, old.ncomp(), ngf);

  const int ncf = nU.ncomp();

  // DISTRIBUTED parent (de-replicated coarse): par.fab only holds the LOCAL boxes, so

  // mf_find_box would return -1 for a coarse cell owned by a REMOTE rank and the patch

  // would stay uninitialized there. We bring the needed parent regions onto a LOCAL

  // child-coarsen grid (coarsen of the fine BoxArray) via parallel_copy, then interpolate from

  // it. Replicated parent: par is entirely local, direct read via mf_find_box.

  const bool par_replicated = (pk != 0) || coarse_replicated;

  MultiFab parloc;

  if (!par_replicated) {

    parloc = MultiFab(coarsen(nU.box_array(), kAmrRefRatio), nU.dmap(), par.ncomp(), 0);

    parallel_copy(parloc, par);

    // parallel_copy launches async kernels under Cuda and, at np=1, returns WITHOUT a fence: without this

    // fence the read of parloc below would read device memory not yet written -> NaN.

    device_fence();

  }

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

    Array4 a = nU.fab(li).array();

    const Box2D nb = nU.box(li);

    if (par_replicated) {

      for (int j = nb.lo[1]; j <= nb.hi[1]; ++j)  // 1) interp from the parent (local)

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

          const int pb =

              mf_find_box(par, coarsen_index(i, kAmrRefRatio), coarsen_index(j, kAmrRefRatio));

          if (pb < 0)

            continue;

          const ConstArray4 pp = par.fab(pb).const_array();

          for (int k = 0; k < ncf; ++k)

            a(i, j, k) = pp(coarsen_index(i, kAmrRefRatio), coarsen_index(j, kAmrRefRatio), k);

        }

    } else {

      const ConstArray4 pp = parloc.fab(li).const_array();  // local child-coarsen grid

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

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

          for (int k = 0; k < ncf; ++k)

            a(i, j, k) = pp(coarsen_index(i, kAmrRefRatio), coarsen_index(j, kAmrRefRatio), k);

    }

    for (int ol = 0; ol < old.local_size(); ++ol) {  // 2) carry-over of the fine data

      const ConstArray4 o = old.fab(ol).const_array();

      const Box2D inter = nb.intersect(old.box(ol));

      if (inter.empty())

        continue;

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

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

          for (int k = 0; k < ncf; ++k)

            a(i, j, k) = o(i, j, k);

    }

  }

  return nU;

}


template <class Crit>


void amr_regrid_finest(std::vector<AmrLevelMP>& L, std::vector<MultiFab>& aux, const Box2D& dom,

                       Crit crit, int grow, int margin, int aux_ncomp = kAuxBaseComps,

                       bool coarse_replicated = true) {

  const int nlev = static_cast<int>(L.size());

  if (nlev < 2)

    return;

  const int fk = nlev - 1, pk = fk - 1;  // fine and its parent

  const int PNX = dom.nx() << pk, PNY = dom.ny() << pk;

  const Box2D pdom = Box2D::from_extents(PNX, PNY);

  TagBox tags = tag_cells(L[pk].U, pdom, crit);

  TagBox grown = grow_tags(tags, grow, pdom);

  // (1) Compute the fine layout (tags -> grow [already done] -> all_reduce_or -> clustering -> clamp).

  auto [fb, dmap] =

      regrid_compute_fine_layout(std::move(grown), pdom, pk, margin, coarse_replicated);

  if (fb.size() == 0)

    return;  // nothing to refine: keep the current grid

  // The new patches INHERIT the ghost width of the level being replaced (not a frozen 1): a

  // level rebuilt in 2nd-order MUSCL (Minmod / VanLeer) carries 2 ghosts, which the regrid must

  // preserve, otherwise the reconstruction would read out of bounds after re-refinement.

  const int ngf = L[fk].U.n_grow();

  // (2) Re-grid the field U on this layout (parent interp + fine carry-over): SAME body as before,

  // so the single-block path stays BIT-IDENTICAL (chaining of (1) then (2) on a single block).

  L[fk].U = regrid_field_on_layout(fb, dmap, L[pk].U, L[fk].U, pk, ngf, coarse_replicated);

  aux[fk] = MultiFab(L[fk].U.box_array(), L[fk].U.dmap(), aux_ncomp, 1);  // stable address

  L[fk].aux = &aux[fk];

}


}  // namespace pops

amr_reflux_mf.hpp
Umbrella for the AMR MultiFab stack: includes the numerics/time sub-headers in dependency order (flux...

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

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::DistributionMapping
Owning MPI rank of each box, indexed by GLOBAL box index (parallel to a BoxArray).
Definition distribution_mapping.hpp:19

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::ncomp
int ncomp() const
Number of components.
Definition multifab.hpp:60

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

cluster.hpp
Berger-Rigoutsos clustering: tag grid -> small set of boxes covering the tagged cells.

comm.hpp
Parallel seam: minimal MPI abstraction (rank/size + collectives) with serial fallback.

distribution_mapping.hpp
DistributionMapping: maps each box (by global index) to its owning MPI rank.

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

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

pops
Definition amr_hierarchy.hpp:29

pops::mf_find_box
int mf_find_box(const MultiFab &mf, int I, int J)
Definition amr_subcycling.hpp:231

pops::kAuxBaseComps
constexpr int kAuxBaseComps
Definition state.hpp:147

pops::n_ranks
int n_ranks()
Definition comm.hpp:139

pops::all_reduce_or_inplace
void all_reduce_or_inplace(char *, int)
Definition comm.hpp:156

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::regrid_compute_fine_layout
std::pair< BoxArray, DistributionMapping > regrid_compute_fine_layout(TagBox grown, const Box2D &pdom, int pk, int margin, bool coarse_replicated=true)
Compute the fine layout (BoxArray + DistributionMapping) of a Berger-Rigoutsos regrid from the grown ...
Definition amr_regrid_coupler.hpp:53

pops::device_fence
void device_fence()
Device barrier: waits for in-flight kernels to finish before a HOST access to unified memory.
Definition kokkos_env.hpp:43

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

pops::amr_regrid_finest
void amr_regrid_finest(std::vector< AmrLevelMP > &L, std::vector< MultiFab > &aux, const Box2D &dom, Crit crit, int grow, int margin, int aux_ncomp=kAuxBaseComps, bool coarse_replicated=true)
Regrid the finest level (L.back()) by Berger-Rigoutsos on the criterion crit applied to the parent: r...
Definition amr_regrid_coupler.hpp:151

pops::parallel_copy
void parallel_copy(MultiFab &dst, const MultiFab &src)
Copies the valid regions that OVERLAP from src to dst (same indices, no shift).
Definition refinement.hpp:50

pops::berger_rigoutsos
std::vector< Box2D > berger_rigoutsos(const TagBox &tags, const ClusterParams &p={})
Cluster a TagBox into boxes covering the tagged cells (Berger-Rigoutsos), then final chop.
Definition cluster.hpp:185

pops::grow_tags
TagBox grow_tags(const TagBox &in, int n, const Box2D &domain)
Grows the tags by n cells (square neighborhood), staying within the domain.
Definition regrid.hpp:58

pops::coarsen_index
POPS_HD int coarsen_index(int a, int r)
Index of the coarse cell containing the fine cell a (FLOOR division by r, handles a < 0).
Definition refinement.hpp:34

pops::regrid_field_on_layout
MultiFab regrid_field_on_layout(const BoxArray &fb, const DistributionMapping &dmap, const MultiFab &par, const MultiFab &old, int pk, int ngf, bool coarse_replicated=true)
Rebuild ONE fine MultiFab on the IMPOSED layout fb / dmap (the same one for all blocks in multi-block...
Definition amr_regrid_coupler.hpp:86

pops::tag_cells
TagBox tag_cells(const MultiFab &mf, const Box2D &domain, Crit crit)
Marks the valid cells where the predicate is true, on a TagBox covering the domain.
Definition regrid.hpp:40

refinement.hpp
AMR inter-level transfer operators (integer ratio r) + parallel_copy.

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

regrid.hpp
Dynamic regrid: (re)builds a fine level from the tagging of a coarse level.

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::empty
bool empty() const
true if the box contains no cell (hi < lo in one direction).
Definition box2d.hpp:58

pops::Box2D::from_extents
static Box2D from_extents(int nx, int ny)
Box [0, nx-1] x [0, ny-1] covering nx*ny cells from the index origin.
Definition box2d.hpp:42

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::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::Box2D::intersect
Box2D intersect(const Box2D &o) const
Intersection of the two boxes (possibly empty: hi < lo if they do not overlap).
Definition box2d.hpp:95

pops::ClusterParams
Berger-Rigoutsos clustering parameters (configuration object).
Definition cluster.hpp:31

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

pops::TagBox
Dense grid of 0/1 markers over a box, input to Berger-Rigoutsos clustering.
Definition tag_box.hpp:35

pops::TagBox::t
std::vector< char > t
Definition tag_box.hpp:37

tag_box.hpp
TagBox: dense grid of markers (0/1) over a region, input to Berger-Rigoutsos clustering.

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