ImpactX
Triangle.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_TRIANGLE
11 #define IMPACTX_DISTRIBUTION_TRIANGLE
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 Triangle
26  {
39  amrex::ParticleReal const sigx, amrex::ParticleReal const sigy,
40  amrex::ParticleReal const sigt, amrex::ParticleReal const sigpx,
41  amrex::ParticleReal const sigpy, amrex::ParticleReal const sigpt,
42  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
43  amrex::ParticleReal const mutpt=0.0
44  )
45  : m_sigmaX(sigx), m_sigmaY(sigy), m_sigmaT(sigt), m_sigmaPx(sigpx), m_sigmaPy(sigpy),
46  m_sigmaPt(sigpt), m_muxpx(muxpx), m_muypy(muypy), m_mutpt(mutpt)
47  {
48  }
49 
57  void initialize ([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const & ref)
58  {
59  }
60 
65  void
67  {
68  }
69 
82  amrex::ParticleReal & x,
83  amrex::ParticleReal & y,
84  amrex::ParticleReal & t,
85  amrex::ParticleReal & px,
86  amrex::ParticleReal & py,
87  amrex::ParticleReal & pt,
88  amrex::RandomEngine const& engine
89  ) const
90  {
91 
92  using namespace amrex::literals;
94 
95  amrex::ParticleReal ln1, norm, u0, u1, u2;
96  amrex::ParticleReal g1, g2, g3, g4, g5;
97  amrex::ParticleReal d, root, a1, a2;
98 
99  // Sample the t coordinate for a ramped triangular profile (unit
100  // variance):
101  u0 = amrex::Random(engine);
102  t = sqrt(2_prt)*(2_prt-3_prt*sqrt(u0));
103 
104  // Generate five standard normal random variables using Box-Muller:
105  u1 = amrex::Random(engine);
106  u2 = amrex::Random(engine);
107  ln1 = sqrt(-2_prt*log(u1));
108  g1 = ln1*cos(2_prt*pi*u2);
109  g2 = ln1*sin(2_prt*pi*u2);
110  u1 = amrex::Random(engine);
111  u2 = amrex::Random(engine);
112  ln1 = sqrt(-2_prt*log(u1));
113  g3 = ln1*cos(2_prt*pi*u2);
114  g4 = ln1*sin(2_prt*pi*u2);
115  u1 = amrex::Random(engine);
116  u2 = amrex::Random(engine);
117  ln1 = sqrt(-2_prt*log(u1));
118  g5 = ln1*cos(2_prt*pi*u2);
119 
120  // Use one of these normal random variables for pt:
121  pt = g5;
122 
123  // Normalize the rest to produce uniform samples on the unit sphere:
124  norm = sqrt(g1*g1+g2*g2+g3*g3+g4*g4);
125  g1 /= norm;
126  g2 /= norm;
127  g3 /= norm;
128  g4 /= norm;
129 
130  // Scale to produce uniform samples in a 4D ball (unit variance):
131  d = 4_prt; // unit ball dimension
132  u1 = amrex::Random(engine); // uniform sample
133  u2 = sqrt(d+2_prt)*pow(u1,1_prt/d);
134  x = g1*u2;
135  y = g2*u2;
136  px = g3*u2;
137  py = g4*u2;
138 
139  // Transform to produce the desired second moments/correlations:
140  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
141  a1 = m_sigmaX*x/root;
142  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
143  x = a1;
144  px = a2;
145  root = sqrt(1.0_prt-m_muypy*m_muypy);
146  a1 = m_sigmaY*y/root;
147  a2 = m_sigmaPy*(-m_muypy*y/root+py);
148  y = a1;
149  py = a2;
150  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
151  a1 = m_sigmaT*t/root;
152  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
153  t = a1;
154  pt = a2;
155  }
156  private:
157  amrex::ParticleReal m_sigmaX, m_sigmaY, m_sigmaT;
158  amrex::ParticleReal m_sigmaPx, m_sigmaPy, m_sigmaPt;
159  amrex::ParticleReal m_muxpx, m_muypy, m_mutpt;
160  };
161 
162 } // namespace impactx::distribution
163 
164 #endif // IMPACTX_DISTRIBUTION_TRIANGLE
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
static constexpr amrex::Real pi
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:27
@ t
fixed t as the independent variable
Definition: ReferenceParticle.H:30
Definition: Triangle.H:26
amrex::ParticleReal m_sigmaY
Definition: Triangle.H:157
amrex::ParticleReal m_sigmaPt
Definition: Triangle.H:158
amrex::ParticleReal m_sigmaPy
Definition: Triangle.H:158
amrex::ParticleReal m_mutpt
Definition: Triangle.H:159
amrex::ParticleReal m_muypy
Definition: Triangle.H:159
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Triangle.H:158
amrex::ParticleReal m_sigmaT
Definition: Triangle.H:157
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: Triangle.H:81
Triangle(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: Triangle.H:38
amrex::ParticleReal m_sigmaX
Definition: Triangle.H:157
void initialize([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const &ref)
Definition: Triangle.H:57
void finalize()
Definition: Triangle.H:66
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Triangle.H:159