WarpX
DefaultInitialization.H
Go to the documentation of this file.
1 /* Copyright 2019-2020 Andrew Myers, Axel Huebl,
2  * Maxence Thevenet
3  *
4  * This file is part of WarpX.
5  *
6  * License: BSD-3-Clause-LBNL
7  */
8 #ifndef WARPX_DEFAULTINITIALIZATION_H_
9 #define WARPX_DEFAULTINITIALIZATION_H_
10 
11 #include <WarpX.H>
12 #ifdef WARPX_QED
15 #endif
16 
17 #include <AMReX_GpuContainers.H>
18 #include <AMReX_REAL.H>
19 
20 #include <cmath>
21 #include <map>
22 #include <string>
23 
40 
45 static std::map<std::string, InitializationPolicy> initialization_policies = {
50 #ifdef WARPX_DIM_RZ
51  {"theta", InitializationPolicy::Zero},
52 #endif
53 
54 #ifdef WARPX_QED
55  {"opticalDepthBW", InitializationPolicy::RandomExp},
56  {"opticalDepthQSR", InitializationPolicy::RandomExp}
57 #endif
58 
59 };
60 
62 amrex::ParticleReal initializeRealValue (const InitializationPolicy policy, amrex::RandomEngine const& engine) noexcept
63 {
64  switch (policy) {
65  case InitializationPolicy::Zero : return 0.0;
66  case InitializationPolicy::One : return 1.0;
68  return -std::log(amrex::Random(engine));
69  }
70  default : {
71  amrex::Abort("Initialization Policy not recognized");
72  return 1.0;
73  }
74  }
75 }
76 
78 int initializeIntValue (const InitializationPolicy policy) noexcept
79 {
80  switch (policy) {
81  case InitializationPolicy::Zero : return 0;
82  case InitializationPolicy::One : return 1;
83  default : {
84  amrex::Abort("Initialization Policy not recognized");
85  return 1;
86  }
87  }
88 }
89 
90 namespace ParticleCreation {
91 
117 template <typename PTile>
119  const int n_external_attr_real,
120  const int n_external_attr_int,
121  const std::vector<std::string>& user_real_attribs,
122  const std::vector<std::string>& user_int_attribs,
123  const std::map<std::string, int>& particle_comps,
124  const std::map<std::string, int>& particle_icomps,
125  const std::vector<amrex::Parser*>& user_real_attrib_parser,
126  const std::vector<amrex::Parser*>& user_int_attrib_parser,
127 #ifdef WARPX_QED
128  const bool do_qed_comps,
129  BreitWheelerEngine* p_bw_engine,
130  QuantumSynchrotronEngine* p_qs_engine,
131 #endif
132  const int ionization_initial_level,
133  int start, int stop)
134 {
135  using namespace amrex::literals;
136 
137  // Preparing data needed for user defined attributes
138  const auto n_user_real_attribs = static_cast<int>(user_real_attribs.size());
139  const auto n_user_int_attribs = static_cast<int>(user_int_attribs.size());
140  const auto get_position = GetParticlePosition<PIdx>(ptile);
141  const auto soa = ptile.getParticleTileData();
142  const amrex::ParticleReal* AMREX_RESTRICT ux = soa.m_rdata[PIdx::ux];
143  const amrex::ParticleReal* AMREX_RESTRICT uy = soa.m_rdata[PIdx::uy];
144  const amrex::ParticleReal* AMREX_RESTRICT uz = soa.m_rdata[PIdx::uz];
145  constexpr int lev = 0;
146  const amrex::Real t = WarpX::GetInstance().gett_new(lev);
147 
148  // Initialize the last NumRuntimeRealComps() - n_external_attr_real runtime real attributes
149  for (int j = PIdx::nattribs + n_external_attr_real; j < ptile.NumRealComps() ; ++j)
150  {
151  auto attr_ptr = ptile.GetStructOfArrays().GetRealData(j).data();
152 #ifdef WARPX_QED
153  // Current runtime comp is quantum synchrotron optical depth
154  if (particle_comps.find("opticalDepthQSR") != particle_comps.end() &&
155  particle_comps.at("opticalDepthQSR") == j)
156  {
157  if (!do_qed_comps) { continue; }
158  const QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt =
159  p_qs_engine->build_optical_depth_functor();
160  // If the particle tile was allocated in a memory pool that can run on GPU, launch GPU kernel
161  if constexpr (amrex::RunOnGpu<typename PTile::template AllocatorType<amrex::Real>>::value) {
162  amrex::ParallelForRNG(stop - start,
163  [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept {
164  const int ip = i + start;
165  attr_ptr[ip] = quantum_sync_get_opt(engine);
166  });
167  // Otherwise (e.g. particle tile allocated in pinned memory), run on CPU
168  } else {
169  for (int ip = start; ip < stop; ++ip) {
170  attr_ptr[ip] = quantum_sync_get_opt(amrex::RandomEngine{});
171  }
172  }
173  }
174 
175  // Current runtime comp is Breit-Wheeler optical depth
176  if (particle_comps.find("opticalDepthBW") != particle_comps.end() &&
177  particle_comps.at("opticalDepthBW") == j)
178  {
179  if (!do_qed_comps) { continue; }
180  const BreitWheelerGetOpticalDepth breit_wheeler_get_opt =
181  p_bw_engine->build_optical_depth_functor();;
182  // If the particle tile was allocated in a memory pool that can run on GPU, launch GPU kernel
183  if constexpr (amrex::RunOnGpu<typename PTile::template AllocatorType<amrex::Real>>::value) {
184  amrex::ParallelForRNG(stop - start,
185  [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept {
186  const int ip = i + start;
187  attr_ptr[ip] = breit_wheeler_get_opt(engine);
188  });
189  // Otherwise (e.g. particle tile allocated in pinned memory), run on CPU
190  } else {
191  for (int ip = start; ip < stop; ++ip) {
192  attr_ptr[ip] = breit_wheeler_get_opt(amrex::RandomEngine{});
193  }
194  }
195  }
196 #endif
197 
198  for (int ia = 0; ia < n_user_real_attribs; ++ia)
199  {
200  // Current runtime comp is ia-th user defined attribute
201  if (particle_comps.find(user_real_attribs[ia]) != particle_comps.end() &&
202  particle_comps.at(user_real_attribs[ia]) == j)
203  {
204  const amrex::ParserExecutor<7> user_real_attrib_parserexec =
205  user_real_attrib_parser[ia]->compile<7>();
206  // If the particle tile was allocated in a memory pool that can run on GPU, launch GPU kernel
207  if constexpr (amrex::RunOnGpu<typename PTile::template AllocatorType<amrex::Real>>::value) {
208  amrex::ParallelFor(stop - start,
209  [=] AMREX_GPU_DEVICE (int i) noexcept {
210  const int ip = i + start;
211  amrex::ParticleReal xp, yp, zp;
212  get_position(ip, xp, yp, zp);
213  attr_ptr[ip] = user_real_attrib_parserexec(xp, yp, zp,
214  ux[ip], uy[ip], uz[ip], t);
215  });
216  // Otherwise (e.g. particle tile allocated in pinned memory), run on CPU
217  } else {
218  for (int ip = start; ip < stop; ++ip) {
219  amrex::ParticleReal xp, yp, zp;
220  get_position(ip, xp, yp, zp);
221  attr_ptr[ip] = user_real_attrib_parserexec(xp, yp, zp,
222  ux[ip], uy[ip], uz[ip], t);
223  }
224  }
225  }
226  }
227  }
228 
229  // Initialize the last NumRuntimeIntComps() - n_external_attr_int runtime int attributes
230  for (int j = n_external_attr_int; j < ptile.NumIntComps() ; ++j)
231  {
232  auto attr_ptr = ptile.GetStructOfArrays().GetIntData(j).data();
233 
234  // Current runtime comp is ionization level
235  if (particle_icomps.find("ionizationLevel") != particle_icomps.end() &&
236  particle_icomps.at("ionizationLevel") == j)
237  {
238  if constexpr (amrex::RunOnGpu<typename PTile::template AllocatorType<int>>::value) {
239  amrex::ParallelFor(stop - start,
240  [=] AMREX_GPU_DEVICE (int i) noexcept {
241  const int ip = i + start;
242  attr_ptr[ip] = ionization_initial_level;
243  });
244  } else {
245  for (int ip = start; ip < stop; ++ip) {
246  attr_ptr[ip] = ionization_initial_level;
247  }
248  }
249  }
250 
251  for (int ia = 0; ia < n_user_int_attribs; ++ia)
252  {
253  // Current runtime comp is ia-th user defined attribute
254  if (particle_icomps.find(user_int_attribs[ia]) != particle_icomps.end() &&
255  particle_icomps.at(user_int_attribs[ia]) == j)
256  {
257  const amrex::ParserExecutor<7> user_int_attrib_parserexec =
258  user_int_attrib_parser[ia]->compile<7>();
259  if constexpr (amrex::RunOnGpu<typename PTile::template AllocatorType<int>>::value) {
260  amrex::ParallelFor(stop - start,
261  [=] AMREX_GPU_DEVICE (int i) noexcept {
262  const int ip = i + start;
263  amrex::ParticleReal xp, yp, zp;
264  get_position(ip, xp, yp, zp);
265  attr_ptr[ip] = static_cast<int>(
266  user_int_attrib_parserexec(xp, yp, zp, ux[ip], uy[ip], uz[ip], t));
267  });
268  } else {
269  for (int ip = start; ip < stop; ++ip) {
270  amrex::ParticleReal xp, yp, zp;
271  get_position(ip, xp, yp, zp);
272  attr_ptr[ip] = static_cast<int>(
273  user_int_attrib_parserexec(xp, yp, zp, ux[ip], uy[ip], uz[ip], t));
274  }
275  }
276  }
277  }
278  }
279 }
280 
281 }
282 
283 #endif //WARPX_DEFAULTINITIALIZATION_H_
#define AMREX_FORCE_INLINE
#define AMREX_RESTRICT
#define AMREX_GPU_DEVICE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int initializeIntValue(const InitializationPolicy policy) noexcept
Definition: DefaultInitialization.H:78
InitializationPolicy
This set of initialization policies describes what happens when we need to create a new particle due ...
Definition: DefaultInitialization.H:39
static std::map< std::string, InitializationPolicy > initialization_policies
This map sets the initialization policy for each particle component used in WarpX.
Definition: DefaultInitialization.H:45
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal initializeRealValue(const InitializationPolicy policy, amrex::RandomEngine const &engine) noexcept
Definition: DefaultInitialization.H:62
Definition: BreitWheelerEngineWrapper.H:294
BreitWheelerGetOpticalDepth build_optical_depth_functor() const
Definition: BreitWheelerEngineWrapper.cpp:39
Definition: BreitWheelerEngineWrapper.H:78
Definition: QuantumSyncEngineWrapper.H:273
QuantumSynchrotronGetOpticalDepth build_optical_depth_functor()
Definition: QuantumSyncEngineWrapper.cpp:39
Definition: QuantumSyncEngineWrapper.H:76
amrex::Vector< amrex::Real > gett_new() const
Definition: WarpX.H:848
static WarpX & GetInstance()
Definition: WarpX.cpp:239
Definition: DefaultInitialization.H:90
void DefaultInitializeRuntimeAttributes(PTile &ptile, const int n_external_attr_real, const int n_external_attr_int, const std::vector< std::string > &user_real_attribs, const std::vector< std::string > &user_int_attribs, const std::map< std::string, int > &particle_comps, const std::map< std::string, int > &particle_icomps, const std::vector< amrex::Parser * > &user_real_attrib_parser, const std::vector< amrex::Parser * > &user_int_attrib_parser, const bool do_qed_comps, BreitWheelerEngine *p_bw_engine, QuantumSynchrotronEngine *p_qs_engine, const int ionization_initial_level, int start, int stop)
Default initialize runtime attributes in a tile. This routine does not initialize the first n_externa...
Definition: DefaultInitialization.H:118
Real Random()
void Abort(const std::string &msg)
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
AMREX_ATTRIBUTE_FLATTEN_FOR void ParallelForRNG(T n, L const &f) noexcept
i
Definition: check_interp_points_and_weights.py:174
default
Definition: run_alltests.py:113
value
Definition: updateAMReX.py:141
@ nattribs
number of compile-time attributes
Definition: NamedComponentParticleContainer.H:38
@ uz
Definition: NamedComponentParticleContainer.H:34
@ uy
Definition: NamedComponentParticleContainer.H:34
@ ux
Definition: NamedComponentParticleContainer.H:34