include/pops/coupling/schur/core/schur_condensation.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
• schur
• core

schur_condensation.hpp

Go to the documentation of this file.

1#pragma once

2

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

4#include <pops/core/state/variables.hpp>

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

6#include <pops/mesh/geometry/geometry.hpp>

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

8#include <pops/mesh/boundary/physical_bc.hpp>

9#include <pops/numerics/elliptic/poisson/poisson_operator.hpp> // apply_laplacian (Lap phi^n)

10#include <pops/numerics/linalg/lorentz_eliminator.hpp> // closed 2x2 B^{-1}

11

12#include <stdexcept>

13

55

56namespace pops {

57

58namespace detail {

59

67struct SchurOperatorCoeffKernel {

68 ConstArray4 s;

69 ConstArray4 bz;

70 Array4 ex, ey;

71 Array4 axy, ayx;

72 Real c;

73 Real th_dt;

74 int c_rho;

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

76 const Real rho = s(i, j, c_rho);

77 // Closed B^{-1} (2x2 rotation-dilation): only w = theta dt B_z enters binv. We reconstruct

78 // the eliminator via (theta=th_dt, dt=1, B_z) to obtain w = th_dt * B_z (same binv).

79 const LorentzEliminator le(th_dt, Real(1), bz(i, j, 0));

80 const Real cr = c * rho; // c rho: common factor of the 4 entries of A - I

81 ex(i, j, 0) = Real(1) + cr * le.binv_11();

82 ey(i, j, 0) = Real(1) + cr * le.binv_22();

83 axy(i, j, 0) = cr * le.binv_12();

84 ayx(i, j, 0) = cr * le.binv_21();

85 }

86};

87

91struct SchurExplicitFluxKernel {

92 ConstArray4 s;

93 ConstArray4 bz;

94 Array4 fx, fy;

95 Real th_dt;

96 int c_mx, c_my;

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

98 const LorentzEliminator le(th_dt, Real(1), bz(i, j, 0));

99 Real Fx, Fy;

100 le.apply_Binv(s(i, j, c_mx), s(i, j, c_my), Fx, Fy); // B^{-1} (mx, my) = rho B^{-1} v

101 fx(i, j, 0) = Fx;

102 fy(i, j, 0) = Fy;

103 }

104};

105

112struct SchurRhsAssembleKernel {

113 ConstArray4 neg_lap;

114 ConstArray4 fx, fy;

115 Array4 rhs;

116 Real g;

117 Real half_idx, half_idy;

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

119 const Real divF = (fx(i + 1, j, 0) - fx(i - 1, j, 0)) * half_idx +

120 (fy(i, j + 1, 0) - fy(i, j - 1, 0)) * half_idy;

121 rhs(i, j, 0) = neg_lap(i, j, 0) - g * divF;

122 }

123};

124

126struct NegateKernel {

127 ConstArray4 src;

128 Array4 dst;

129 POPS_HD void operator()(int i, int j) const { dst(i, j, 0) = -src(i, j, 0); }

130};

131

132} // namespace detail

133

141struct SchurCondensationOperator {

142 MultiFab eps_x;

143 MultiFab eps_y;

144 MultiFab a_xy;

145 MultiFab a_yx;

146 MultiFab rhs;

147};

148

