include/pops/coupling/source/coupled_source_program.hpp Source FileΒΆ

adc_cpp: include/pops/coupling/source/coupled_source_program.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
coupled_source_program.hpp
Go to the documentation of this file.
1
13
14#pragma once
15
16#include <pops/core/foundation/types.hpp>
17#include <pops/mesh/storage/fab2d.hpp> // Array4 (POD device-copyable)
18
19#include <cassert>
20#include <cmath>
21
22namespace pops {
23
26enum class CsOp : int {
27 PushReg = 0, // push r[arg] (input or constant)
28 Add = 1, // b = pop; a = pop; push a + b
29 Sub = 2, // b = pop; a = pop; push a - b
30 Mul = 3, // b = pop; a = pop; push a * b
31 Div = 4, // b = pop; a = pop; push a / b
32 Neg = 5, // a = pop; push -a
33 Pow = 6, // b = pop; a = pop; push pow(a, b)
34 Sqrt = 7, // a = pop; push sqrt(a)
35};
36
37// Fixed capacities (bound the stack and the number of registers in the device kernel, where dynamic
38// allocation is forbidden). Generous for realistic coupling formulas; if exceeded, the Python API
39// raises an EXPLICIT error before reaching the device.
40inline constexpr int kCsMaxReg = 32; // inputs + constants
41inline constexpr int kCsMaxStack = 32; // postfix stack depth
42inline constexpr int kCsMaxProg = 256; // length of a program (opcodes)
43inline constexpr int kCsMaxTerms = 16; // source terms (one per .add)
44
47struct CsProgram {
48 int len = 0;
49 int op[kCsMaxProg] = {};
50 int arg[kCsMaxProg] = {};
51
60 POPS_HD Real eval(const Real* reg) const {
61 Real st[kCsMaxStack];
62 int sp = 0;
63 for (int k = 0; k < len; ++k) {
64 switch (static_cast<CsOp>(op[k])) {
65 case CsOp::PushReg:
66 if (sp < kCsMaxStack)
67 st[sp++] = reg[arg[k]];
68 break;
69 case CsOp::Add: {
70 assert(sp >= 2);
71 if (sp < 2)
72 break;
73 Real b = st[--sp];
74 Real a = st[--sp];
75 st[sp++] = a + b;
76 } break;
77 case CsOp::Sub: {
78 assert(sp >= 2);
79 if (sp < 2)
80 break;
81 Real b = st[--sp];
82 Real a = st[--sp];
83 st[sp++] = a - b;
84 } break;
85 case CsOp::Mul: {
86 assert(sp >= 2);
87 if (sp < 2)
88 break;
89 Real b = st[--sp];
90 Real a = st[--sp];
91 st[sp++] = a * b;
92 } break;
93 case CsOp::Div: {
94 assert(sp >= 2);
95 if (sp < 2)
96 break;
97 Real b = st[--sp];
98 Real a = st[--sp];
99 st[sp++] = a / b;
100 } break;
101 case CsOp::Neg: {
102 assert(sp >= 1);
103 if (sp < 1)
104 break;
105 Real a = st[--sp];
106 st[sp++] = -a;
107 } break;
108 case CsOp::Pow: {
109 assert(sp >= 2);
110 if (sp < 2)
111 break;
112 Real b = st[--sp];
113 Real a = st[--sp];
114 st[sp++] = std::pow(a, b);
115 } break;
116 case CsOp::Sqrt: {
117 assert(sp >= 1);
118 if (sp < 1)
119 break;
120 Real a = st[--sp];
121 st[sp++] = std::sqrt(a);
122 } break;
123 }
124 }
125 return sp > 0 ? st[sp - 1] : Real(0);
126 }
127};
128
133struct CoupledSourceKernel {
134 Array4 in[kCsMaxReg]; // input fields (one per (block, role) read); only the first n_in are valid
135 int in_comp[kCsMaxReg];
136 int n_in = 0;
137
138 Real consts[kCsMaxReg]; // constants (parameters), loaded into r[n_in ..]
139 int n_const = 0;
140
141 Array4 out[kCsMaxTerms]; // target of each term (may alias an input: same fab)
142 int out_comp[kCsMaxTerms];
143 CsProgram prog[kCsMaxTerms];
144 int n_terms = 0;
145
146 Real dt = Real(0);
147
148 POPS_HD void operator()(int i, int j) const {
149 Real reg[kCsMaxReg];
150 for (int c = 0; c < n_in; ++c)
151 reg[c] = in[c](i, j, in_comp[c]);
152 for (int c = 0; c < n_const; ++c)
153 reg[n_in + c] = consts[c];
154 // We evaluate ALL terms on the state AT THE START of the step (frozen reg), then write: a term writing
155 // a target that is also an input does not perturb the evaluation of the other terms (explicit additive
156 // splitting consistent, order of the .add irrelevant to the result at 1st order).
157 Real sval[kCsMaxTerms];
158 for (int t = 0; t < n_terms; ++t)
159 sval[t] = prog[t].eval(reg);
160 for (int t = 0; t < n_terms; ++t)
161 out[t](i, j, out_comp[t]) += dt * sval[t];
162 }
163};
164
165// MAX REDUCTION functor of a PER-CELL COUPLED FREQUENCY (CoupledSource.frequency with an Expr,
166// refinement of the declared CONSTANT frequency). mu(i, j) = prog.eval(cell registers): same
167// vocabulary as the source terms ((block, role) input fields + .param() constants). Modeled on
168// CoupledSourceKernel but READ-ONLY (no writes); it reduces the MAX instead of writing outputs.
169// Captures the PODs by VALUE (input Array4, program, constants) -> device-clean (nvcc/Kokkos);
170// the signature (i, j, Real& acc) is the one expected by reduce_max_cell. Registers loaded
171// EXACTLY like the source kernel: r[0 .. n_in-1] = inputs, r[n_in ..] = constants -> the register
172// indices of the program (emitted on the Python side against the SAME table) are consistent. The step
173// bound that derives from it is dt <= cfl / max(mu) (global max aggregated by step_cfl, all_reduce_max
174// over all ranks).
175struct CoupledFreqKernel {
176 Array4 in[kCsMaxReg]; // input fields (one per (block, role) read); only the first n_in are valid
177 int in_comp[kCsMaxReg];
178 int n_in = 0;
179
180 Real consts[kCsMaxReg]; // constants (parameters), loaded into r[n_in ..]
181 int n_const = 0;
182
183 CsProgram prog; // postfix program of the frequency (final top = mu of the cell)
184
185 POPS_HD void operator()(int i, int j, Real& acc) const {
186 Real reg[kCsMaxReg];
187 for (int c = 0; c < n_in; ++c)
188 reg[c] = in[c](i, j, in_comp[c]);
189 for (int c = 0; c < n_const; ++c)
190 reg[n_in + c] = consts[c];
191 const Real mu = prog.eval(reg);
192 if (mu > acc)
193 acc = mu;
194 }
195};
196
197} // namespace pops
fab2d.hpp
Fab2D: single-grid data on a Box2D (in-house equivalent of AMReX's FArrayBox); Array4 / ConstArray4: ...
pops
Definition amr_hierarchy.hpp:29
pops::kCsMaxStack
constexpr int kCsMaxStack
Definition coupled_source_program.hpp:41
pops::kCsMaxTerms
constexpr int kCsMaxTerms
Definition coupled_source_program.hpp:43
pops::Real
double Real
Definition types.hpp:30
pops::kCsMaxReg
constexpr int kCsMaxReg
Definition coupled_source_program.hpp:40
pops::CsOp
CsOp
Opcodes of the postfix stack machine.
Definition coupled_source_program.hpp:26
pops::kCsMaxProg
constexpr int kCsMaxProg
Definition coupled_source_program.hpp:42
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::CoupledFreqKernel
Definition coupled_source_program.hpp:175
pops::CoupledFreqKernel::operator()
POPS_HD void operator()(int i, int j, Real &acc) const
Definition coupled_source_program.hpp:185
pops::CoupledFreqKernel::prog
CsProgram prog
Definition coupled_source_program.hpp:183
pops::CoupledFreqKernel::consts
Real consts[kCsMaxReg]
Definition coupled_source_program.hpp:180
pops::CoupledFreqKernel::in_comp
int in_comp[kCsMaxReg]
Definition coupled_source_program.hpp:177
pops::CoupledFreqKernel::n_const
int n_const
Definition coupled_source_program.hpp:181
pops::CoupledFreqKernel::n_in
int n_in
Definition coupled_source_program.hpp:178
pops::CoupledFreqKernel::in
Array4 in[kCsMaxReg]
Definition coupled_source_program.hpp:176
pops::CoupledSourceKernel
Device functor applying ONE coupled source over a box: captures the PODs by VALUE (input/output Array...
Definition coupled_source_program.hpp:133
pops::CoupledSourceKernel::out
Array4 out[kCsMaxTerms]
Definition coupled_source_program.hpp:141
pops::CoupledSourceKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition coupled_source_program.hpp:148
pops::CoupledSourceKernel::out_comp
int out_comp[kCsMaxTerms]
Definition coupled_source_program.hpp:142
pops::CoupledSourceKernel::n_terms
int n_terms
Definition coupled_source_program.hpp:144
pops::CoupledSourceKernel::dt
Real dt
Definition coupled_source_program.hpp:146
pops::CoupledSourceKernel::consts
Real consts[kCsMaxReg]
Definition coupled_source_program.hpp:138
pops::CoupledSourceKernel::in
Array4 in[kCsMaxReg]
Definition coupled_source_program.hpp:134
pops::CoupledSourceKernel::n_const
int n_const
Definition coupled_source_program.hpp:139
pops::CoupledSourceKernel::n_in
int n_in
Definition coupled_source_program.hpp:136
pops::CoupledSourceKernel::in_comp
int in_comp[kCsMaxReg]
Definition coupled_source_program.hpp:135
pops::CoupledSourceKernel::prog
CsProgram prog[kCsMaxTerms]
Definition coupled_source_program.hpp:143
pops::CsProgram
Fixed-capacity postfix program (POD device-copyable): len opcodes, arg read only by PushReg (register...
Definition coupled_source_program.hpp:47
pops::CsProgram::eval
POPS_HD Real eval(const Real *reg) const
Evaluate the program on the registers reg (loaded for the cell).
Definition coupled_source_program.hpp:60
pops::CsProgram::len
int len
Definition coupled_source_program.hpp:48
pops::CsProgram::arg
int arg[kCsMaxProg]
Definition coupled_source_program.hpp:50
pops::CsProgram::op
int op[kCsMaxProg]
Definition coupled_source_program.hpp:49
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