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 
13 #include <ablastr/constant.H>
14 
15 #include <AMReX_Random.H>
16 #include <AMReX_REAL.H>
17 
18 #include <cmath>
19 
20 
21 namespace impactx::distribution
22 {
23  struct Triangle
24  {
37  amrex::ParticleReal const sigx, amrex::ParticleReal const sigy,
38  amrex::ParticleReal const sigt, amrex::ParticleReal const sigpx,
39  amrex::ParticleReal const sigpy, amrex::ParticleReal const sigpt,
40  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
41  amrex::ParticleReal const mutpt=0.0
42  )
43  : m_sigmaX(sigx), m_sigmaY(sigy), m_sigmaT(sigt), m_sigmaPx(sigpx), m_sigmaPy(sigpy),
44  m_sigmaPt(sigpt), m_muxpx(muxpx), m_muypy(muypy), m_mutpt(mutpt)
45  {
46  }
47 
60  amrex::ParticleReal & x,
61  amrex::ParticleReal & y,
62  amrex::ParticleReal & t,
63  amrex::ParticleReal & px,
64  amrex::ParticleReal & py,
65  amrex::ParticleReal & pt,
66  amrex::RandomEngine const& engine
67  ) const
68  {
69 
70  using namespace amrex::literals;
72 
73  amrex::ParticleReal ln1, norm, u0, u1, u2;
74  amrex::ParticleReal g1, g2, g3, g4, g5;
75  amrex::ParticleReal d, root, a1, a2;
76 
77  // Sample the t coordinate for a ramped triangular profile (unit
78  // variance):
79  u0 = amrex::Random(engine);
80  t = sqrt(2_prt)*(2_prt-3_prt*sqrt(u0));
81 
82  // Generate five standard normal random variables using Box-Muller:
83  u1 = amrex::Random(engine);
84  u2 = amrex::Random(engine);
85  ln1 = sqrt(-2_prt*log(u1));
86  g1 = ln1*cos(2_prt*pi*u2);
87  g2 = ln1*sin(2_prt*pi*u2);
88  u1 = amrex::Random(engine);
89  u2 = amrex::Random(engine);
90  ln1 = sqrt(-2_prt*log(u1));
91  g3 = ln1*cos(2_prt*pi*u2);
92  g4 = ln1*sin(2_prt*pi*u2);
93  u1 = amrex::Random(engine);
94  u2 = amrex::Random(engine);
95  ln1 = sqrt(-2_prt*log(u1));
96  g5 = ln1*cos(2_prt*pi*u2);
97 
98  // Use one of these normal random variables for pt:
99  pt = g5;
100 
101  // Normalize the rest to produce uniform samples on the unit sphere:
102  norm = sqrt(g1*g1+g2*g2+g3*g3+g4*g4);
103  g1 /= norm;
104  g2 /= norm;
105  g3 /= norm;
106  g4 /= norm;
107 
108  // Scale to produce uniform samples in a 4D ball (unit variance):
109  d = 4_prt; // unit ball dimension
110  u1 = amrex::Random(engine); // uniform sample
111  u2 = sqrt(d+2_prt)*pow(u1,1_prt/d);
112  x = g1*u2;
113  y = g2*u2;
114  px = g3*u2;
115  py = g4*u2;
116 
117  // Transform to produce the desired second moments/correlations:
118  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
119  a1 = m_sigmaX*x/root;
120  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
121  x = a1;
122  px = a2;
123  root = sqrt(1.0_prt-m_muypy*m_muypy);
124  a1 = m_sigmaY*y/root;
125  a2 = m_sigmaPy*(-m_muypy*y/root+py);
126  y = a1;
127  py = a2;
128  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
129  a1 = m_sigmaT*t/root;
130  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
131  t = a1;
132  pt = a2;
133  }
134  private:
135  amrex::ParticleReal m_sigmaX, m_sigmaY, m_sigmaT;
136  amrex::ParticleReal m_sigmaPx, m_sigmaPy, m_sigmaPt;
137  amrex::ParticleReal m_muxpx, m_muypy, m_mutpt;
138  };
139 
140 } // namespace impactx::distribution
141 
142 #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:26
Definition: Triangle.H:24
amrex::ParticleReal m_sigmaY
Definition: Triangle.H:135
amrex::ParticleReal m_sigmaPt
Definition: Triangle.H:136
amrex::ParticleReal m_sigmaPy
Definition: Triangle.H:136
amrex::ParticleReal m_mutpt
Definition: Triangle.H:137
amrex::ParticleReal m_muypy
Definition: Triangle.H:137
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Triangle.H:136
amrex::ParticleReal m_sigmaT
Definition: Triangle.H:135
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:59
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:36
amrex::ParticleReal m_sigmaX
Definition: Triangle.H:135
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Triangle.H:137