ImpactX
RFCavity.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: Chad Mitchell, Axel Huebl
8  * License: BSD-3-Clause-LBNL
9  */
10 #ifndef IMPACTX_RFCAVITY_H
11 #define IMPACTX_RFCAVITY_H
12 
15 #include "mixin/beamoptic.H"
16 #include "mixin/thick.H"
17 
18 #include <ablastr/constant.H>
19 
20 #include <AMReX.H>
21 #include <AMReX_Array.H>
22 #include <AMReX_Extension.H>
23 #include <AMReX_REAL.H>
24 
25 #include <array>
26 #include <cmath>
27 #include <stdexcept>
28 #include <tuple>
29 #include <vector>
30 
31 
32 namespace impactx
33 {
44  {
46  0.1644024074311037,
47  -0.1324009958969339,
48  4.3443060026047219e-002,
49  8.5602654094946495e-002,
50  -0.2433578169042885,
51  0.5297150596779437,
52  0.7164884680963959,
53  -5.2579522442877296e-003,
54  -5.5025369142193678e-002,
55  4.6845673335028933e-002,
56  -2.3279346335638568e-002,
57  4.0800777539657775e-003,
58  4.1378326533752169e-003,
59  -2.5040533340490805e-003,
60  -4.0654981400000964e-003,
61  9.6630592067498289e-003,
62  -8.5275895985990214e-003,
63  -5.8078747006425020e-002,
64  -2.4044337836660403e-002,
65  1.0968240064697212e-002,
66  -3.4461179858301418e-003,
67  -8.1201564869443749e-004,
68  2.1438992904959380e-003,
69  -1.4997753525697276e-003,
70  1.8685171825676386e-004
71  };
72 
74  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
75  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
76  0, 0, 0
77  };
78  };
79 
80 namespace data
81 {
85  {
86  amrex::ParticleReal m_escale;
87  amrex::ParticleReal m_freq;
88  amrex::ParticleReal m_phase;
89  int m_mapsteps;
90 
91  int m_ncoef = 0;
92  amrex::ParticleReal* m_cos_data = nullptr;
93  amrex::ParticleReal* m_sin_data = nullptr;
94 
96  amrex::ParticleReal escale,
97  amrex::ParticleReal freq,
98  amrex::ParticleReal phase,
99  int mapsteps = 1
100  ) : m_escale(escale), m_freq(freq), m_phase(phase), m_mapsteps(mapsteps)
101  {}
102 
104  RFCavity_device_copyable& operator=(RFCavity_device_copyable const &) = default;
106  RFCavity_device_copyable& operator=(RFCavity_device_copyable &&) = default;
107  };
108 } // namespace data
109 
110  struct RFCavity
111  : public elements::BeamOptic<RFCavity>,
112  public elements::Thick,
114  {
115  static constexpr auto name = "RFCavity";
117 
132  amrex::ParticleReal ds,
133  amrex::ParticleReal escale,
134  amrex::ParticleReal freq,
135  amrex::ParticleReal phase,
136  std::vector<amrex::ParticleReal> cos_coef,
137  std::vector<amrex::ParticleReal> sin_coef,
138  int mapsteps = 1,
139  int nslice = 1
140  )
141  : Thick(ds, nslice),
142  RFCavity_device_copyable(escale, freq, phase, mapsteps)
143  {
144  m_ncoef = cos_coef.size();
145  if (m_ncoef != int(sin_coef.size()))
146  throw std::runtime_error("RFCavity: cos and sin coefficients must have same length!");
147 
148  m_cos_coef.resize(m_ncoef);
149  m_sin_coef.resize(m_ncoef);
151  cos_coef.begin(), cos_coef.end(),
152  m_cos_coef.begin());
154  sin_coef.begin(), sin_coef.end(),
155  m_sin_coef.begin());
157 
158  m_cos_data = m_cos_coef.data();
159  m_sin_data = m_sin_coef.data();
160  }
161 
162  // copy and move constructors
163  RFCavity(RFCavity const & other)
164  : Thick(other.m_ds, other.m_nslice),
165  RFCavity_device_copyable(other.m_escale, other.m_freq, other.m_phase, other.m_mapsteps)
166  {
167 #if !AMREX_DEVICE_COMPILE
168  // copy the data container if we copy the host element
169  m_cos_coef = other.m_cos_coef;
170  m_sin_coef = other.m_sin_coef;
172 #endif
173  m_ncoef = m_cos_coef.size();
174  m_cos_data = m_cos_coef.data();
175  m_sin_data = m_sin_coef.data();
176  }
177  RFCavity& operator=(RFCavity const& other)
178  {
179  if (this == &other)
180  return *this;
181 
182  Thick::operator=(other);
183  RFCavity_device_copyable::operator=(other);
184 #if !AMREX_DEVICE_COMPILE
185  // copy the data container if we copy the host element
186  m_cos_coef = other.m_cos_coef;
187  m_sin_coef = other.m_sin_coef;
189 #endif
190  m_ncoef = m_cos_coef.size();
191  m_cos_data = m_cos_coef.data();
192  m_sin_data = m_sin_coef.data();
193  return *this;
194  }
195  RFCavity(RFCavity && other) = default;
196  RFCavity& operator=(RFCavity && other) = default;
197 
199  ~RFCavity() = default;
200 
202  using BeamOptic::operator();
203 
214  void operator() (
215  PType& AMREX_RESTRICT p,
216  amrex::ParticleReal & AMREX_RESTRICT px,
217  amrex::ParticleReal & AMREX_RESTRICT py,
218  amrex::ParticleReal & AMREX_RESTRICT pt,
219  [[maybe_unused]] RefPart const & refpart
220  ) const
221  {
222  using namespace amrex::literals; // for _rt and _prt
223 
224  // access AoS data such as positions and cpu/id
225  amrex::ParticleReal const x = p.pos(0);
226  amrex::ParticleReal const y = p.pos(1);
227  amrex::ParticleReal const t = p.pos(2);
228 
229  // initialize output values of momenta
230  amrex::ParticleReal pxout = px;
231  amrex::ParticleReal pyout = py;
232  amrex::ParticleReal ptout = pt;
233 
234  // get the linear map
236 
237  // symplectic linear map for the RF cavity is computed using the
238  // Hamiltonian formalism as described in:
239  // https://uspas.fnal.gov/materials/09UNM/ComputationalMethods.pdf.
240  // R denotes the transfer matrix in the basis (x,px,y,py,t,pt),
241  // so that, e.g., R(3,4) = dyf/dpyi.
242 
243  // push particles using the linear map
244  p.pos(0) = R(1,1)*x + R(1,2)*px + R(1,3)*y
245  + R(1,4)*py + R(1,5)*t + R(1,6)*pt;
246  pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y
247  + R(2,4)*py + R(2,5)*t + R(2,6)*pt;
248  p.pos(1) = R(3,1)*x + R(3,2)*px + R(3,3)*y
249  + R(3,4)*py + R(3,5)*t + R(3,6)*pt;
250  pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y
251  + R(4,4)*py + R(4,5)*t + R(4,6)*pt;
252  p.pos(2) = R(5,1)*x + R(5,2)*px + R(5,3)*y
253  + R(5,4)*py + R(5,5)*t + R(5,6)*pt;
254  ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y
255  + R(6,4)*py + R(6,5)*t + R(6,6)*pt;
256 
257  // assign updated momenta
258  px = pxout;
259  py = pyout;
260  pt = ptout;
261  }
262 
268  void operator() (RefPart & AMREX_RESTRICT refpart) const
269  {
270  using namespace amrex::literals; // for _rt and _prt
271 
272  // assign input reference particle values
273  amrex::ParticleReal const x = refpart.x;
274  amrex::ParticleReal const px = refpart.px;
275  amrex::ParticleReal const y = refpart.y;
276  amrex::ParticleReal const py = refpart.py;
277  amrex::ParticleReal const z = refpart.z;
278  amrex::ParticleReal const pz = refpart.pz;
279  amrex::ParticleReal const pt = refpart.pt;
280  amrex::ParticleReal const s = refpart.s;
281  amrex::ParticleReal const sedge = refpart.sedge;
282 
283  // initialize linear map (deviation) values
284  for (int i=1; i<7; i++) {
285  for (int j=1; j<7; j++) {
286  if (i == j)
287  refpart.map(i, j) = 1.0_prt;
288  else
289  refpart.map(i, j) = 0.0_prt;
290  }
291  }
292 
293  // length of the current slice
294  amrex::ParticleReal const slice_ds = m_ds / nslice();
295 
296  // compute intial value of beta*gamma
297  amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
298 
299  // call integrator to advance (t,pt)
300  amrex::ParticleReal const zin = s - sedge;
301  amrex::ParticleReal const zout = zin + slice_ds;
302  int const nsteps = m_mapsteps;
303 
304  integrators::symp2_integrate_split3(refpart,zin,zout,nsteps,*this);
305  amrex::ParticleReal const ptf = refpart.pt;
306 
307  // advance position (x,y,z)
308  refpart.x = x + slice_ds*px/bgi;
309  refpart.y = y + slice_ds*py/bgi;
310  refpart.z = z + slice_ds*pz/bgi;
311 
312  // compute final value of beta*gamma
313  amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
314 
315  // advance momentum (px,py,pz)
316  refpart.px = px*bgf/bgi;
317  refpart.py = py*bgf/bgi;
318  refpart.pz = pz*bgf/bgi;
319 
320  // convert linear map from dynamic to static units
321  amrex::ParticleReal scale_in = 1.0_prt;
322  amrex::ParticleReal scale_fin = 1.0_prt;
323 
324  for (int i=1; i<7; i++) {
325  for (int j=1; j<7; j++) {
326  if( i % 2 == 0)
327  scale_fin = bgf;
328  else
329  scale_fin = 1.0_prt;
330  if( j % 2 == 0)
331  scale_in = bgi;
332  else
333  scale_in = 1.0_prt;
334  refpart.map(i, j) = refpart.map(i, j) * scale_in / scale_fin;
335  }
336  }
337 
338  // advance integrated path length
339  refpart.s = s + slice_ds;
340  }
341 
348  std::tuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>
350  RF_Efield (amrex::ParticleReal const zeval) const
351  {
352  using namespace amrex::literals; // for _rt and _prt
353 
354  // specify constants
356  amrex::ParticleReal const zlen = m_ds;
357  amrex::ParticleReal const zmid = zlen / 2.0_prt;
358 
359  // compute on-axis electric field (z is relative to cavity midpoint)
360  amrex::ParticleReal efield = 0.0;
361  amrex::ParticleReal efieldp = 0.0;
362  amrex::ParticleReal efieldpp = 0.0;
363  amrex::ParticleReal efieldint = 0.0;
364  amrex::ParticleReal const z = zeval - zmid;
365 
366  if (abs(z)<=zmid)
367  {
368  efield = 0.5_prt*m_cos_data[0];
369  efieldint = z*efield;
370  for (int j=1; j < m_ncoef; ++j)
371  {
372  efield = efield + m_cos_data[j]*cos(j*2*pi*z/zlen) +
373  m_sin_data[j]*sin(j*2*pi*z/zlen);
374  efieldp = efieldp-j*2*pi*m_cos_data[j]*sin(j*2*pi*z/zlen)/zlen +
375  j*2*pi*m_sin_data[j]*cos(j*2*pi*z/zlen)/zlen;
376  efieldpp = efieldpp- pow(j*2*pi*m_cos_data[j]/zlen,2) *cos(j*2*pi*z/zlen) -
377  pow(j*2*pi*m_sin_data[j]/zlen,2) *sin(j*2*pi*z/zlen);
378  efieldint = efieldint + zlen*m_cos_data[j]*sin(j*2*pi*z/zlen)/(j*2*pi) -
379  zlen*m_sin_data[j]*cos(j*2*pi*z/zlen)/(j*2*pi);
380  }
381  }
382  return std::make_tuple(efield, efieldp, efieldint);
383  }
384 
394  void map3 (amrex::ParticleReal const tau,
395  RefPart & refpart,
396  [[maybe_unused]] amrex::ParticleReal & zeval) const
397  {
398  using namespace amrex::literals; // for _rt and _prt
399 
400  // push the reference particle
401  amrex::ParticleReal const t = refpart.t;
402  amrex::ParticleReal const pt = refpart.pt;
403 
404  if (pt < -1.0_prt) {
405  refpart.t = t + tau/sqrt(1.0_prt - pow(pt, -2));
406  refpart.pt = pt;
407  }
408  else {
409  refpart.t = t;
410  refpart.pt = pt;
411  }
412 
413  // push the linear map equations
415  amrex::ParticleReal const betgam = refpart.beta_gamma();
416 
417  refpart.map(5,5) = R(5,5) + tau*R(6,5)/pow(betgam,3);
418  refpart.map(5,6) = R(5,6) + tau*R(6,6)/pow(betgam,3);
419  }
420 
430  void map2 (amrex::ParticleReal const tau,
431  RefPart & refpart,
432  amrex::ParticleReal & zeval) const
433  {
434  using namespace amrex::literals; // for _rt and _prt
435 
436  amrex::ParticleReal const t = refpart.t;
437  amrex::ParticleReal const pt = refpart.pt;
438 
439  // Define parameters and intermediate constants
442  amrex::ParticleReal const k = (2.0_prt*pi/c)*m_freq;
443  amrex::ParticleReal const phi = m_phase*(pi/180.0_prt);
444  amrex::ParticleReal const E0 = m_escale;
445 
446  // push the reference particle
447  auto [ez, ezp, ezint] = RF_Efield(zeval);
448  amrex::ignore_unused(ez, ezint);
449 
450  refpart.t = t;
451  refpart.pt = pt;
452 
453  // push the linear map equations
455  amrex::ParticleReal const s = tau/refpart.beta_gamma();
456  amrex::ParticleReal const L = E0*ezp*sin(k*t+phi)/(2.0_prt*k);
457 
458  refpart.map(1,1) = (1.0_prt-s*L)*R(1,1) + s*R(2,1);
459  refpart.map(1,2) = (1.0_prt-s*L)*R(1,2) + s*R(2,2);
460  refpart.map(2,1) = -s*pow(L,2)*R(1,1) + (1.0_prt+s*L)*R(2,1);
461  refpart.map(2,2) = -s*pow(L,2)*R(1,2) + (1.0_prt+s*L)*R(2,2);
462 
463  refpart.map(3,3) = (1.0_prt-s*L)*R(3,3) + s*R(4,3);
464  refpart.map(3,4) = (1.0_prt-s*L)*R(3,4) + s*R(4,4);
465  refpart.map(4,3) = -s*pow(L,2)*R(3,3) + (1.0_prt+s*L)*R(4,3);
466  refpart.map(4,4) = -s*pow(L,2)*R(3,4) + (1.0_prt+s*L)*R(4,4);
467  }
468 
478  void map1 (amrex::ParticleReal const tau,
479  RefPart & refpart,
480  amrex::ParticleReal & zeval) const
481  {
482  using namespace amrex::literals; // for _rt and _prt
483 
484  amrex::ParticleReal const t = refpart.t;
485  amrex::ParticleReal const pt = refpart.pt;
486  amrex::ParticleReal const z = zeval;
487 
488  // Define parameters and intermediate constants
491  amrex::ParticleReal const k = (2.0_prt*pi/c)*m_freq;
492  amrex::ParticleReal const phi = m_phase*(pi/180.0_prt);
493  amrex::ParticleReal const E0 = m_escale;
494 
495  // push the reference particle
496  auto [ez, ezp, ezint] = RF_Efield(z);
498  zeval = z + tau;
499  auto [ezf, ezpf, ezintf] = RF_Efield(zeval);
501 
502  refpart.t = t;
503  refpart.pt = pt - E0*(ezintf-ezint)*cos(k*t+phi);
504 
505  // push the linear map equations
507  amrex::ParticleReal const M = E0*(ezintf-ezint)*k*sin(k*t+phi);
508  amrex::ParticleReal const L = E0*(ezpf-ezp)*sin(k*t+phi)/(2.0_prt*k)+M/2.0_prt;
509 
510  refpart.map(2,1) = L*R(1,1) + R(2,1);
511  refpart.map(2,2) = L*R(1,2) + R(2,2);
512 
513  refpart.map(4,3) = L*R(3,3) + R(4,3);
514  refpart.map(4,4) = L*R(3,4) + R(4,4);
515 
516  refpart.map(6,5) = M*R(5,5) + R(6,5);
517  refpart.map(6,6) = M*R(5,6) + R(6,6);
518  }
519 
520  private:
521  // we cannot copy these to device with a memcpy when we copy the element class
524  };
525 
526 } // namespace impactx
527 
528 #endif // IMPACTX_RFCAVITY_H
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void symp2_integrate_split3(RefPart &refpart, amrex::ParticleReal const zin, amrex::ParticleReal const zout, int const nsteps, T_Element const &element)
Definition: Integrators.H:53
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_cos_coef
Definition: RFCavity.H:522
constexpr std::enable_if_t< std::is_floating_point< T >::value, T > pi()
data
static constexpr auto c
Definition: ImpactX.cpp:31
static constexpr HostToDevice hostToDevice
Definition: beamoptic.H:134
c
amrex::ParticleReal pt
energy deviation, normalized by rest energy
Definition: ReferenceParticle.H:39
RFCavity & operator=(RFCavity const &other)
Definition: RFCavity.H:177
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMREX_GPU_HOST RFCavity(amrex::ParticleReal ds, amrex::ParticleReal escale, amrex::ParticleReal freq, amrex::ParticleReal phase, std::vector< amrex::ParticleReal > cos_coef, std::vector< amrex::ParticleReal > sin_coef, int mapsteps=1, int nslice=1)
Definition: RFCavity.H:131
#define AMREX_FORCE_INLINE
amrex::ParticleReal m_freq
scaling factor for RF electric field
Definition: RFCavity.H:87
nslice
Definition: __init__.py:32
amrex::Vector< amrex::ParticleReal > default_sin_coef
Definition: RFCavity.H:73
amrex::ParticleReal m_escale
Definition: RFCavity.H:86
#define AMREX_GPU_HOST_DEVICE
amrex::ParticleReal t
clock time * c in meters
Definition: ReferenceParticle.H:35
int nsteps
i
Definition: thick.H:23
amrex::Array2D< amrex::ParticleReal, 1, 6, 1, 6 > map
linearized map
Definition: ReferenceParticle.H:44
int ez
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Definition: ReferenceParticle.H:29
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, const T &a_y) noexcept
RFCavity(RFCavity const &other)
Definition: RFCavity.H:163
string name
amrex::Vector< amrex::ParticleReal > default_cos_coef
Definition: RFCavity.H:45
s
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal beta_gamma() const
Definition: ReferenceParticle.H:79
std::tuple< amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal > AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RF_Efield(amrex::ParticleReal const zeval) const
Definition: RFCavity.H:350
static constexpr amrex::Real pi
Definition: RFCavity.H:43
Definition: RFCavity.H:110
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map2(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: RFCavity.H:430
amrex::ParticleReal m_phase
RF frequency in Hz.
Definition: RFCavity.H:88
RFCavity_device_copyable(amrex::ParticleReal escale, amrex::ParticleReal freq, amrex::ParticleReal phase, int mapsteps=1)
non-owning pointer to device sine coefficients
Definition: RFCavity.H:95
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
#define AMREX_GPU_HOST
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map3(amrex::ParticleReal const tau, RefPart &refpart, [[maybe_unused]] amrex::ParticleReal &zeval) const
Definition: RFCavity.H:394
void synchronize() noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map1(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: RFCavity.H:478
int m_mapsteps
RF driven phase in deg.
Definition: RFCavity.H:89
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_sin_coef
cosine coefficients in Fourier expansion of on-axis electric field Ez
Definition: RFCavity.H:523