10 #ifndef IMPACTX_SOFTSOL_H 11 #define IMPACTX_SOFTSOL_H 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
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,
104 amrex::ParticleReal* m_cos_data =
nullptr;
105 amrex::ParticleReal* m_sin_data =
nullptr;
108 amrex::ParticleReal bscale,
110 ) : m_bscale(bscale), m_mapsteps(mapsteps)
126 static constexpr
auto name =
"SoftSolenoid";
141 amrex::ParticleReal ds,
142 amrex::ParticleReal bscale,
143 std::vector<amrex::ParticleReal> cos_coef,
144 std::vector<amrex::ParticleReal> sin_coef,
149 SoftSolenoid_device_copyable(bscale, mapsteps)
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!");
155 m_cos_coef.resize(m_ncoef);
156 m_sin_coef.resize(m_ncoef);
158 cos_coef.begin(), cos_coef.end(),
161 sin_coef.begin(), sin_coef.end(),
166 m_cos_data = m_cos_coef.data();
167 m_sin_data = m_sin_coef.data();
172 : Thick(other.m_ds, other.m_nslice),
173 SoftSolenoid_device_copyable(other.m_bscale, other.m_mapsteps)
175 #if !AMREX_DEVICE_COMPILE 181 m_ncoef = m_cos_coef.size();
182 m_cos_data = m_cos_coef.data();
183 m_sin_data = m_sin_coef.data();
190 Thick::operator=(other);
191 SoftSolenoid_device_copyable::operator=(other);
192 #if !AMREX_DEVICE_COMPILE 198 m_ncoef = m_cos_coef.size();
199 m_cos_data = m_cos_coef.data();
200 m_sin_data = m_sin_coef.data();
210 using BeamOptic::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
230 using namespace amrex::literals;
233 amrex::ParticleReal
const x = p.pos(0);
234 amrex::ParticleReal
const y = p.pos(1);
235 amrex::ParticleReal
const t = p.pos(2);
238 amrex::ParticleReal pxout = px;
239 amrex::ParticleReal pyout = py;
240 amrex::ParticleReal ptout = pt;
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;
276 void operator() (
RefPart & AMREX_RESTRICT refpart)
const 278 using namespace amrex::literals;
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;
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;
300 amrex::ParticleReal
const slice_ds = m_ds /
nslice();
303 amrex::ParticleReal
const bgi =
sqrt(
pow(pt, 2) - 1.0_prt);
306 amrex::ParticleReal
const zin = s - sedge;
307 amrex::ParticleReal
const zout = zin + slice_ds;
308 int const nsteps = m_mapsteps;
311 amrex::ParticleReal
const ptf = refpart.pt;
323 refpart.x = x + slice_ds*px/bgi;
324 refpart.y = y + slice_ds*py/bgi;
325 refpart.z = z + slice_ds*pz/bgi;
328 amrex::ParticleReal
const bgf =
sqrt(
pow(ptf, 2) - 1.0_prt);
331 refpart.px = px*bgf/bgi;
332 refpart.py = py*bgf/bgi;
333 refpart.pz = pz*bgf/bgi;
336 refpart.s = s + slice_ds;
345 std::tuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>
349 using namespace amrex::literals;
353 amrex::ParticleReal
const zlen = m_ds;
354 amrex::ParticleReal
const zmid = zlen / 2.0_prt;
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;
364 bfield = 0.5_prt*m_cos_data[0];
365 bfieldint = z*bfield;
366 for (
int j=1; j < m_ncoef; ++j)
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);
376 return std::make_tuple(bfield, bfieldp, bfieldint);
388 void map1 (amrex::ParticleReal
const tau,
390 [[maybe_unused]] amrex::ParticleReal & zeval)
const 392 using namespace amrex::literals;
395 amrex::ParticleReal
const t = refpart.
t;
396 amrex::ParticleReal
const pt = refpart.
pt;
397 amrex::ParticleReal
const z = zeval;
400 refpart.
t = t + tau/
sqrt(1.0_prt -
pow(pt, -2));
412 amrex::ParticleReal
const betgam = refpart.
beta_gamma();
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);
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);
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);
438 void map2 (amrex::ParticleReal
const tau,
440 amrex::ParticleReal & zeval)
const 442 using namespace amrex::literals;
444 amrex::ParticleReal
const t = refpart.
t;
445 amrex::ParticleReal
const pt = refpart.
pt;
448 amrex::ParticleReal
const B0 = m_bscale;
451 auto [bz, bzp, bzint] = Sol_Bfield(zeval);
459 amrex::ParticleReal
const alpha = B0*bz/2.0_prt;
460 amrex::ParticleReal
const alpha2 =
pow(alpha,2);
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);
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);
483 void map3 (amrex::ParticleReal
const tau,
485 amrex::ParticleReal & zeval)
const 487 using namespace amrex::literals;
489 amrex::ParticleReal
const t = refpart.
t;
490 amrex::ParticleReal
const pt = refpart.
pt;
491 amrex::ParticleReal
const z = zeval;
494 amrex::ParticleReal
const B0 = m_bscale;
497 auto [bz, bzp, bzint] = Sol_Bfield(z);
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);
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;
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;
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;
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;
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()
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
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
Particle< NStructReal, NStructInt > ParticleType
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void map3(amrex::ParticleReal const tau, RefPart &refpart, amrex::ParticleReal &zeval) const
Definition: SoftSol.H:483
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
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