ImpactX
Waterbag.H
Go to the documentation of this file.
1 /* Copyright 2022-2023 The Regents of the University of California, through Lawrence
2  * Berkeley National Laboratory (subject to receipt of any required
3  * approvals from the U.S. Dept. of Energy). All rights reserved.
4  *
5  * This file is part of ImpactX.
6  *
7  * Authors: Ji Qiang, Axel Huebl
8  * License: BSD-3-Clause-LBNL
9  */
10 #ifndef IMPACTX_DISTRIBUTION_WATERBAG
11 #define IMPACTX_DISTRIBUTION_WATERBAG
12 
13 #include <AMReX_Random.H>
14 #include <AMReX_REAL.H>
15 
16 #include <cmath>
17 
18 
19 namespace impactx::distribution
20 {
21  struct Waterbag
22  {
32  Waterbag(amrex::ParticleReal const sigx, amrex::ParticleReal const sigy,
33  amrex::ParticleReal const sigt,amrex::ParticleReal const sigpx,
34  amrex::ParticleReal const sigpy,amrex::ParticleReal const sigpt,
35  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
36  amrex::ParticleReal const mutpt=0.0
37  )
38  : m_sigmaX(sigx),m_sigmaY(sigy),m_sigmaT(sigt),m_sigmaPx(sigpx),m_sigmaPy(sigpy),
39  m_sigmaPt(sigpt),m_muxpx(muxpx),m_muypy(muypy),m_mutpt(mutpt)
40  {
41  }
42 
55  amrex::ParticleReal & x,
56  amrex::ParticleReal & y,
57  amrex::ParticleReal & t,
58  amrex::ParticleReal & px,
59  amrex::ParticleReal & py,
60  amrex::ParticleReal & pt,
61  amrex::RandomEngine const& engine) const {
62 
63  using namespace amrex::literals;
64 
65  amrex::ParticleReal ln1,norm,u1,u2;
66  amrex::ParticleReal g1,g2,g3,g4,g5,g6;
67  amrex::ParticleReal root,a1,a2;
68 
69  constexpr amrex::ParticleReal pi = 3.14159265358979_prt;
70 
71  // Generate six standard normal random variables using Box-Muller:
72  u1 = amrex::Random(engine);
73  u2 = amrex::Random(engine);
74  ln1 = sqrt(-2_prt*log(u1));
75  g1 = ln1*cos(2_prt*pi*u2);
76  g2 = ln1*sin(2_prt*pi*u2);
77  u1 = amrex::Random(engine);
78  u2 = amrex::Random(engine);
79  ln1 = sqrt(-2_prt*log(u1));
80  g3 = ln1*cos(2_prt*pi*u2);
81  g4 = ln1*sin(2_prt*pi*u2);
82  u1 = amrex::Random(engine);
83  u2 = amrex::Random(engine);
84  ln1 = sqrt(-2_prt*log(u1));
85  g5 = ln1*cos(2_prt*pi*u2);
86  g6 = ln1*sin(2_prt*pi*u2);
87 
88  // Normalize to produce uniform samples on the unit sphere:
89  norm = sqrt(g1*g1+g2*g2+g3*g3+g4*g4+g5*g5+g6*g6);
90  g1 /= norm;
91  g2 /= norm;
92  g3 /= norm;
93  g4 /= norm;
94  g5 /= norm;
95  g6 /= norm;
96 
97  // Scale to produce uniform samples in a ball (unit variance):
98  u1 = amrex::Random(engine);
99  u2 = sqrt(8_prt)*pow(u1,1_prt/6);
100  x = g1*u2;
101  y = g2*u2;
102  t = g3*u2;
103  px = g4*u2;
104  py = g5*u2;
105  pt = g6*u2;
106 
107  // Transform to produce the desired second moments/correlations:
108  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
109  a1 = m_sigmaX*x/root;
110  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
111  x = a1;
112  px = a2;
113  root = sqrt(1.0_prt-m_muypy*m_muypy);
114  a1 = m_sigmaY*y/root;
115  a2 = m_sigmaPy*(-m_muypy*y/root+py);
116  y = a1;
117  py = a2;
118  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
119  a1 = m_sigmaT*t/root;
120  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
121  t = a1;
122  pt = a2;
123  }
124  private:
125  amrex::ParticleReal m_sigmaX,m_sigmaY,m_sigmaT;
126  amrex::ParticleReal m_sigmaPx,m_sigmaPy,m_sigmaPt;
127  amrex::ParticleReal m_muxpx,m_muypy,m_mutpt;
128  };
129 
130 } // namespace impactx::distribution
131 
132 #endif // IMPACTX_DISTRIBUTION_WATERBAG
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
constexpr std::enable_if_t< std::is_floating_point< T >::value, T > pi()
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > log(const GpuComplex< T > &a_z) noexcept
Real Random()
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T norm(const GpuComplex< T > &a_z) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, const T &a_y) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Definition: All.H:26
Definition: Waterbag.H:22
amrex::ParticleReal m_muypy
Definition: Waterbag.H:127
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Waterbag.H:126
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Waterbag.H:127
amrex::ParticleReal m_mutpt
Definition: Waterbag.H:127
amrex::ParticleReal m_sigmaPt
Definition: Waterbag.H:126
amrex::ParticleReal m_sigmaT
Definition: Waterbag.H:125
amrex::ParticleReal m_sigmaY
Definition: Waterbag.H:125
amrex::ParticleReal m_sigmaX
Definition: Waterbag.H:125
Waterbag(amrex::ParticleReal const sigx, amrex::ParticleReal const sigy, amrex::ParticleReal const sigt, amrex::ParticleReal const sigpx, amrex::ParticleReal const sigpy, amrex::ParticleReal const sigpt, amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0, amrex::ParticleReal const mutpt=0.0)
Definition: Waterbag.H:32
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()(amrex::ParticleReal &x, amrex::ParticleReal &y, amrex::ParticleReal &t, amrex::ParticleReal &px, amrex::ParticleReal &py, amrex::ParticleReal &pt, amrex::RandomEngine const &engine) const
Definition: Waterbag.H:54
amrex::ParticleReal m_sigmaPy
Definition: Waterbag.H:126