include/pops/numerics/time/reference/amr_level.hpp Source FileΒΆ

adc_cpp: include/pops/numerics/time/reference/amr_level.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
amr_level.hpp
Go to the documentation of this file.
1#pragma once
2#include <pops/numerics/time/amr/reflux/amr_flux_helpers.hpp>
3#include <pops/amr/hierarchy/refinement_ratio.hpp>
4
25
26namespace pops {
27
28static_assert(kAmrRefRatio == 2, "ratio-2-structural kernels below assume kAmrRefRatio == 2");
29
30namespace detail { // single-box MF oracle: Fab2D -> MF -> MP
31
32// One level of the MultiFab hierarchy. aux held elsewhere (pointer). rC* = region
33// (coords of THIS level) refined by the child; valid if has_fine.
41
42// One conservative 2-level step. Uc: coarse (periodic domain, ghosts for the
43// Limiter). Uf: fine (refined box). auxc/auxf: (phi, grad phi) prescribed, ghosts
44// filled. dt = coarse step; the fine one does r=2 substeps of dt/2 then reflux.
45template <class Limiter = NoSlope, class NumericalFlux = RusanovFlux, class Model>
46void amr_step_2level_mf(const Model& m, MultiFab& Uc, const Box2D& dom, Real dxc, Real dyc,
47 MultiFab& Uf, int CI0, int CI1, int CJ0, int CJ1, const MultiFab& auxc,
48 const MultiFab& auxf, Real dt) {
49 const int r = kAmrRefRatio, nc = Uc.ncomp();
50 const Real dxf = dxc / kAmrRefRatio, dyf = dyc / kAmrRefRatio, dtf = dt / r;
51 const int nJ = CJ1 - CJ0 + 1, nI = CI1 - CI0 + 1;
52 MultiFab Uc_old = Uc; // coarse state at time t (temporal interp of fine ghosts)
53
54 // --- coarse fluxes at the 4 faces of the fine region (before update) ---
55 fill_boundary(Uc, dom, Periodicity{true, true});
56 MultiFab fxc(BoxArray(std::vector<Box2D>{xface_box(Uc.box(0))}), Uc.dmap(), nc, 0);
57 MultiFab fyc(BoxArray(std::vector<Box2D>{yface_box(Uc.box(0))}), Uc.dmap(), nc, 0);
58 compute_face_fluxes<Limiter, NumericalFlux>(m, Uc, auxc, fxc, fyc, dxc, dyc);
59 std::vector<Real> cL(nJ * nc), cR(nJ * nc), cB(nI * nc), cT(nI * nc);
60 {
61 device_fence();
62 const ConstArray4 FX = fxc.fab(0).const_array(), FY = fyc.fab(0).const_array();
63 for (int J = CJ0; J <= CJ1; ++J)
64 for (int k = 0; k < nc; ++k) {
65 cL[(J - CJ0) * nc + k] = FX(CI0, J, k);
66 cR[(J - CJ0) * nc + k] = FX(CI1 + 1, J, k);
67 }
68 for (int I = CI0; I <= CI1; ++I)
69 for (int k = 0; k < nc; ++k) {
70 cB[(I - CI0) * nc + k] = FY(I, CJ0, k);
71 cT[(I - CI0) * nc + k] = FY(I, CJ1 + 1, k);
72 }
73 }
74 mf_advance_faces(Uc, fxc, fyc, dxc, dyc, dt); // Uc -> state t+dt
75 mf_apply_source(m, Uc, auxc, dt); // source S(U,aux) at the substep
76
77 // --- fine subcycling: r substeps, accumulation of fine fluxes (x dtf) ---
78 std::vector<Real> fL(nJ * nc, 0), fR(nJ * nc, 0), fB(nI * nc, 0), fT(nI * nc, 0);
79 MultiFab fxf(BoxArray(std::vector<Box2D>{xface_box(Uf.box(0))}), Uf.dmap(), nc, 0);
80 MultiFab fyf(BoxArray(std::vector<Box2D>{yface_box(Uf.box(0))}), Uf.dmap(), nc, 0);
81 for (int s = 0; s < r; ++s) {
82 mf_fill_fine_ghosts_t(Uf, Uc_old, Uc, Real(s) / r);
83 compute_face_fluxes<Limiter, NumericalFlux>(m, Uf, auxf, fxf, fyf, dxf, dyf);
84 device_fence();
85 const ConstArray4 FX = fxf.fab(0).const_array(), FY = fyf.fab(0).const_array();
86 for (int J = CJ0; J <= CJ1; ++J)
87 for (int k = 0; k < nc; ++k) {
88 fL[(J - CJ0) * nc + k] +=
89 Real(0.5) * (FX(2 * CI0, 2 * J, k) + FX(2 * CI0, 2 * J + 1, k)) * dtf;
90 fR[(J - CJ0) * nc + k] +=
91 Real(0.5) * (FX(2 * CI1 + 2, 2 * J, k) + FX(2 * CI1 + 2, 2 * J + 1, k)) * dtf;
92 }
93 for (int I = CI0; I <= CI1; ++I)
94 for (int k = 0; k < nc; ++k) {
95 fB[(I - CI0) * nc + k] +=
96 Real(0.5) * (FY(2 * I, 2 * CJ0, k) + FY(2 * I + 1, 2 * CJ0, k)) * dtf;
97 fT[(I - CI0) * nc + k] +=
98 Real(0.5) * (FY(2 * I, 2 * CJ1 + 2, k) + FY(2 * I + 1, 2 * CJ1 + 2, k)) * dtf;
99 }
100 mf_advance_faces(Uf, fxf, fyf, dxf, dyf, dtf);
101 mf_apply_source(m, Uf, auxf, dtf); // source S(U,aux) at the substep
102 }
103
104 mf_average_down(Uf, Uc, CI0, CI1, CJ0, CJ1); // sync of covered cells
105
106 // --- reflux: coarse flux (x dt) replaced by sum of fine fluxes (x dtf) ---
107 device_fence();
108 Array4 c = Uc.fab(0).array();
109 for (int J = CJ0; J <= CJ1; ++J)
110 for (int k = 0; k < nc; ++k) {
111 c(CI0 - 1, J, k) -= (fL[(J - CJ0) * nc + k] - cL[(J - CJ0) * nc + k] * dt) / dxc;
112 c(CI1 + 1, J, k) += (fR[(J - CJ0) * nc + k] - cR[(J - CJ0) * nc + k] * dt) / dxc;
113 }
114 for (int I = CI0; I <= CI1; ++I)
115 for (int k = 0; k < nc; ++k) {
116 c(I, CJ0 - 1, k) -= (fB[(I - CI0) * nc + k] - cB[(I - CI0) * nc + k] * dt) / dyc;
117 c(I, CJ1 + 1, k) += (fT[(I - CI0) * nc + k] - cT[(I - CI0) * nc + k] * dt) / dyc;
118 }
119}
120
121// --- N-level recursion (MultiFab counterpart of amr_multilevel.hpp) ---
122
123// Recursively advances level lev by dt (Berger-Oliger subcycling r=2 + reflux).
124// pOld/pNew = parent states bounding the step; preg* = parent register (fine fluxes of
125// THIS level), null if lev==0. Generic (Limiter, NumericalFlux, N comp), GPU seam.
126template <class Limiter = NoSlope, class NumericalFlux = RusanovFlux, class Model>
127void subcycle_level_mf(const Model& m, std::vector<AmrLevelMF>& L, int lev, Real dt,
128 const Box2D& dom, const MultiFab* pOld, const MultiFab* pNew, Real frac,
129 std::vector<Real>* pregL, std::vector<Real>* pregR, std::vector<Real>* pregB,
130 std::vector<Real>* pregT) {
131 const int r = kAmrRefRatio;
132 AmrLevelMF& lv = L[lev];
133 const int nc = lv.U.ncomp();
134
135 if (lev == 0)
136 fill_boundary(lv.U, dom, Periodicity{true, true});
137 else
138 mf_fill_fine_ghosts_t(lv.U, *pOld, *pNew, frac);
139
140 MultiFab fx(BoxArray(std::vector<Box2D>{xface_box(lv.U.box(0))}), lv.U.dmap(), nc, 0);
141 MultiFab fy(BoxArray(std::vector<Box2D>{yface_box(lv.U.box(0))}), lv.U.dmap(), nc, 0);
142 compute_face_fluxes<Limiter, NumericalFlux>(m, lv.U, *lv.aux, fx, fy, lv.dx, lv.dy);
143
144 if (lev > 0) { // contribution to the parent register (fine fluxes x dt)
145 device_fence();
146 const ConstArray4 FX = fx.fab(0).const_array(), FY = fy.fab(0).const_array();
147 const AmrLevelMF& par = L[lev - 1];
148 const int pI0 = par.rCI0, pI1 = par.rCI1, pJ0 = par.rCJ0, pJ1 = par.rCJ1;
149 for (int J = pJ0; J <= pJ1; ++J)
150 for (int k = 0; k < nc; ++k) {
151 (*pregL)[(J - pJ0) * nc + k] +=
152 Real(0.5) * (FX(2 * pI0, 2 * J, k) + FX(2 * pI0, 2 * J + 1, k)) * dt;
153 (*pregR)[(J - pJ0) * nc + k] +=
154 Real(0.5) * (FX(2 * pI1 + 2, 2 * J, k) + FX(2 * pI1 + 2, 2 * J + 1, k)) * dt;
155 }
156 for (int I = pI0; I <= pI1; ++I)
157 for (int k = 0; k < nc; ++k) {
158 (*pregB)[(I - pI0) * nc + k] +=
159 Real(0.5) * (FY(2 * I, 2 * pJ0, k) + FY(2 * I + 1, 2 * pJ0, k)) * dt;
160 (*pregT)[(I - pI0) * nc + k] +=
161 Real(0.5) * (FY(2 * I, 2 * pJ1 + 2, k) + FY(2 * I + 1, 2 * pJ1 + 2, k)) * dt;
162 }
163 }
164
165 if (!lv.has_fine) {
166 mf_advance_faces(lv.U, fx, fy, lv.dx, lv.dy, dt); // leaf
167 mf_apply_source(m, lv.U, *lv.aux, dt); // source S(U,aux) at the substep
168 return;
169 }
170
171 // level with a child: local register + coarse flux saved (without dt)
172 const int cI0 = lv.rCI0, cI1 = lv.rCI1, cJ0 = lv.rCJ0, cJ1 = lv.rCJ1;
173 const int nI = cI1 - cI0 + 1, nJ = cJ1 - cJ0 + 1;
174 std::vector<Real> cL(nJ * nc), cR(nJ * nc), cB(nI * nc), cT(nI * nc);
175 {
176 device_fence();
177 const ConstArray4 FX = fx.fab(0).const_array(), FY = fy.fab(0).const_array();
178 for (int J = cJ0; J <= cJ1; ++J)
179 for (int k = 0; k < nc; ++k) {
180 cL[(J - cJ0) * nc + k] = FX(cI0, J, k);
181 cR[(J - cJ0) * nc + k] = FX(cI1 + 1, J, k);
182 }
183 for (int I = cI0; I <= cI1; ++I)
184 for (int k = 0; k < nc; ++k) {
185 cB[(I - cI0) * nc + k] = FY(I, cJ0, k);
186 cT[(I - cI0) * nc + k] = FY(I, cJ1 + 1, k);
187 }
188 }
189 std::vector<Real> fL(nJ * nc, 0), fR(nJ * nc, 0), fB(nI * nc, 0), fT(nI * nc, 0);
190
191 MultiFab U_old = lv.U; // state t (temporal interp of the child)
192 mf_advance_faces(lv.U, fx, fy, lv.dx, lv.dy, dt); // lv.U -> t+dt
193 mf_apply_source(m, lv.U, *lv.aux, dt); // source S(U,aux) at the substep
194
195 for (int s = 0; s < r; ++s)
196 subcycle_level_mf<Limiter, NumericalFlux>(m, L, lev + 1, dt / r, dom, &U_old, &lv.U,
197 Real(s) / r, &fL, &fR, &fB, &fT);
198
199 mf_average_down(L[lev + 1].U, lv.U, cI0, cI1, cJ0, cJ1);
200
201 device_fence();
202 Array4 c = lv.U.fab(0).array();
203 for (int J = cJ0; J <= cJ1; ++J)
204 for (int k = 0; k < nc; ++k) {
205 c(cI0 - 1, J, k) -= (fL[(J - cJ0) * nc + k] - cL[(J - cJ0) * nc + k] * dt) / lv.dx;
206 c(cI1 + 1, J, k) += (fR[(J - cJ0) * nc + k] - cR[(J - cJ0) * nc + k] * dt) / lv.dx;
207 }
208 for (int I = cI0; I <= cI1; ++I)
209 for (int k = 0; k < nc; ++k) {
210 c(I, cJ0 - 1, k) -= (fB[(I - cI0) * nc + k] - cB[(I - cI0) * nc + k] * dt) / lv.dy;
211 c(I, cJ1 + 1, k) += (fT[(I - cI0) * nc + k] - cT[(I - cI0) * nc + k] * dt) / lv.dy;
212 }
213}
214
215// Driver: one dt step of the full hierarchy (level 0 = coarse).
216template <class Limiter = NoSlope, class NumericalFlux = RusanovFlux, class Model>
217void amr_step_multilevel_mf(const Model& m, std::vector<AmrLevelMF>& L, const Box2D& dom, Real dt) {
218 subcycle_level_mf<Limiter, NumericalFlux>(m, L, 0, dt, dom, nullptr, nullptr, Real(0), nullptr,
219 nullptr, nullptr, nullptr);
220}
221
222} // namespace detail
223
224} // namespace pops
amr_flux_helpers.hpp
Basic MultiFab building blocks of an AMR step: AmrTimeMethod enum, device-clean functors and advance ...
pops::BoxArray
Ordered list of boxes tiling a level.
Definition box_array.hpp:22
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::dmap
const DistributionMapping & dmap() const
GLOBAL distribution (owner rank per box).
Definition multifab.hpp:58
pops::MultiFab::box
const Box2D & box(int li) const
VALID box of local fab li.
Definition multifab.hpp:71
pops::detail::amr_step_2level_mf
void amr_step_2level_mf(const Model &m, MultiFab &Uc, const Box2D &dom, Real dxc, Real dyc, MultiFab &Uf, int CI0, int CI1, int CJ0, int CJ1, const MultiFab &auxc, const MultiFab &auxf, Real dt)
Definition amr_level.hpp:46
pops::detail::amr_step_multilevel_mf
void amr_step_multilevel_mf(const Model &m, std::vector< AmrLevelMF > &L, const Box2D &dom, Real dt)
Definition amr_level.hpp:217
pops::detail::subcycle_level_mf
void subcycle_level_mf(const Model &m, std::vector< AmrLevelMF > &L, int lev, Real dt, const Box2D &dom, const MultiFab *pOld, const MultiFab *pNew, Real frac, std::vector< Real > *pregL, std::vector< Real > *pregR, std::vector< Real > *pregB, std::vector< Real > *pregT)
Definition amr_level.hpp:127
pops
Definition amr_hierarchy.hpp:29
pops::mf_apply_source
void mf_apply_source(const Model &m, MultiFab &U, const MultiFab &aux, Real dt)
Definition amr_flux_helpers.hpp:128
pops::Real
double Real
Definition types.hpp:30
pops::mf_fill_fine_ghosts_t
void mf_fill_fine_ghosts_t(MultiFab &Uf, const MultiFab &Uc_old, const MultiFab &Uc_new, Real frac, Real pos_floor=Real(0), int pos_comp=0)
Definition amr_flux_helpers.hpp:215
pops::mf_average_down
void mf_average_down(const MultiFab &Uf, MultiFab &Uc, int CI0, int CI1, int CJ0, int CJ1)
Definition amr_flux_helpers.hpp:181
pops::xface_box
Box2D xface_box(const Box2D &v)
xface_box / yface_box: face boxes normal to x (resp.
Definition face_flux.hpp:147
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::fill_boundary
void fill_boundary(MultiFab &mf, const Box2D &domain, Periodicity per={})
BLOCKING halo exchange: begin then end immediately (no overlap).
Definition fill_boundary.hpp:333
pops::yface_box
Box2D yface_box(const Box2D &v)
Definition face_flux.hpp:150
pops::kAmrRefRatio
constexpr int kAmrRefRatio
The native AMR refinement ratio between two consecutive levels.
Definition refinement_ratio.hpp:30
pops::mf_advance_faces
void mf_advance_faces(MultiFab &U, const MultiFab &Fx, const MultiFab &Fy, Real dx, Real dy, Real dt)
Definition amr_flux_helpers.hpp:96
refinement_ratio.hpp
Single source of truth for the native AMR refinement ratio.
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::Periodicity
Per-direction periodicity: halo wrapping in x and/or y during the exchange (false = open edge,...
Definition fill_boundary.hpp:37
pops::detail::AmrLevelMF
Definition amr_level.hpp:34
pops::detail::AmrLevelMF::dy
Real dy
Definition amr_level.hpp:37
pops::detail::AmrLevelMF::has_fine
bool has_fine
Definition amr_level.hpp:39
pops::detail::AmrLevelMF::rCI0
int rCI0
Definition amr_level.hpp:38
pops::detail::AmrLevelMF::rCI1
int rCI1
Definition amr_level.hpp:38
pops::detail::AmrLevelMF::U
MultiFab U
Definition amr_level.hpp:35
pops::detail::AmrLevelMF::rCJ1
int rCJ1
Definition amr_level.hpp:38
pops::detail::AmrLevelMF::dx
Real dx
Definition amr_level.hpp:37
pops::detail::AmrLevelMF::aux
const MultiFab * aux
Definition amr_level.hpp:36
pops::detail::AmrLevelMF::rCJ0
int rCJ0
Definition amr_level.hpp:38
Generated by doxygen 1.9.8