include/pops/mesh/layout/refinement.hpp Source File

adc_cpp 0.3.0 Model-free C++23 core for coupled hyperbolic-elliptic systems on adaptive (AMR) meshes, with MPI and GPU (Kokkos) backends

• include
• pops
• mesh
• layout

refinement.hpp

Go to the documentation of this file.

1

12

13#pragma once

14#include <cassert>

15#include <cstdint>

16#include <cstdlib>

17#include <cstdio>

18

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

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

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

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

23#include <pops/mesh/boundary/fill_boundary.hpp> // detail::copy_shifted

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

25

26#include <utility>

27#include <vector>

28

29namespace pops {

30

34POPS_HD inline int coarsen_index(int a, int r) {

35 return floor_div(a, r);

36}

37

39inline BoxArray coarsen(const BoxArray& ba, int r) {

40 std::vector<Box2D> b;

41 b.reserve(ba.size());

42 for (int i = 0; i < ba.size(); ++i)

43 b.push_back(ba[i].coarsen(r));

44 return BoxArray{std::move(b)};

45}

46

50inline void parallel_copy(MultiFab& dst, const MultiFab& src) {

51 const int nc = std::min(dst.ncomp(), src.ncomp());

52 const BoxArray& sba = src.box_array();

53 const BoxArray& dba = dst.box_array();

54 const BoxHash shash(sba, suggest_bin(sba));

55

56 // --- local copies (dst local AND src local) ---

57 for (int ld = 0; ld < dst.local_size(); ++ld) {

58 Fab2D& D = dst.fab(ld);

59 const Box2D vd = D.box();

60 for (int gs : shash.query(vd)) {

61 const int ls = src.local_index_of(gs);

62 if (ls < 0)

63 continue; // src not local -> MPI below

64 const Box2D region = vd.intersect(sba[gs]);

65 if (region.empty())

66 continue;

67 detail::copy_shifted(D, src.fab(ls), region, 0, 0, nc);

68 }

69 }

70

71#ifdef POPS_HAS_MPI

72 if (n_ranks() <= 1)

73 return;

74 const int me = my_rank(), np = n_ranks();

75 const DistributionMapping& sdm = src.dmap();

76 const DistributionMapping& ddm = dst.dmap();

77

78 struct Job {

79 int gs, gd;

80 Box2D region;

81 };

82 std::vector<std::vector<Job>> send(np), recv(np);

83 for (int gd = 0; gd < dba.size(); ++gd) {

84 const int od = ddm[gd];

85 const Box2D vd = dba[gd];

86 for (int gs : shash.query(vd)) {

87 const int os = sdm[gs];

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

89 continue; // does not concern us

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

91 continue; // local, already done

92 const Box2D region = vd.intersect(sba[gs]);

93 if (region.empty())

94 continue;

95 if (os == me)

96 send[od].push_back({gs, gd, region}); // I own the src

97 else

98 recv[os].push_back({gs, gd, region}); // I own the dst

99 }

100 }

101

102 auto bufsz = [&](const std::vector<Job>& js) {

103 std::int64_t n = 0;

104 for (const auto& j : js)

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

106 return n;

107 };

108 device_fence(); // GPU: the local (device) copies precede the HOST pack

109 std::vector<std::vector<Real>> sbuf(np), rbuf(np);

110 std::vector<MPI_Request> reqs;

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

112 if (!send[r].empty()) {

113 sbuf[r].resize(bufsz(send[r]));

114 std::int64_t k = 0;

115 for (const auto& j : send[r]) {

116 const ConstArray4 s = src.fab(src.local_index_of(j.gs)).const_array();

117 for (int c = 0; c < nc; ++c)

118 for (int jj = j.region.lo[1]; jj <= j.region.hi[1]; ++jj)

119 for (int ii = j.region.lo[0]; ii <= j.region.hi[0]; ++ii)

120 sbuf[r][k++] = s(ii, jj, c);

121 }

122 reqs.emplace_back();

123 MPI_Isend(sbuf[r].data(), static_cast<int>(sbuf[r].size()), MPI_DOUBLE, r, 1, MPI_COMM_WORLD,

124 &reqs.back());

125 }

126 if (!recv[r].empty()) {

127 rbuf[r].resize(bufsz(recv[r]));

128 reqs.emplace_back();

129 MPI_Irecv(rbuf[r].data(), static_cast<int>(rbuf[r].size()), MPI_DOUBLE, r, 1, MPI_COMM_WORLD,

130 &reqs.back());

131 }

132 }

133 if (!reqs.empty())

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

135

136 device_fence(); // GPU: barrier before the HOST write of the received cells

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

138 std::int64_t k = 0;

139 for (const auto& j : recv[r]) {

140 Array4 d = dst.fab(dst.local_index_of(j.gd)).array();

141 for (int c = 0; c < nc; ++c)

142 for (int jj = j.region.lo[1]; jj <= j.region.hi[1]; ++jj)

143 for (int ii = j.region.lo[0]; ii <= j.region.hi[0]; ++ii)

144 d(ii, jj, c) = rbuf[r][k++];

145 }

146 }

