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 
14 
15 #include <ablastr/constant.H>
16 
17 #include <AMReX_Random.H>
18 #include <AMReX_REAL.H>
19 
20 #include <cmath>
21 
22 
23 namespace impactx::distribution
24 {
25  struct Waterbag
26  {
37  amrex::ParticleReal lambdax,
38  amrex::ParticleReal lambday,
39  amrex::ParticleReal lambdat,
40  amrex::ParticleReal lambdapx,
41  amrex::ParticleReal lambdapy,
42  amrex::ParticleReal lambdapt,
43  amrex::ParticleReal muxpx=0.0,
44  amrex::ParticleReal muypy=0.0,
45  amrex::ParticleReal mutpt=0.0
46  )
47  : m_lambdaX(lambdax), m_lambdaY(lambday), m_lambdaT(lambdat), m_lambdaPx(lambdapx), m_lambdaPy(lambdapy),
48  m_lambdaPt(lambdapt), m_muxpx(muxpx), m_muypy(muypy), m_mutpt(mutpt)
49  {
50  }
51 
59  void initialize ([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const & ref)
60  {
61  }
62 
67  void
69  {
70  }
71 
84  amrex::ParticleReal & AMREX_RESTRICT x,
85  amrex::ParticleReal & AMREX_RESTRICT y,
86  amrex::ParticleReal & AMREX_RESTRICT t,
87  amrex::ParticleReal & AMREX_RESTRICT px,
88  amrex::ParticleReal & AMREX_RESTRICT py,
89  amrex::ParticleReal & AMREX_RESTRICT pt,
90  amrex::RandomEngine const & engine
91  ) const
92  {
93  using namespace amrex::literals;
95 
96  amrex::ParticleReal ln1,norm,u1,u2;
97  amrex::ParticleReal g1,g2,g3,g4,g5,g6;
98  amrex::ParticleReal root,a1,a2;
99 
100  // Generate six standard normal random variables using Box-Muller:
101  u1 = amrex::Random(engine);
102  u2 = amrex::Random(engine);
103  ln1 = std::sqrt(-2_prt*std::log(u1));
104  g1 = ln1*std::cos(2_prt*pi*u2);
105  g2 = ln1*std::sin(2_prt*pi*u2);
106  u1 = amrex::Random(engine);
107  u2 = amrex::Random(engine);
108  ln1 = std::sqrt(-2_prt*std::log(u1));
109  g3 = ln1*std::cos(2_prt*pi*u2);
110  g4 = ln1*std::sin(2_prt*pi*u2);
111  u1 = amrex::Random(engine);
112  u2 = amrex::Random(engine);
113  ln1 = std::sqrt(-2_prt*std::log(u1));
114  g5 = ln1*std::cos(2_prt*pi*u2);
115  g6 = ln1*std::sin(2_prt*pi*u2);
116 
117  // Normalize to produce uniform samples on the unit sphere:
118  norm = std::sqrt(g1*g1+g2*g2+g3*g3+g4*g4+g5*g5+g6*g6);
119  g1 /= norm;
120  g2 /= norm;
121  g3 /= norm;
122  g4 /= norm;
123  g5 /= norm;
124  g6 /= norm;
125 
126  // Scale to produce uniform samples in a ball (unit variance):
127  u1 = amrex::Random(engine);
128  u2 = std::sqrt(8_prt)*std::pow(u1,1_prt/6_prt);
129  x = g1*u2;
130  y = g2*u2;
131  t = g3*u2;
132  px = g4*u2;
133  py = g5*u2;
134  pt = g6*u2;
135 
136  // Transform to produce the desired second moments/correlations:
137  root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
138  a1 = m_lambdaX * x / root;
139  a2 = m_lambdaPx * (-m_muxpx * x / root + px);
140  x = a1;
141  px = a2;
142  root = std::sqrt(1.0_prt-m_muypy*m_muypy);
143  a1 = m_lambdaY * y / root;
144  a2 = m_lambdaPy * (-m_muypy * y / root + py);
145  y = a1;
146  py = a2;
147  root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
148  a1 = m_lambdaT * t / root;
149  a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
150  t = a1;
151  pt = a2;
152  }
153 
154  private:
155  amrex::ParticleReal m_lambdaX,m_lambdaY,m_lambdaT;
156  amrex::ParticleReal m_lambdaPx,m_lambdaPy,m_lambdaPt;
157  amrex::ParticleReal m_muxpx,m_muypy,m_mutpt;
158  };
159 
160 } // namespace impactx::distribution
161 
162 #endif // IMPACTX_DISTRIBUTION_WATERBAG
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
static constexpr amrex::Real pi
constexpr std::enable_if_t< std::is_floating_point_v< T >, T > pi()
Real Random()
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T norm(const GpuComplex< T > &a_z) noexcept
Definition: All.H:27
@ t
fixed t as the independent variable
Definition: ReferenceParticle.H:30
Definition: Waterbag.H:26
amrex::ParticleReal m_lambdaX
Definition: Waterbag.H:155
amrex::ParticleReal m_muypy
Definition: Waterbag.H:157
amrex::ParticleReal m_muxpx
related momentum axis intercepts of the phase space ellipse
Definition: Waterbag.H:157
amrex::ParticleReal m_mutpt
Definition: Waterbag.H:157
void finalize()
Definition: Waterbag.H:68
Waterbag(amrex::ParticleReal lambdax, amrex::ParticleReal lambday, amrex::ParticleReal lambdat, amrex::ParticleReal lambdapx, amrex::ParticleReal lambdapy, amrex::ParticleReal lambdapt, amrex::ParticleReal muxpx=0.0, amrex::ParticleReal muypy=0.0, amrex::ParticleReal mutpt=0.0)
Definition: Waterbag.H:36
amrex::ParticleReal m_lambdaPt
Definition: Waterbag.H:156
amrex::ParticleReal m_lambdaT
Definition: Waterbag.H:155
amrex::ParticleReal m_lambdaPx
related position axis intercepts (length) of the phase space ellipse
Definition: Waterbag.H:156
amrex::ParticleReal m_lambdaPy
Definition: Waterbag.H:156
amrex::ParticleReal m_lambdaY
Definition: Waterbag.H:155
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()(amrex::ParticleReal &AMREX_RESTRICT x, amrex::ParticleReal &AMREX_RESTRICT y, amrex::ParticleReal &AMREX_RESTRICT t, amrex::ParticleReal &AMREX_RESTRICT px, amrex::ParticleReal &AMREX_RESTRICT py, amrex::ParticleReal &AMREX_RESTRICT pt, amrex::RandomEngine const &engine) const
Definition: Waterbag.H:83
void initialize([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const &ref)
Definition: Waterbag.H:59