PoPS/doxygen/fill__boundary_8hpp_source.html

#pragma once


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

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

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

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

#include <pops/mesh/boundary/halo_schedule.hpp>

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

#include <pops/parallel/comm.hpp>


#include <cstdint>

#include <memory>

#include <utility>

#include <vector>


namespace pops {


struct Periodicity {

  bool x = false;

  bool y = false;

};


namespace detail {


// NAMED FUNCTORS (not POPS_HD lambdas) for the halo-exchange kernels. Same reasons as the rest of the

// elliptic/mesh path (#93, recipe #64): fill_boundary is first instantiated from the MG V-cycle pulled

// from an external TU; an extended lambda there stalls device kernel emission under nvcc (-O Release

// without -g). Strictly identical body -> bit-identical CPU and device.


struct CopyShiftedKernel {

  Array4 d;

  ConstArray4 s;

  int sx, sy, c;

  POPS_HD void operator()(int i, int j) const { d(i, j, c) = s(i - sx, j - sy, c); }

};


// dst(i, j, c) = src(i - sx, j - sy, c) for (i, j) in region.


inline void copy_shifted(Fab2D& dst, const Fab2D& src, const Box2D& region, int sx, int sy,

                         int ncomp) {

  Array4 d = dst.array();

  ConstArray4 s = src.const_array();

  for (int c = 0; c < ncomp; ++c)

    for_each_cell(region, CopyShiftedKernel{d, s, sx, sy, c});

}


// Pack of a send job: sb[b0 + c*rsz + off] = s(i - sx, jc - sy, c), off = (jc-lo1)*rnx + (i-lo0).


struct PackKernel {

  Real* sb;

  ConstArray4 s;

  std::int64_t b0, rsz;

  int lo0, lo1, rnx, sx, sy, ncl;


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

    const std::int64_t off = static_cast<std::int64_t>(jc - lo1) * rnx + (i - lo0);

    for (int c = 0; c < ncl; ++c)

      sb[b0 + static_cast<std::int64_t>(c) * rsz + off] = s(i - sx, jc - sy, c);

  }


};


// Unpack of a receive job: d(i, jc, c) = rb[b0 + c*rsz + off], off = (jc-lo1)*rnx + (i-lo0).


struct UnpackKernel {

  const Real* rb;

  Array4 d;

  std::int64_t b0, rsz;

  int lo0, lo1, rnx, ncl;


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

    const std::int64_t off = static_cast<std::int64_t>(jc - lo1) * rnx + (i - lo0);

    for (int c = 0; c < ncl; ++c)

      d(i, jc, c) = rb[b0 + static_cast<std::int64_t>(c) * rsz + off];

  }


};


// Enumerates the halo schedule for (mf layout, per, domain): the BoxHash build + the local (dst AND

// src local) and, under MPI with n_ranks()>1, the global (cross-rank send/recv) job lists. This is

// the per-call work that ADC-260 hoists out of fill_boundary_begin; it runs ONCE per distinct

// (layout, Periodicity, domain) and is then replayed from the cache. Jobs are produced in the SAME

// deterministic order as the legacy inline loops (local: li x shifts x sorted gB; global: gF x

// shifts x sorted gB), so the packed buffers stay bit-identical and the per-rank send/recv lists

// stay aligned. Bumps the build counter (cache-engagement test hook).


inline void build_halo_schedule(const MultiFab& mf, const Box2D& domain, Periodicity per,

                                HaloSchedule& sched) {

  ++halo_schedule_build_counter();

  const int ng = mf.n_grow();

  const int Lx = domain.nx();

  const int Ly = domain.ny();

  const BoxArray& ba = mf.box_array();


  std::vector<int> sxv = {0};

  if (per.x) {

    sxv.push_back(Lx);

    sxv.push_back(-Lx);

  }

  std::vector<int> syv = {0};

  if (per.y) {

    syv.push_back(Ly);

    syv.push_back(-Ly);

  }

  std::vector<std::pair<int, int>> shifts;

  for (int sx : sxv)

    for (int sy : syv)

      shifts.push_back({sx, sy});


  // spatial hash: restricts the neighbor-box search (see box_hash.hpp).

  const BoxHash hash(ba, suggest_bin(ba));


  // --- local jobs (local dst AND local src) ---

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

    const int gF = mf.global_index(li);

    const Box2D gbox = mf.fab(li).box().grow(ng);

    for (auto [sx, sy] : shifts) {

      const Box2D Q = gbox.shift(0, -sx).shift(1, -sy);

      for (int gB : hash.query(Q)) {

        if (gB == gF && sx == 0 && sy == 0)

          continue;  // self, without shift

        const int srcLocal = mf.local_index_of(gB);

        if (srcLocal < 0)

          continue;  // non-local src -> MPI below

        const Box2D region = gbox.intersect(ba[gB].shift(0, sx).shift(1, sy));

        if (region.empty())

          continue;

        sched.local.push_back({gB, gF, sx, sy, region});

      }

    }

  }


