ImpactX
RFCavity.H
Go to the documentation of this file.
1 /* Copyright 2022 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 
16 #include <ablastr/constant.H>
17 
18 #include <AMReX.H>
19 #include <AMReX_Array.H>
20 #include <AMReX_Extension.H>
21 #include <AMReX_REAL.H>
22 
23 #include <array>
24 #include <cmath>
25 #include <tuple>
26 
27 
28 namespace impactx
29 {
30  struct RFCavity
31  {
32  static constexpr auto name = "RFCavity";
34 
46  amrex::ParticleReal const ds,
47  amrex::ParticleReal const escale,
48  amrex::ParticleReal const freq,
49  amrex::ParticleReal const phase,
50  int const mapsteps,
51  int const nslice
52  )
53  : m_ds(ds), m_escale(escale), m_freq(freq), m_phase(phase),
54  m_mapsteps(mapsteps), m_nslice(nslice)
55  {
56  }
57 
69  PType& AMREX_RESTRICT p,
70  amrex::ParticleReal & AMREX_RESTRICT px,
71  amrex::ParticleReal & AMREX_RESTRICT py,
72  amrex::ParticleReal & AMREX_RESTRICT pt,
73  [[maybe_unused]] RefPart const & refpart
74  ) const
75  {
76  using namespace amrex::literals; // for _rt and _prt
77 
78  // access AoS data such as positions and cpu/id
79  amrex::ParticleReal const x = p.pos(0);
80  amrex::ParticleReal const y = p.pos(1);
81  amrex::ParticleReal const t = p.pos(2);
82 
83  // initialize output values of momenta
84  amrex::ParticleReal pxout = px;
85  amrex::ParticleReal pyout = py;
86  amrex::ParticleReal ptout = pt;
87 
88  // get the linear map
90 
91  // symplectic linear map for the RF cavity is computed using the
92  // Hamiltonian formalism as described in:
93  // https://uspas.fnal.gov/materials/09UNM/ComputationalMethods.pdf.
94  // R denotes the transfer matrix in the basis (x,px,y,py,t,pt),
95  // so that, e.g., R(3,4) = dyf/dpyi.
96 
97  // push particles using the linear map
98  p.pos(0) = R(1,1)*x + R(1,2)*px + R(1,3)*y
99  + R(1,4)*py + R(1,5)*t + R(1,6)*pt;
100  pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y
101  + R(2,4)*py + R(2,5)*t + R(2,6)*pt;
102  p.pos(1) = R(3,1)*x + R(3,2)*px + R(3,3)*y
103  + R(3,4)*py + R(3,5)*t + R(3,6)*pt;
104  pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y
105  + R(4,4)*py + R(4,5)*t + R(4,6)*pt;
106  p.pos(2) = R(5,1)*x + R(5,2)*px + R(5,3)*y
107  + R(5,4)*py + R(5,5)*t + R(5,6)*pt;
108  ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y
109  + R(6,4)*py + R(6,5)*t + R(6,6)*pt;
110 
111  // assign updated momenta
112  px = pxout;
113  py = pyout;
114  pt = ptout;
115  }
116 
122  void operator() (RefPart & AMREX_RESTRICT refpart) const
123  {
124  using namespace amrex::literals; // for _rt and _prt
125 
126  // assign input reference particle values
127  amrex::ParticleReal const x = refpart.x;
128  amrex::ParticleReal const px = refpart.px;
129  amrex::ParticleReal const y = refpart.y;
130  amrex::ParticleReal const py = refpart.py;
131  amrex::ParticleReal const z = refpart.z;
132  amrex::ParticleReal const pz = refpart.pz;
133  amrex::ParticleReal const pt = refpart.pt;
134  amrex::ParticleReal const s = refpart.s;
135  amrex::ParticleReal const sedge = refpart.sedge;
136 
137  // initialize linear map (deviation) values
138  for (int i=1; i<7; i++) {
139  for (int j=1; j<7; j++) {
140  if (i == j)
141  refpart.map(i, j) = 1.0_prt;
142  else
143  refpart.map(i, j) = 0.0_prt;
144  }
145  }
146 
147  // length of the current slice
148  amrex::ParticleReal const slice_ds = m_ds / nslice();
149 
150  // compute intial value of beta*gamma
151  amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
152 
153  // call integrator to advance (t,pt)
154  amrex::ParticleReal const zin = s - sedge;
155  amrex::ParticleReal const zout = zin + slice_ds;
156  int const nsteps = m_mapsteps;
157 
158  integrators::symp2_integrate_split3(refpart,zin,zout,nsteps,*this);
159  amrex::ParticleReal const ptf = refpart.pt;
160 
161  // advance position (x,y,z)
162  refpart.x = x + slice_ds*px/bgi;
163  refpart.y = y + slice_ds*py/bgi;
164  refpart.z = z + slice_ds*pz/bgi;
165 
166  // compute final value of beta*gamma
167  amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
168 
169  // advance momentum (px,py,pz)
170  refpart.px = px*bgf/bgi;
171  refpart.py = py*bgf/bgi;
172  refpart.pz = pz*bgf/bgi;
173 
174  // convert linear map from dynamic to static units
175  amrex::ParticleReal scale_in = 1.0_prt;
176  amrex::ParticleReal scale_fin = 1.0_prt;
177 
178  for (int i=1; i<7; i++) {
179  for (int j=1; j<7; j++) {
180  if( i % 2 == 0)
181  scale_fin = bgf;
182  else
183  scale_fin = 1.0_prt;
184  if( j % 2 == 0)
185  scale_in = bgi;
186  else
187  scale_in = 1.0_prt;
188  refpart.map(i, j) = refpart.map(i, j) * scale_in / scale_fin;
189  }
190  }
191 
192  // advance integrated path length
193  refpart.s = s + slice_ds;
194  }
195 
202  std::tuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>
204  RF_Efield (amrex::ParticleReal const zeval) const
205  {
206  using namespace amrex::literals; // for _rt and _prt
207 
208  // specify constants
210  constexpr amrex::ParticleReal zlen = 1.31879807_prt;
211  constexpr amrex::ParticleReal zmid = zlen / 2.0_prt;
212  constexpr int ncoef = 25;
213 
214  // specify Fourier coefficients:
215  // Fourier coefficients for the on-axis longitudinal
216  // electric field Ez of the 9-cell TESLA superconducting
217  // cavity (default cavity model):
218  // https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.3.092001.
219  // May be superseded by user-provided input values (to do).
220  constexpr std::array<amrex::ParticleReal, ncoef> cos_coef = {
221  0.1644024074311037,
222  -0.1324009958969339,
223  4.3443060026047219e-002,
224  8.5602654094946495e-002,
225  -0.2433578169042885,
226  0.5297150596779437,
227  0.7164884680963959,
228  -5.2579522442877296e-003,
229  -5.5025369142193678e-002,
230  4.6845673335028933e-002,
231  -2.3279346335638568e-002,
232  4.0800777539657775e-003,
233  4.1378326533752169e-003,
234  -2.5040533340490805e-003,
235  -4.0654981400000964e-003,
236  9.6630592067498289e-003,
237  -8.5275895985990214e-003,
238  -5.8078747006425020e-002,
239  -2.4044337836660403e-002,
240  1.0968240064697212e-002,
241  -3.4461179858301418e-003,
242  -8.1201564869443749e-004,
243  2.1438992904959380e-003,
244  -1.4997753525697276e-003,
245  1.8685171825676386e-004
246  };
247  constexpr std::array<amrex::ParticleReal, ncoef> sin_coef = {}; // all zero
248 
249  // compute on-axis electric field (z is relative to cavity midpoint)
250  amrex::ParticleReal efield = 0.0;
251  amrex::ParticleReal efieldp = 0.0;
252  amrex::ParticleReal efieldpp = 0.0;
253  amrex::ParticleReal efieldint = 0.0;
254  amrex::ParticleReal const z = zeval - zmid;
255 
256  if (abs(z)<=zmid)
257  {
258  efield = 0.5_prt*cos_coef[0];
259  efieldint = z*efield;
260  for (int j=1; j < ncoef; ++j)
261  {
262  efield = efield + cos_coef[j]*cos(j*2*pi*z/zlen) +
263  sin_coef[j]*sin(j*2*pi*z/zlen);
264  efieldp = efieldp-j*2*pi*cos_coef[j]*sin(j*2*pi*z/zlen)/zlen +
265  j*2*pi*sin_coef[j]*cos(j*2*pi*z/zlen)/zlen;
266  efieldpp = efieldpp- pow(j*2*pi*cos_coef[j]/zlen,2) *cos(j*2*pi*z/zlen) -
267  pow(j*2*pi*sin_coef[j]/zlen,2) *sin(j*2*pi*z/zlen);
268  efieldint = efieldint + zlen*cos_coef[j]*sin(j*2*pi*z/zlen)/(j*2*pi) -
269  zlen*sin_coef[j]*cos(j*2*pi*z/zlen)/(j*2*pi);
270  }
271  }
272  return std::make_tuple(efield, efieldp, efieldint);
273  }
274 
284  void map3 (amrex::ParticleReal const tau,
285  RefPart & refpart,
286  [[maybe_unused]] amrex::ParticleReal & zeval) const
287  {
288  using namespace amrex::literals; // for _rt and _prt
289 
290  // push the reference particle
291  amrex::ParticleReal const t = refpart.t;
292  amrex::ParticleReal const pt = refpart.pt;
293 
294  if (pt < -1.0_prt) {
295  refpart.t = t + tau/sqrt(1.0_prt - pow(pt, -2));
296  refpart.pt = pt;
297  }
298  else {
299  refpart.t = t;
300  refpart.pt = pt;
301  }
302 
303  // push the linear map equations
305  amrex::ParticleReal const betgam = refpart.beta_gamma();
306 
307  refpart.map(5,5) = R(5,5) + tau*R(6,5)/pow(betgam,3);
308  refpart.map(5,6) = R(5,6) + tau*R(6,6)/pow(betgam,3);
309  }
310 
320  void map2 (amrex::ParticleReal const tau,
321  RefPart & refpart,
322  amrex::ParticleReal & zeval) const
323  {
324  using namespace amrex::literals; // for _rt and _prt
325 
326  amrex::ParticleReal const t = refpart.t;
327  amrex::ParticleReal const pt = refpart.pt;
328 
329  // Define parameters and intermediate constants
332  amrex::ParticleReal const k = (2.0_prt*pi/c)*m_freq;
333  amrex::ParticleReal const phi = m_phase*(pi/180.0_prt);
334  amrex::ParticleReal const E0 = m_escale;
335 
336  // push the reference particle
337  auto [ez, ezp, ezint] = RF_Efield(zeval);
338  amrex::ignore_unused(ez, ezint);
339 
340  refpart.t = t;
341  refpart.pt = pt;
342 
343  // push the linear map equations
345  amrex::ParticleReal const s = tau/refpart.beta_gamma();
346  amrex::ParticleReal const L = E0*ezp*sin(k*t+phi)/(2.0_prt*k);
347 
348  refpart.map(1,1) = (1.0_prt-s*L)*R(1,1) + s*R(2,1);
349  refpart.map(1,2) = (1.0_prt-s*L)*R(1,2) + s*R(2,2);
350  refpart.map(2,1) = -s*pow(L,2)*R(1,1) + (1.0_prt+s*L)*R(2,1);
351  refpart.map(2,2) = -s*pow(L,2)*R(1,2) + (1.0_prt+s*L)*R(2,2);
352 
353  refpart.map(3,3) = (1.0_prt-s*L)*R(3,3) + s*R(4,3);
354  refpart.map(3,4) = (1.0_prt-s*L)*R(3,4) + s*R(4,4);
355  refpart.map(4,3) = -s*pow(L,2)*R(3,3) + (1.0_prt+s*L)*R(4,3);
356  refpart.map(4,4) = -s*pow(L,2)*R(3,4) + (1.0_prt+s*L)*R(4,4);
357  }
358 
368  void map1 (amrex::ParticleReal const tau,
369  RefPart & refpart,
370  amrex::ParticleReal & zeval) const
371  {
372  using namespace amrex::literals; // for _rt and _prt
373 
374  amrex::ParticleReal const t = refpart.t;
375  amrex::ParticleReal const pt = refpart.pt;
376  amrex::ParticleReal const z = zeval;
377 
378  // Define parameters and intermediate constants
381  amrex::ParticleReal const k = (2.0_prt*pi/c)*m_freq;
382  amrex::ParticleReal const phi = m_phase*(pi/180.0_prt);
383  amrex::ParticleReal const E0 = m_escale;
384 
385  // push the reference particle
386  auto [ez, ezp, ezint] = RF_Efield(z);
388  zeval = z + tau;
389  auto [ezf, ezpf, ezintf] = RF_Efield(zeval);
391 
392  refpart.t = t;
393  refpart.pt = pt - E0*(ezintf-ezint)*cos(k*t+phi);
394 
395  // push the linear map equations
397  amrex::ParticleReal const M = E0*(ezintf-ezint)*k*sin(k*t+phi);
398  amrex::ParticleReal const L = E0*(ezpf-ezp)*sin(k*t+phi)/(2.0_prt*k)+M/2.0_prt;
399 
400  refpart.map(2,1) = L*R(1,1) + R(2,1);
401  refpart.map(2,2) = L*R(1,2) + R(2,2);
402 
403  refpart.map(4,3) = L*R(3,3) + R(4,3);
404  refpart.map(4,4) = L*R(3,4) + R(4,4);
405 
406  refpart.map(6,5) = M*R(5,5) + R(6,5);
407  refpart.map(6,6) = M*R(5,6) + R(6,6);
408  }
409 
415  int nslice () const
416  {
417  return m_nslice;
418  }
419 
425  amrex::ParticleReal ds () const
426  {
427  return m_ds;
428  }
429 
430  private:
431  amrex::ParticleReal m_ds;
432  amrex::ParticleReal m_escale;
433  amrex::ParticleReal m_freq;
434  amrex::ParticleReal m_phase;
436  int m_nslice;
437  };
438 
439 } // namespace impactx
440 
441 #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
constexpr std::enable_if_t< std::is_floating_point< T >::value, T > pi()
int m_mapsteps
RF driven phase in deg.
Definition: RFCavity.H:435
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()(PType &AMREX_RESTRICT p, amrex::ParticleReal &AMREX_RESTRICT px, amrex::ParticleReal &AMREX_RESTRICT py, amrex::ParticleReal &AMREX_RESTRICT pt, [[maybe_unused]] RefPart const &refpart) const
Definition: RFCavity.H:68
static constexpr auto c
Definition: ImpactX.cpp:31
RFCavity(amrex::ParticleReal const ds, amrex::ParticleReal const escale, amrex::ParticleReal const freq, amrex::ParticleReal const phase, int const mapsteps, int const nslice)
Definition: RFCavity.H:45
amrex::ParticleReal m_phase
RF frequency in Hz.
Definition: RFCavity.H:434
amrex::ParticleReal m_ds
Definition: RFCavity.H:431
c
amrex::ParticleReal pt
energy deviation, normalized by rest energy
Definition: ReferenceParticle.H:39
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
amrex::ParticleReal t
clock time * c in meters
Definition: ReferenceParticle.H:35
int nsteps
i
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int nslice() const
Definition: RFCavity.H:415
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
int m_nslice
number of map integration steps per slice
Definition: RFCavity.H:436
amrex::ParticleReal m_escale
segment length in m
Definition: RFCavity.H:432
static constexpr auto name
Definition: RFCavity.H:32
s
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal ds() const
Definition: RFCavity.H:425
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:204
static constexpr amrex::Real pi
Definition: RFCavity.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map2(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: RFCavity.H:320
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map3(amrex::ParticleReal const tau, RefPart &refpart, [[maybe_unused]] amrex::ParticleReal &zeval) const
Definition: RFCavity.H:284
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map1(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: RFCavity.H:368
amrex::ParticleReal m_freq
scaling factor for RF electric field
Definition: RFCavity.H:433