include/pops/mesh/storage/mf_arith.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
• storage

mf_arith.hpp

Go to the documentation of this file.

1

11

12#pragma once

13

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

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

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

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

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

19#include <pops/parallel/comm.hpp> // all_reduce_sum: COLLECTIVE dot product (Krylov under MPI)

20

21#include <algorithm>

22#include <limits>

23

24namespace pops {

25

26namespace detail {

27// NAMED FUNCTORS (not POPS_HD lambdas) for the MultiFab arithmetic kernels. Same recipe as

28// the block path (#64): these operations are first-instantiated from the MG V-cycle, itself pulled

29// from an external TU (native harness/loader); an extended lambda at this spot makes nvcc stumble on

30// device kernel emission (null kernel-stub -> Cuda segfault in -O Release without -g, #93). Body

31// strictly identical to the old lambdas -> bit-identical on CPU and device.

32struct SaxpyKernel {

33 Array4 Y;

34 ConstArray4 X;

35 Real a;

36 int c;

37 POPS_HD void operator()(int i, int j) const { Y(i, j, c) += a * X(i, j, c); }

38};

39

40struct LincombKernel {

41 Array4 Z;

42 ConstArray4 X, Y;

43 Real a, b;

44 int c;

45 POPS_HD void operator()(int i, int j) const { Z(i, j, c) = a * X(i, j, c) + b * Y(i, j, c); }

46};

47

48// Reducer |f(i,j,comp)| -> max, passed DIRECTLY to reduce_max_cell (no wrapping extended

49// lambda, unlike for_each_cell_reduce_max). This is the device-clean path documented

50// in for_each.hpp. Reducer signature (i, j, Real& acc); same Kokkos::Max / same sequential host

51// loop -> bit-identical to the old norm_inf (max and fabs without rounding).

52struct NormInfKernel {

53 ConstArray4 a;

54 int comp;

55 POPS_HD void operator()(int i, int j, Real& acc) const {

56 const Real v = a(i, j, comp);

57 const Real av = v < 0 ? -v : v;

58 if (av > acc)

59 acc = av;

60 }

61};

62

63// Reducer x(i,j,comp) * y(i,j,comp) -> sum, passed DIRECTLY to reduce_sum_cell (no wrapping

64// extended lambda). Device-clean NAMED functor (same recipe as NormInfKernel) for the Krylov

65// solver dot product, pulled from an external TU. Reducer signature (i, j, Real& acc).

66struct DotKernel {

67 ConstArray4 x, y;

68 int comp;

69 POPS_HD void operator()(int i, int j, Real& acc) const { acc += x(i, j, comp) * y(i, j, comp); }

70};

71

72// Reducer f(i,j,comp) -> sum / signed max / signed min over one component. Same device-clean named

73// functor recipe as DotKernel / NormInfKernel (the compiled time Program reductions are first

74// instantiated from a generated problem.so, an external TU). MaxKernel/MinKernel are SIGNED (no fabs,

75// unlike NormInfKernel): they reduce the value itself, the contract of P.max / P.min.

76struct SumKernel {

77 ConstArray4 a;

78 int comp;

79 POPS_HD void operator()(int i, int j, Real& acc) const { acc += a(i, j, comp); }

80};

81struct MaxKernel {

82 ConstArray4 a;

83 int comp;

84 POPS_HD void operator()(int i, int j, Real& acc) const {

85 const Real v = a(i, j, comp);

86 if (v > acc)

87 acc = v;

88 }

89};

90struct MinKernel {

91 ConstArray4 a;

92 int comp;

93 POPS_HD void operator()(int i, int j, Real& acc) const {

94 const Real v = a(i, j, comp);

95 if (v < acc)

96 acc = v;

97 }

98};

99} // namespace detail

100

102inline void saxpy(MultiFab& y, Real a, const MultiFab& x) {

103 const int nc = y.ncomp();

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

105 Array4 Y = y.fab(li).array();

106 const ConstArray4 X = x.fab(li).const_array();

107 const Box2D b = y.fab(li).box();

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

109 for_each_cell(b, detail::SaxpyKernel{Y, X, a, c});

110 }

111}

112

113// Infinity norm over the valid cells of one component. Each local fab is

114// reduced by for_each_cell_reduce_max over |f(i,j,comp)| (true Kokkos reduction,

115// Kokkos::Max), aggregated by host max over the fabs.

116//

117// No more device_fence() up front: under Kokkos parallel_reduce is blocking and

118// absorbs the barrier. EXACT everywhere: max and fabs are without rounding and max

119// is associative/commutative in IEEE754, so bit-identical to the old norm_inf

120// regardless of backend (the reduction order changes no bit).

123inline Real norm_inf(const MultiFab& mf, int comp = 0) {

124 Real m = 0;

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

126 const ConstArray4 a = mf.fab(li).const_array();

127 m = std::max(m, reduce_max_cell(mf.box(li), detail::NormInfKernel{a, comp}));

128 }

129 return m; // MPI all-reduce max later (iso-behavior, not added here)

130}

131