#ifdef POPS_HAS_MPI

  const int np = n_ranks();

  if (np > 1) {

    const int me = my_rank();

    const DistributionMapping& dm = mf.dmap();

    sched.send.assign(np, {});

    sched.recv.assign(np, {});

    // deterministic global enumeration: (dst gF) x shifts x hash candidates (gB sorted). Identical

    // on all ranks -> aligned send/recv lists.

    for (int gF = 0; gF < ba.size(); ++gF) {

      const int od = dm[gF];

      const Box2D gbox = ba[gF].grow(ng);

      for (auto [sx, sy] : shifts) {

        const Box2D Q = gbox.shift(0, -sx).shift(1, -sy);

        for (int gB : hash.query(Q)) {

          if (gB == gF && sx == 0 && sy == 0)

            continue;

          const int os = dm[gB];

          if (od != me && os != me)

            continue;

          if (od == me && os == me)

            continue;

          const Box2D region = gbox.intersect(ba[gB].shift(0, sx).shift(1, sy));

          if (region.empty())

            continue;

          if (os == me)

            sched.send[od].push_back({gB, gF, sx, sy, region});

          else

            sched.recv[os].push_back({gB, gF, sx, sy, region});

        }

      }

    }

  }

#endif

}


// Returns the cached schedule for (mf layout, per, domain), building and memoizing it on first use.


inline std::shared_ptr<const HaloSchedule> get_halo_schedule(const MultiFab& mf,

                                                             const Box2D& domain, Periodicity per) {

  HaloScheduleCache& cache = mf.halo_cache();

  if (std::shared_ptr<const HaloSchedule> hit = cache.find(per.x, per.y, domain))

    return hit;

  std::shared_ptr<HaloSchedule> s = cache.add();

  s->per_x = per.x;

  s->per_y = per.y;

  s->domain = domain;

  build_halo_schedule(mf, domain, per, *s);

  return s;

}


}  // namespace detail


struct HaloExchange {

  std::shared_ptr<const HaloSchedule> sched;  // replayed plan (null if mf has no ghost)

#ifdef POPS_HAS_MPI

  // Buffers in PINNED HOST memory (comm_allocator = Kokkos::SharedHostPinnedSpace under Kokkos,

  // std::allocator otherwise), NOT managed. The pack/unpack in for_each (device under Kokkos)

  // writes/reads directly into them since pinned host is device-accessible; BUT the pointer passed to

  // MPI is seen as HOST (cuPointerGetAttribute = HOST), so a CUDA-aware MPI (BTL smcuda) does NOT

  // attempt CUDA IPC on it. A managed/UVM pointer, on the other hand, triggered IPC, which DEADLOCKS

  // between two GPUs isolated by cgroup (srun --gpus-per-task=1: each rank sees only its GPU as device

  // 0, cuIpcOpenMemHandle of the peer's buffer impossible). See core/allocator.hpp (comm_allocator).

  std::vector<std::vector<Real, comm_allocator<Real>>> sbuf, rbuf;  // alive until end

  std::vector<MPI_Request> reqs;

  int nc = 0;

#endif

};


