WarpX
InjectorMomentum.H
Go to the documentation of this file.
1 /* Copyright 2019-2020 Andrew Myers, Axel Huebl, Cameron Yang,
2  * Maxence Thevenet, Weiqun Zhang
3  *
4  * This file is part of WarpX.
5  *
6  * License: BSD-3-Clause-LBNL
7  */
8 #ifndef WARPX_INJECTOR_MOMENTUM_H_
9 #define WARPX_INJECTOR_MOMENTUM_H_
10 
11 #include "GetTemperature.H"
12 #include "GetVelocity.H"
13 #include "TemperatureProperties.H"
14 #include "VelocityProperties.H"
16 #include "Utils/TextMsg.H"
17 #include "Utils/WarpXConst.H"
18 
19 #include <AMReX.H>
20 #include <AMReX_Config.H>
21 #include <AMReX_Dim3.H>
22 #include <AMReX_GpuQualifiers.H>
23 #include <AMReX_REAL.H>
24 #include <AMReX_Parser.H>
25 #include <AMReX_Random.H>
26 
27 #include <AMReX_BaseFwd.H>
28 
29 #include <cmath>
30 #include <string>
31 
32 // struct whose getMomentum returns constant momentum.
34 {
35  InjectorMomentumConstant (amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) noexcept
36  : m_ux(a_ux), m_uy(a_uy), m_uz(a_uz) {}
37 
38  [[nodiscard]]
41  getMomentum (amrex::Real, amrex::Real, amrex::Real,
42  amrex::RandomEngine const&) const noexcept
43  {
44  return amrex::XDim3{m_ux,m_uy,m_uz};
45  }
46 
47  [[nodiscard]]
50  getBulkMomentum (amrex::Real, amrex::Real, amrex::Real) const noexcept
51  {
52  return amrex::XDim3{m_ux,m_uy,m_uz};
53  }
54 
55 private:
56  amrex::Real m_ux, m_uy, m_uz;
57 };
58 
59 // struct whose getMomentum returns momentum for 1 particle, from random
60 // gaussian distribution.
62 {
63  InjectorMomentumGaussian (amrex::Real a_ux_m, amrex::Real a_uy_m,
64  amrex::Real a_uz_m, amrex::Real a_ux_th,
65  amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept
66  : m_ux_m(a_ux_m), m_uy_m(a_uy_m), m_uz_m(a_uz_m),
67  m_ux_th(a_ux_th), m_uy_th(a_uy_th), m_uz_th(a_uz_th)
68  {}
69 
70  [[nodiscard]]
73  getMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/,
74  amrex::RandomEngine const& engine) const noexcept
75  {
79  }
80 
81  [[nodiscard]]
84  getBulkMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/) const noexcept
85  {
86  return amrex::XDim3{m_ux_m, m_uy_m, m_uz_m};
87  }
88 
89 private:
90  amrex::Real m_ux_m, m_uy_m, m_uz_m;
91  amrex::Real m_ux_th, m_uy_th, m_uz_th;
92 };
93 
94 // struct whose getMomentum returns momentum for 1 particle, from random
95 // gaussian flux distribution in the specified direction.
96 // Along the normal axis, the distribution is v*Gaussian,
97 // with the sign set by flux_direction.
99 {
100  InjectorMomentumGaussianFlux (amrex::Real a_ux_m, amrex::Real a_uy_m,
101  amrex::Real a_uz_m, amrex::Real a_ux_th,
102  amrex::Real a_uy_th, amrex::Real a_uz_th,
103  int a_flux_normal_axis, int a_flux_direction) noexcept
104  : m_ux_m(a_ux_m), m_uy_m(a_uy_m), m_uz_m(a_uz_m),
105  m_ux_th(a_ux_th), m_uy_th(a_uy_th), m_uz_th(a_uz_th),
106  m_flux_normal_axis(a_flux_normal_axis),
107  m_flux_direction(a_flux_direction)
108  {
109  }
110 
111  [[nodiscard]]
114  getMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/,
115  amrex::RandomEngine const& engine) const noexcept
116  {
117  using namespace amrex::literals;
118 
119  // Generate the distribution in the direction of the flux
120  amrex::Real u_m = 0, u_th = 0;
121  if (m_flux_normal_axis == 0) {
122  u_m = m_ux_m;
123  u_th = m_ux_th;
124  } else if (m_flux_normal_axis == 1) {
125  u_m = m_uy_m;
126  u_th = m_uy_th;
127  } else if (m_flux_normal_axis == 2) {
128  u_m = m_uz_m;
129  u_th = m_uz_th;
130  }
131  amrex::Real u = generateGaussianFluxDist(u_m, u_th, engine);
132  if (m_flux_direction < 0) { u = -u; }
133 
134  // Note: Here, in RZ geometry, the variables `ux` and `uy` actually
135  // correspond to the radial and azimuthal component of the momentum
136  // (and e.g.`m_flux_normal_axis==1` corresponds to v*Gaussian along theta)
137  amrex::Real const ux = (m_flux_normal_axis == 0 ? u : amrex::RandomNormal(m_ux_m, m_ux_th, engine));
138  amrex::Real const uy = (m_flux_normal_axis == 1 ? u : amrex::RandomNormal(m_uy_m, m_uy_th, engine));
139  amrex::Real const uz = (m_flux_normal_axis == 2 ? u : amrex::RandomNormal(m_uz_m, m_uz_th, engine));
140  return amrex::XDim3{ux, uy, uz};
141  }
142 
143  [[nodiscard]]
146  getBulkMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/) const noexcept
147  {
148  return amrex::XDim3{m_ux_m, m_uy_m, m_uz_m};
149  }
150 
151 private:
152  amrex::Real m_ux_m, m_uy_m, m_uz_m;
153  amrex::Real m_ux_th, m_uy_th, m_uz_th;
156 };
157 
158 
159 // struct whose getMomentum returns momentum for 1 particle, from random
160 // uniform distribution u_min < u < u_max
162 {
163  InjectorMomentumUniform (amrex::Real a_ux_min, amrex::Real a_uy_min,
164  amrex::Real a_uz_min, amrex::Real a_ux_max,
165  amrex::Real a_uy_max, amrex::Real a_uz_max) noexcept
166  : m_ux_min(a_ux_min), m_uy_min(a_uy_min), m_uz_min(a_uz_min),
167  m_ux_max(a_ux_max), m_uy_max(a_uy_max), m_uz_max(a_uz_max),
168  m_Dux(m_ux_max - m_ux_min),
169  m_Duy(m_uy_max - m_uy_min),
170  m_Duz(m_uz_max - m_uz_min),
171  m_ux_h(amrex::Real(0.5) * (m_ux_max + m_ux_min)),
172  m_uy_h(amrex::Real(0.5) * (m_uy_max + m_uy_min)),
173  m_uz_h(amrex::Real(0.5) * (m_uz_max + m_uz_min))
174  {}
175 
176  [[nodiscard]]
179  getMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/,
180  amrex::RandomEngine const& engine) const noexcept
181  {
182  return amrex::XDim3{m_ux_min + amrex::Random(engine) * m_Dux,
183  m_uy_min + amrex::Random(engine) * m_Duy,
184  m_uz_min + amrex::Random(engine) * m_Duz};
185  }
186 
187  [[nodiscard]]
190  getBulkMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/) const noexcept
191  {
192  return amrex::XDim3{m_ux_h, m_uy_h, m_uz_h};
193  }
194 
195 private:
196  amrex::Real m_ux_min, m_uy_min, m_uz_min;
197  amrex::Real m_ux_max, m_uy_max, m_uz_max;
198  amrex::Real m_Dux, m_Duy, m_Duz;
199  amrex::Real m_ux_h, m_uy_h, m_uz_h;
200 };
201 
202 // struct whose getMomentum returns momentum for 1 particle with relativistic
203 // drift velocity beta, from the Maxwell-Boltzmann distribution.
205 {
206  // Constructor whose inputs are:
207  // a reference to the initial temperature container t,
208  // a reference to the initial velocity container b
210  : velocity(b), temperature(t)
211  {}
212 
213  [[nodiscard]]
216  getMomentum (amrex::Real const x, amrex::Real const y, amrex::Real const z,
217  amrex::RandomEngine const& engine) const noexcept
218  {
219  using namespace amrex::literals;
220 
221  // Calculate the local temperature and check if it's too
222  // high for Boltzmann or less than zero
223  amrex::Real const theta = temperature(x,y,z);
224  if (theta < 0._rt) {
225  amrex::Abort("Negative temperature parameter theta encountered, which is not allowed");
226  }
227  // Calculate local velocity and abort if |beta|>=1
228  amrex::Real const beta = velocity(x,y,z);
229  if (beta <= -1._rt || beta >= 1._rt) {
230  amrex::Abort("beta = v/c magnitude greater than or equal to 1");
231  }
232  // Calculate the value of vave from the local temperature
233  amrex::Real const vave = std::sqrt(theta);
234  int const dir = velocity.direction();
235 
236  amrex::Real u[3];
237  u[dir] = amrex::RandomNormal(0._rt, vave, engine);
238  u[(dir+1)%3] = amrex::RandomNormal(0._rt, vave, engine);
239  u[(dir+2)%3] = amrex::RandomNormal(0._rt, vave, engine);
240  amrex::Real const gamma = std::sqrt(1._rt + u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
241 
242  // The following condition is equation 32 in Zenitani 2015
243  // (Phys. Plasmas 22, 042116) , called the flipping method. It
244  // transforms the integral: d3x' -> d3x where d3x' is the volume
245  // element for positions in the boosted frame. The particle positions
246  // and densities can be initialized in the simulation frame.
247  // The flipping method can transform any symmetric distribution from one
248  // reference frame to another moving at a relative velocity of beta.
249  // An equivalent alternative to this method native to WarpX would be to
250  // initialize the particle positions and densities in the frame moving
251  // at speed beta, and then perform a Lorentz transform on the positions
252  // and MB sampled velocities to the simulation frame.
253  if(-beta*u[dir]/gamma > amrex::Random(engine))
254  {
255  u[dir] = -u[dir];
256  }
257  // This Lorentz transform is equation 17 in Zenitani.
258  // It transforms the integral d3u' -> d3u
259  // where d3u' is the volume element for momentum in the boosted frame.
260  u[dir] = 1._rt/std::sqrt(1._rt-beta*beta)*(u[dir]+gamma*beta);
261  // Note that if beta = 0 then the flipping method and Lorentz transform
262  // have no effect on the u[dir] direction.
263  return amrex::XDim3 {u[0],u[1],u[2]};
264  }
265 
266  [[nodiscard]]
269  getBulkMomentum (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
270  {
271  using namespace amrex::literals;
272  amrex::Real u[3];
273  for (auto& el : u) { el = 0.0_rt; }
274  const amrex::Real beta = velocity(x,y,z);
275  int const dir = velocity.direction();
276  const auto gamma = 1._rt/std::sqrt(1._rt-beta*beta);
277  u[dir] = gamma*beta;
278  return amrex::XDim3 {u[0],u[1],u[2]};
279  }
280 
281 private:
284 };
285 
286 // struct whose getMomentum returns momentum for 1 particle with relativistic
287 // drift velocity beta, from the Maxwell-Juttner distribution. Method is from
288 // Zenitani 2015 (Phys. Plasmas 22, 042116).
290 {
291  // Constructor whose inputs are:
292  // a reference to the initial temperature container t,
293  // a reference to the initial velocity container b
295  : velocity(b), temperature(t)
296  {}
297 
298  [[nodiscard]]
301  getMomentum (amrex::Real const x, amrex::Real const y, amrex::Real const z,
302  amrex::RandomEngine const& engine) const noexcept
303  {
304  using namespace amrex::literals;
305  // Sobol method for sampling MJ Speeds,
306  // from Zenitani 2015 (Phys. Plasmas 22, 042116).
307  amrex::Real x1, x2, gamma;
308  amrex::Real u [3];
309  amrex::Real const theta = temperature(x,y,z);
310  // Check if temperature is too low to do sampling method. Abort for now,
311  // in future should implement an alternate method e.g. inverse transform
312  if (theta < 0.1_rt) {
313  amrex::Abort("Temeprature parameter theta is less than minimum 0.1 allowed for Maxwell-Juttner");
314  }
315  // Calculate local velocity and abort if |beta|>=1
316  amrex::Real const beta = velocity(x,y,z);
317  if (beta <= -1._rt || beta >= 1._rt) {
318  amrex::Abort("beta = v/c magnitude greater than or equal to 1");
319  }
320  int const dir = velocity.direction();
321  x1 = static_cast<amrex::Real>(0._rt);
322  gamma = static_cast<amrex::Real>(0._rt);
323  u[dir] = static_cast<amrex::Real>(0._rt);
324  // This condition is equation 10 in Zenitani,
325  // though x1 is defined differently.
326  while(u[dir]-gamma <= x1)
327  {
328  u[dir] = -theta*
329  std::log(amrex::Random(engine)*amrex::Random(engine)*amrex::Random(engine));
330  gamma = std::sqrt(1._rt+u[dir]*u[dir]);
331  x1 = theta*std::log(amrex::Random(engine));
332  }
333  // The following code samples a random unit vector
334  // and multiplies the result by speed u[dir].
335  x1 = amrex::Random(engine);
336  x2 = amrex::Random(engine);
337  // Direction dir is an input parameter that sets the boost direction:
338  // 'x' -> d = 0, 'y' -> d = 1, 'z' -> d = 2.
339  u[(dir+1)%3] = 2._rt*u[dir]*std::sqrt(x1*(1._rt-x1))*std::sin(2._rt*MathConst::pi*x2);
340  u[(dir+2)%3] = 2._rt*u[dir]*std::sqrt(x1*(1._rt-x1))*std::cos(2._rt*MathConst::pi*x2);
341  // The value of dir is the boost direction to be transformed.
342  u[dir] = u[dir]*(2._rt*x1-1._rt);
343  x1 = amrex::Random(engine);
344  // The following condition is equation 32 in Zenitani, called
345  // The flipping method. It transforms the integral: d3x' -> d3x
346  // where d3x' is the volume element for positions in the boosted frame.
347  // The particle positions and densities can be initialized in the
348  // simulation frame with this method.
349  // The flipping method can similarly transform any
350  // symmetric distribution from one reference frame to another moving at
351  // a relative velocity of beta.
352  // An equivalent alternative to this method native to WarpX
353  // would be to initialize the particle positions and densities in the
354  // frame moving at speed beta, and then perform a Lorentz transform
355  // on their positions and MJ sampled velocities to the simulation frame.
356  if(-beta*u[dir]/gamma>x1)
357  {
358  u[dir] = -u[dir];
359  }
360  // This Lorentz transform is equation 17 in Zenitani.
361  // It transforms the integral d3u' -> d3u
362  // where d3u' is the volume element for momentum in the boosted frame.
363  u[dir] = 1._rt/std::sqrt(1._rt-beta*beta)*(u[dir]+gamma*beta);
364  // Note that if beta = 0 then the flipping method and Lorentz transform
365  // have no effect on the u[dir] direction.
366  return amrex::XDim3 {u[0],u[1],u[2]};
367  }
368 
369  [[nodiscard]]
372  getBulkMomentum (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
373  {
374  using namespace amrex::literals;
375  amrex::Real u[3];
376  for (auto& el : u) { el = 0.0_rt; }
377  amrex::Real const beta = velocity(x,y,z);
378  int const dir = velocity.direction();
379  auto const gamma = 1._rt/std::sqrt(1._rt-beta*beta);
380  u[dir] = gamma*beta;
381  return amrex::XDim3 {u[0],u[1],u[2]};
382  }
383 
384 private:
387 };
388 
397 {
398  InjectorMomentumRadialExpansion (amrex::Real a_u_over_r) noexcept
399  : u_over_r(a_u_over_r)
400  {}
401 
402  [[nodiscard]]
405  getMomentum (amrex::Real x, amrex::Real y, amrex::Real z,
406  amrex::RandomEngine const&) const noexcept
407  {
408  return {x*u_over_r, y*u_over_r, z*u_over_r};
409  }
410 
411  [[nodiscard]]
414  getBulkMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
415  {
416  return {x*u_over_r, y*u_over_r, z*u_over_r};
417  }
418 
419 private:
420  amrex::Real u_over_r;
421 };
422 
423 // struct whose getMomentum returns local momentum computed from parser.
425 {
427  amrex::ParserExecutor<3> const& a_uy_parser,
428  amrex::ParserExecutor<3> const& a_uz_parser) noexcept
429  : m_ux_parser(a_ux_parser), m_uy_parser(a_uy_parser),
430  m_uz_parser(a_uz_parser) {}
431 
432  [[nodiscard]]
435  getMomentum (amrex::Real x, amrex::Real y, amrex::Real z,
436  amrex::RandomEngine const&) const noexcept
437  {
438  return amrex::XDim3{m_ux_parser(x,y,z),m_uy_parser(x,y,z),m_uz_parser(x,y,z)};
439  }
440 
441  [[nodiscard]]
444  getBulkMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
445  {
446  return amrex::XDim3{m_ux_parser(x,y,z),m_uy_parser(x,y,z),m_uz_parser(x,y,z)};
447  }
448 
449  amrex::ParserExecutor<3> m_ux_parser, m_uy_parser, m_uz_parser;
450 };
451 
452 // struct whose getMomentum returns local momentum and thermal spread computed from parser.
454 {
456  amrex::ParserExecutor<3> const& a_uy_m_parser,
457  amrex::ParserExecutor<3> const& a_uz_m_parser,
458  amrex::ParserExecutor<3> const& a_ux_th_parser,
459  amrex::ParserExecutor<3> const& a_uy_th_parser,
460  amrex::ParserExecutor<3> const& a_uz_th_parser) noexcept
461  : m_ux_m_parser(a_ux_m_parser), m_uy_m_parser(a_uy_m_parser), m_uz_m_parser(a_uz_m_parser),
462  m_ux_th_parser(a_ux_th_parser), m_uy_th_parser(a_uy_th_parser), m_uz_th_parser(a_uz_th_parser) {}
463 
464  [[nodiscard]]
467  getMomentum (amrex::Real x, amrex::Real y, amrex::Real z,
468  amrex::RandomEngine const& engine) const noexcept
469  {
470  amrex::Real const ux_m = m_ux_m_parser(x,y,z);
471  amrex::Real const uy_m = m_uy_m_parser(x,y,z);
472  amrex::Real const uz_m = m_uz_m_parser(x,y,z);
473  amrex::Real const ux_th = m_ux_th_parser(x,y,z);
474  amrex::Real const uy_th = m_uy_th_parser(x,y,z);
475  amrex::Real const uz_th = m_uz_th_parser(x,y,z);
476  return amrex::XDim3{amrex::RandomNormal(ux_m, ux_th, engine),
477  amrex::RandomNormal(uy_m, uy_th, engine),
478  amrex::RandomNormal(uz_m, uz_th, engine)};
479  }
480 
481  [[nodiscard]]
484  getBulkMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
485  {
486  return amrex::XDim3{m_ux_m_parser(x,y,z), m_uy_m_parser(x,y,z), m_uz_m_parser(x,y,z)};
487  }
488 
489  amrex::ParserExecutor<3> m_ux_m_parser, m_uy_m_parser, m_uz_m_parser;
490  amrex::ParserExecutor<3> m_ux_th_parser, m_uy_th_parser, m_uz_th_parser;
491 };
492 
493 // Base struct for momentum injector.
494 // InjectorMomentum contains a union (called Object) that holds any one
495 // instance of:
496 // - InjectorMomentumConstant : to generate constant density;
497 // - InjectorMomentumGaussian : to generate gaussian distribution;
498 // - InjectorMomentumGaussianFlux : to generate v*gaussian distribution;
499 // - InjectorMomentumRadialExpansion: to generate radial expansion;
500 // - InjectorMomentumParser : to generate momentum from parser;
501 // - InjectorMomentumGaussianParser : to generate momentum and thermal spread from parser;
502 // The choice is made at runtime, depending in the constructor called.
503 // This mimics virtual functions.
505 {
506  // This constructor stores a InjectorMomentumConstant in union object.
508  amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz)
509  : type(Type::constant),
510  object(t, a_ux, a_uy, a_uz)
511  { }
512 
513  // This constructor stores a InjectorMomentumParser in union object.
515  amrex::ParserExecutor<3> const& a_ux_parser,
516  amrex::ParserExecutor<3> const& a_uy_parser,
517  amrex::ParserExecutor<3> const& a_uz_parser)
518  : type(Type::parser),
519  object(t, a_ux_parser, a_uy_parser, a_uz_parser)
520  { }
521 
522  // This constructor stores a InjectorMomentumGaussianParser in union object.
524  amrex::ParserExecutor<3> const& a_ux_m_parser,
525  amrex::ParserExecutor<3> const& a_uy_m_parser,
526  amrex::ParserExecutor<3> const& a_uz_m_parser,
527  amrex::ParserExecutor<3> const& a_ux_th_parser,
528  amrex::ParserExecutor<3> const& a_uy_th_parser,
529  amrex::ParserExecutor<3> const& a_uz_th_parser)
530  : type(Type::gaussianparser),
531  object(t, a_ux_m_parser, a_uy_m_parser, a_uz_m_parser,
532  a_ux_th_parser, a_uy_th_parser, a_uz_th_parser)
533  { }
534 
535  // This constructor stores a InjectorMomentumGaussian in union object.
537  amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m,
538  amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th)
539  : type(Type::gaussian),
540  object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th)
541  { }
542 
543  // This constructor stores a InjectorMomentumGaussianFlux in union object.
545  amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m,
546  amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th,
547  int a_flux_normal_axis, int a_flux_direction)
548  : type(Type::gaussianflux),
549  object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th,a_flux_normal_axis,a_flux_direction)
550  { }
551 
552  // This constructor stores a InjectorMomentumUniform in union object.
554  amrex::Real a_ux_min, amrex::Real a_uy_min, amrex::Real a_uz_min,
555  amrex::Real a_ux_max, amrex::Real a_uy_max, amrex::Real a_uz_max)
556  : type(Type::uniform),
557  object(t,a_ux_min,a_uy_min,a_uz_min,a_ux_max,a_uy_max,a_uz_max)
558  { }
559 
561  GetTemperature const& temperature, GetVelocity const& velocity)
562  : type(Type::boltzmann),
563  object(t, temperature, velocity)
564  { }
565 
566  // This constructor stores a InjectorMomentumJuttner in union object.
568  GetTemperature const& temperature, GetVelocity const& velocity)
569  : type(Type::juttner),
570  object(t, temperature, velocity)
571  { }
572 
573  // This constructor stores a InjectorMomentumRadialExpansion in union object.
575  amrex::Real u_over_r)
576  : type(Type::radial_expansion),
577  object(t, u_over_r)
578  { }
579 
580  // Explicitly prevent the compiler from generating copy constructors
581  // and copy assignment operators.
584  void operator= (InjectorMomentum const&) = delete;
585  void operator= (InjectorMomentum &&) = delete;
586 
587  // Default destructor
588  ~InjectorMomentum() = default;
589 
590  void clear ();
591 
592  // call getMomentum from the object stored in the union
593  // (the union is called Object, and the instance is called object).
594  [[nodiscard]]
597  getMomentum (amrex::Real x, amrex::Real y, amrex::Real z,
598  amrex::RandomEngine const& engine) const noexcept
599  {
600  switch (type)
601  {
602  case Type::parser:
603  {
604  return object.parser.getMomentum(x,y,z,engine);
605  }
606  case Type::gaussian:
607  {
608  return object.gaussian.getMomentum(x,y,z,engine);
609  }
610  case Type::gaussianparser:
611  {
612  return object.gaussianparser.getMomentum(x,y,z,engine);
613  }
614  case Type::gaussianflux:
615  {
616  return object.gaussianflux.getMomentum(x,y,z,engine);
617  }
618  case Type::uniform:
619  {
620  return object.uniform.getMomentum(x,y,z,engine);
621  }
622  case Type::boltzmann:
623  {
624  return object.boltzmann.getMomentum(x,y,z,engine);
625  }
626  case Type::juttner:
627  {
628  return object.juttner.getMomentum(x,y,z,engine);
629  }
630  case Type::constant:
631  {
632  return object.constant.getMomentum(x,y,z,engine);
633  }
634  case Type::radial_expansion:
635  {
636  return object.radial_expansion.getMomentum(x,y,z,engine);
637  }
638  default:
639  {
640  amrex::Abort("InjectorMomentum: unknown type");
641  return {0.0,0.0,0.0};
642  }
643  }
644  }
645 
646  // call getBulkMomentum from the object stored in the union
647  // (the union is called Object, and the instance is called object).
648  [[nodiscard]]
651  getBulkMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
652  {
653  switch (type)
654  {
655  case Type::parser:
656  {
657  return object.parser.getBulkMomentum(x,y,z);
658  }
659  case Type::gaussian:
660  {
661  return object.gaussian.getBulkMomentum(x,y,z);
662  }
663  case Type::gaussianparser:
664  {
665  return object.gaussianparser.getBulkMomentum(x,y,z);
666  }
667  case Type::gaussianflux:
668  {
669  return object.gaussianflux.getBulkMomentum(x,y,z);
670  }
671  case Type::uniform:
672  {
673  return object.uniform.getBulkMomentum(x,y,z);
674  }
675  case Type::boltzmann:
676  {
677  return object.boltzmann.getBulkMomentum(x,y,z);
678  }
679  case Type::juttner:
680  {
681  return object.juttner.getBulkMomentum(x,y,z);
682  }
683  case Type::constant:
684  {
685  return object.constant.getBulkMomentum(x,y,z);
686  }
687  case Type::radial_expansion:
688  {
689  return object.radial_expansion.getBulkMomentum(x,y,z);
690  }
691  default:
692  {
693  amrex::Abort("InjectorMomentum: unknown type");
694  return {0.0,0.0,0.0};
695  }
696  }
697  }
698 
699  enum struct Type { constant, gaussian, gaussianflux, uniform, boltzmann, juttner, radial_expansion, parser, gaussianparser };
701 
702 private:
703 
704  // An instance of union Object constructs and stores any one of
705  // the objects declared (constant or gaussian or
706  // radial_expansion or parser).
707  union Object {
709  amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) noexcept
710  : constant(a_ux,a_uy,a_uz) {}
712  amrex::Real a_ux_m, amrex::Real a_uy_m,
713  amrex::Real a_uz_m, amrex::Real a_ux_th,
714  amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept
715  : gaussian(a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th) {}
717  amrex::Real a_ux_m, amrex::Real a_uy_m,
718  amrex::Real a_uz_m, amrex::Real a_ux_th,
719  amrex::Real a_uy_th, amrex::Real a_uz_th,
720  int a_flux_normal_axis, int a_flux_direction) noexcept
721  : gaussianflux(a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th,a_flux_normal_axis,a_flux_direction) {}
723  amrex::Real a_ux_min, amrex::Real a_uy_min,
724  amrex::Real a_uz_min, amrex::Real a_ux_max,
725  amrex::Real a_uy_max, amrex::Real a_uz_max) noexcept
726  : uniform(a_ux_min,a_uy_min,a_uz_min,a_ux_max,a_uy_max,a_uz_max) {}
728  GetTemperature const& t, GetVelocity const& b) noexcept
729  : boltzmann(t,b) {}
731  GetTemperature const& t, GetVelocity const& b) noexcept
732  : juttner(t,b) {}
734  amrex::Real u_over_r) noexcept
735  : radial_expansion(u_over_r) {}
737  amrex::ParserExecutor<3> const& a_ux_parser,
738  amrex::ParserExecutor<3> const& a_uy_parser,
739  amrex::ParserExecutor<3> const& a_uz_parser) noexcept
740  : parser(a_ux_parser, a_uy_parser, a_uz_parser) {}
742  amrex::ParserExecutor<3> const& a_ux_m_parser,
743  amrex::ParserExecutor<3> const& a_uy_m_parser,
744  amrex::ParserExecutor<3> const& a_uz_m_parser,
745  amrex::ParserExecutor<3> const& a_ux_th_parser,
746  amrex::ParserExecutor<3> const& a_uy_th_parser,
747  amrex::ParserExecutor<3> const& a_uz_th_parser) noexcept
748  : gaussianparser(a_ux_m_parser, a_uy_m_parser, a_uz_m_parser,
749  a_ux_th_parser, a_uy_th_parser, a_uz_th_parser) {}
759  };
761 };
762 
763 // In order for InjectorMomentum to be trivially copyable, its destructor
764 // must be trivial. So we have to rely on a custom deleter for unique_ptr.
766  void operator () (InjectorMomentum* p) const {
767  if (p) {
768  p->clear();
769  delete p;
770  }
771  }
772 };
773 
774 #endif //WARPX_INJECTOR_MOMENTUM_H_
#define AMREX_GPU_HOST_DEVICE
Real Random()
Real RandomNormal(Real mean, Real stddev)
void Abort(const std::string &msg)
list x1
Definition: plot_particle_path.py:130
type
Definition: run_alltests_1node.py:72
beta
Definition: stencil.py:434
int gamma
boosted frame
Definition: stencil.py:431
parser
Definition: stencil.py:411
Get temperature at a point on the grid.
Definition: GetTemperature.H:23
Definition: GetVelocity.H:21
Definition: InjectorMomentum.H:205
GetTemperature temperature
Definition: InjectorMomentum.H:283
GetVelocity velocity
Definition: InjectorMomentum.H:282
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real const x, amrex::Real const y, amrex::Real const z, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:216
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
Definition: InjectorMomentum.H:269
InjectorMomentumBoltzmann(GetTemperature const &t, GetVelocity const &b) noexcept
Definition: InjectorMomentum.H:209
Definition: InjectorMomentum.H:34
amrex::Real m_uy
Definition: InjectorMomentum.H:56
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real, amrex::Real, amrex::Real) const noexcept
Definition: InjectorMomentum.H:50
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real, amrex::Real, amrex::Real, amrex::RandomEngine const &) const noexcept
Definition: InjectorMomentum.H:41
amrex::Real m_ux
Definition: InjectorMomentum.H:56
InjectorMomentumConstant(amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) noexcept
Definition: InjectorMomentum.H:35
amrex::Real m_uz
Definition: InjectorMomentum.H:56
Definition: InjectorMomentum.H:765
Definition: InjectorMomentum.H:99
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real, amrex::Real, amrex::Real) const noexcept
Definition: InjectorMomentum.H:146
amrex::Real m_uy_m
Definition: InjectorMomentum.H:152
int m_flux_normal_axis
Definition: InjectorMomentum.H:154
amrex::Real m_ux_m
Definition: InjectorMomentum.H:152
amrex::Real m_uz_th
Definition: InjectorMomentum.H:153
amrex::Real m_ux_th
Definition: InjectorMomentum.H:153
amrex::Real m_uy_th
Definition: InjectorMomentum.H:153
InjectorMomentumGaussianFlux(amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th, int a_flux_normal_axis, int a_flux_direction) noexcept
Definition: InjectorMomentum.H:100
amrex::Real m_uz_m
Definition: InjectorMomentum.H:152
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real, amrex::Real, amrex::Real, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:114
int m_flux_direction
Definition: InjectorMomentum.H:155
Definition: InjectorMomentum.H:62
amrex::Real m_uz_m
Definition: InjectorMomentum.H:90
amrex::Real m_ux_m
Definition: InjectorMomentum.H:90
amrex::Real m_ux_th
Definition: InjectorMomentum.H:91
amrex::Real m_uy_th
Definition: InjectorMomentum.H:91
amrex::Real m_uy_m
Definition: InjectorMomentum.H:90
InjectorMomentumGaussian(amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept
Definition: InjectorMomentum.H:63
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real, amrex::Real, amrex::Real, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:73
amrex::Real m_uz_th
Definition: InjectorMomentum.H:91
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real, amrex::Real, amrex::Real) const noexcept
Definition: InjectorMomentum.H:84
Definition: InjectorMomentum.H:454
InjectorMomentumGaussianParser(amrex::ParserExecutor< 3 > const &a_ux_m_parser, amrex::ParserExecutor< 3 > const &a_uy_m_parser, amrex::ParserExecutor< 3 > const &a_uz_m_parser, amrex::ParserExecutor< 3 > const &a_ux_th_parser, amrex::ParserExecutor< 3 > const &a_uy_th_parser, amrex::ParserExecutor< 3 > const &a_uz_th_parser) noexcept
Definition: InjectorMomentum.H:455
amrex::ParserExecutor< 3 > m_ux_th_parser
Definition: InjectorMomentum.H:490
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
Definition: InjectorMomentum.H:484
amrex::ParserExecutor< 3 > m_ux_m_parser
Definition: InjectorMomentum.H:489
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real x, amrex::Real y, amrex::Real z, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:467
Definition: InjectorMomentum.H:505
InjectorMomentum(InjectorMomentumRadialExpansion *t, amrex::Real u_over_r)
Definition: InjectorMomentum.H:574
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real x, amrex::Real y, amrex::Real z, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:597
Type type
Definition: InjectorMomentum.H:700
InjectorMomentum(InjectorMomentum &&)=delete
~InjectorMomentum()=default
InjectorMomentum(InjectorMomentumParser *t, amrex::ParserExecutor< 3 > const &a_ux_parser, amrex::ParserExecutor< 3 > const &a_uy_parser, amrex::ParserExecutor< 3 > const &a_uz_parser)
Definition: InjectorMomentum.H:514
void clear()
Definition: InjectorMomentum.cpp:12
InjectorMomentum(InjectorMomentumConstant *t, amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz)
Definition: InjectorMomentum.H:507
InjectorMomentum(InjectorMomentum const &)=delete
InjectorMomentum(InjectorMomentumGaussian *t, amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th)
Definition: InjectorMomentum.H:536
InjectorMomentum(InjectorMomentumUniform *t, amrex::Real a_ux_min, amrex::Real a_uy_min, amrex::Real a_uz_min, amrex::Real a_ux_max, amrex::Real a_uy_max, amrex::Real a_uz_max)
Definition: InjectorMomentum.H:553
Type
Definition: InjectorMomentum.H:699
Object object
Definition: InjectorMomentum.H:760
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
Definition: InjectorMomentum.H:651
InjectorMomentum(InjectorMomentumGaussianFlux *t, amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th, int a_flux_normal_axis, int a_flux_direction)
Definition: InjectorMomentum.H:544
InjectorMomentum(InjectorMomentumBoltzmann *t, GetTemperature const &temperature, GetVelocity const &velocity)
Definition: InjectorMomentum.H:560
InjectorMomentum(InjectorMomentumJuttner *t, GetTemperature const &temperature, GetVelocity const &velocity)
Definition: InjectorMomentum.H:567
InjectorMomentum(InjectorMomentumGaussianParser *t, amrex::ParserExecutor< 3 > const &a_ux_m_parser, amrex::ParserExecutor< 3 > const &a_uy_m_parser, amrex::ParserExecutor< 3 > const &a_uz_m_parser, amrex::ParserExecutor< 3 > const &a_ux_th_parser, amrex::ParserExecutor< 3 > const &a_uy_th_parser, amrex::ParserExecutor< 3 > const &a_uz_th_parser)
Definition: InjectorMomentum.H:523
Definition: InjectorMomentum.H:290
GetVelocity velocity
Definition: InjectorMomentum.H:385
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
Definition: InjectorMomentum.H:372
GetTemperature temperature
Definition: InjectorMomentum.H:386
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real const x, amrex::Real const y, amrex::Real const z, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:301
InjectorMomentumJuttner(GetTemperature const &t, GetVelocity const &b) noexcept
Definition: InjectorMomentum.H:294
Definition: InjectorMomentum.H:425
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real x, amrex::Real y, amrex::Real z, amrex::RandomEngine const &) const noexcept
Definition: InjectorMomentum.H:435
InjectorMomentumParser(amrex::ParserExecutor< 3 > const &a_ux_parser, amrex::ParserExecutor< 3 > const &a_uy_parser, amrex::ParserExecutor< 3 > const &a_uz_parser) noexcept
Definition: InjectorMomentum.H:426
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
Definition: InjectorMomentum.H:444
amrex::ParserExecutor< 3 > m_ux_parser
Definition: InjectorMomentum.H:449
struct whose getMomentum returns momentum for 1 particle, for radial expansion.
Definition: InjectorMomentum.H:397
InjectorMomentumRadialExpansion(amrex::Real a_u_over_r) noexcept
Definition: InjectorMomentum.H:398
amrex::Real u_over_r
Definition: InjectorMomentum.H:420
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real x, amrex::Real y, amrex::Real z, amrex::RandomEngine const &) const noexcept
Definition: InjectorMomentum.H:405
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
Definition: InjectorMomentum.H:414
Definition: InjectorMomentum.H:162
amrex::Real m_Dux
Definition: InjectorMomentum.H:198
amrex::Real m_ux_h
Definition: InjectorMomentum.H:199
amrex::Real m_ux_max
Definition: InjectorMomentum.H:197
AMREX_GPU_HOST_DEVICE amrex::XDim3 getBulkMomentum(amrex::Real, amrex::Real, amrex::Real) const noexcept
Definition: InjectorMomentum.H:190
AMREX_GPU_HOST_DEVICE amrex::XDim3 getMomentum(amrex::Real, amrex::Real, amrex::Real, amrex::RandomEngine const &engine) const noexcept
Definition: InjectorMomentum.H:179
InjectorMomentumUniform(amrex::Real a_ux_min, amrex::Real a_uy_min, amrex::Real a_uz_min, amrex::Real a_ux_max, amrex::Real a_uy_max, amrex::Real a_uz_max) noexcept
Definition: InjectorMomentum.H:163
amrex::Real m_ux_min
Definition: InjectorMomentum.H:196
Definition: InjectorMomentum.H:707
InjectorMomentumGaussian gaussian
Definition: InjectorMomentum.H:751
InjectorMomentumUniform uniform
Definition: InjectorMomentum.H:753
Object(InjectorMomentumGaussianFlux *, amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th, int a_flux_normal_axis, int a_flux_direction) noexcept
Definition: InjectorMomentum.H:716
Object(InjectorMomentumBoltzmann *, GetTemperature const &t, GetVelocity const &b) noexcept
Definition: InjectorMomentum.H:727
Object(InjectorMomentumParser *, amrex::ParserExecutor< 3 > const &a_ux_parser, amrex::ParserExecutor< 3 > const &a_uy_parser, amrex::ParserExecutor< 3 > const &a_uz_parser) noexcept
Definition: InjectorMomentum.H:736
InjectorMomentumGaussianParser gaussianparser
Definition: InjectorMomentum.H:758
InjectorMomentumJuttner juttner
Definition: InjectorMomentum.H:755
Object(InjectorMomentumUniform *, amrex::Real a_ux_min, amrex::Real a_uy_min, amrex::Real a_uz_min, amrex::Real a_ux_max, amrex::Real a_uy_max, amrex::Real a_uz_max) noexcept
Definition: InjectorMomentum.H:722
InjectorMomentumRadialExpansion radial_expansion
Definition: InjectorMomentum.H:756
InjectorMomentumBoltzmann boltzmann
Definition: InjectorMomentum.H:754
InjectorMomentumGaussianFlux gaussianflux
Definition: InjectorMomentum.H:752
Object(InjectorMomentumJuttner *, GetTemperature const &t, GetVelocity const &b) noexcept
Definition: InjectorMomentum.H:730
Object(InjectorMomentumRadialExpansion *, amrex::Real u_over_r) noexcept
Definition: InjectorMomentum.H:733
Object(InjectorMomentumGaussian *, amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept
Definition: InjectorMomentum.H:711
InjectorMomentumParser parser
Definition: InjectorMomentum.H:757
InjectorMomentumConstant constant
Definition: InjectorMomentum.H:750
Object(InjectorMomentumGaussianParser *, amrex::ParserExecutor< 3 > const &a_ux_m_parser, amrex::ParserExecutor< 3 > const &a_uy_m_parser, amrex::ParserExecutor< 3 > const &a_uz_m_parser, amrex::ParserExecutor< 3 > const &a_ux_th_parser, amrex::ParserExecutor< 3 > const &a_uy_th_parser, amrex::ParserExecutor< 3 > const &a_uz_th_parser) noexcept
Definition: InjectorMomentum.H:741
Object(InjectorMomentumConstant *, amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) noexcept
Definition: InjectorMomentum.H:708