133inline void lincomb(MultiFab& z, Real a, const MultiFab& x, Real b, const MultiFab& y) {

134 const int nc = z.ncomp();

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

136 Array4 Z = z.fab(li).array();

137 const ConstArray4 X = x.fab(li).const_array();

138 const ConstArray4 Y = y.fab(li).const_array();

139 const Box2D bb = z.fab(li).box();

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

141 for_each_cell(bb, detail::LincombKernel{Z, X, Y, a, b, c});

142 }

143}

144

145// Dot product sum_cells x . y over the VALID cells of component comp, reduced over all

146// ranks (all-reduce). Building block of Krylov solvers (BiCGStab: rho, alpha, omega, betas). Each

147// local fab is reduced by reduce_sum_cell (true Kokkos reduction, Kokkos::Sum), the local fabs

148// aggregated by host sum, then all_reduce_sum aggregates the ranks.

149//

150// COLLECTIVE, MANDATORY UNDER MPI: all_reduce_sum is called on EVERY rank, including a rank

151// WITH NO box (local_size()==0, which then contributes 0 to the local sum). Without this call on all

152// ranks, MPI_Allreduce deadlocks (desynchronized collective); the Krylov solver must therefore

153// NEVER short-circuit dot() on an empty rank. In serial all_reduce_sum is the identity.

154//

155// FP NOTE (like sum()): Kokkos::Sum re-associates the sum per tile, so dot is not bit-identical

156// to a lexicographic sum (deterministic/idempotent nonetheless, all Kokkos spaces). Under MPI, the all-reduce

157// returns the SAME value to all ranks (MPI_SUM over one same set of local contributions), so the

158// Krylov stopping criterion triggers at the SAME iteration everywhere (no desynchronization).

163inline Real dot(const MultiFab& x, const MultiFab& y, int comp = 0) {

164 Real s = 0;

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

166 const ConstArray4 X = x.fab(li).const_array();

167 const ConstArray4 Y = y.fab(li).const_array();

168 s += reduce_sum_cell(x.box(li), detail::DotKernel{X, Y, comp});

169 }

170 return static_cast<Real>(all_reduce_sum(static_cast<double>(s)));

171}

172

182inline Real dot_all(const MultiFab& x, const MultiFab& y) {

183 const int nc = x.ncomp();

184 Real s = 0;

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

186 const ConstArray4 X = x.fab(li).const_array();

187 const ConstArray4 Y = y.fab(li).const_array();

188 const Box2D b = x.box(li);

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

190 s += reduce_sum_cell(b, detail::DotKernel{X, Y, c});

191 }

192 return static_cast<Real>(all_reduce_sum(static_cast<double>(s)));

193}

194

198inline Real reduce_sum(const MultiFab& mf, int comp = 0) {

199 Real s = 0;

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

201 const ConstArray4 a = mf.fab(li).const_array();

202 s += reduce_sum_cell(mf.box(li), detail::SumKernel{a, comp});

203 }

204 return static_cast<Real>(all_reduce_sum(static_cast<double>(s)));

205}

206

211inline Real reduce_max(const MultiFab& mf, int comp = 0) {

212 Real m = -std::numeric_limits<Real>::infinity();

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

214 const ConstArray4 a = mf.fab(li).const_array();

215 m = std::max(m, reduce_max_cell(mf.box(li), detail::MaxKernel{a, comp}));

216 }

217 return static_cast<Real>(all_reduce_max(static_cast<double>(m)));

218}

219

223inline Real reduce_min(const MultiFab& mf, int comp = 0) {

224 Real m = std::numeric_limits<Real>::infinity();

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

226 const ConstArray4 a = mf.fab(li).const_array();

227 m = std::min(m, reduce_min_cell(mf.box(li), detail::MinKernel{a, comp}));

228 }

229 return static_cast<Real>(all_reduce_min(static_cast<double>(m)));

230}

231

232} // namespace pops

box2d.hpp

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

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

int ncomp() const

Number of components.

Definition multifab.hpp:60

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

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

multifab.hpp

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

pops

Definition amr_hierarchy.hpp:29

pops::dot_all

Real dot_all(const MultiFab &x, const MultiFab &y)