inline HaloExchange fill_boundary_begin(MultiFab& mf, const Box2D& domain, Periodicity per = {}) {

  HaloExchange h;

  const int ng = mf.n_grow();

  if (ng == 0)

    return h;

  const int nc = mf.ncomp();

  // memoized schedule (BoxHash + enumeration) for this (layout, Periodicity, domain).

  const std::shared_ptr<const HaloSchedule> sched = detail::get_halo_schedule(mf, domain, per);

  h.sched = sched;


  // --- local copies (local dst AND local src), replayed from the cached plan ---

  for (const HaloJob& j : sched->local) {

    Fab2D& dst = mf.fab(mf.local_index_of(j.dst));

    const Fab2D& src = mf.fab(mf.local_index_of(j.src));

    detail::copy_shifted(dst, src, j.region, j.sx, j.sy, nc);

  }


#ifdef POPS_HAS_MPI

  if (n_ranks() <= 1)

    return h;

  const int np = n_ranks();

  h.nc = nc;


  auto buf_size = [&](const std::vector<HaloJob>& js) {

    std::int64_t n = 0;

    for (const auto& j : js)

      n += j.region.num_cells() * nc;

    return n;

  };

  h.sbuf.assign(np, {});

  h.rbuf.assign(np, {});

  // device PACK (for_each, parallel under Kokkos) into the pinned host buffers. Per-job layout:

  // c-major then (jj, ii), IDENTICAL to the old k++ order -> buffer bit-identical to the host path

  // (the CPU MPI ctests stay bit-identical at np=1/2/4). The peer rank enumerates in the same order,

  // so sbuf[A->B] and rbuf[B<-A] align without negotiating sizes.

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

    const std::vector<HaloJob>& send_r = sched->send[r];

    if (send_r.empty())

      continue;

    h.sbuf[r].resize(buf_size(send_r));

    Real* sb = h.sbuf[r].data();

    std::int64_t base = 0;

    for (const auto& jb : send_r) {

      const ConstArray4 s = mf.fab(mf.local_index_of(jb.src)).const_array();

      const int lo0 = jb.region.lo[0], lo1 = jb.region.lo[1], rnx = jb.region.nx();

      const std::int64_t rsz = static_cast<std::int64_t>(rnx) * jb.region.ny();

      const int sx = jb.sx, sy = jb.sy, ncl = nc;

      const std::int64_t b0 = base;

      for_each_cell(jb.region, detail::PackKernel{sb, s, b0, rsz, lo0, lo1, rnx, sx, sy, ncl});

      base += rsz * nc;

    }

  }

  for (int r = 0; r < np; ++r)  // allocate the receive buffers

    if (!sched->recv[r].empty())

      h.rbuf[r].resize(buf_size(sched->recv[r]));

  device_fence();  // the pack kernels (and the local copies) must finish before MPI reads sbuf

  for (

      int r = 0; r < np;

      ++r) {  // non-blocking posting; MPI receives PINNED HOST pointers (seen HOST, no GPUDirect/CUDA IPC)

    if (!h.sbuf[r].empty()) {

      h.reqs.emplace_back();

      MPI_Isend(h.sbuf[r].data(), static_cast<int>(h.sbuf[r].size()), MPI_DOUBLE, r, 0,

                MPI_COMM_WORLD, &h.reqs.back());

    }

    if (!h.rbuf[r].empty()) {

      h.reqs.emplace_back();

      MPI_Irecv(h.rbuf[r].data(), static_cast<int>(h.rbuf[r].size()), MPI_DOUBLE, r, 0,

                MPI_COMM_WORLD, &h.reqs.back());

    }

  }

#endif

  return h;

}


inline void fill_boundary_end(MultiFab& mf, HaloExchange& h) {

#ifdef POPS_HAS_MPI

  if (h.reqs.empty())

    return;

  MPI_Waitall(static_cast<int>(h.reqs.size()), h.reqs.data(), MPI_STATUSES_IGNORE);

  // device UNPACK (for_each) from the received PINNED HOST buffers. Waitall guarantees the transfer is

  // complete; the kernel launched next reads the pinned host (device-accessible, coherent). Replayed

  // from the SAME cached recv list begin used (h.sched), so base offsets match the sender's layout.

  const HaloSchedule& sched = *h.sched;

  for (std::size_t r = 0; r < sched.recv.size(); ++r) {

    if (h.rbuf[r].empty())

      continue;

    const Real* rb = h.rbuf[r].data();

    std::int64_t base = 0;

    for (const auto& jb : sched.recv[r]) {

      Array4 d = mf.fab(mf.local_index_of(jb.dst)).array();

      const int lo0 = jb.region.lo[0], lo1 = jb.region.lo[1], rnx = jb.region.nx();

      const std::int64_t rsz = static_cast<std::int64_t>(rnx) * jb.region.ny();

      const int ncl = h.nc;

      const std::int64_t b0 = base;

      for_each_cell(jb.region, detail::UnpackKernel{rb, d, b0, rsz, lo0, lo1, rnx, ncl});

      base += rsz * ncl;

    }

  }

  // The unpack kernels above are ASYNC (device) and read h.rbuf; comm_allocator (pinned host) frees

  // IMMEDIATELY on destruction of h (no deferred free like the ManagedArena). So we drain the device

  // BEFORE the pinned buffers are freed -> no use-after-free.

  device_fence();

#else

  (void)mf;

  (void)h;

#endif

}