147#endif

148}

149

150namespace detail {

151

152// NAMED FUNCTORS (rather than POPS_HD lambdas) for the AMR transfer kernels. Same reasons as the

153// rest of the mesh/elliptic path (cf. fill_boundary.hpp): refinement is first instantiated from

154// the MG V-cycle pulled in by an external TU; an extended lambda there trips up device kernel

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

156

158struct AverageDownKernel {

159 ConstArray4 F;

160 Array4 C;

161 Real inv;

162 int r, c;

163 POPS_HD void operator()(int I, int J) const {

164 Real s = 0;

165 for (int b = 0; b < r; ++b)

166 for (int a = 0; a < r; ++a)

167 s += F(r * I + a, r * J + b, c);

168 C(I, J, c) = s * inv;

169 }

170};

171

172} // namespace detail

173

177// PROVIDED-BUFFER variant: @p cfine is the "fine coarsen" grid (layout coarsen(fine.box_array(), r),

178// dmap = fine.dmap(), >= min(ncomp) components, 0 ghost) ALLOCATED by the caller and reused on each

179// call (hot path of the MG V-cycle: avoids one MultiFab allocation per restriction). Computation

180// STRICTLY identical to the allocating variant below.

181inline void average_down(const MultiFab& fine, MultiFab& coarse, int r, MultiFab& cfine) {

182 const int nc = std::min(fine.ncomp(), coarse.ncomp());

183 assert(cfine.box_array().size() == fine.box_array().size() && cfine.ncomp() >= nc &&

184 "average_down(scratch): cfine must be coarsen(fine, r) on the fine dmap");

185 const Real inv = Real(1) / (r * r);

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

187 const ConstArray4 F = fine.fab(li).const_array();

188 Array4 C = cfine.fab(li).array();

189 const Box2D cb = cfine.fab(li).box();

190 for (int c = 0; c < nc; ++c)

191 for_each_cell(cb, detail::AverageDownKernel{F, C, inv, r, c});

192 }

193 parallel_copy(coarse, cfine);

194}

195inline void average_down(const MultiFab& fine, MultiFab& coarse, int r) {

196 MultiFab cfine(coarsen(fine.box_array(), r), fine.dmap(), fine.ncomp(), 0);

197 average_down(fine, coarse, r, cfine);

198}

199

200namespace detail {

201

203struct InterpolateKernel {

204 Array4 F;

205 ConstArray4 C;

206 int r, c;

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

208 F(i, j, c) = C(coarsen_index(i, r), coarsen_index(j, r), c);

209 }

210};

211

212} // namespace detail

213

217// PROVIDED-BUFFER variant: @p cfine is the "fine coarsen" grid (same layout contract as

218// average_down above) allocated by the caller and reused (hot path of the MG V-cycle: avoids one

219// allocation per prolongation). Computation STRICTLY identical to the allocating variant.

220inline void interpolate(const MultiFab& coarse, MultiFab& fine, int r, MultiFab& cfine) {

221 const int nc = std::min(fine.ncomp(), coarse.ncomp());

222 assert(cfine.box_array().size() == fine.box_array().size() && cfine.ncomp() >= nc &&

223 "interpolate(scratch): cfine must be coarsen(fine, r) on the fine dmap");

224 parallel_copy(cfine, coarse); // bring the coarse values onto the fine-coarsen grid

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

226 Array4 F = fine.fab(li).array();

227 const ConstArray4 C = cfine.fab(li).const_array();

228 const Box2D fb = fine.fab(li).box();

229 for (int c = 0; c < nc; ++c)

230 for_each_cell(fb, detail::InterpolateKernel{F, C, r, c});

231 }