153class ElectrostaticLorentzCondensation {

154 public:

160 ElectrostaticLorentzCondensation(const VariableSet& vars, Real alpha, Real theta, Real dt)

161 : c_rho_(vars.index_of(VariableRole::Density)),

162 c_mx_(vars.index_of(VariableRole::MomentumX)),

163 c_my_(vars.index_of(VariableRole::MomentumY)),

164 alpha_(alpha),

165 theta_(theta),

166 dt_(dt) {

167 if (c_rho_ < 0 || c_mx_ < 0 || c_my_ < 0)

168 throw std::runtime_error(

169 "ElectrostaticLorentzCondensation: the fluid block must expose the roles Density, "

170 "MomentumX and MomentumY (VariableSet.roles populated).");

171 }

172

173 Real c_coeff() const { return theta_ * theta_ * dt_ * dt_ * alpha_; }

174 int density_comp() const { return c_rho_; }

175 int momentum_x_comp() const { return c_mx_; }

176 int momentum_y_comp() const { return c_my_; }

177

185 void assemble_operator(const MultiFab& state, const MultiFab& bz, const Geometry& geom,

186 const BCRec& bc, MultiFab& eps_x, MultiFab& eps_y, MultiFab& a_xy,

187 MultiFab& a_yx) const {

188 const Real c = c_coeff();

189 const Real th_dt = theta_ * dt_;

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

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

192 const ConstArray4 b = bz.fab(li).const_array();

193 Array4 ex = eps_x.fab(li).array();

194 Array4 ey = eps_y.fab(li).array();

195 Array4 fxy = a_xy.fab(li).array();

196 Array4 fyx = a_yx.fab(li).array();

197 for_each_cell(eps_x.box(li),

198 detail::SchurOperatorCoeffKernel{s, b, ex, ey, fxy, fyx, c, th_dt, c_rho_});

199 }

200 // coefficient ghosts: the operator face mean reads the neighbor at +-1 (box boundary AND

201 // physical boundary). We extend by zero-gradient on the physical boundaries (Foextrap), periodic otherwise

202 // (eps_bc of GeometricMG). a_xy/a_yx: face arithmetic mean -> same extension.

203 const BCRec ebc = coeff_bc(bc);

204 fill_ghosts(eps_x, geom.domain, ebc);

205 fill_ghosts(eps_y, geom.domain, ebc);

206 fill_ghosts(a_xy, geom.domain, ebc);

207 fill_ghosts(a_yx, geom.domain, ebc);

208 }

209

216 void assemble_rhs(MultiFab& phi_n, const MultiFab& state, const MultiFab& bz,

217 const Geometry& geom, const BCRec& bc, MultiFab& rhs) const {

218 const Real th_dt = theta_ * dt_;

219 const Real g =

220 theta_ * dt_ * alpha_; // theta dt alpha (coefficient of the div(rho B^{-1} v) term)

221 const BoxArray& ba = rhs.box_array();

222 const DistributionMapping& dm = rhs.dmap();

223

224 // 1) -Lap phi^n: apply_laplacian WITHOUT coefficient = div(grad phi^n) (canonical Laplacian), negated.

225 // We first fill the ghosts of phi^n (physical BC) so that the boundary stencil is correct.

226 device_fence();

227 fill_ghosts(phi_n, geom.domain, bc);

228 MultiFab lap(ba, dm, 1, 0);

229 apply_laplacian(phi_n, geom, lap); // lap = Lap phi^n (A=I, kappa=0)

230 MultiFab neg_lap(ba, dm, 1, 0);

231 for (int li = 0; li < neg_lap.local_size(); ++li)

232 for_each_cell(neg_lap.box(li),

233 detail::NegateKernel{lap.fab(li).const_array(), neg_lap.fab(li).array()});

234

235 // 2) EXPLICIT flux F = rho B^{-1} v^n = B^{-1} (mx, my), at the center (1 ghost for the centered div).

236 MultiFab fx(ba, dm, 1, 1), fy(ba, dm, 1, 1);

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

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

239 const ConstArray4 b = bz.fab(li).const_array();

240 for_each_cell(fx.box(li),

241 detail::SchurExplicitFluxKernel{s, b, fx.fab(li).array(), fy.fab(li).array(),

242 th_dt, c_mx_, c_my_});

243 }

244 // ghosts of F: the centered divergence reads Fx(i+-1), Fy(j+-1) (box boundary AND physical boundary).

245 const BCRec ebc = coeff_bc(bc);

246 fill_ghosts(fx, geom.domain, ebc);

247 fill_ghosts(fy, geom.domain, ebc);

248

249 // 3) rhs = -Lap phi^n - g div F (second-order centered divergence).

250 const Real half_idx = Real(1) / (Real(2) * geom.dx());

251 const Real half_idy = Real(1) / (Real(2) * geom.dy());

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

253 for_each_cell(rhs.box(li),

254 detail::SchurRhsAssembleKernel{

255 neg_lap.fab(li).const_array(), fx.fab(li).const_array(),

256 fy.fab(li).const_array(), rhs.fab(li).array(), g, half_idx, half_idy});

257 }

258

262 SchurCondensationOperator assemble(MultiFab& phi_n, const MultiFab& state, const MultiFab& bz,

263 const Geometry& geom, const BCRec& bc) const {

264 const BoxArray& ba = state.box_array();

265 const DistributionMapping& dm = state.dmap();

266 SchurCondensationOperator op{MultiFab(ba, dm, 1, 1), MultiFab(ba, dm, 1, 1),

267 MultiFab(ba, dm, 1, 1), MultiFab(ba, dm, 1, 1),

268 MultiFab(ba, dm, 1, 0)};

269 assemble_operator(state, bz, geom, bc, op.eps_x, op.eps_y, op.a_xy, op.a_yx);

270 assemble_rhs(phi_n, state, bz, geom, bc, op.rhs);

271 return op;

272 }

