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 lambdax, amrex::ParticleReal const lambday,
40  amrex::ParticleReal const lambdat, amrex::ParticleReal const lambdapx,
41  amrex::ParticleReal const lambdapy, amrex::ParticleReal const lambdapt,
42  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
43  amrex::ParticleReal const mutpt=0.0
44  )
45  : m_lambdaX(lambdax), m_lambdaY(lambday), m_lambdaT(lambdat), m_lambdaPx(lambdapx), m_lambdaPy(lambdapy),
46  m_lambdaPt(lambdapt), 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_lambdaX * x / root;
142  a2 = m_lambdaPx * (-m_muxpx * x / root + px);
143  x = a1;
144  px = a2;
145  root = sqrt(1.0_prt-m_muypy*m_muypy);
146  a1 = m_lambdaY * y / root;
147  a2 = m_lambdaPy * (-m_muypy * y / root + py);
148  y = a1;
149  py = a2;
150  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
151  a1 = m_lambdaT * t / root;
152  a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
153  t = a1;
154  pt = a2;
155  }
156  private:
157  amrex::ParticleReal m_lambdaX, m_lambdaY, m_lambdaT;
158  amrex::ParticleReal m_lambdaPx, m_lambdaPy, m_lambdaPt;
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_v< T >, 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_lambdaPy
Definition: Triangle.H:158
amrex::ParticleReal m_lambdaPx
related position axis intercepts (length) of the phase space ellipse
Definition: Triangle.H:158
amrex::ParticleReal m_lambdaT
Definition: Triangle.H:157
amrex::ParticleReal m_lambdaX
Definition: Triangle.H:157
amrex::ParticleReal m_mutpt
Definition: Triangle.H:159
amrex::ParticleReal m_lambdaY
Definition: Triangle.H:157
amrex::ParticleReal m_muypy
Definition: Triangle.H:159
amrex::ParticleReal m_lambdaPt
Definition: Triangle.H:158
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 lambdax, amrex::ParticleReal const lambday, amrex::ParticleReal const lambdat, amrex::ParticleReal const lambdapx, amrex::ParticleReal const lambdapy, amrex::ParticleReal const lambdapt, amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0, amrex::ParticleReal const mutpt=0.0)
Definition: Triangle.H:38
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
related momentum axis intercepts of the phase space ellipse
Definition: Triangle.H:159