inline void fill_boundary(MultiFab& mf, const Box2D& domain, Periodicity per = {}) {

  HaloExchange h = fill_boundary_begin(mf, domain, per);

  fill_boundary_end(mf, h);

}


}  // namespace pops

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

box_hash.hpp
BoxHash: spatial hash for fast lookup of boxes intersecting a query.

pops::BoxArray
Ordered list of boxes tiling a level.
Definition box_array.hpp:22

pops::BoxHash
Spatial index of a BoxArray's boxes via a bin grid.
Definition box_hash.hpp:26

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

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::HaloScheduleCache
Small per-MultiFab cache of halo schedules, one entry per distinct (Periodicity, domain).
Definition halo_schedule.hpp:57

pops::HaloScheduleCache::add
std::shared_ptr< HaloSchedule > add()
Appends a fresh, empty schedule and returns it for the caller to populate.
Definition halo_schedule.hpp:70

pops::HaloScheduleCache::find
std::shared_ptr< const HaloSchedule > find(bool px, bool py, const Box2D &dom) const
Existing schedule for (px, py, dom), or nullptr if none is cached yet.
Definition halo_schedule.hpp:60

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::local_index_of
int local_index_of(int global) const
LOCAL index of the global box global, or -1 if it is not owned by this rank.
Definition multifab.hpp:75

pops::MultiFab::ncomp
int ncomp() const
Number of components.
Definition multifab.hpp:60

pops::MultiFab::halo_cache
HaloScheduleCache & halo_cache() const
Internal (ADC-260): memoized halo-exchange schedule used by fill_boundary.
Definition multifab.hpp:100

pops::MultiFab::n_grow
int n_grow() const
Number of ghost layers.
Definition multifab.hpp:62

pops::MultiFab::box_array
const BoxArray & box_array() const
GLOBAL decomposition of the level (all boxes, all ranks).
Definition multifab.hpp:56

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

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

halo_schedule.hpp
HaloSchedule: memoized intra-level halo-exchange plan for fill_boundary (ADC-260).

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

pops::detail::copy_shifted
void copy_shifted(Fab2D &dst, const Fab2D &src, const Box2D &region, int sx, int sy, int ncomp)
Definition fill_boundary.hpp:56

pops::detail::get_halo_schedule
std::shared_ptr< const HaloSchedule > get_halo_schedule(const MultiFab &mf, const Box2D &domain, Periodicity per)
Definition fill_boundary.hpp:180

pops::detail::build_halo_schedule
void build_halo_schedule(const MultiFab &mf, const Box2D &domain, Periodicity per, HaloSchedule &sched)
Definition fill_boundary.hpp:97

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::n_ranks
int n_ranks()
Definition comm.hpp:139

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::fill_boundary
void fill_boundary(MultiFab &mf, const Box2D &domain, Periodicity per={})
BLOCKING halo exchange: begin then end immediately (no overlap).
Definition fill_boundary.hpp:333

pops::fill_boundary_end
void fill_boundary_end(MultiFab &mf, HaloExchange &h)
Phase 2 (blocking): MPI_Waitall on the transfers posted by begin, then unpacks the received buffers i...
Definition fill_boundary.hpp:297

pops::my_rank
int my_rank()
Definition comm.hpp:136

pops::fill_boundary_begin
HaloExchange fill_boundary_begin(MultiFab &mf, const Box2D &domain, Periodicity per={})
Phase 1 (non-blocking): does the LOCAL halo copies and posts the Isend/Irecv of the distant halos.
Definition fill_boundary.hpp:220

