ImpactX
SoftSol.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_SOFTSOL_H
11 #define IMPACTX_SOFTSOL_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 {
47  {
49  0.350807812299706,
50  0.323554693720069,
51  0.260320578919415,
52  0.182848575294969,
53  0.106921016050403,
54  4.409581845710694E-002,
55  -9.416427163897508E-006,
56  -2.459452716865687E-002,
57  -3.272762575737291E-002,
58  -2.936414401076162E-002,
59  -1.995780078926890E-002,
60  -9.102893342953847E-003,
61  -2.456410658713271E-006,
62  5.788233017324325E-003,
63  8.040408292420691E-003,
64  7.480064552867431E-003,
65  5.230254569468851E-003,
66  2.447614547094685E-003,
67  -1.095525090532255E-006,
68  -1.614586867387170E-003,
69  -2.281365457438345E-003,
70  -2.148709081338292E-003,
71  -1.522541739363011E-003,
72  -7.185505862719508E-004,
73  -6.171194824600157E-007,
74  4.842109305036943E-004,
75  6.874508102002901E-004,
76  6.535550288205728E-004,
77  4.648795813759210E-004,
78  2.216564722797528E-004,
79  -4.100982995210341E-007,
80  -1.499332112463395E-004,
81  -2.151538438342482E-004,
82  -2.044590946652016E-004,
83  -1.468242784844341E-004
84  };
85 
87  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
88  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
89  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
90  0, 0, 0, 0, 0
91  };
92  };
93 
94 namespace data
95 {
99  {
100  amrex::ParticleReal m_bscale;
102 
103  int m_ncoef = 0;
104  amrex::ParticleReal* m_cos_data = nullptr;
105  amrex::ParticleReal* m_sin_data = nullptr;
106 
108  amrex::ParticleReal bscale,
109  int mapsteps = 1
110  ) : m_bscale(bscale), m_mapsteps(mapsteps)
111  {}
112 
114  SoftSolenoid_device_copyable& operator=(SoftSolenoid_device_copyable const &) = default;
117 
118  };
119 } // namespace data
120 
122  : public elements::BeamOptic<SoftSolenoid>,
123  public elements::Thick,
125  {
126  static constexpr auto name = "SoftSolenoid";
128 
141  amrex::ParticleReal ds,
142  amrex::ParticleReal bscale,
143  std::vector<amrex::ParticleReal> cos_coef,
144  std::vector<amrex::ParticleReal> sin_coef,
145  int mapsteps = 1,
146  int nslice = 1
147  )
148  : Thick(ds, nslice),
149  SoftSolenoid_device_copyable(bscale, mapsteps)
150  {
151  m_ncoef = cos_coef.size();
152  if (m_ncoef != int(sin_coef.size()))
153  throw std::runtime_error("SoftSolenoid: cos and sin coefficients must have same length!");
154 
155  m_cos_coef.resize(m_ncoef);
156  m_sin_coef.resize(m_ncoef);
158  cos_coef.begin(), cos_coef.end(),
159  m_cos_coef.begin());
161  sin_coef.begin(), sin_coef.end(),
162  m_sin_coef.begin());
164 
165  // low-level objects we can use on device
166  m_cos_data = m_cos_coef.data();
167  m_sin_data = m_sin_coef.data();
168  }
169 
170  // copy and move constructors
171  SoftSolenoid(SoftSolenoid const & other)
172  : Thick(other.m_ds, other.m_nslice),
173  SoftSolenoid_device_copyable(other.m_bscale, other.m_mapsteps)
174  {
175 #if !AMREX_DEVICE_COMPILE
176  // copy the data container if we copy the host element
177  m_cos_coef = other.m_cos_coef;
178  m_sin_coef = other.m_sin_coef;
180 #endif
181  m_ncoef = m_cos_coef.size();
182  m_cos_data = m_cos_coef.data();
183  m_sin_data = m_sin_coef.data();
184  }
186  {
187  if (this == &other)
188  return *this;
189 
190  Thick::operator=(other);
191  SoftSolenoid_device_copyable::operator=(other);
192 #if !AMREX_DEVICE_COMPILE
193  // copy the data container if we copy the host element
194  m_cos_coef = other.m_cos_coef;
195  m_sin_coef = other.m_sin_coef;
197 #endif
198  m_ncoef = m_cos_coef.size();
199  m_cos_data = m_cos_coef.data();
200  m_sin_data = m_sin_coef.data();
201  return *this;
202  }
203  SoftSolenoid(SoftSolenoid && other) = default;
204  SoftSolenoid& operator=(SoftSolenoid && other) = default;
205 
207  ~SoftSolenoid() = default;
208 
210  using BeamOptic::operator();
211 
222  void operator() (
223  PType& AMREX_RESTRICT p,
224  amrex::ParticleReal & AMREX_RESTRICT px,
225  amrex::ParticleReal & AMREX_RESTRICT py,
226  amrex::ParticleReal & AMREX_RESTRICT pt,
227  [[maybe_unused]] RefPart const & refpart
228  ) const
229  {
230  using namespace amrex::literals; // for _rt and _prt
231 
232  // access AoS data such as positions and cpu/id
233  amrex::ParticleReal const x = p.pos(0);
234  amrex::ParticleReal const y = p.pos(1);
235  amrex::ParticleReal const t = p.pos(2);
236 
237  // initialize output values of momenta
238  amrex::ParticleReal pxout = px;
239  amrex::ParticleReal pyout = py;
240  amrex::ParticleReal ptout = pt;
241 
242  // get the linear map
244 
245  // symplectic linear map for a solenoid is computed using the
246  // Hamiltonian formalism as described in:
247  // https://uspas.fnal.gov/materials/09UNM/ComputationalMethods.pdf.
248  // R denotes the transfer matrix in the basis (x,px,y,py,t,pt),
249  // so that, e.g., R(3,4) = dyf/dpyi.
250 
251  // push particles using the linear map
252  p.pos(0) = R(1,1)*x + R(1,2)*px + R(1,3)*y
253  + R(1,4)*py + R(1,5)*t + R(1,6)*pt;
254  pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y
255  + R(2,4)*py + R(2,5)*t + R(2,6)*pt;
256  p.pos(1) = R(3,1)*x + R(3,2)*px + R(3,3)*y
257  + R(3,4)*py + R(3,5)*t + R(3,6)*pt;
258  pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y
259  + R(4,4)*py + R(4,5)*t + R(4,6)*pt;
260  p.pos(2) = R(5,1)*x + R(5,2)*px + R(5,3)*y
261  + R(5,4)*py + R(5,5)*t + R(5,6)*pt;
262  ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y
263  + R(6,4)*py + R(6,5)*t + R(6,6)*pt;
264 
265  // assign updated momenta
266  px = pxout;
267  py = pyout;
268  pt = ptout;
269  }
270 
276  void operator() (RefPart & AMREX_RESTRICT refpart) const
277  {
278  using namespace amrex::literals; // for _rt and _prt
279 
280  // assign input reference particle values
281  amrex::ParticleReal const x = refpart.x;
282  amrex::ParticleReal const px = refpart.px;
283  amrex::ParticleReal const y = refpart.y;
284  amrex::ParticleReal const py = refpart.py;
285  amrex::ParticleReal const z = refpart.z;
286  amrex::ParticleReal const pz = refpart.pz;
287  amrex::ParticleReal const pt = refpart.pt;
288  amrex::ParticleReal const s = refpart.s;
289  amrex::ParticleReal const sedge = refpart.sedge;
290 
291  // initialize linear map (deviation) values
292  for (int i=1; i<7; i++) {
293  for (int j=1; j<7; j++) {
294  auto const default_value = (i == j) ? 1.0_prt : 0.0_prt;
295  refpart.map(i, j) = default_value;
296  }
297  }
298 
299  // length of the current slice
300  amrex::ParticleReal const slice_ds = m_ds / nslice();
301 
302  // compute intial value of beta*gamma
303  amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
304 
305  // call integrator to advance (t,pt)
306  amrex::ParticleReal const zin = s - sedge;
307  amrex::ParticleReal const zout = zin + slice_ds;
308  int const nsteps = m_mapsteps;
309 
310  integrators::symp2_integrate_split3(refpart,zin,zout,nsteps,*this);
311  amrex::ParticleReal const ptf = refpart.pt;
312 
313  /* print computed linear map:
314  for(int i=1; i<7; ++i){
315  for(int j=1; j<7; ++j){
316  amrex::PrintToFile("SolMap.txt") << i << " " <<
317  j << " " << refpart.map(i,j) << "\n";
318  }
319  }
320  */
321 
322  // advance position (x,y,z)
323  refpart.x = x + slice_ds*px/bgi;
324  refpart.y = y + slice_ds*py/bgi;
325  refpart.z = z + slice_ds*pz/bgi;
326 
327  // compute final value of beta*gamma
328  amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
329 
330  // advance momentum (px,py,pz)
331  refpart.px = px*bgf/bgi;
332  refpart.py = py*bgf/bgi;
333  refpart.pz = pz*bgf/bgi;
334 
335  // advance integrated path length
336  refpart.s = s + slice_ds;
337  }
338 
345  std::tuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>
347  Sol_Bfield (amrex::ParticleReal const zeval) const
348  {
349  using namespace amrex::literals; // for _rt and _prt
350 
351  // specify constants
353  amrex::ParticleReal const zlen = m_ds;
354  amrex::ParticleReal const zmid = zlen / 2.0_prt;
355 
356  // compute on-axis magnetic field (z is relative to solenoid midpoint)
357  amrex::ParticleReal bfield = 0.0;
358  amrex::ParticleReal bfieldp = 0.0;
359  amrex::ParticleReal bfieldint = 0.0;
360  amrex::ParticleReal const z = zeval - zmid;
361 
362  if (abs(z)<=zmid)
363  {
364  bfield = 0.5_prt*m_cos_data[0];
365  bfieldint = z*bfield;
366  for (int j=1; j < m_ncoef; ++j)
367  {
368  bfield = bfield + m_cos_data[j]*cos(j*2*pi*z/zlen) +
369  m_sin_data[j]*sin(j*2*pi*z/zlen);
370  bfieldp = bfieldp-j*2*pi*m_cos_data[j]*sin(j*2*pi*z/zlen)/zlen +
371  j*2*pi*m_sin_data[j]*cos(j*2*pi*z/zlen)/zlen;
372  bfieldint = bfieldint + zlen*m_cos_data[j]*sin(j*2*pi*z/zlen)/(j*2*pi) -
373  zlen*m_sin_data[j]*cos(j*2*pi*z/zlen)/(j*2*pi);
374  }
375  }
376  return std::make_tuple(bfield, bfieldp, bfieldint);
377  }
378 
388  void map1 (amrex::ParticleReal const tau,
389  RefPart & refpart,
390  [[maybe_unused]] amrex::ParticleReal & zeval) const
391  {
392  using namespace amrex::literals; // for _rt and _prt
393 
394  // push the reference particle
395  amrex::ParticleReal const t = refpart.t;
396  amrex::ParticleReal const pt = refpart.pt;
397  amrex::ParticleReal const z = zeval;
398 
399  if (pt < -1.0_prt) {
400  refpart.t = t + tau/sqrt(1.0_prt - pow(pt, -2));
401  refpart.pt = pt;
402  }
403  else {
404  refpart.t = t;
405  refpart.pt = pt;
406  }
407 
408  zeval = z + tau;
409 
410  // push the linear map equations
412  amrex::ParticleReal const betgam = refpart.beta_gamma();
413 
414  refpart.map(1,1) = R(1,1) + tau*R(2,1);
415  refpart.map(1,2) = R(1,2) + tau*R(2,2);
416  refpart.map(1,3) = R(1,3) + tau*R(2,3);
417  refpart.map(1,4) = R(1,4) + tau*R(2,4);
418 
419  refpart.map(3,1) = R(3,1) + tau*R(4,1);
420  refpart.map(3,2) = R(3,2) + tau*R(4,2);
421  refpart.map(3,3) = R(3,3) + tau*R(4,3);
422  refpart.map(3,4) = R(3,4) + tau*R(4,4);
423 
424  refpart.map(5,5) = R(5,5) + tau*R(6,5)/pow(betgam,2);
425  refpart.map(5,6) = R(5,6) + tau*R(6,6)/pow(betgam,2);
426 
427  }
428 
438  void map2 (amrex::ParticleReal const tau,
439  RefPart & refpart,
440  amrex::ParticleReal & zeval) const
441  {
442  using namespace amrex::literals; // for _rt and _prt
443 
444  amrex::ParticleReal const t = refpart.t;
445  amrex::ParticleReal const pt = refpart.pt;
446 
447  // Define parameters and intermediate constants
448  amrex::ParticleReal const B0 = m_bscale;
449 
450  // push the reference particle
451  auto [bz, bzp, bzint] = Sol_Bfield(zeval);
452  amrex::ignore_unused(bzp, bzint);
453 
454  refpart.t = t;
455  refpart.pt = pt;
456 
457  // push the linear map equations
459  amrex::ParticleReal const alpha = B0*bz/2.0_prt;
460  amrex::ParticleReal const alpha2 = pow(alpha,2);
461 
462  refpart.map(2,1) = R(2,1) - tau*alpha2*R(1,1);
463  refpart.map(2,2) = R(2,2) - tau*alpha2*R(1,2);
464  refpart.map(2,3) = R(2,3) - tau*alpha2*R(1,3);
465  refpart.map(2,4) = R(2,4) - tau*alpha2*R(1,4);
466 
467  refpart.map(4,1) = R(4,1) - tau*alpha2*R(3,1);
468  refpart.map(4,2) = R(4,2) - tau*alpha2*R(3,2);
469  refpart.map(4,3) = R(4,3) - tau*alpha2*R(3,3);
470  refpart.map(4,4) = R(4,4) - tau*alpha2*R(3,4);
471 
472  }
473 
483  void map3 (amrex::ParticleReal const tau,
484  RefPart & refpart,
485  amrex::ParticleReal & zeval) const
486  {
487  using namespace amrex::literals; // for _rt and _prt
488 
489  amrex::ParticleReal const t = refpart.t;
490  amrex::ParticleReal const pt = refpart.pt;
491  amrex::ParticleReal const z = zeval;
492 
493  // Define parameters and intermediate constants
494  amrex::ParticleReal const B0 = m_bscale;
495 
496  // push the reference particle
497  auto [bz, bzp, bzint] = Sol_Bfield(z);
498  amrex::ignore_unused(bzp, bzint);
499 
500  refpart.t = t;
501  refpart.pt = pt;
502 
503  // push the linear map equations
505  amrex::ParticleReal const theta = tau*B0*bz/2.0_prt;
506  amrex::ParticleReal const cs = cos(theta);
507  amrex::ParticleReal const sn = sin(theta);
508 
509  refpart.map(1,1) = R(1,1)*cs + R(3,1)*sn;
510  refpart.map(1,2) = R(1,2)*cs + R(3,2)*sn;
511  refpart.map(1,3) = R(1,3)*cs + R(3,3)*sn;
512  refpart.map(1,4) = R(1,4)*cs + R(3,4)*sn;
513 
514  refpart.map(2,1) = R(2,1)*cs + R(4,1)*sn;
515  refpart.map(2,2) = R(2,2)*cs + R(4,2)*sn;
516  refpart.map(2,3) = R(2,3)*cs + R(4,3)*sn;
517  refpart.map(2,4) = R(2,4)*cs + R(4,4)*sn;
518 
519  refpart.map(3,1) = R(3,1)*cs - R(1,1)*sn;
520  refpart.map(3,2) = R(3,2)*cs - R(1,2)*sn;
521  refpart.map(3,3) = R(3,3)*cs - R(1,3)*sn;
522  refpart.map(3,4) = R(3,4)*cs - R(1,4)*sn;
523 
524  refpart.map(4,1) = R(4,1)*cs - R(2,1)*sn;
525  refpart.map(4,2) = R(4,2)*cs - R(2,2)*sn;
526  refpart.map(4,3) = R(4,3)*cs - R(2,3)*sn;
527  refpart.map(4,4) = R(4,4)*cs - R(2,4)*sn;
528 
529  }
530 
531  private:
532  // we cannot copy these to device with a memcpy when we copy the element class
535  };
536 
537 } // namespace impactx
538 
539 #endif // IMPACTX_SOFTSOL_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()
data
std::tuple< amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal > AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Sol_Bfield(amrex::ParticleReal const zeval) const
Definition: SoftSol.H:347
int m_mapsteps
scaling factor for solenoid Bz field
Definition: SoftSol.H:101
Definition: ImpactX.cpp:31
SoftSolenoid_device_copyable(amrex::ParticleReal bscale, int mapsteps=1)
non-owning pointer to device sine coefficients
Definition: SoftSol.H:107
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_cos_coef
Definition: SoftSol.H:533
SoftSolenoid(SoftSolenoid const &other)
Definition: SoftSol.H:171
static constexpr HostToDevice hostToDevice
amrex::Vector< amrex::ParticleReal > default_sin_coef
Definition: SoftSol.H:86
Definition: beamoptic.H:134
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
nslice
Definition: __init__.py:32
#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_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map1(amrex::ParticleReal const tau, RefPart &refpart, [[maybe_unused]] amrex::ParticleReal &zeval) const
Definition: SoftSol.H:388
amrex::Array2D< amrex::ParticleReal, 1, 6, 1, 6 > map
linearized map
Definition: ReferenceParticle.H:44
amrex::ParticleReal m_bscale
Definition: SoftSol.H:100
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map2(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: SoftSol.H:438
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_sin_coef
cosine coefficients in Fourier expansion of on-axis magnetic field Bz
Definition: SoftSol.H:534
Definition: ReferenceParticle.H:29
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, const T &a_y) noexcept
SoftSolenoid & operator=(SoftSolenoid const &other)
Definition: SoftSol.H:185
Definition: SoftSol.H:46
string name
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map3(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: SoftSol.H:483
s
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal beta_gamma() const
Definition: ReferenceParticle.H:79
static constexpr amrex::Real pi
AMREX_GPU_HOST SoftSolenoid(amrex::ParticleReal ds, amrex::ParticleReal bscale, std::vector< amrex::ParticleReal > cos_coef, std::vector< amrex::ParticleReal > sin_coef, int mapsteps=1, int nslice=1)
Definition: SoftSol.H:140
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
#define AMREX_GPU_HOST
void synchronize() noexcept
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
amrex::Vector< amrex::ParticleReal > default_cos_coef
Definition: SoftSol.H:48
Definition: SoftSol.H:121