ImpactX
Thermal.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, Chad Mitchell
8  * License: BSD-3-Clause-LBNL
9  */
10 #ifndef IMPACTX_DISTRIBUTION_THERMAL
11 #define IMPACTX_DISTRIBUTION_THERMAL
12 
14 
15 #include <ablastr/constant.H>
16 
17 #include <AMReX_GpuContainers.H>
18 #include <AMReX_Math.H>
19 #include <AMReX_Print.H>
20 #include <AMReX_Random.H>
21 #include <AMReX_REAL.H>
22 
23 #include <cmath>
24 #include <cstdint>
25 #include <random>
26 #include <utility>
27 #include <vector>
28 
29 
30 namespace impactx
31 {
32 namespace distribution
33 {
34  struct ThermalData
35  {
36  static constexpr amrex::ParticleReal tolerance = 1.0e-3;
37  static constexpr amrex::ParticleReal rin = 1.0e-10;
38  static constexpr amrex::ParticleReal rout = 10.0;
39  static constexpr int nsteps = 2000;
40 
41  amrex::ParticleReal m_f1;
42  amrex::ParticleReal m_f2;
43  amrex::ParticleReal m_phi1;
44  amrex::ParticleReal m_phi2;
45  amrex::ParticleReal m_p1;
46  amrex::ParticleReal m_p2;
47  amrex::ParticleReal m_rmin;
48  amrex::ParticleReal m_rmax;
49  int m_nbins;
50  amrex::ParticleReal* m_cdf1 = nullptr;
51  amrex::ParticleReal* m_cdf2 = nullptr;
52  amrex::ParticleReal m_Cintensity;
53  amrex::ParticleReal m_bg;
54  amrex::ParticleReal m_k;
55  amrex::ParticleReal m_T1;
56  amrex::ParticleReal m_T2;
57  amrex::ParticleReal m_w;
58 
59  // GPU data
62 
64  amrex::ParticleReal kin,
65  amrex::ParticleReal T1in,
66  amrex::ParticleReal T2in,
67  amrex::ParticleReal p1in,
68  amrex::ParticleReal p2in,
69  amrex::ParticleReal win
70  )
71  {
72  m_k = kin;
73  m_T1 = T1in;
74  m_T2 = T2in;
75  m_p1 = p1in;
76  m_p2 = p2in;
77  m_w = win;
78  }
79 
85  void
86  generate_radial_dist (amrex::ParticleReal bunch_charge, RefPart const & refpart)
87  {
88  using namespace amrex::literals; // for _rt and _prt
89  using namespace ablastr::constant::math; // for pi
90 
91  // Get relativistic beta*gamma
92  amrex::ParticleReal bg = refpart.beta_gamma();
93  m_bg = bg;
94 
95  // Get reference particle rest energy (eV) and charge (q_e)
96  amrex::ParticleReal Erest = refpart.mass_MeV()*1.0e6;
97  amrex::ParticleReal q_e = refpart.charge_qe();
98 
99  // Set space charge intensity
100  m_Cintensity = q_e*bunch_charge/(pow(bg,2)*Erest*ablastr::constant::SI::ep0);
101 
102  // Set minimum and maximum radius
103  amrex::ParticleReal r_scale = matched_scale_radius();
104  amrex::ParticleReal rmin = rin*r_scale;
105  amrex::ParticleReal rmax = rout*r_scale;
106  // amrex::PrintToFile("equilibrium_params.out") << r_scale << " " << data.Cintensity << "\n";
107 
108  // Scale the parameters p1 and p2
109  amrex::ParticleReal rt2pi = sqrt(2.0_prt*pi);
110  amrex::ParticleReal p_scale = pow(r_scale*rt2pi,-3);
111  m_p1 = m_p1*p_scale;
112  m_p2 = m_p2*p_scale;
113 
114  // store integration parameters
115  m_nbins = nsteps;
116  m_rmin = rmin;
117  m_rmax = rmax;
118 
119  // allocate CDFs
120  std::vector<amrex::ParticleReal> cdf1(m_nbins+1);
121  std::vector<amrex::ParticleReal> cdf2(m_nbins+1);
122  m_cdf1 = cdf1.data();
123  m_cdf2 = cdf2.data();
124 
125  // set initial conditions
126  m_f1 = 0.0_prt;
127  m_f2 = 0.0_prt;
128  m_phi1 = 0.0_prt;
129  m_phi2 = 0.0_prt;
130  m_cdf1[0] = 0.0_prt;
131  m_cdf2[0] = 0.0_prt;
132 
133  // integrate ODEs for the radial density
134  integrate(rmin,rmax,nsteps);
135 
136  // a search over normalization parameters p1, p2 can be added here
137 
138  // rescale cdf to ensure the exact range [0,1]
139  for (int n = 0; n < m_nbins; ++n) {
140  m_cdf1[n] = m_cdf1[n] / m_cdf1[m_nbins];
141  m_cdf2[n] = m_cdf2[n] / m_cdf2[m_nbins];
142  }
143 
144  // copy CFD data to device
148  cdf1.begin(), cdf1.end(),
149  m_d_cdf1.begin());
151  cdf2.begin(), cdf2.end(),
152  m_d_cdf2.begin());
154  }
155 
156  amrex::ParticleReal
158  {
159  using namespace amrex::literals; // for _rt and _prt
160  using namespace ablastr::constant::math; // for pi
161 
162  amrex::ParticleReal k = m_k;
163  amrex::ParticleReal kT = (1.0_prt - m_w) * m_T1 + m_w * m_T2;
164  amrex::ParticleReal a = m_Cintensity/(4.0_prt*pi*5.0_prt*sqrt(5.0_prt));
165  amrex::ParticleReal rscale = sqrt(kT + pow(a*k,2.0/3.0))/k;
166 
167  return rscale;
168  }
169 
170  void
172  amrex::ParticleReal in,
173  amrex::ParticleReal out,
174  int steps
175  )
176  {
177  using namespace amrex::literals; // for _rt and _prt
178 
179  // initialize numerical integration parameters
180  amrex::ParticleReal const tau = (out-in)/steps;
181  amrex::ParticleReal const half = tau/2.0_prt;
182 
183  // initialize the value of the independent variable
184  amrex::ParticleReal reval = in;
185 
186  // loop over integration steps
187  for (int j=0; j < steps; ++j)
188  {
189  // for a second-order Strang splitting
190  map1(half, reval);
191  map2(tau, reval);
192  map1(half, reval);
193  // store tabulated CDF
194  m_cdf1[j+1] = m_f1;
195  m_cdf2[j+1] = m_f2;
196  // write cumulative density to file for debugging
197  /* comment in for debugging
198  amrex::PrintToFile("density1.out") << reval << " " << data.f1 << "\n";
199  amrex::PrintToFile("density2.out") << reval << " " << data.f2 << "\n";
200  amrex::PrintToFile("phi1.out") << reval << " " << data.phi1 << "\n";
201  amrex::PrintToFile("phi2.out") << reval << " " << data.phi2 << "\n";
202  */
203  }
204 
205  }
206 
207  void
209  amrex::ParticleReal const tau,
210  amrex::ParticleReal & reval
211  )
212  {
213  using namespace amrex::literals; // for _rt and _prt
214  using namespace ablastr::constant::math; // for pi
215 
216  amrex::ParticleReal const f1 = m_f1;
217  amrex::ParticleReal const f2 = m_f2;
218  amrex::ParticleReal const phi1 = m_phi1;
219  amrex::ParticleReal const phi2 = m_phi2;
220  amrex::ParticleReal const r = reval;
221 
222  // Apply map to update phi1, phi2, and r:
223  reval = r + tau;
224  m_phi1 = phi1 + f1/(4.0_prt*pi*reval) - f1/(4.0_prt*pi*r);
225  m_phi2 = phi2 + f2/(4.0_prt*pi*reval) - f2/(4.0_prt*pi*r);;
226  m_f1 = f1;
227  m_f2 = f2;
228  }
229 
230  void
232  amrex::ParticleReal const tau,
233  amrex::ParticleReal & reval
234  )
235  {
236  using namespace amrex::literals; // for _rt and _prt
237  using namespace ablastr::constant::math; // for pi
238 
239  amrex::ParticleReal const f1 = m_f1;
240  amrex::ParticleReal const f2 = m_f2;
241  amrex::ParticleReal const phi1 = m_phi1;
242  amrex::ParticleReal const phi2 = m_phi2;
243  amrex::ParticleReal const r = reval;
244  amrex::ParticleReal const k = m_k;
245  amrex::ParticleReal const w = m_w;
246  amrex::ParticleReal const T1 = m_T1;
247  amrex::ParticleReal const T2 = m_T2;
248 
249  // Define space charge intensity parameters
250  amrex::ParticleReal const c1 = (1.0_prt-w) * m_Cintensity;
251  amrex::ParticleReal const c2 = w * m_Cintensity;
252 
253  // Define intermediate quantities
254  amrex::ParticleReal potential = 0.0_prt;
255  potential = pow(k*r,2.0)/2.0_prt + c1*phi1 + c2*phi2;
256  amrex::ParticleReal Pdensity1 = m_p1*exp(-potential/T1);
257  amrex::ParticleReal Pdensity2 = m_p2*exp(-potential/T2);
258  // amrex::ParticleReal Pdensity_tot = (1.0_prt-w)*Pdensity1 + w*Pdensity2;
259  // amrex::PrintToFile("Pdensity.out") << reval << " " << Pdensity_tot << "\n";
260 
261  // Apply map to update f1 and f2:
262  m_phi1 = phi1;
263  m_phi2 = phi2;
264  m_f1 = f1 + tau*4.0_prt*pi*pow(r,2.0)*Pdensity1;
265  m_f2 = f2 + tau*4.0_prt*pi*pow(r,2.0)*Pdensity2;
266  reval = r;
267  }
268  };
269 
270  struct Thermal
271  {
286  amrex::ParticleReal k,
287  amrex::ParticleReal kT,
288  amrex::ParticleReal kT_halo,
289  amrex::ParticleReal normalize,
290  amrex::ParticleReal normalize_halo,
291  amrex::ParticleReal halo
292  )
293  : m_k(k),
294  m_T1(kT),
295  m_T2(kT_halo),
296  m_normalize(normalize),
297  m_normalize_halo(normalize_halo),
298  m_halo(halo)
299  {
300  }
301 
309  void initialize (amrex::ParticleReal bunch_charge, RefPart const & ref)
310  {
311  // Generate the struct "data" containing the radial CDF:
312  ThermalData data(m_k, m_T1, m_T2, m_normalize, m_normalize_halo, m_halo);
313  data.generate_radial_dist(bunch_charge, ref);
314 
315  // share data
316  m_w = data.m_w;
317  m_T1 = data.m_T1;
318  m_T2 = data.m_T2;
319  m_rmin = data.m_rmin;
320  m_rmax = data.m_rmax;
321  m_nbins = data.m_nbins;
322  m_bg = data.m_bg;
323  m_cdf1 = data.m_d_cdf1.data();
324  m_cdf2 = data.m_d_cdf2.data();
325  }
326 
329  void
331  {
332  // deallocate
333  ThermalData::m_d_cdf1.clear();
334  ThermalData::m_d_cdf2.clear();
335  }
336 
348  void operator() (
349  amrex::ParticleReal & AMREX_RESTRICT x,
350  amrex::ParticleReal & AMREX_RESTRICT y,
351  amrex::ParticleReal & AMREX_RESTRICT t,
352  amrex::ParticleReal & AMREX_RESTRICT px,
353  amrex::ParticleReal & AMREX_RESTRICT py,
354  amrex::ParticleReal & AMREX_RESTRICT pt,
355  amrex::RandomEngine const & engine
356  ) const
357  {
358 
359  using namespace amrex::literals; // for _rt and _prt
360  using namespace ablastr::constant::math; // for pi
361 
362  amrex::ParticleReal ln1,norm,u1,u2,uhalo;
363  amrex::ParticleReal g1,g2,g3,g4,g5,g6;
364  amrex::ParticleReal z,pz;
365 
366  // Use a Bernoulli random variable to select between core and halo:
367  // If u < w, the particle is in the halo population.
368  // If u >= w, the particle is in the core population.
369  uhalo = amrex::Random(engine);
370 
371  // Generate six standard normal random variables using Box-Muller:
372  u1 = amrex::Random(engine);
373  u2 = amrex::Random(engine);
374  ln1 = sqrt(-2_prt*log(u1));
375  g1 = ln1*cos(2_prt*pi*u2);
376  g2 = ln1*sin(2_prt*pi*u2);
377  u1 = amrex::Random(engine);
378  u2 = amrex::Random(engine);
379  ln1 = sqrt(-2_prt*log(u1));
380  g3 = ln1*cos(2_prt*pi*u2);
381  g4 = ln1*sin(2_prt*pi*u2);
382  u1 = amrex::Random(engine);
383  u2 = amrex::Random(engine);
384  ln1 = sqrt(-2_prt*log(u1));
385  g5 = ln1*cos(2_prt*pi*u2);
386  g6 = ln1*sin(2_prt*pi*u2);
387 
388  // Scale the last three variables to produce the momenta:
389  amrex::ParticleReal kT = (uhalo > m_w) ? m_T1 : m_T2; // select core or halo value
390  px = sqrt(kT)*g4;
391  py = sqrt(kT)*g5;
392  pz = sqrt(kT)*g6;
393 
394  // Normalize the first three variables to produce uniform samples on the unit 3-sphere:
395  norm = sqrt(g1*g1+g2*g2+g3*g3);
396  g1 /= norm;
397  g2 /= norm;
398  g3 /= norm;
399 
400  // Collect the radial CDF returned by generate_radial_dist:
401 
402  // select core or halo CDF
403  amrex::ParticleReal const * cdf = (uhalo > m_w) ? m_cdf1 : m_cdf2;
404 
405  // Generate a radial coordinate from the CDF
406  amrex::ParticleReal u = amrex::Random(engine);
407  amrex::ParticleReal const * ptr = amrex::lower_bound(cdf, cdf + m_nbins + 1, u);
408  int const off = amrex::max(0, int(ptr - cdf - 1));
409  amrex::ParticleReal tv = (u - cdf[off]) / (cdf[off + 1] - cdf[off]);
410  amrex::ParticleReal xv = (off + tv) / amrex::ParticleReal(m_nbins);
411  amrex::ParticleReal r = m_rmin + (m_rmax - m_rmin) * xv;
412 
413  // Scale to produce samples with the correct radial density:
414  x = g1*r;
415  y = g2*r;
416  z = g3*r;
417 
418  // Transform longitudinal variables into the laboratory frame:
419  t = -z / m_bg;
420  pt = -pz * m_bg;
421  }
422 
423  private:
424  amrex::ParticleReal m_k;
425  amrex::ParticleReal m_T1, m_T2;
426  amrex::ParticleReal m_normalize, m_normalize_halo;
427  amrex::ParticleReal m_halo;
428  amrex::ParticleReal m_rmin;
429  amrex::ParticleReal m_rmax;
430  int m_nbins;
431  amrex::ParticleReal m_bg;
432  amrex::ParticleReal m_w;
433 
434  // radial profile data
435  amrex::ParticleReal const * m_cdf1 = nullptr;
436  amrex::ParticleReal const * m_cdf2 = nullptr;
437  };
438 
439 } // namespace distribution
440 } // namespace impactx
441 
442 #endif // IMPACTX_DISTRIBUTION_THERMAL
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
static constexpr auto ep0
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
static constexpr HostToDevice hostToDevice
void streamSynchronize() noexcept
constexpr std::enable_if_t< std::is_floating_point< T >::value, T > pi()
real(kind=amrex_real), parameter half
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > log(const GpuComplex< T > &a_z) noexcept
Real Random()
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
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 > exp(const GpuComplex< T > &a_z) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
AMREX_GPU_HOST_DEVICE ItType lower_bound(ItType first, ItType last, const ValType &val)
Definition: ImpactX.cpp:33
@ t
fixed t as the independent variable
data
int n
tuple w
Definition: ReferenceParticle.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal beta_gamma() const
Definition: ReferenceParticle.H:79
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal mass_MeV() const
Definition: ReferenceParticle.H:94
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal charge_qe() const
Definition: ReferenceParticle.H:186
Definition: Thermal.H:35
void generate_radial_dist(amrex::ParticleReal bunch_charge, RefPart const &refpart)
Definition: Thermal.H:86
amrex::ParticleReal m_bg
reference value of relativistic beta*gamma
Definition: Thermal.H:53
static amrex::Gpu::DeviceVector< amrex::ParticleReal > m_d_cdf1
Definition: Thermal.H:60
amrex::ParticleReal m_phi2
potential generated by second population
Definition: Thermal.H:44
void map1(amrex::ParticleReal const tau, amrex::ParticleReal &reval)
Definition: Thermal.H:208
amrex::ParticleReal m_phi1
potential generated by first population
Definition: Thermal.H:43
ThermalData(amrex::ParticleReal kin, amrex::ParticleReal T1in, amrex::ParticleReal T2in, amrex::ParticleReal p1in, amrex::ParticleReal p2in, amrex::ParticleReal win)
Definition: Thermal.H:63
amrex::ParticleReal m_f1
cumulative distribution of first population
Definition: Thermal.H:41
void map2(amrex::ParticleReal const tau, amrex::ParticleReal &reval)
Definition: Thermal.H:231
amrex::ParticleReal * m_cdf2
tabulated cumulative distribution (second)
Definition: Thermal.H:51
amrex::ParticleReal * m_cdf1
tabulated cumulative distribution (first)
Definition: Thermal.H:50
amrex::ParticleReal m_Cintensity
space charge intensity parameter
Definition: Thermal.H:52
static amrex::Gpu::DeviceVector< amrex::ParticleReal > m_d_cdf2
Definition: Thermal.H:61
amrex::ParticleReal m_rmin
minimum r value for tabulated cdf
Definition: Thermal.H:47
static constexpr amrex::ParticleReal rin
initial r value for numerical integration
Definition: Thermal.H:37
void integrate(amrex::ParticleReal in, amrex::ParticleReal out, int steps)
Definition: Thermal.H:171
amrex::ParticleReal m_rmax
maximum r value for tabulated cdf
Definition: Thermal.H:48
static constexpr amrex::ParticleReal tolerance
tolerance for matching condition
Definition: Thermal.H:36
static constexpr int nsteps
number of radial steps for numerical integration
Definition: Thermal.H:39
amrex::ParticleReal matched_scale_radius()
Definition: Thermal.H:157
amrex::ParticleReal m_k
linear focusing strength (1/meters)
Definition: Thermal.H:54
amrex::ParticleReal m_f2
cumulative distribution of second population
Definition: Thermal.H:42
amrex::ParticleReal m_p2
normalization constant of second population
Definition: Thermal.H:46
int m_nbins
number of radial bins for tabulated cdf
Definition: Thermal.H:49
amrex::ParticleReal m_p1
normalization constant of first population
Definition: Thermal.H:45
amrex::ParticleReal m_T1
temperature k*T of the primary (core) population
Definition: Thermal.H:55
static constexpr amrex::ParticleReal rout
final r value for numerical integration
Definition: Thermal.H:38
amrex::ParticleReal m_T2
temperature k*T of the secondary (halo) population
Definition: Thermal.H:56
amrex::ParticleReal m_w
weight of the secondary (halo) population
Definition: Thermal.H:57
Definition: Thermal.H:271
void finalize()
Definition: Thermal.H:330
Thermal(amrex::ParticleReal k, amrex::ParticleReal kT, amrex::ParticleReal kT_halo, amrex::ParticleReal normalize, amrex::ParticleReal normalize_halo, amrex::ParticleReal halo)
Definition: Thermal.H:285
amrex::ParticleReal m_bg
reference value of relativistic beta*gamma
Definition: Thermal.H:431
int m_nbins
number of radial bins for tabulated cdf
Definition: Thermal.H:430
amrex::ParticleReal m_T1
linear focusing strength (1/meters)
Definition: Thermal.H:425
void initialize(amrex::ParticleReal bunch_charge, RefPart const &ref)
Definition: Thermal.H:309
amrex::ParticleReal m_rmin
relative weight of halo population
Definition: Thermal.H:428
amrex::ParticleReal m_normalize
temperature of each particle population
Definition: Thermal.H:426
amrex::ParticleReal m_w
weight of the secondary (halo) population
Definition: Thermal.H:432
amrex::ParticleReal m_halo
normalization constant of first/second population
Definition: Thermal.H:427
amrex::ParticleReal m_rmax
maximum r value for tabulated cdf
Definition: Thermal.H:429
amrex::ParticleReal m_k
Definition: Thermal.H:424