FULL-component dot Sum_{cells, c} x(.,.,c) * y(.,.,c) over ALL components, reduced over ALL ranks (al...

Definition mf_arith.hpp:182

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

Real dot(const MultiFab &x, const MultiFab &y, int comp=0)

Dot product Sum_cells x.y over component comp, reduced over ALL ranks (all_reduce).

Definition mf_arith.hpp:163

pops::reduce_min_cell

Real reduce_min_cell(const Box2D &b, F f)

Definition for_each.hpp:253

pops::all_reduce_sum

double all_reduce_sum(double x)

Definition comm.hpp:143

pops::saxpy

void saxpy(MultiFab &y, Real a, const MultiFab &x)

y <- y + a x over ALL components of the valid cells. Identical layouts required.

Definition mf_arith.hpp:102

pops::reduce_min

Real reduce_min(const MultiFab &mf, int comp=0)

Signed minimum min_cells f(.,.,comp) over component comp, reduced over ALL ranks (all_reduce_min) – t...

Definition mf_arith.hpp:223

pops::reduce_sum_cell

Real reduce_sum_cell(const Box2D &b, F f)

SUM reduction with a REDUCING FUNCTOR: f receives (i, j, Real& acc) and accumulates,...

Definition for_each.hpp:274

pops::reduce_max_cell

Real reduce_max_cell(const Box2D &b, F f)

MAX reduction with a REDUCING FUNCTOR: f receives (i, j, Real& acc) and updates acc,...

Definition for_each.hpp:240

pops::all_reduce_max

double all_reduce_max(double x)

Definition comm.hpp:146

pops::all_reduce_min

double all_reduce_min(double x)

Definition comm.hpp:149

pops::lincomb

void lincomb(MultiFab &z, Real a, const MultiFab &x, Real b, const MultiFab &y)

z <- a x + b y over ALL components of the valid cells. Identical layouts; aliasing safe.

Definition mf_arith.hpp:133

pops::reduce_max

Real reduce_max(const MultiFab &mf, int comp=0)

Signed maximum max_cells f(.,.,comp) over component comp, reduced over ALL ranks (all_reduce_max) – t...

Definition mf_arith.hpp:211

pops::reduce_sum

Real reduce_sum(const MultiFab &mf, int comp=0)

Sum Sum_cells f(.,.,comp) over component comp, reduced over ALL ranks (all_reduce_sum) – the compiled...

Definition mf_arith.hpp:198

pops::norm_inf

Real norm_inf(const MultiFab &mf, int comp=0)

Infinity norm max |f(.,.,comp)| over the valid cells (LOCAL, without MPI all_reduce).

Definition mf_arith.hpp:123

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

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

Definition fab2d.hpp:44

pops::detail::DotKernel

Definition mf_arith.hpp:66

pops::detail::DotKernel::x

ConstArray4 x

Definition mf_arith.hpp:67

pops::detail::DotKernel::comp

int comp

Definition mf_arith.hpp:68

pops::detail::DotKernel::y

ConstArray4 y

Definition mf_arith.hpp:67

pops::detail::DotKernel::operator()

POPS_HD void operator()(int i, int j, Real &acc) const

Definition mf_arith.hpp:69

pops::detail::LincombKernel

Definition mf_arith.hpp:40

pops::detail::LincombKernel::b

Real b

Definition mf_arith.hpp:43

pops::detail::LincombKernel::c

int c

Definition mf_arith.hpp:44

pops::detail::LincombKernel::operator()

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

Definition mf_arith.hpp:45

pops::detail::LincombKernel::Z

Array4 Z

Definition mf_arith.hpp:41

pops::detail::LincombKernel::Y

ConstArray4 Y

Definition mf_arith.hpp:42

pops::detail::LincombKernel::a

Real a

Definition mf_arith.hpp:43

pops::detail::LincombKernel::X

ConstArray4 X

Definition mf_arith.hpp:42

pops::detail::MaxKernel

Definition mf_arith.hpp:81

pops::detail::MaxKernel::a

ConstArray4 a

Definition mf_arith.hpp:82

pops::detail::MaxKernel::comp

int comp

Definition mf_arith.hpp:83

pops::detail::MaxKernel::operator()

POPS_HD void operator()(int i, int j, Real &acc) const

Definition mf_arith.hpp:84

pops::detail::MinKernel

Definition mf_arith.hpp:90

pops::detail::MinKernel::a

ConstArray4 a

Definition mf_arith.hpp:91

pops::detail::MinKernel::comp

int comp

Definition mf_arith.hpp:92

pops::detail::MinKernel::operator()

POPS_HD void operator()(int i, int j, Real &acc) const

Definition mf_arith.hpp:93

pops::detail::NormInfKernel

Definition mf_arith.hpp:52

pops::detail::NormInfKernel::a

ConstArray4 a

Definition mf_arith.hpp:53

pops::detail::NormInfKernel::comp

int comp

Definition mf_arith.hpp:54

pops::detail::NormInfKernel::operator()

POPS_HD void operator()(int i, int j, Real &acc) const

Definition mf_arith.hpp:55

pops::detail::SaxpyKernel

Definition mf_arith.hpp:32

pops::detail::SaxpyKernel::c

int c

Definition mf_arith.hpp:36

pops::detail::SaxpyKernel::a

Real a

Definition mf_arith.hpp:35

pops::detail::SaxpyKernel::X

ConstArray4 X

Definition mf_arith.hpp:34

pops::detail::SaxpyKernel::operator()

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

Definition mf_arith.hpp:37

pops::detail::SaxpyKernel::Y

Array4 Y

Definition mf_arith.hpp:33

pops::detail::SumKernel

Definition mf_arith.hpp:76

pops::detail::SumKernel::operator()

POPS_HD void operator()(int i, int j, Real &acc) const

Definition mf_arith.hpp:79

pops::detail::SumKernel::a

ConstArray4 a

Definition mf_arith.hpp:77

pops::detail::SumKernel::comp

int comp

Definition mf_arith.hpp:78

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