273

274 private:

278 static BCRec coeff_bc(const BCRec& bc) {

279 auto fo = [](BCType t) { return t == BCType::Periodic ? t : BCType::Foextrap; };

280 BCRec b;

281 b.xlo = fo(bc.xlo);

282 b.xhi = fo(bc.xhi);

283 b.ylo = fo(bc.ylo);

284 b.yhi = fo(bc.yhi);

285 return b;

286 }

287

288 int c_rho_, c_mx_, c_my_;

289 Real alpha_, theta_, dt_;

290};

291

292} // namespace pops

pops::BoxArray

Ordered list of boxes tiling a level.

Definition box_array.hpp:22

pops::DistributionMapping

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

Definition distribution_mapping.hpp:19

pops::ElectrostaticLorentzCondensation

GENERIC builder of the condensed source stage for the electrostatic + Lorentz source (kind="electrost...

Definition schur_condensation.hpp:153

pops::ElectrostaticLorentzCondensation::assemble_operator

void assemble_operator(const MultiFab &state, const MultiFab &bz, const Geometry &geom, const BCRec &bc, MultiFab &eps_x, MultiFab &eps_y, MultiFab &a_xy, MultiFab &a_yx) const

Assembles ONLY the coefficients of the tensor operator A_op = I + c rho B^{-1} into MultiFab eps_x/ep...

Definition schur_condensation.hpp:185

pops::ElectrostaticLorentzCondensation::density_comp

int density_comp() const

Definition schur_condensation.hpp:174

pops::ElectrostaticLorentzCondensation::assemble

SchurCondensationOperator assemble(MultiFab &phi_n, const MultiFab &state, const MultiFab &bz, const Geometry &geom, const BCRec &bc) const

COMPLETE assembly (operator + RHS) into one SchurCondensationOperator object allocated on the layout ...

Definition schur_condensation.hpp:262

pops::ElectrostaticLorentzCondensation::c_coeff

Real c_coeff() const

c = theta^2 dt^2 alpha

Definition schur_condensation.hpp:173

pops::ElectrostaticLorentzCondensation::momentum_y_comp

int momentum_y_comp() const

Definition schur_condensation.hpp:176

pops::ElectrostaticLorentzCondensation::ElectrostaticLorentzCondensation

ElectrostaticLorentzCondensation(const VariableSet &vars, Real alpha, Real theta, Real dt)

vars: descriptor of the fluid block; MUST expose the roles Density / MomentumX / MomentumY.

Definition schur_condensation.hpp:160

pops::ElectrostaticLorentzCondensation::assemble_rhs

void assemble_rhs(MultiFab &phi_n, const MultiFab &state, const MultiFab &bz, const Geometry &geom, const BCRec &bc, MultiFab &rhs) const

Assembles ONLY the condensed right-hand side -Lap phi^n - theta dt alpha div(rho B^{-1} v^n).

Definition schur_condensation.hpp:216

pops::ElectrostaticLorentzCondensation::momentum_x_comp

int momentum_x_comp() const

Definition schur_condensation.hpp:175

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

for_each.hpp

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

geometry.hpp

Geometry: index-space (Box2D) <-> Cartesian physical-space mapping; PolarGeometry: SIBLING for a glob...

lorentz_eliminator.hpp

LorentzEliminator: 2x2 operator B of the Schur scheme for implicit velocity elimination.

multifab.hpp

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

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

void apply_laplacian(const MultiFab &phi, const Geometry &geom, MultiFab &lap, const MultiFab *coef=nullptr, const MultiFab *eps=nullptr, const MultiFab *kappa=nullptr, const MultiFab *eps_y=nullptr, const MultiFab *a_xy=nullptr, const MultiFab *a_yx=nullptr)

Definition poisson_operator.hpp:191

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

VariableRole

PHYSICAL role of a component.

Definition variables.hpp:27

pops::VariableRole::Density

@ Density

pops::VariableRole::MomentumX

@ MomentumX

pops::VariableRole::MomentumY

@ MomentumY

pops::BCType

BCType

Boundary condition type for a face: Periodic (handled by fill_boundary), Foextrap (zero gradient,...

Definition physical_bc.hpp:25

pops::BCType::Foextrap

@ Foextrap

pops::BCType::Periodic

@ Periodic

pops::fill_ghosts

void fill_ghosts(MultiFab &mf, const Box2D &domain, const BCRec &bc)

COMPLETE ghost filling: fill_boundary (interior + periodic, periodicity deduced from bc) THEN fill_ph...

Definition physical_bc.hpp:227

physical_bc.hpp

PHYSICAL boundary conditions at the domain edge (BCType, BCRec, fill_physical_bc, fill_ghosts).

poisson_operator.hpp

Free functions of the elliptic operator: apply_laplacian (matvec), poisson_residual (residual),...

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

Boundary conditions for the FOUR faces of the domain (type + associated Dirichlet value).

Definition physical_bc.hpp:29

pops::BCRec::yhi

BCType yhi

Definition physical_bc.hpp:31

pops::BCRec::xlo

BCType xlo

Definition physical_bc.hpp:30

pops::BCRec::ylo

BCType ylo

Definition physical_bc.hpp:31

pops::BCRec::xhi

BCType xhi

Definition physical_bc.hpp:30

pops::ConstArray4

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

Definition fab2d.hpp:44

pops::Geometry

Cartesian geometry of a level: index domain + physical bounds [xlo, xhi] x [ylo, yhi].

Definition geometry.hpp:20

pops::Geometry::dy

POPS_HD Real dy() const

Grid spacing in y (= (yhi - ylo) / domain.ny()). POPS_HD.

Definition geometry.hpp:33

pops::Geometry::dx

POPS_HD Real dx() const

Grid spacing in x (= (xhi - xlo) / domain.nx()). POPS_HD.

Definition geometry.hpp:31

pops::Geometry::domain

Box2D domain

Definition geometry.hpp:21

pops::LorentzEliminator

LorentzEliminator: operator B = [[1,-w],[w,1]] and its analytic inverse.

Definition lorentz_eliminator.hpp:55

pops::LorentzEliminator::binv_11

POPS_HD Real binv_11() const

Definition lorentz_eliminator.hpp:81

pops::LorentzEliminator::binv_21

POPS_HD Real binv_21() const

Definition lorentz_eliminator.hpp:83

pops::LorentzEliminator::binv_12

POPS_HD Real binv_12() const

Definition lorentz_eliminator.hpp:82

pops::LorentzEliminator::apply_Binv

POPS_HD void apply_Binv(Real vx, Real vy, Real &vxp, Real &vyp) const

apply_Binv: applies B^{-1} = (1/det)*[[1,w],[-w,1]] to (vx, vy), writes (vxp, vyp)....

Definition lorentz_eliminator.hpp:73

pops::LorentzEliminator::binv_22

POPS_HD Real binv_22() const

Definition lorentz_eliminator.hpp:84

pops::SchurCondensationOperator

Result of the Schur assembly: the coefficient MultiFab of the tensor operator A_op and the condensed ...

Definition schur_condensation.hpp:141

pops::SchurCondensationOperator::a_xy

MultiFab a_xy

cross A_xy = c rho binv_12

Definition schur_condensation.hpp:144

pops::SchurCondensationOperator::a_yx

MultiFab a_yx

cross A_yx = c rho binv_21

Definition schur_condensation.hpp:145

pops::SchurCondensationOperator::eps_y

MultiFab eps_y

diagonal A_yy = 1 + c rho binv_22

Definition schur_condensation.hpp:143

pops::SchurCondensationOperator::eps_x

MultiFab eps_x

diagonal A_xx = 1 + c rho binv_11

Definition schur_condensation.hpp:142

pops::SchurCondensationOperator::rhs

MultiFab rhs

condensed right-hand side -Lap phi^n - theta dt alpha div(rho B^{-1} v^n)

Definition schur_condensation.hpp:146

pops::VariableSet

A model's variable set: kind (cons/prim), names, size, canonical roles (optional, parallel to names; ...

Definition variables.hpp:58

pops::detail::NegateKernel

neg(i,j) = -src(i,j) (negation of component 0). Device-clean NAMED functor.

Definition schur_condensation.hpp:126

pops::detail::NegateKernel::src

ConstArray4 src

Definition schur_condensation.hpp:127

pops::detail::NegateKernel::operator()

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

Definition schur_condensation.hpp:129

pops::detail::NegateKernel::dst

Array4 dst

Definition schur_condensation.hpp:128

pops::detail::SchurExplicitFluxKernel

EXPLICIT flux F = rho B^{-1} v^n (v = (mx,my)/rho) ASSEMBLED per cell, WRITTEN at the center into fx/...

Definition schur_condensation.hpp:91

pops::detail::SchurExplicitFluxKernel::th_dt

Real th_dt

theta * dt (w = th_dt * B_z)

Definition schur_condensation.hpp:95

pops::detail::SchurExplicitFluxKernel::fx

Array4 fx

Definition schur_condensation.hpp:94

pops::detail::SchurExplicitFluxKernel::operator()

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

Definition schur_condensation.hpp:97

pops::detail::SchurExplicitFluxKernel::c_mx

int c_mx

Definition schur_condensation.hpp:96

pops::detail::SchurExplicitFluxKernel::bz

ConstArray4 bz

B_z field at the center.

Definition schur_condensation.hpp:93

pops::detail::SchurExplicitFluxKernel::c_my

int c_my

MomentumX / MomentumY components.

Definition schur_condensation.hpp:96

pops::detail::SchurExplicitFluxKernel::fy

Array4 fy

output: components of the flux F = B^{-1} (mx, my)

Definition schur_condensation.hpp:94

pops::detail::SchurExplicitFluxKernel::s

ConstArray4 s

fluid state (mx, my read at components c_mx, c_my)

Definition schur_condensation.hpp:92

pops::detail::SchurOperatorCoeffKernel

Coefficients of the tensor operator A_op = I + c rho B^{-1} ASSEMBLED per cell from the fluid state a...

Definition schur_condensation.hpp:67

pops::detail::SchurOperatorCoeffKernel::s

ConstArray4 s

fluid state

Definition schur_condensation.hpp:68

pops::detail::SchurOperatorCoeffKernel::operator()

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

Definition schur_condensation.hpp:75

pops::detail::SchurOperatorCoeffKernel::c_rho

int c_rho

Density component.

Definition schur_condensation.hpp:74

pops::detail::SchurOperatorCoeffKernel::ey

Array4 ey

output: eps_x, eps_y (diagonal of A)

Definition schur_condensation.hpp:70

pops::detail::SchurOperatorCoeffKernel::axy

Array4 axy

Definition schur_condensation.hpp:71

pops::detail::SchurOperatorCoeffKernel::ex

Array4 ex

Definition schur_condensation.hpp:70

pops::detail::SchurOperatorCoeffKernel::th_dt

Real th_dt

theta * dt (for LorentzEliminator: w = th_dt * B_z, and binv depends only on w)

Definition schur_condensation.hpp:73

pops::detail::SchurOperatorCoeffKernel::ayx

Array4 ayx

output: cross terms a_xy, a_yx

Definition schur_condensation.hpp:71

pops::detail::SchurOperatorCoeffKernel::c

Real c

c = theta^2 dt^2 alpha

Definition schur_condensation.hpp:72

pops::detail::SchurOperatorCoeffKernel::bz

ConstArray4 bz

B_z field at the center.

Definition schur_condensation.hpp:69

pops::detail::SchurRhsAssembleKernel

rhs(i,j) = lap(i,j) (= -Lap phi^n, already negated by the caller) - g * div F, second-order centered ...

Definition schur_condensation.hpp:112

pops::detail::SchurRhsAssembleKernel::neg_lap

ConstArray4 neg_lap

-Lap phi^n (already negated)

Definition schur_condensation.hpp:113

pops::detail::SchurRhsAssembleKernel::half_idy

Real half_idy

1/(2 dx), 1/(2 dy)

Definition schur_condensation.hpp:117

pops::detail::SchurRhsAssembleKernel::g

Real g

theta dt alpha

Definition schur_condensation.hpp:116

pops::detail::SchurRhsAssembleKernel::operator()

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

Definition schur_condensation.hpp:118

pops::detail::SchurRhsAssembleKernel::half_idx

Real half_idx

Definition schur_condensation.hpp:117

pops::detail::SchurRhsAssembleKernel::fy

ConstArray4 fy

flux F at the center (ghosts filled)

Definition schur_condensation.hpp:114

pops::detail::SchurRhsAssembleKernel::rhs

Array4 rhs

output: condensed right-hand side

Definition schur_condensation.hpp:115

pops::detail::SchurRhsAssembleKernel::fx

ConstArray4 fx

Definition schur_condensation.hpp:114

types.hpp

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

POPS_HD

#define POPS_HD

Definition types.hpp:25

variables.hpp

Descriptor of a model's variables (Vars).

---

Generated by  1.9.8