REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_PC_Init.cpp
Go to the documentation of this file.
1#ifdef REMORA_USE_PARTICLES
2
3#include <REMORA_PC.H>
4#include <AMReX_ParmParse.H>
5
6using namespace amrex;
7
8void REMORAPC::readInputs ()
9{
10 BL_PROFILE("REMORAPC::readInputs");
11
12 ParmParse pp("remora."+m_name);
13
14 m_initialization_type = REMORAParticleInitializations::init_box_uniform;
15 pp.queryAdd("initial_distribution_type", m_initialization_type);
16
17 if (m_initialization_type == REMORAParticleInitializations::init_box_uniform)
18 {
19 Vector<Real> particle_box_lo(AMREX_SPACEDIM);
20 Vector<Real> particle_box_hi(AMREX_SPACEDIM);
21
22 // Defaults
23 for (int i = 0; i < AMREX_SPACEDIM; i++) { particle_box_lo[i] = Geom(0).ProbLo(i); }
24 for (int i = 0; i < AMREX_SPACEDIM; i++) { particle_box_hi[i] = Geom(0).ProbHi(i); }
25
26 pp.queryAdd("particle_box_lo", particle_box_lo, AMREX_SPACEDIM);
27 AMREX_ASSERT(particle_box_lo.size() == AMREX_SPACEDIM);
28
29 pp.queryAdd("particle_box_hi", particle_box_hi, AMREX_SPACEDIM);
30 AMREX_ASSERT(particle_box_hi.size() == AMREX_SPACEDIM);
31
32 m_particle_box.setLo(particle_box_lo);
33 m_particle_box.setHi(particle_box_hi);
34
35 // We default to placing the particles randomly within each cell,
36 // but can override this for regression testing
37 place_randomly_in_cells = true;
38 pp.queryAdd("place_randomly_in_cells", place_randomly_in_cells);
39 }
40
41 m_ppc_init = 1;
42 pp.queryAdd("initial_particles_per_cell", m_ppc_init);
43
44 m_advect_w_flow = (m_name == REMORAParticleNames::tracers ? true : false);
45 pp.queryAdd("advect_with_flow", m_advect_w_flow);
46
47 return;
48}
49
50/*! Initialize particles in domain */
51void REMORAPC::InitializeParticles (const std::unique_ptr<MultiFab>& a_height_ptr)
52{
53 BL_PROFILE("REMORAPC::initializeParticles");
54
55 if (m_initialization_type == REMORAParticleInitializations::init_box_uniform) {
56 initializeParticlesUniformDistributionInBox( a_height_ptr , m_particle_box );
57 } else {
58 Print() << "Error: " << m_initialization_type
59 << " is not a valid initialization for "
60 << m_name << " particle species.\n";
61 Error("See error message!");
62 }
63 return;
64}
65
66/*! Uniform distribution: the number of particles per grid cell is specified
67 * by "initial_particles_per_cell", and they are randomly distributed. */
68void REMORAPC::initializeParticlesUniformDistributionInBox (const std::unique_ptr<MultiFab>& a_height_ptr,
69 const RealBox& particle_init_domain)
70{
71 BL_PROFILE("REMORAPC::initializeParticlesUniformDistributionInBox");
72
73 const int lev = 0;
74 const auto dx = Geom(lev).CellSizeArray();
75 const auto plo = Geom(lev).ProbLoArray();
76
77 int particles_per_cell = m_ppc_init;
78
79 iMultiFab num_particles( ParticleBoxArray(lev),
80 ParticleDistributionMap(lev),
81 1, 0 );
82 num_particles.setVal(0);
83#ifdef _OPENMP
84#pragma omp parallel if (Gpu::notInLaunchRegion())
85#endif
86 for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
87 const Box& tile_box = mfi.tilebox();
88 auto num_particles_arr = num_particles[mfi].array();
89 if (a_height_ptr) {
90 const auto height_arr = (*a_height_ptr)[mfi].array();
91 ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
92 {
93 Real x = plo[0] + (i + Real(0.5))*dx[0];
94 Real y = plo[1] + (j + Real(0.5))*dx[1];
95 Real z = Real(0.125) * (height_arr(i,j ,k ) + height_arr(i+1,j ,k ) +
96 height_arr(i,j+1,k ) + height_arr(i+1,j+1,k ) +
97 height_arr(i,j ,k+1) + height_arr(i+1,j ,k+1) +
98 height_arr(i,j+1,k+1) + height_arr(i+1,j+1,k+1) );
99 if (particle_init_domain.contains(RealVect(x,y,z))) {
100 num_particles_arr(i,j,k) = particles_per_cell;
101 }
102 });
103
104 } else {
105 ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
106 {
107 Real x = plo[0] + (i + Real(0.5))*dx[0];
108 Real y = plo[1] + (j + Real(0.5))*dx[1];
109 Real z = plo[2] + (k + Real(0.5))*dx[2];
110 if (particle_init_domain.contains(RealVect(x,y,z))) {
111 num_particles_arr(i,j,k) = particles_per_cell;
112 }
113 });
114 }
115 }
116
117 iMultiFab offsets( ParticleBoxArray(lev),
118 ParticleDistributionMap(lev),
119 1, 0 );
120 offsets.setVal(0);
121
122 // Define the tiles in serial; the map of tiles is not safe to grow
123 // concurrently in the OpenMP region below. Same as AMReX's addParticles.
124 for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
125 DefineAndReturnParticleTile(lev, mfi);
126 }
127
128#ifdef _OPENMP
129#pragma omp parallel if (Gpu::notInLaunchRegion())
130#endif
131 for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
132 const Box& tile_box = mfi.tilebox();
133
134 // Keep this loop per *tile*: the prefix sum, the offsets it writes, the resize
135 // and the fill must all describe tile_box. Grid-wide quantities here (scanning
136 // the whole num_particles fab, say) agree only while
137 // ParticleContainerBase::do_tiling is false; with tiling on they oversize every
138 // tile and race on the shared offsets fab.
139 const auto num_particles_arr = num_particles[mfi].const_array();
140 auto offset_arr = offsets[mfi].array();
141
142 // atOffset walks tile_box with i fastest, so with one tile per grid the scan
143 // sees exactly the fab's own layout.
144 const int np = Scan::PrefixSum<int>( static_cast<int>(tile_box.numPts()),
145 [=] AMREX_GPU_DEVICE (int idx) -> int {
146 return num_particles_arr(tile_box.atOffset(idx));
147 },
148 [=] AMREX_GPU_DEVICE (int idx, int const &x) {
149 offset_arr(tile_box.atOffset(idx)) = x;
150 },
151 Scan::Type::exclusive,
152 Scan::retSum );
153
154 // already defined in the serial pass above
155 auto& particle_tile = ParticlesAt(lev, mfi);
156 particle_tile.resize(np);
157
158 // Nothing to place: a tile (or, with tiling off, a whole grid) can lie entirely
159 // outside particle_init_domain. Bail before taking &aos[0] on an empty tile.
160 if (np == 0) { continue; }
161
162 auto aos = &particle_tile.GetArrayOfStructs()[0];
163 auto& soa = particle_tile.GetStructOfArrays();
164 auto* vx_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::vx).data();
165 auto* vy_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::vy).data();
166 auto* vz_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::vz).data();
167 auto* mass_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::mass).data();
168
169 auto my_proc = ParallelDescriptor::MyProc();
170 Long pid;
171#ifdef _OPENMP
172#pragma omp critical (remora_particle_nextid)
173#endif
174 {
175 pid = ParticleType::NextID();
176 ParticleType::NextID(pid+np);
177 }
178 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( static_cast<Long>(pid + np) < LastParticleID,
179 "Error: overflow on particle id numbers!" );
180
181 if (a_height_ptr && place_randomly_in_cells) {
182
183 const auto height_arr = (*a_height_ptr)[mfi].array();
184
185 ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k,
186 const RandomEngine& rnd_engine) noexcept
187 {
188 int start = offset_arr(i,j,k);
189 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
190 Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)};
191 Real v[3] = {zero, zero, zero};
192
193 Real x = plo[0] + (i + r[0])*dx[0];
194 Real y = plo[1] + (j + r[1])*dx[1];
195
196 Real sx[] = {one - r[0], r[0]};
197 Real sy[] = {one - r[1], r[1]};
198
199 Real height_at_pxy_lo = zero;
200 for (int ii = 0; ii < 2; ++ii) {
201 for (int jj = 0; jj < 2; ++jj) {
202 height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k);
203 }
204 }
205 Real height_at_pxy_hi = zero;
206 for (int ii = 0; ii < 2; ++ii) {
207 for (int jj = 0; jj < 2; ++jj) {
208 height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1);
209 }
210 }
211
212 Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo);
213
214 auto& p = aos[n];
215 p.id() = pid + n;
216 p.cpu() = my_proc;
217
218 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
219
220 p.idata(REMORAParticlesIntIdxAoS::k) = k;
221
222 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
223
224 mass_ptr[n] = Real(1.0e-6);
225 }
226 });
227
228 } else if (a_height_ptr && !place_randomly_in_cells) {
229
230 const auto height_arr = (*a_height_ptr)[mfi].array();
231
232 ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
233 {
234 int start = offset_arr(i,j,k);
235 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
236 Real r[3] = {Real(0.3), Real(0.7), Real(0.25)};
237 Real v[3] = {zero, zero, zero};
238
239 Real x = plo[0] + (i + r[0])*dx[0];
240 Real y = plo[1] + (j + r[1])*dx[1];
241
242 Real sx[] = { one - r[0], r[0]};
243 Real sy[] = { one - r[1], r[1]};
244
245 Real height_at_pxy_lo = zero;
246 for (int ii = 0; ii < 2; ++ii) {
247 for (int jj = 0; jj < 2; ++jj) {
248 height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k);
249 }
250 }
251 Real height_at_pxy_hi = zero;
252 for (int ii = 0; ii < 2; ++ii) {
253 for (int jj = 0; jj < 2; ++jj) {
254 height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1);
255 }
256 }
257
258 Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo);
259
260 auto& p = aos[n];
261 p.id() = pid + n;
262 p.cpu() = my_proc;
263
264 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
265
266 p.idata(REMORAParticlesIntIdxAoS::k) = k;
267
268 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
269
270 mass_ptr[n] = Real(1.0e-6);
271 }
272 });
273
274 } else if (!a_height_ptr && place_randomly_in_cells) {
275
276 ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k,
277 const RandomEngine& rnd_engine) noexcept
278 {
279 int start = offset_arr(i,j,k);
280 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
281 Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)};
282 Real v[3] = {zero, zero, zero};
283
284 Real x = plo[0] + (i + r[0])*dx[0];
285 Real y = plo[1] + (j + r[1])*dx[1];
286 Real z = plo[2] + (k + r[2])*dx[2];
287
288 auto& p = aos[n];
289 p.id() = pid + n;
290 p.cpu() = my_proc;
291
292 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
293
294 p.idata(REMORAParticlesIntIdxAoS::k) = k;
295
296 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
297
298 mass_ptr[n] = Real(1.0e-6);
299 }
300 });
301
302 } else { // if (!a_height_ptr && !place_randomly_in_cells) {
303
304 ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
305 {
306 int start = offset_arr(i,j,k);
307 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
308 Real r[3] = {Real(0.3), Real(0.7), Real(0.25)};
309 Real v[3] = {zero, zero, zero};
310
311 Real x = plo[0] + (i + r[0])*dx[0];
312 Real y = plo[1] + (j + r[1])*dx[1];
313 Real z = plo[2] + (k + r[2])*dx[2];
314
315 auto& p = aos[n];
316 p.id() = pid + n;
317 p.cpu() = my_proc;
318
319 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
320
321 p.idata(REMORAParticlesIntIdxAoS::k) = k;
322
323 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
324
325 mass_ptr[n] = Real(1.0e-6);
326 }
327 });
328 }
329 }
330
331 return;
332}
333
334#endif
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
Coord::y
@ y
Coord::x
@ x
Coord::z
@ z
setVal
mf_h setVal(geomdata.ProbHi(2))
dx
const Real dx
Definition REMORA_InitAnalyticMasks_DogboneAnalytic.H:6
REMORA_PC.H
REMORAPC::m_initialization_type
std::string m_initialization_type
Definition REMORA_PC.H:216
REMORAPC::m_ppc_init
int m_ppc_init
Definition REMORA_PC.H:217
REMORAPC::m_particle_box
amrex::RealBox m_particle_box
Definition REMORA_PC.H:212
REMORAPC::m_advect_w_flow
bool m_advect_w_flow
Definition REMORA_PC.H:210
REMORAPC::initializeParticlesUniformDistributionInBox
void initializeParticlesUniformDistributionInBox(const std::unique_ptr< amrex::MultiFab > &a_ptr, const amrex::RealBox &particle_box)
Definition REMORA_PC_Init.cpp:68
REMORAPC::InitializeParticles
virtual void InitializeParticles(const std::unique_ptr< amrex::MultiFab > &a_ptr=nullptr)
Definition REMORA_PC_Init.cpp:51
REMORAPC::readInputs
virtual void readInputs()
Definition REMORA_PC_Init.cpp:8
REMORAPC::place_randomly_in_cells
bool place_randomly_in_cells
Definition REMORA_PC.H:224
REMORAPC::m_name
std::string m_name
Definition REMORA_PC.H:214
REMORAParticleInitializations::init_box_uniform
const std::string init_box_uniform
Definition REMORA_PC.H:46
REMORAParticleNames::tracers
const std::string tracers
Definition REMORA_PC.H:51