include/pops/numerics/spatial/primitives/wave_speed.hpp Source File

adc_cpp: include/pops/numerics/spatial/primitives/wave_speed.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
wave_speed.hpp
Go to the documentation of this file.
1
15
16#pragma once
17
18#include <pops/core/state/state.hpp>
19#include <pops/core/foundation/types.hpp>
20#include <pops/mesh/storage/fab2d.hpp>
21#include <pops/mesh/execution/for_each.hpp> // reduce_max_cell, reduce_min_cell
22#include <pops/mesh/storage/multifab.hpp>
23#include <pops/numerics/spatial/primitives/state_access.hpp> // load_state, load_aux, aux_comps
24#include <pops/parallel/comm.hpp> // all_reduce_max, all_reduce_min
25
26#include <algorithm>
27#include <limits>
28
29namespace pops {
30
31namespace detail {
37template <class Model>
38struct MaxWaveSpeedKernel {
39 Model model;
40 ConstArray4 u, a;
41 POPS_HD void operator()(int i, int j, Real& acc) const {
42 const auto s = load_state<Model>(u, i, j);
43 const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);
44 const Real wx = model.max_wave_speed(s, ax, 0);
45 const Real wy = model.max_wave_speed(s, ax, 1);
46 const Real w = wx > wy ? wx : wy;
47 if (w > acc)
48 acc = w;
49 }
50};
51} // namespace detail
52
59//
60// COLLECTIVE UNDER MPI: we aggregate via all_reduce_max over ALL ranks (same convention as
61// AmrCouplerMp::max_wave_speed and GeometricMG::current_residual). Without this all-reduce, each
62// rank only sees the max of ITS boxes: step_cfl / step_adaptive then choose a DIFFERENT dt per
63// rank (the rank whose local max is lower takes too large a step) and the simulation diverges or
64// desynchronizes the ranks. In serial all_reduce_max is the identity (behavior unchanged).
65template <class Model>
66inline Real max_wave_speed_mf(const Model& model, const MultiFab& U, const MultiFab& aux) {
67 Real m = 0;
68 for (int li = 0; li < U.local_size(); ++li) {
69 const ConstArray4 u = U.fab(li).const_array();
70 const ConstArray4 a = aux.fab(li).const_array();
71 m = std::max(m, reduce_max_cell(U.box(li), detail::MaxWaveSpeedKernel<Model>{model, u, a}));
72 }
73 return static_cast<Real>(all_reduce_max(static_cast<double>(m)));
74}
75
76namespace detail {
82template <class Model>
83struct WaveSpeedMatchKernel {
84 Model model;
85 ConstArray4 u, a;
86 Real target;
87 Real nx; // encoding stride (nx of the DOMAIN, global indices)
88 POPS_HD void operator()(int i, int j, Real& acc) const {
89 const auto s = load_state<Model>(u, i, j);
90 const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);
91 const Real wx = model.max_wave_speed(s, ax, 0);
92 const Real wy = model.max_wave_speed(s, ax, 1);
93 const Real w = wx > wy ? wx : wy;
94 if (w == target) {
95 const Real idx = static_cast<Real>(j) * nx + static_cast<Real>(i);
96 if (idx < acc)
97 acc = idx;
98 }
99 }
100};
101} // namespace detail
102
108template <class Model>
109inline void max_wave_speed_hotspot_mf(const Model& model, const MultiFab& U, const MultiFab& aux,
110 int nx, Real& w_out, int& i_out, int& j_out) {
111 const Real w = max_wave_speed_mf(model, U, aux);
112 Real best = std::numeric_limits<Real>::infinity();
113 for (int li = 0; li < U.local_size(); ++li) {
114 const ConstArray4 u = U.fab(li).const_array();
115 const ConstArray4 a = aux.fab(li).const_array();
116 best = std::min(best, reduce_min_cell(U.box(li), detail::WaveSpeedMatchKernel<Model>{
117 model, u, a, w, static_cast<Real>(nx)}));
118 }
119 best = static_cast<Real>(all_reduce_min(static_cast<double>(best)));
120 w_out = w;
121 // identity of Kokkos::Min = max_real (finite): a rank/box without a cell equaling the max
122 // leaves this value -> we only decode if a REAL index was encoded.
123 if (best >= Real(0) && best < std::numeric_limits<Real>::max() * Real(0.5)) {
124 const long long idx = static_cast<long long>(best);
125 i_out = static_cast<int>(idx % nx);
126 j_out = static_cast<int>(idx / nx);
127 } else { // empty domain / degenerate state: no cell (w may be 0)
128 i_out = -1;
129 j_out = -1;
130 }
131}
132
133// ============================================================================
134// OPTIONAL STEP-BOUND REDUCTIONS (audit 2026-06, step_cfl effort).
135// Counterparts of max_wave_speed_mf for the HasStabilitySpeed / HasSourceFrequency /
136// HasStabilityDt traits (cf. core/physical_model.hpp). Same conventions: reduction via the seam
137// (device under Kokkos), MPI all_reduce (without which each rank would choose a different dt).
138// Instantiated ONLY for a model declaring the trait (if constexpr on the block_builder side):
139// zero codegen, zero cost for a legacy model.
140// ============================================================================
141
142namespace detail {
145template <class Model>
146struct StabilitySpeedKernel {
147 Model model;
148 ConstArray4 u, a;
149 POPS_HD void operator()(int i, int j, Real& acc) const {
150 const auto s = load_state<Model>(u, i, j);
151 const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);
152 const Real wx = model.stability_speed(s, ax, 0);
153 const Real wy = model.stability_speed(s, ax, 1);
154 const Real w = wx > wy ? wx : wy;
155 if (w > acc)
156 acc = w;
157 }
158};
159
161template <class Model>
162struct SourceFrequencyKernel {
163 Model model;
164 ConstArray4 u, a;
165 POPS_HD void operator()(int i, int j, Real& acc) const {
166 const auto s = load_state<Model>(u, i, j);
167 const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);
168 const Real mu = model.source_frequency(s, ax);
169 if (mu > acc)
170 acc = mu;
171 }
172};
173
178template <class Model>
179struct InvStabilityDtKernel {
180 Model model;
181 ConstArray4 u, a;
182 POPS_HD void operator()(int i, int j, Real& acc) const {
183 const auto s = load_state<Model>(u, i, j);
184 const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);
185 const Real db = model.stability_dt(s, ax);
186 if (db > Real(0)) {
187 const Real inv = Real(1) / db;
188 if (inv > acc)
189 acc = inv;
190 }
191 }
192};
193} // namespace detail
194
196template <class Model>
197inline Real max_stability_speed_mf(const Model& model, const MultiFab& U, const MultiFab& aux) {
198 Real m = 0;
199 for (int li = 0; li < U.local_size(); ++li) {
200 const ConstArray4 u = U.fab(li).const_array();
201 const ConstArray4 a = aux.fab(li).const_array();
202 m = std::max(m, reduce_max_cell(U.box(li), detail::StabilitySpeedKernel<Model>{model, u, a}));
203 }
204 return static_cast<Real>(all_reduce_max(static_cast<double>(m)));
205}
206
208template <class Model>
209inline Real max_source_frequency_mf(const Model& model, const MultiFab& U, const MultiFab& aux) {
210 Real m = 0;
211 for (int li = 0; li < U.local_size(); ++li) {
212 const ConstArray4 u = U.fab(li).const_array();
213 const ConstArray4 a = aux.fab(li).const_array();
214 m = std::max(m, reduce_max_cell(U.box(li), detail::SourceFrequencyKernel<Model>{model, u, a}));
215 }
216 return static_cast<Real>(all_reduce_max(static_cast<double>(m)));
217}
218
221template <class Model>
222inline Real min_stability_dt_mf(const Model& model, const MultiFab& U, const MultiFab& aux) {
223 Real inv = 0;
224 for (int li = 0; li < U.local_size(); ++li) {
225 const ConstArray4 u = U.fab(li).const_array();
226 const ConstArray4 a = aux.fab(li).const_array();
227 inv =
228 std::max(inv, reduce_max_cell(U.box(li), detail::InvStabilityDtKernel<Model>{model, u, a}));
229 }
230 inv = static_cast<Real>(all_reduce_max(static_cast<double>(inv)));
231 return inv > Real(0) ? Real(1) / inv : Real(0);
232}
233
234// ============================================================================
235// HLL WAVE SPEED CACHE (OPT-IN -- default path strictly untouched)
236// ============================================================================
237// The HLL flux bounds each face by the signal speeds of BOTH adjacent cells (Davis estimates, cf.
238// hll_speeds). The default path recalls model.wave_speeds per face: for a model whose wave speeds are
239// expensive (moment hierarchy + factorizations at each call), wave_speeds is recomputed several times
240// per cell and per RK stage. This OPT-IN path evaluates wave_speeds ONCE per cell and per direction in
241// a scratch (4 components lo_x/hi_x/lo_y/hi_y), then assembles the residual by reading the scratch: the
242// face speed at i-1/2 becomes min/max of the signed speeds of the two neighbor cells (union over the
243// neighbors, as in the per-cell reference). The scratch is CONSUMED by assemble_rhs_hll_cached
244// (cartesian_operator.hpp).
245
246namespace detail {
250template <class Model>
251struct WaveSpeedCacheKernel {
252 Model model;
253 ConstArray4 u, a;
254 Array4 ws;
255 POPS_HD void operator()(int i, int j) const {
256 const auto s = load_state<Model>(u, i, j);
257 const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);
258 Real lox, hix, loy, hiy;
259 model.wave_speeds(s, ax, 0, lox, hix);
260 model.wave_speeds(s, ax, 1, loy, hiy);
261 ws(i, j, 0) = lox;
262 ws(i, j, 1) = hix;
263 ws(i, j, 2) = loy;
264 ws(i, j, 3) = hiy;
265 }
266};
267} // namespace detail
268
276template <class Model>
277inline void fill_wave_speed_cache(const Model& model, const MultiFab& U, const MultiFab& aux,
278 MultiFab& cache) {
279 for (int li = 0; li < U.local_size(); ++li) {
280 const ConstArray4 u = U.fab(li).const_array();
281 const ConstArray4 a = aux.fab(li).const_array();
282 Array4 ws = cache.fab(li).array();
283 for_each_cell(U.box(li).grow(1), detail::WaveSpeedCacheKernel<Model>{model, u, a, ws});
284 }
285}
286
287} // 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::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::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::max_wave_speed_hotspot_mf
void max_wave_speed_hotspot_mf(const Model &model, const MultiFab &U, const MultiFab &aux, int nx, Real &w_out, int &i_out, int &j_out)
dt_hotspot diagnostic (ADC-182): the cell (GLOBAL indices) that dominates the block's transport CFL b...
Definition wave_speed.hpp:109
pops::Real
double Real
Definition types.hpp:30
pops::reduce_min_cell
Real reduce_min_cell(const Box2D &b, F f)
Definition for_each.hpp:253
pops::min_stability_dt_mf
Real min_stability_dt_mf(const Model &model, const MultiFab &U, const MultiFab &aux)
Global min of the declared admissible step (HasStabilityDt trait), via max(1/dt) (cf.
Definition wave_speed.hpp:222
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::max_source_frequency_mf
Real max_source_frequency_mf(const Model &model, const MultiFab &U, const MultiFab &aux)
Global max of the source frequency (HasSourceFrequency trait). 0 if the source does not constrain.
Definition wave_speed.hpp:209
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::fill_wave_speed_cache
void fill_wave_speed_cache(const Model &model, const MultiFab &U, const MultiFab &aux, MultiFab &cache)
fill_wave_speed_cache: fills the per-cell wave speed scratch (lo_x, hi_x, lo_y, hi_y).
Definition wave_speed.hpp:277
pops::max_stability_speed_mf
Real max_stability_speed_mf(const Model &model, const MultiFab &U, const MultiFab &aux)
Global max of the STABILITY speed (HasStabilitySpeed trait) – counterpart of max_wave_speed_mf.
Definition wave_speed.hpp:197
pops::max_wave_speed_mf
Real max_wave_speed_mf(const Model &model, const MultiFab &U, const MultiFab &aux)
max_wave_speed_mf: global max of the wave speed over the whole MultiFab (CFL).
Definition wave_speed.hpp:66
state.hpp
Pointwise types of the physics layer: StateVec<N> (conserved state) and Aux (auxiliary fields from th...
state_access.hpp
Model/state/aux access layer of the Cartesian spatial operator.
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::Aux
POINTWISE auxiliary fields shared with the physics: single coupling channel.
Definition state.hpp:123
pops::Box2D::grow
Box2D grow(int n) const
Grows the box by n cells in ALL directions (uniform ghost layer).
Definition box2d.hpp:69
pops::ConstArray4
READ-only handle (const counterpart of Array4): same layout and same contract (POD device-copyable,...
Definition fab2d.hpp:44
pops::detail::InvStabilityDtKernel
InvStabilityDtKernel: max over cells of 1/model.stability_dt.
Definition wave_speed.hpp:179
pops::detail::InvStabilityDtKernel::a
ConstArray4 a
Definition wave_speed.hpp:181
pops::detail::InvStabilityDtKernel::operator()
POPS_HD void operator()(int i, int j, Real &acc) const
Definition wave_speed.hpp:182
pops::detail::InvStabilityDtKernel::model
Model model
Definition wave_speed.hpp:180
pops::detail::InvStabilityDtKernel::u
ConstArray4 u
Definition wave_speed.hpp:181
pops::detail::MaxWaveSpeedKernel
MaxWaveSpeedKernel<Model>: device reduction functor for max_wave_speed_mf.
Definition wave_speed.hpp:38
pops::detail::MaxWaveSpeedKernel::a
ConstArray4 a
Definition wave_speed.hpp:40
pops::detail::MaxWaveSpeedKernel::u
ConstArray4 u
Definition wave_speed.hpp:40
pops::detail::MaxWaveSpeedKernel::model
Model model
Definition wave_speed.hpp:39
pops::detail::MaxWaveSpeedKernel::operator()
POPS_HD void operator()(int i, int j, Real &acc) const
Definition wave_speed.hpp:41
pops::detail::SourceFrequencyKernel
SourceFrequencyKernel: max over cells of model.source_frequency (mu >= 0, 1/s).
Definition wave_speed.hpp:162
pops::detail::SourceFrequencyKernel::a
ConstArray4 a
Definition wave_speed.hpp:164
pops::detail::SourceFrequencyKernel::model
Model model
Definition wave_speed.hpp:163
pops::detail::SourceFrequencyKernel::u
ConstArray4 u
Definition wave_speed.hpp:164
pops::detail::SourceFrequencyKernel::operator()
POPS_HD void operator()(int i, int j, Real &acc) const
Definition wave_speed.hpp:165
pops::detail::StabilitySpeedKernel
StabilitySpeedKernel: max over cells/directions of model.stability_speed (replaces MaxWaveSpeedKernel...
Definition wave_speed.hpp:146
pops::detail::StabilitySpeedKernel::operator()
POPS_HD void operator()(int i, int j, Real &acc) const
Definition wave_speed.hpp:149
pops::detail::StabilitySpeedKernel::model
Model model
Definition wave_speed.hpp:147
pops::detail::StabilitySpeedKernel::u
ConstArray4 u
Definition wave_speed.hpp:148
pops::detail::StabilitySpeedKernel::a
ConstArray4 a
Definition wave_speed.hpp:148
pops::detail::WaveSpeedCacheKernel
WaveSpeedCacheKernel: evaluates model.wave_speeds per cell in both directions and stores (lo_x,...
Definition wave_speed.hpp:251
pops::detail::WaveSpeedCacheKernel::ws
Array4 ws
Definition wave_speed.hpp:254
pops::detail::WaveSpeedCacheKernel::a
ConstArray4 a
Definition wave_speed.hpp:253
pops::detail::WaveSpeedCacheKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition wave_speed.hpp:255
pops::detail::WaveSpeedCacheKernel::model
Model model
Definition wave_speed.hpp:252
pops::detail::WaveSpeedCacheKernel::u
ConstArray4 u
Definition wave_speed.hpp:253
pops::detail::WaveSpeedMatchKernel
Locates the cell DOMINATING the CFL (dt_hotspot diagnostic, ADC-182): EQUALITY scan of the recomputed...
Definition wave_speed.hpp:83
pops::detail::WaveSpeedMatchKernel::target
Real target
Definition wave_speed.hpp:86
pops::detail::WaveSpeedMatchKernel::a
ConstArray4 a
Definition wave_speed.hpp:85
pops::detail::WaveSpeedMatchKernel::operator()
POPS_HD void operator()(int i, int j, Real &acc) const
Definition wave_speed.hpp:88
pops::detail::WaveSpeedMatchKernel::nx
Real nx
Definition wave_speed.hpp:87
pops::detail::WaveSpeedMatchKernel::model
Model model
Definition wave_speed.hpp:84
pops::detail::WaveSpeedMatchKernel::u
ConstArray4 u
Definition wave_speed.hpp:85
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