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 int np = 0;
135 {
136 int ncell = num_particles[mfi].numPts();
137 const int* in = num_particles[mfi].dataPtr();
138 int* out = offsets[mfi].dataPtr();
139 np = Scan::PrefixSum<int>( ncell,
140 [=] AMREX_GPU_DEVICE (int i) -> int { return in[i]; },
141 [=] AMREX_GPU_DEVICE (int i, int const &x) { out[i] = x; },
142 Scan::Type::exclusive,
143 Scan::retSum );
144 }
145 auto offset_arr = offsets[mfi].array();
146
147 // already defined in the serial pass above
148 auto& particle_tile = ParticlesAt(lev, mfi);
149 particle_tile.resize(np);
150 auto aos = &particle_tile.GetArrayOfStructs()[0];
151 auto& soa = particle_tile.GetStructOfArrays();
152 auto* vx_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::vx).data();
153 auto* vy_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::vy).data();
154 auto* vz_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::vz).data();
155 auto* mass_ptr = soa.GetRealData(REMORAParticlesRealIdxSoA::mass).data();
156
157 const auto num_particles_arr = num_particles[mfi].array();
158
159 auto my_proc = ParallelDescriptor::MyProc();
160 Long pid;
161#ifdef _OPENMP
162#pragma omp critical (remora_particle_nextid)
163#endif
164 {
165 pid = ParticleType::NextID();
166 ParticleType::NextID(pid+np);
167 }
168 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( static_cast<Long>(pid + np) < LastParticleID,
169 "Error: overflow on particle id numbers!" );
170
171 if (a_height_ptr && place_randomly_in_cells) {
172
173 const auto height_arr = (*a_height_ptr)[mfi].array();
174
175 ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k,
176 const RandomEngine& rnd_engine) noexcept
177 {
178 int start = offset_arr(i,j,k);
179 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
180 Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)};
181 Real v[3] = {zero, zero, zero};
182
183 Real x = plo[0] + (i + r[0])*dx[0];
184 Real y = plo[1] + (j + r[1])*dx[1];
185
186 Real sx[] = {one - r[0], r[0]};
187 Real sy[] = {one - r[1], r[1]};
188
189 Real height_at_pxy_lo = zero;
190 for (int ii = 0; ii < 2; ++ii) {
191 for (int jj = 0; jj < 2; ++jj) {
192 height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k);
193 }
194 }
195 Real height_at_pxy_hi = zero;
196 for (int ii = 0; ii < 2; ++ii) {
197 for (int jj = 0; jj < 2; ++jj) {
198 height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1);
199 }
200 }
201
202 Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo);
203
204 auto& p = aos[n];
205 p.id() = pid + n;
206 p.cpu() = my_proc;
207
208 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
209
210 p.idata(REMORAParticlesIntIdxAoS::k) = k;
211
212 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
213
214 mass_ptr[n] = Real(1.0e-6);
215 }
216 });
217
218 } else if (a_height_ptr && !place_randomly_in_cells) {
219
220 const auto height_arr = (*a_height_ptr)[mfi].array();
221
222 ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
223 {
224 int start = offset_arr(i,j,k);
225 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
226 Real r[3] = {Real(0.3), Real(0.7), Real(0.25)};
227 Real v[3] = {zero, zero, zero};
228
229 Real x = plo[0] + (i + r[0])*dx[0];
230 Real y = plo[1] + (j + r[1])*dx[1];
231
232 Real sx[] = { one - r[0], r[0]};
233 Real sy[] = { one - r[1], r[1]};
234
235 Real height_at_pxy_lo = zero;
236 for (int ii = 0; ii < 2; ++ii) {
237 for (int jj = 0; jj < 2; ++jj) {
238 height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k);
239 }
240 }
241 Real height_at_pxy_hi = zero;
242 for (int ii = 0; ii < 2; ++ii) {
243 for (int jj = 0; jj < 2; ++jj) {
244 height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1);
245 }
246 }
247
248 Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo);
249
250 auto& p = aos[n];
251 p.id() = pid + n;
252 p.cpu() = my_proc;
253
254 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
255
256 p.idata(REMORAParticlesIntIdxAoS::k) = k;
257
258 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
259
260 mass_ptr[n] = Real(1.0e-6);
261 }
262 });
263
264 } else if (!a_height_ptr && place_randomly_in_cells) {
265
266 ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k,
267 const RandomEngine& rnd_engine) noexcept
268 {
269 int start = offset_arr(i,j,k);
270 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
271 Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)};
272 Real v[3] = {zero, zero, zero};
273
274 Real x = plo[0] + (i + r[0])*dx[0];
275 Real y = plo[1] + (j + r[1])*dx[1];
276 Real z = plo[2] + (k + r[2])*dx[2];
277
278 auto& p = aos[n];
279 p.id() = pid + n;
280 p.cpu() = my_proc;
281
282 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
283
284 p.idata(REMORAParticlesIntIdxAoS::k) = k;
285
286 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
287
288 mass_ptr[n] = Real(1.0e-6);
289 }
290 });
291
292 } else { // if (!a_height_ptr && !place_randomly_in_cells) {
293
294 ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
295 {
296 int start = offset_arr(i,j,k);
297 for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
298 Real r[3] = {Real(0.3), Real(0.7), Real(0.25)};
299 Real v[3] = {zero, zero, zero};
300
301 Real x = plo[0] + (i + r[0])*dx[0];
302 Real y = plo[1] + (j + r[1])*dx[1];
303 Real z = plo[2] + (k + r[2])*dx[2];
304
305 auto& p = aos[n];
306 p.id() = pid + n;
307 p.cpu() = my_proc;
308
309 p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;
310
311 p.idata(REMORAParticlesIntIdxAoS::k) = k;
312
313 vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];
314
315 mass_ptr[n] = Real(1.0e-6);
316 }
317 });
318 }
319 }
320
321 return;
322}
323
324#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 auto dx
Definition REMORA_InitAnalyticSMFlux_DoubleGyre.H:3
REMORA_PC.H
mod_param::n
integer, dimension(ngrids) n
Definition REMORA_fennel_roms.F:79