232}

233inline void interpolate(const MultiFab& coarse, MultiFab& fine, int r) {

234 MultiFab cfine(coarsen(fine.box_array(), r), fine.dmap(), fine.ncomp(), 0);

235 interpolate(coarse, fine, r, cfine);

236}

237

238} // namespace pops

box_array.hpp

BoxArray: the set of boxes tiling a level (disjoint, covering).

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::BoxArray::size

int size() const

Number of boxes in the tiling.

Definition box_array.hpp:42

pops::BoxHash

Spatial index of a BoxArray's boxes via a bin grid.

Definition box_hash.hpp:26

pops::BoxHash::query

std::vector< int > query(const Box2D &q) const

Indices (SORTED, without duplicates) of the boxes that may intersect q: guaranteed SUPERSET (no inter...

Definition box_hash.hpp:41

pops::DistributionMapping

Owning MPI rank of each box, indexed by GLOBAL box index (parallel to a BoxArray).

Definition distribution_mapping.hpp:19

pops::Fab2D

Single-grid data on a Box2D: VALID box + ng ghost layers, ncomp components, component-slow layout.

Definition fab2d.hpp:59

pops::Fab2D::box

const Box2D & box() const

VALID box (without ghosts).

Definition fab2d.hpp:74

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

fill_boundary.hpp

fill_boundary: INTRA-level halo exchange (fills ghosts from neighbors).

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

Definition amr_hierarchy.hpp:29

pops::average_down

void average_down(const MultiFab &fine, MultiFab &coarse, int r, MultiFab &cfine)

CONSERVATIVE average fine -> coarse (ratio r): coarse(I, J) = average of the r^2 fine cells of the bl...

Definition refinement.hpp:181

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::floor_div

POPS_HD constexpr int floor_div(int a, int b)

Integer division of a by b rounded down (handles a < 0 AND b < 0). POPS_HD constexpr (kernels).

Definition box2d.hpp:28

pops::Real

double Real

Definition types.hpp:30

pops::n_ranks

int n_ranks()

Definition comm.hpp:139

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::interpolate

void interpolate(const MultiFab &coarse, MultiFab &fine, int r, MultiFab &cfine)

Interpolation coarse -> fine (ratio r) by piecewise-CONSTANT injection: each fine cell (including the...

Definition refinement.hpp:220

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::my_rank

int my_rank()

Definition comm.hpp:136

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::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::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::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::Box2D::num_cells

std::int64_t num_cells() const

Total number of cells (nx*ny, floored at 0 per direction): 0 if the box is empty.

Definition box2d.hpp:54

pops::ConstArray4

READ-only handle (const counterpart of Array4): same layout and same contract (POD device-copyable,...

Definition fab2d.hpp:44

pops::detail::AverageDownKernel

CONSERVATIVE average of an r x r block: C(I, J, c) = (sum of the r^2 fine cells) * inv.

Definition refinement.hpp:158

pops::detail::AverageDownKernel::c

int c

Definition refinement.hpp:162

pops::detail::AverageDownKernel::operator()

POPS_HD void operator()(int I, int J) const

Definition refinement.hpp:163

pops::detail::AverageDownKernel::C

Array4 C

Definition refinement.hpp:160

pops::detail::AverageDownKernel::r

int r

Definition refinement.hpp:162

pops::detail::AverageDownKernel::F

ConstArray4 F

Definition refinement.hpp:159

pops::detail::AverageDownKernel::inv

Real inv

Definition refinement.hpp:161

pops::detail::InterpolateKernel

Piecewise-CONSTANT injection: F(i, j, c) receives the value of its covering coarse cell.

Definition refinement.hpp:203

pops::detail::InterpolateKernel::c

int c

Definition refinement.hpp:206

pops::detail::InterpolateKernel::C

ConstArray4 C

Definition refinement.hpp:205

pops::detail::InterpolateKernel::operator()

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

Definition refinement.hpp:207

pops::detail::InterpolateKernel::r

int r

Definition refinement.hpp:206

pops::detail::InterpolateKernel::F

Array4 F

Definition refinement.hpp:204

types.hpp

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

POPS_HD

#define POPS_HD

Definition types.hpp:25

---

Generated by  1.9.8