pops::suggest_bin
int suggest_bin(const BoxArray &ba)
Recommended bin size for a BoxArray: the largest box extent (at least 1), so that neighboring boxes f...
Definition box_hash.hpp:72

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::grow
Box2D grow(int n) const
Grows the box by n cells in ALL directions (uniform ghost layer).
Definition box2d.hpp:69

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::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::shift
Box2D shift(int d, int s) const
Translates the box by s cells in direction d (lo and hi shifted by the same s).
Definition box2d.hpp:78

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

pops::HaloExchange
Opaque state of an in-flight halo exchange, returned by fill_boundary_begin and consumed by fill_boun...
Definition fill_boundary.hpp:201

pops::HaloExchange::sched
std::shared_ptr< const HaloSchedule > sched
Definition fill_boundary.hpp:202

pops::HaloSchedule
Memoized schedule for ONE (Periodicity, domain) over a fixed layout.
Definition halo_schedule.hpp:40

pops::HaloSchedule::recv
std::vector< std::vector< HaloJob > > recv
Definition halo_schedule.hpp:46

pops::Periodicity
Per-direction periodicity: halo wrapping in x and/or y during the exchange (false = open edge,...
Definition fill_boundary.hpp:37

pops::Periodicity::x
bool x
Definition fill_boundary.hpp:38

pops::Periodicity::y
bool y
Definition fill_boundary.hpp:39

pops::detail::CopyShiftedKernel
Definition fill_boundary.hpp:48

pops::detail::CopyShiftedKernel::sy
int sy
Definition fill_boundary.hpp:51

pops::detail::CopyShiftedKernel::d
Array4 d
Definition fill_boundary.hpp:49

pops::detail::CopyShiftedKernel::c
int c
Definition fill_boundary.hpp:51

pops::detail::CopyShiftedKernel::s
ConstArray4 s
Definition fill_boundary.hpp:50

pops::detail::CopyShiftedKernel::sx
int sx
Definition fill_boundary.hpp:51

pops::detail::CopyShiftedKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition fill_boundary.hpp:52

pops::detail::PackKernel
Definition fill_boundary.hpp:65

pops::detail::PackKernel::s
ConstArray4 s
Definition fill_boundary.hpp:67

pops::detail::PackKernel::rnx
int rnx
Definition fill_boundary.hpp:69

pops::detail::PackKernel::ncl
int ncl
Definition fill_boundary.hpp:69

pops::detail::PackKernel::lo1
int lo1
Definition fill_boundary.hpp:69

pops::detail::PackKernel::b0
std::int64_t b0
Definition fill_boundary.hpp:68

pops::detail::PackKernel::rsz
std::int64_t rsz
Definition fill_boundary.hpp:68

pops::detail::PackKernel::sy
int sy
Definition fill_boundary.hpp:69

pops::detail::PackKernel::lo0
int lo0
Definition fill_boundary.hpp:69

pops::detail::PackKernel::operator()
POPS_HD void operator()(int i, int jc) const
Definition fill_boundary.hpp:70

pops::detail::PackKernel::sx
int sx
Definition fill_boundary.hpp:69

pops::detail::PackKernel::sb
Real * sb
Definition fill_boundary.hpp:66

pops::detail::UnpackKernel
Definition fill_boundary.hpp:78

pops::detail::UnpackKernel::b0
std::int64_t b0
Definition fill_boundary.hpp:81

pops::detail::UnpackKernel::operator()
POPS_HD void operator()(int i, int jc) const
Definition fill_boundary.hpp:83

pops::detail::UnpackKernel::lo1
int lo1
Definition fill_boundary.hpp:82

pops::detail::UnpackKernel::rb
const Real * rb
Definition fill_boundary.hpp:79

pops::detail::UnpackKernel::rnx
int rnx
Definition fill_boundary.hpp:82

pops::detail::UnpackKernel::rsz
std::int64_t rsz
Definition fill_boundary.hpp:81

pops::detail::UnpackKernel::lo0
int lo0
Definition fill_boundary.hpp:82

pops::detail::UnpackKernel::ncl
int ncl
Definition fill_boundary.hpp:82

pops::detail::UnpackKernel::d
Array4 d
Definition fill_boundary.hpp:80

POPS_HD
#define POPS_HD
Definition types.hpp:25