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