From d2476c8843de693a192fd16e4899741502665134 Mon Sep 17 00:00:00 2001 From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com> Date: Mon, 27 Jun 2022 07:56:51 -0700 Subject: [PATCH] Added operator to push reference particle in global coordinates to all elements (#149) * Added math to push the reference particle in global coordinates to all existing elements. * Elements: Pass by Reference * Pass the reference particle by reference, so we can manipulate it in elements. * Remove true "no-op" implementations with a comment. * Add Doxygen comments. * Push the Reference Particle Co-authored-by: Axel Huebl --- src/particles/ImpactXParticleContainer.H | 4 +-- src/particles/ImpactXParticleContainer.cpp | 4 +-- src/particles/Push.cpp | 12 ++++--- src/particles/elements/ConstF.H | 31 ++++++++++++++++++ src/particles/elements/DipEdge.H | 11 +++++++ src/particles/elements/Drift.H | 34 +++++++++++++++++++- src/particles/elements/Multipole.H | 11 +++++++ src/particles/elements/NonlinearLens.H | 11 +++++++ src/particles/elements/Quad.H | 31 ++++++++++++++++++ src/particles/elements/Sbend.H | 37 ++++++++++++++++++++++ src/particles/elements/ShortRF.H | 12 +++++++ 11 files changed, 188 insertions(+), 10 deletions(-) diff --git a/src/particles/ImpactXParticleContainer.H b/src/particles/ImpactXParticleContainer.H index 6066c9c1a..ca38f13d3 100644 --- a/src/particles/ImpactXParticleContainer.H +++ b/src/particles/ImpactXParticleContainer.H @@ -144,8 +144,8 @@ namespace impactx * * @returns refpart */ - RefPart - GetRefParticle () const; + RefPart & + GetRefParticle (); /** Get particle shape */ diff --git a/src/particles/ImpactXParticleContainer.cpp b/src/particles/ImpactXParticleContainer.cpp index 055005ce5..c245efe1c 100644 --- a/src/particles/ImpactXParticleContainer.cpp +++ b/src/particles/ImpactXParticleContainer.cpp @@ -123,8 +123,8 @@ namespace impactx m_refpart = refpart; } - RefPart - ImpactXParticleContainer::GetRefParticle () const + RefPart & + ImpactXParticleContainer::GetRefParticle () { return m_refpart; } diff --git a/src/particles/Push.cpp b/src/particles/Push.cpp index 9707ec679..8cf9f1713 100644 --- a/src/particles/Push.cpp +++ b/src/particles/Push.cpp @@ -126,18 +126,20 @@ namespace detail // ... // preparing to access reference particle data: RefPart - RefPart ref_part; - ref_part = pc.GetRefParticle(); + RefPart & ref_part = pc.GetRefParticle(); // loop over all beamline elements for (auto & element_variant : lattice) { // here we just access the element by its respective type - std::visit([=](auto&& element) { + std::visit([=, &ref_part](auto&& element) { + // push beam particles relative to reference particle detail::PushSingleParticle const pushSingleParticle( element, aos_ptr, part_px, part_py, part_pt, ref_part); - - // loop over particles in the box + // loop over beam particles in the box amrex::ParallelFor(np, pushSingleParticle); + + // push reference particle in global coordinates + element(ref_part); }, element_variant); }; // end loop over all beamline elements } // end loop over all particle boxes diff --git a/src/particles/elements/ConstF.H b/src/particles/elements/ConstF.H index ded5519c3..30ad7b76d 100644 --- a/src/particles/elements/ConstF.H +++ b/src/particles/elements/ConstF.H @@ -88,6 +88,37 @@ namespace impactx } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + RefPart & AMREX_RESTRICT refpart) const { + + using namespace amrex::literals; // for _rt and _prt + + // assign input reference particle values + amrex::ParticleReal const x = refpart.x; + amrex::ParticleReal const px = refpart.px; + amrex::ParticleReal const y = refpart.y; + amrex::ParticleReal const py = refpart.py; + amrex::ParticleReal const z = refpart.z; + amrex::ParticleReal const pz = refpart.pz; + amrex::ParticleReal const t = refpart.t; + amrex::ParticleReal const pt = refpart.pt; + + // assign intermediate parameter + amrex::ParticleReal const s = m_ds/sqrt(pow(pt,2)-1.0_prt); + + // advance position and momentum (straight element) + refpart.x = x + s*px; + refpart.y = y + s*py; + refpart.z = z + s*pz; + refpart.t = t - s*pt; + + } + private: amrex::ParticleReal m_ds; //! segment length in m amrex::ParticleReal m_kx; //! focusing x strength in 1/m diff --git a/src/particles/elements/DipEdge.H b/src/particles/elements/DipEdge.H index a9a48ccde..949fa2eb7 100644 --- a/src/particles/elements/DipEdge.H +++ b/src/particles/elements/DipEdge.H @@ -84,6 +84,17 @@ namespace impactx py = py + R43*y; } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + [[maybe_unused]] RefPart & AMREX_RESTRICT refpart) const { + + // nothing to do: this is a zero-length element + } + private: amrex::ParticleReal m_psi; //! pole face angle in rad amrex::ParticleReal m_rc; //! bend radius in m diff --git a/src/particles/elements/Drift.H b/src/particles/elements/Drift.H index 1f773f7a0..c1ad7b99b 100644 --- a/src/particles/elements/Drift.H +++ b/src/particles/elements/Drift.H @@ -56,7 +56,7 @@ namespace impactx amrex::ParticleReal const y = p.pos(1); amrex::ParticleReal const t = p.pos(2); - // intialize output values of momenta + // initialize output values of momenta amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; @@ -77,6 +77,38 @@ namespace impactx px = pxout; py = pyout; pt = ptout; + + } + + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + RefPart & AMREX_RESTRICT refpart) const { + + using namespace amrex::literals; // for _rt and _prt + + // assign input reference particle values + amrex::ParticleReal const x = refpart.x; + amrex::ParticleReal const px = refpart.px; + amrex::ParticleReal const y = refpart.y; + amrex::ParticleReal const py = refpart.py; + amrex::ParticleReal const z = refpart.z; + amrex::ParticleReal const pz = refpart.pz; + amrex::ParticleReal const t = refpart.t; + amrex::ParticleReal const pt = refpart.pt; + + // assign intermediate parameter + amrex::ParticleReal const s = m_ds/sqrt(pow(pt,2)-1.0_prt); + + // advance position and momentum (drift) + refpart.x = x + s*px; + refpart.y = y + s*py; + refpart.z = z + s*pz; + refpart.t = t - s*pt; + } private: diff --git a/src/particles/elements/Multipole.H b/src/particles/elements/Multipole.H index cfb108989..d5868f5a7 100644 --- a/src/particles/elements/Multipole.H +++ b/src/particles/elements/Multipole.H @@ -108,6 +108,17 @@ namespace impactx } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + [[maybe_unused]] RefPart & AMREX_RESTRICT refpart) const { + + // nothing to do: this is a zero-length element + } + private: int m_multipole; //! multipole index int m_mfactorial; //! factorial of multipole index diff --git a/src/particles/elements/NonlinearLens.H b/src/particles/elements/NonlinearLens.H index 526067c86..89638affd 100644 --- a/src/particles/elements/NonlinearLens.H +++ b/src/particles/elements/NonlinearLens.H @@ -117,6 +117,17 @@ namespace impactx } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + [[maybe_unused]] RefPart & AMREX_RESTRICT refpart) const { + + // nothing to do: this is a zero-length element + } + private: amrex::ParticleReal m_knll; //! integrated strength of the nonlinear lens (m) amrex::ParticleReal m_cnll; //! distance of singularities from the origin (m) diff --git a/src/particles/elements/Quad.H b/src/particles/elements/Quad.H index 64516c904..7cc7febd6 100644 --- a/src/particles/elements/Quad.H +++ b/src/particles/elements/Quad.H @@ -101,6 +101,37 @@ namespace impactx } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + RefPart & AMREX_RESTRICT refpart) const { + + using namespace amrex::literals; // for _rt and _prt + + // assign input reference particle values + amrex::ParticleReal const x = refpart.x; + amrex::ParticleReal const px = refpart.px; + amrex::ParticleReal const y = refpart.y; + amrex::ParticleReal const py = refpart.py; + amrex::ParticleReal const z = refpart.z; + amrex::ParticleReal const pz = refpart.pz; + amrex::ParticleReal const t = refpart.t; + amrex::ParticleReal const pt = refpart.pt; + + // assign intermediate parameter + amrex::ParticleReal const s = m_ds/sqrt(pow(pt,2)-1.0_prt); + + // advance position and momentum (straight element) + refpart.x = x + s*px; + refpart.y = y + s*py; + refpart.z = z + s*pz; + refpart.t = t - s*pt; + + } + private: amrex::ParticleReal m_ds; //! segment length in m amrex::ParticleReal m_k; //! quadrupole strength in 1/m diff --git a/src/particles/elements/Sbend.H b/src/particles/elements/Sbend.H index a1d3fb2a3..f0ae78563 100644 --- a/src/particles/elements/Sbend.H +++ b/src/particles/elements/Sbend.H @@ -91,6 +91,43 @@ namespace impactx } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + RefPart & AMREX_RESTRICT refpart) const { + + using namespace amrex::literals; // for _rt and _prt + + // assign input reference particle values + amrex::ParticleReal const x = refpart.x; + amrex::ParticleReal const px = refpart.px; + amrex::ParticleReal const y = refpart.y; + amrex::ParticleReal const py = refpart.py; + amrex::ParticleReal const z = refpart.z; + amrex::ParticleReal const pz = refpart.pz; + amrex::ParticleReal const t = refpart.t; + amrex::ParticleReal const pt = refpart.pt; + + // assign intermediate parameter + amrex::ParticleReal const theta = m_ds/m_rc; + amrex::ParticleReal const B = sqrt(pow(pt,2)-1.0_prt)/m_rc; + + // advance position and momentum (bend) + refpart.px = px*cos(theta) - pz*sin(theta); + refpart.py = py; + refpart.pz = pz*cos(theta) + px*sin(theta); + refpart.pt = pt; + + refpart.x = x + (refpart.pz - pz)/B; + refpart.y = y + (theta/B)*py; + refpart.z = z - (refpart.px - px)/B; + refpart.t = t - (theta/B)*pt; + + } + private: amrex::ParticleReal m_ds; //! segment length in m amrex::ParticleReal m_rc; //! bend radius in m diff --git a/src/particles/elements/ShortRF.H b/src/particles/elements/ShortRF.H index 615195f96..c980484cc 100644 --- a/src/particles/elements/ShortRF.H +++ b/src/particles/elements/ShortRF.H @@ -85,6 +85,18 @@ namespace impactx } + /** This pushes the reference particle. + * + * @param[in,out] refpart reference particle + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() ( + [[maybe_unused]] RefPart & AMREX_RESTRICT refpart) const { + + // nothing to do: this is a zero-length element + + } + private: amrex::ParticleReal m_V; //! normalized (max) RF voltage drop. amrex::ParticleReal m_k; //! RF wavenumber in 1/m.