include/pops/coupling/base/elliptic_rhs.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
• coupling
• base

elliptic_rhs.hpp

Go to the documentation of this file.

1#pragma once

2

3#include <pops/core/model/coupled_system.hpp>

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

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

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

7#include <pops/numerics/spatial_operator.hpp> // load_state

8

9#include <cstddef>

10#include <stdexcept>

11#include <vector>

12

22

23namespace pops {

24

25namespace detail {

26

28template <class Model>

29struct SingleModelEllipticRhsKernel {

30 Model m;

31 ConstArray4 s;

32 Array4 f;

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

34 f(i, j, 0) = m.elliptic_rhs(load_state<Model>(s, i, j));

35 }

36};

37

38} // namespace detail

39

42template <class Model>

43struct SingleModelEllipticRhs {

44 Model model;

45

47 void operator()(const MultiFab& state, MultiFab& rhs) const {

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

49 const ConstArray4 s = state.fab(li).const_array();

50 Array4 f = rhs.fab(li).array();

51 const Box2D v = rhs.box(li);

52 const Model m = model;

53 for_each_cell(v, detail::SingleModelEllipticRhsKernel<Model>{m, s, f});

54 }

55 }

56};

57

58namespace detail {

59

61struct TwoFieldChargeDensityRhsKernel {

62 Array4 r;

63 ConstArray4 u0, u1;

64 Real a0, a1;

65 int c0, c1;

66 POPS_HD void operator()(int i, int j) const { r(i, j, 0) = a0 * u0(i, j, c0) + a1 * u1(i, j, c1); }

67};

68

69} // namespace detail

70

72struct TwoFieldChargeDensityRhs {

73 Real q0 = Real(1);

74 Real q1 = Real(-1);

75 int comp0 = 0;

76 int comp1 = 0;

77

79 void operator()(const MultiFab& U0, const MultiFab& U1, MultiFab& rhs) const {

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

81 const ConstArray4 u0 = U0.fab(li).const_array();

82 const ConstArray4 u1 = U1.fab(li).const_array();

83 Array4 r = rhs.fab(li).array();

84 const Box2D b = rhs.box(li);

85 const Real a0 = q0, a1 = q1;

86 const int c0 = comp0, c1 = comp1;

87 for_each_cell(b, detail::TwoFieldChargeDensityRhsKernel{r, u0, u1, a0, a1, c0, c1});

88 }

89 }

90};

91

94struct TwoBlockChargeDensityRhs {

95 Real q0 = Real(1);

96 Real q1 = Real(-1);

97 int comp0 = 0;

98 int comp1 = 0;

99

101 template <CoupledSystemLike System>

102 void operator()(const System& system, MultiFab& rhs) const {

103 TwoFieldChargeDensityRhs two;

104 two.q0 = q0;

105 two.q1 = q1;

106 two.comp0 = comp0;

107 two.comp1 = comp1;

108 two(system.template block<0>().U(), system.template block<1>().U(), rhs);

109 }

110};

111

116struct SpeciesCharge {

117 Real charge = Real(0);

118 int comp = 0;

119};

120

121namespace detail {

122

124struct AddScaledComponentKernel {

125 Array4 r;

126 ConstArray4 u;

127 Real a;

128 int c;

129 POPS_HD void operator()(int i, int j) const { r(i, j, 0) += a * u(i, j, c); }

130};

131

132} // namespace detail

133

136inline void add_scaled_component(const MultiFab& U, Real q, int comp, MultiFab& rhs) {

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

138 const ConstArray4 u = U.fab(li).const_array();

139 Array4 r = rhs.fab(li).array();

140 const Box2D b = rhs.box(li);

141 const Real a = q;

142 const int c = comp;

143 for_each_cell(b, detail::AddScaledComponentKernel{r, u, a, c});

144 }

145}

146

