WarpX
SampleGaussianFluxDistribution.H
Go to the documentation of this file.
1 /* Copyright 2024 Remi Lehe, Revathi Jambunathan
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_SAMPLE_GAUSSIAN_FLUX_DISTRIBUTION_H
9 #define WARPX_SAMPLE_GAUSSIAN_FLUX_DISTRIBUTION_H
10 
11 #include <AMReX_Random.H>
12 
13 namespace {
21  [[nodiscard]]
24  amrex::Real
25  generateGaussianFluxDist( amrex::Real u_m, amrex::Real u_th, amrex::RandomEngine const& engine ) {
26 
27  using namespace amrex::literals;
28 
29  // Momentum to be returned at the end of this function
30  amrex::Real u = 0._rt;
31 
32  const amrex::Real abs_u_m = std::abs(u_m);
33 
34  if (u_th == 0._rt) {
35  u = u_m; // Trivial case ; avoids division by 0 in the rest of the code below
36  } else if (abs_u_m < 0.6*u_th) {
37  // Mean velocity magnitude is less than thermal velocity
38  // Use the distribution u*exp(-u**2*(1-abs(u_m)/u_th)/(2*u_th**2)) as an approximation
39  // and then use the rejection method to correct it
40  // ( stop rejecting with probability exp(-abs(u_m)/(2*u_th**3)*(u-sign(u_m)*u_th)**2) )
41  // Note that this is the method that is used in the common case u_m=0
42  const amrex::Real umsign = std::copysign(1._rt, u_m);
43  const amrex::Real approx_u_th = u_th/std::sqrt( 1._rt - abs_u_m/u_th );
44  const amrex::Real reject_prefactor = (abs_u_m/u_th)/(2._rt*u_th*u_th); // To save computation
45  bool reject = true;
46  while (reject) {
47  // Generates u according to u*exp(-u**2/(2*approx_u_th**2)),
48  // using the method of the inverse cumulative function
49  amrex::Real xrand = 1._rt - amrex::Random(engine); // ensures urand > 0
50  u = approx_u_th * std::sqrt(2._rt*std::log(1._rt/xrand));
51  // Rejection method
52  xrand = amrex::Random(engine);
53  if (xrand < std::exp(-reject_prefactor*(u - umsign*u_th)*(u - umsign*u_th))) { reject = false; }
54  }
55  } else {
56  // Mean velocity magnitude is greater than thermal velocity
57  // Use the distribution exp(-(u-u_m-u_th**2/abs(u_m))**2/(2*u_th**2)) as an approximation
58  // and then use the rejection method to correct it
59  // ( stop rejecting with probability (u/abs(u_m))*exp(1-(u/abs(u_m))) ; note
60  // that this number is always between 0 and 1 )
61  // Note that in the common case `u_m = 0`, this rejection method
62  // is not used, and the above rejection method is used instead.
63  bool reject = true;
64  const amrex::Real approx_u_m = u_m + u_th*u_th/abs_u_m;
65  const amrex::Real inv_um = 1._rt/abs_u_m; // To save computation
66  while (reject) {
67  // Approximate distribution: normal distribution, where we only retain positive u
68  u = -1._rt;
69  while (u < 0) {
70  u = amrex::RandomNormal(approx_u_m, u_th, engine);
71  }
72  // Rejection method
73  const amrex::Real xrand = amrex::Random(engine);
74  if (xrand < u*inv_um* std::exp(1._rt - u*inv_um)) { reject = false; }
75  }
76  }
77  return u;
78  }
79 }
80 
81 #endif //WARPX_SAMPLE_GAUSSIAN_FLUX_DISTRIBUTION_H
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
Real Random()
Real RandomNormal(Real mean, Real stddev)