include/pops/numerics/elliptic/interface/elliptic_problem.hpp Source FileΒΆ

adc_cpp: include/pops/numerics/elliptic/interface/elliptic_problem.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
elliptic_problem.hpp
Go to the documentation of this file.
1#pragma once
2
24
25#include <pops/core/foundation/types.hpp>
26#include <pops/numerics/elliptic/mg/geometric_mg.hpp> // homogeneous(const BCRec&)
27#include <pops/mesh/layout/box_array.hpp>
28#include <pops/mesh/geometry/geometry.hpp>
29#include <pops/mesh/storage/multifab.hpp>
30#include <pops/mesh/boundary/physical_bc.hpp>
31
32#include <stdexcept>
33#include <utility>
34
35namespace pops {
36
37// (a) Elliptic problem: coefficient, physical BC, nullspace.
38struct EllipticProblem {
39 Real eps = 1; // descriptive: current stencil state (eps = 1)
40 BCRec bc{}; // physical BC already propagated
41 bool nullspace_const = false; // solution defined up to an additive constant
42};
43
44// Homogeneous BCRec associated with the problem: delegates to homogeneous(const BCRec&)
45// already in geometric_mg.hpp (multigrid correction with homogeneous BC).
46inline BCRec homogeneous_bc(const EllipticProblem& p) {
47 return homogeneous(p.bc);
48}
49
50// Additive factory: builds an EllipticSolver (GeometricMG, PoissonFFTSolver)
51// from a named EllipticProblem. Delegates to the existing (geom, ba, bc,
52// active, ...) constructor by extracting problem.bc. eps and nullspace_const are
53// descriptive at this stage (eps = 1 already in the stencil; nullspace handled by
54// the bottom-solve + demean on the MG side, by the k=0 mode set to zero on the FFT side):
55// no existing caller is touched, no numerical value changes. Template
56// to avoid an include cycle (this header already includes geometric_mg.hpp).
57template <class Solver, class... Args>
58inline Solver make_elliptic_solver(const Geometry& geom, const BoxArray& ba,
59 const EllipticProblem& problem, Args&&... args) {
60 // Scientific guard: eps is NOT read by the 5-point stencil
61 // (apply_laplacian / poisson_residual / gs_color write lap without a factor,
62 // so eps = 1). Setting eps != 1 would suggest solving a variable-coefficient
63 // Laplacian with no effect on the values: we forbid this trap instead of
64 // ignoring it silently.
65 if (problem.eps != Real(1))
66 throw std::invalid_argument(
67 "EllipticProblem::eps != 1 unsupported (constant-coefficient Laplacian "
68 "operator)");
69 return Solver(geom, ba, problem.bc, std::forward<Args>(args)...);
70}
71
72// (b) Field post-processing: derivation convention E = -grad phi.
73struct FieldPostProcess {
74 enum class GradSign { Plus, Minus };
75 GradSign sign = GradSign::Plus; // +grad (coupler) or -grad (two_fluid)
76 bool store_phi = true; // phi in component 0 (coupler convention)
77};
78
79namespace detail {
80// Derives aux = (phi, s grad phi) from phi (field_postprocess). Named functor
81// (and not an POPS_HD lambda): same reasons as the elliptic path (#93) -- this kernel is
82// first instantiated from an external TU and an extended lambda breaks the device kernel
83// emission under nvcc. Body identical to the former lambda -> bit-identical.
84struct FieldPostprocessKernel {
85 bool store_phi;
86 Array4 a;
87 ConstArray4 p;
88 int gx;
89 Real s, cx, cy;
90 POPS_HD void operator()(int i, int j) const {
91 if (store_phi)
92 a(i, j, 0) = p(i, j);
93 a(i, j, gx) = s * (p(i + 1, j) - p(i - 1, j)) * cx;
94 a(i, j, gx + 1) = s * (p(i, j + 1) - p(i, j - 1)) * cy;
95 }
96};
97} // namespace detail
98
99// Derives the field from phi per spec. The body reproduces EXACTLY
100// detail::coupler_grad_phi (coupler.hpp): same operation order, same
101// multiplicative factor *cx / *cy, with the sign s = +1 (Plus) or -1 (Minus) being the only
102// degree of freedom. cx, cy remain the centered factors already computed by
103// the caller (1/(2 dx), 1/(2 dy)).
104inline void field_postprocess(const MultiFab& phi, MultiFab& out, Real cx, Real cy,
105 FieldPostProcess spec) {
106 const Real s = (spec.sign == FieldPostProcess::GradSign::Plus) ? Real(1) : Real(-1);
107 const bool store_phi = spec.store_phi;
108 for (int li = 0; li < out.local_size(); ++li) {
109 const ConstArray4 p = phi.fab(li).const_array();
110 Array4 a = out.fab(li).array();
111 const Box2D v = out.box(li);
112 const int gx = store_phi ? 1 : 0; // component offset if phi is stored
113 for_each_cell(v, detail::FieldPostprocessKernel{store_phi, a, p, gx, s, cx, cy});
114 }
115}
116
117} // namespace pops
box_array.hpp
BoxArray: the set of boxes tiling a level (disjoint, covering).
pops::BoxArray
Ordered list of boxes tiling a level.
Definition box_array.hpp:22
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::local_size
int local_size() const
Number of fabs OWNED by this rank (bound on local indices).
Definition multifab.hpp:65
geometric_mg.hpp
GeometricMG: in-house geometric multigrid (V-cycle) for the elliptic operator, Gauss-Seidel smoother ...
geometry.hpp
Geometry: index-space (Box2D) <-> Cartesian physical-space mapping; PolarGeometry: SIBLING for a glob...
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::homogeneous
BCRec homogeneous(const BCRec &b)
Definition geometric_mg.hpp:73
pops::homogeneous_bc
BCRec homogeneous_bc(const EllipticProblem &p)
Definition elliptic_problem.hpp:46
pops::field_postprocess
void field_postprocess(const MultiFab &phi, MultiFab &out, Real cx, Real cy, FieldPostProcess spec)
Definition elliptic_problem.hpp:104
pops::make_elliptic_solver
Solver make_elliptic_solver(const Geometry &geom, const BoxArray &ba, const EllipticProblem &problem, Args &&... args)
Definition elliptic_problem.hpp:58
physical_bc.hpp
PHYSICAL boundary conditions at the domain edge (BCType, BCRec, fill_physical_bc, fill_ghosts).
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::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::EllipticProblem
Definition elliptic_problem.hpp:38
pops::EllipticProblem::eps
Real eps
Definition elliptic_problem.hpp:39
pops::EllipticProblem::nullspace_const
bool nullspace_const
Definition elliptic_problem.hpp:41
pops::EllipticProblem::bc
BCRec bc
Definition elliptic_problem.hpp:40
pops::FieldPostProcess
Definition elliptic_problem.hpp:73
pops::FieldPostProcess::store_phi
bool store_phi
Definition elliptic_problem.hpp:76
pops::FieldPostProcess::GradSign
GradSign
Definition elliptic_problem.hpp:74
pops::FieldPostProcess::sign
GradSign sign
Definition elliptic_problem.hpp:75
pops::Geometry
Cartesian geometry of a level: index domain + physical bounds [xlo, xhi] x [ylo, yhi].
Definition geometry.hpp:20
pops::detail::FieldPostprocessKernel
Definition elliptic_problem.hpp:84
pops::detail::FieldPostprocessKernel::gx
int gx
Definition elliptic_problem.hpp:88
pops::detail::FieldPostprocessKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition elliptic_problem.hpp:90
pops::detail::FieldPostprocessKernel::cx
Real cx
Definition elliptic_problem.hpp:89
pops::detail::FieldPostprocessKernel::p
ConstArray4 p
Definition elliptic_problem.hpp:87
pops::detail::FieldPostprocessKernel::s
Real s
Definition elliptic_problem.hpp:89
pops::detail::FieldPostprocessKernel::cy
Real cy
Definition elliptic_problem.hpp:89
pops::detail::FieldPostprocessKernel::store_phi
bool store_phi
Definition elliptic_problem.hpp:85
pops::detail::FieldPostprocessKernel::a
Array4 a
Definition elliptic_problem.hpp:86
types.hpp
Base scalar types and the POPS_HD macro (host+device portability).
POPS_HD
#define POPS_HD
Definition types.hpp:25
Generated by doxygen 1.9.8