153struct ChargeDensityRhs {

154 std::vector<SpeciesCharge> species;

155

157 template <CoupledSystemLike System>

158 void operator()(const System& system, MultiFab& rhs) const {

159 if (species.size() != System::n_blocks)

160 throw std::runtime_error(

161 "ChargeDensityRhs: exactly one SpeciesCharge per block is required "

162 "(a neutral species is declared with charge = 0)");

163 rhs.set_val(Real(0));

164 std::size_t k = 0;

165 system.for_each_block([&](const auto& block) {

166 const SpeciesCharge sc = species[k++];

167 add_scaled_component(block.U(), sc.charge, sc.comp, rhs);

168 });

169 }

170};

171

172} // namespace pops

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

const Box2D & box(int li) const

VALID box of local fab li.

Definition multifab.hpp:71

pops::MultiFab::set_val

void set_val(Real v)

Fills all cells (valid + ghosts) of every local fab with v.

Definition multifab.hpp:85

pops::MultiFab::local_size

int local_size() const

Number of fabs OWNED by this rank (bound on local indices).

Definition multifab.hpp:65

pops::System

Coupled multi-species system, composed at runtime from generic bricks.

Definition system.hpp:89

pops::System::n_blocks

POPS_EXPORT int n_blocks() const

Number of blocks (species) installed.

coupled_system.hpp

CoupledSystem: heterogeneous collection of equation blocks (multi-species, multi-scheme).

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

void add_scaled_component(const MultiFab &U, Real q, int comp, MultiFab &rhs)

rhs(.,.,0) += q * U(.,.,comp) over the valid cells.

Definition elliptic_rhs.hpp:136

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

spatial_operator.hpp

Cartesian spatial operator: assembles R(U, aux) = -div F + S over the cells of a 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::ChargeDensityRhs

N-species Poisson RHS: f = Sum_s q_s n_s over ALL blocks of the system.

Definition elliptic_rhs.hpp:153

pops::ChargeDensityRhs::operator()

void operator()(const System &system, MultiFab &rhs) const

rhs <- 0 then += q_s n_s for each block. Throws if species.size() != n_blocks.

Definition elliptic_rhs.hpp:158

pops::ChargeDensityRhs::species

std::vector< SpeciesCharge > species

Definition elliptic_rhs.hpp:154

pops::ConstArray4

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

Definition fab2d.hpp:44

pops::SingleModelEllipticRhs

SINGLE-model RHS assembler: rhs(.,.,0) = model.elliptic_rhs(U) over the valid cells.

Definition elliptic_rhs.hpp:43

pops::SingleModelEllipticRhs::operator()

void operator()(const MultiFab &state, MultiFab &rhs) const

rhs <- elliptic_rhs(state) cell by cell (identical layouts between state and rhs).

Definition elliptic_rhs.hpp:47

pops::SingleModelEllipticRhs::model

Model model

Definition elliptic_rhs.hpp:44

pops::SpeciesCharge

Charge (with sign) and density component of a species for the elliptic RHS assembly.

Definition elliptic_rhs.hpp:116

pops::SpeciesCharge::comp

int comp

Definition elliptic_rhs.hpp:118

pops::SpeciesCharge::charge

Real charge

Definition elliptic_rhs.hpp:117

pops::TwoBlockChargeDensityRhs

Two-block RHS: same computation as TwoFieldChargeDensityRhs but reads blocks 0 and 1 of a CoupledSyst...

Definition elliptic_rhs.hpp:94

pops::TwoBlockChargeDensityRhs::operator()

void operator()(const System &system, MultiFab &rhs) const

rhs <- q0 n_0 + q1 n_1 from blocks 0 and 1 of the system.

Definition elliptic_rhs.hpp:102

pops::TwoBlockChargeDensityRhs::q1

Real q1

Definition elliptic_rhs.hpp:96

pops::TwoBlockChargeDensityRhs::comp0

int comp0

Definition elliptic_rhs.hpp:97

pops::TwoBlockChargeDensityRhs::comp1

int comp1

Definition elliptic_rhs.hpp:98

pops::TwoBlockChargeDensityRhs::q0

Real q0

Definition elliptic_rhs.hpp:95

pops::TwoFieldChargeDensityRhs

Two-field RHS: rhs = q0 * U0(.,.,comp0) + q1 * U1(.,.,comp1) (two-species charge density).

Definition elliptic_rhs.hpp:72

pops::TwoFieldChargeDensityRhs::q0

Real q0

Definition elliptic_rhs.hpp:73

pops::TwoFieldChargeDensityRhs::comp0

int comp0

Definition elliptic_rhs.hpp:75

pops::TwoFieldChargeDensityRhs::comp1

int comp1

Definition elliptic_rhs.hpp:76

pops::TwoFieldChargeDensityRhs::q1

Real q1

Definition elliptic_rhs.hpp:74

pops::TwoFieldChargeDensityRhs::operator()

void operator()(const MultiFab &U0, const MultiFab &U1, MultiFab &rhs) const

rhs <- q0 U0 + q1 U1 over the valid cells (identical layouts).

Definition elliptic_rhs.hpp:79

pops::detail::AddScaledComponentKernel

NAMED functor (not an POPS_HD lambda) of the accumulation: r(i,j,0) += a * u(i,j,c).

Definition elliptic_rhs.hpp:124

pops::detail::AddScaledComponentKernel::c

int c

Definition elliptic_rhs.hpp:128

pops::detail::AddScaledComponentKernel::operator()

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

Definition elliptic_rhs.hpp:129

pops::detail::AddScaledComponentKernel::a

Real a

Definition elliptic_rhs.hpp:127

pops::detail::AddScaledComponentKernel::r

Array4 r

Definition elliptic_rhs.hpp:125

pops::detail::AddScaledComponentKernel::u

ConstArray4 u

Definition elliptic_rhs.hpp:126

pops::detail::SingleModelEllipticRhsKernel

NAMED functor (not an POPS_HD lambda) of the single-model RHS: f(i,j,0) = model.elliptic_rhs(U).

Definition elliptic_rhs.hpp:29

pops::detail::SingleModelEllipticRhsKernel::operator()

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

Definition elliptic_rhs.hpp:33

pops::detail::SingleModelEllipticRhsKernel::f

Array4 f

Definition elliptic_rhs.hpp:32

pops::detail::SingleModelEllipticRhsKernel::m

Model m

Definition elliptic_rhs.hpp:30

pops::detail::SingleModelEllipticRhsKernel::s

ConstArray4 s

Definition elliptic_rhs.hpp:31

pops::detail::TwoFieldChargeDensityRhsKernel

NAMED functor (not an POPS_HD lambda) of the two-field RHS: r(i,j,0) = a0 u0 + a1 u1.

Definition elliptic_rhs.hpp:61

pops::detail::TwoFieldChargeDensityRhsKernel::c1

int c1

Definition elliptic_rhs.hpp:65

pops::detail::TwoFieldChargeDensityRhsKernel::u0

ConstArray4 u0

Definition elliptic_rhs.hpp:63

pops::detail::TwoFieldChargeDensityRhsKernel::c0

int c0

Definition elliptic_rhs.hpp:65

pops::detail::TwoFieldChargeDensityRhsKernel::operator()

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

Definition elliptic_rhs.hpp:66

pops::detail::TwoFieldChargeDensityRhsKernel::a0

Real a0

Definition elliptic_rhs.hpp:64

pops::detail::TwoFieldChargeDensityRhsKernel::r

Array4 r

Definition elliptic_rhs.hpp:62

pops::detail::TwoFieldChargeDensityRhsKernel::a1

Real a1

Definition elliptic_rhs.hpp:64

pops::detail::TwoFieldChargeDensityRhsKernel::u1

ConstArray4 u1

Definition elliptic_rhs.hpp:63

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