diff options
| -rw-r--r-- | Source/Particles/PhotonParticleContainer.cpp | 9 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 2 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 176 | ||||
| -rw-r--r-- | Source/Particles/Pusher/CopyParticleAttribs.H | 94 | ||||
| -rw-r--r-- | Source/Particles/Pusher/PushSelector.H | 144 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 6 |
6 files changed, 274 insertions, 157 deletions
diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 90c8a5259..bd80b8e43 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -12,6 +12,7 @@ // Import low-level single-particle kernels #include "Particles/Pusher/UpdatePositionPhoton.H" #include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/Pusher/CopyParticleAttribs.H" #ifdef _OPENMP # include <omp.h> @@ -73,10 +74,9 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type) ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) - { - copy_attribs(pti); - } + auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data); + int do_copy = (WarpX::do_back_transformed_diagnostics && + do_back_transformed_diagnostics && a_dt_type!=DtType::SecondHalf); const auto GetPosition = GetParticlePosition(pti); auto SetPosition = SetParticlePosition(pti); @@ -84,6 +84,7 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type) amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + if (do_copy) copyAttribs(i); ParticleReal x, y, z; GetPosition(i, x, y, z); UpdatePositionPhoton( x, y, z, ux[i], uy[i], uz[i], dt ); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index e1426a482..5e076a2b5 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -171,8 +171,6 @@ public: RealVector& uzp, RealVector& wp ); - void copy_attribs (WarpXParIter& pti); - virtual void PostRestart () final {} void SplitParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2eee6aa79..094a8e556 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -18,17 +18,11 @@ #include "Utils/IonizationEnergiesTable.H" #include "Particles/Gather/FieldGather.H" #include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/Pusher/CopyParticleAttribs.H" +#include "Particles/Pusher/PushSelector.H" #include "Particles/Gather/GetExternalFields.H" - #include "Utils/WarpXAlgorithmSelection.H" -// Import low-level single-particle kernels -#include "Particles/Pusher/UpdatePosition.H" -#include "Particles/Pusher/UpdateMomentumBoris.H" -#include "Particles/Pusher/UpdateMomentumVay.H" -#include "Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H" -#include "Particles/Pusher/UpdateMomentumHigueraCary.H" - #include <AMReX_Print.H> #ifdef WARPX_USE_OPENPMD @@ -1490,10 +1484,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); - if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) - { - copy_attribs(pti); - } + auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data); + int do_copy = (WarpX::do_back_transformed_diagnostics && + do_back_transformed_diagnostics && + (a_dt_type!=DtType::SecondHalf)); int* AMREX_RESTRICT ion_lev = nullptr; if (do_field_ionization){ @@ -1504,111 +1498,28 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) const Real q = this->charge; const Real m = this-> mass; + const auto pusher_algo = WarpX::particle_pusher_algo; + const auto do_crr = do_classical_radiation_reaction; #ifdef WARPX_QED - if(do_classical_radiation_reaction){ - if(m_do_qed_quantum_sync){ - const auto t_chi_max = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - auto chi = QedUtils::chi_lepton(m*ux[i], m*uy[i], m*uz[i], - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i]); - if(chi < t_chi_max){ - UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - } - else{ - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - } - ParticleReal x, y, z; - GetPosition(i, x, y, z); - UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); - SetPosition(i, x, y, z); - } - ); - }else{ - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - ParticleReal x, y, z; - GetPosition(i, x, y, z); - UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); - SetPosition(i, x, y, z); - } - ); - } -#else - if(do_classical_radiation_reaction){ - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Real qp = q; - if (ion_lev){ qp *= ion_lev[i]; } - UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], qp, m, dt); - ParticleReal x, y, z; - GetPosition(i, x, y, z); - UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); - SetPosition(i, x, y, z); - } - ); + const auto do_sync = m_do_qed_quantum_sync; + amrex::Real t_chi_max = 0.0; + if (do_sync) t_chi_max = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; #endif - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Real qp = q; - if (ion_lev){ qp *= ion_lev[i]; } - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], qp, m, dt); - ParticleReal x, y, z; - GetPosition(i, x, y, z); - UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); - SetPosition(i, x, y, z); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Real qp = q; - if (ion_lev){ qp *= ion_lev[i]; } - UpdateMomentumVay( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], qp, m, dt); - ParticleReal x, y, z; - GetPosition(i, x, y, z); - UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); - SetPosition(i, x, y, z); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Real qp = q; - if (ion_lev){ qp *= ion_lev[i]; } - UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], qp, m, dt); - ParticleReal x, y, z; - GetPosition(i, x, y, z); - UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); - SetPosition(i, x, y, z); - } - ); - } else { - amrex::Abort("Unknown particle pusher"); - }; + + amrex::ParallelFor(pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + doParticlePush(GetPosition, SetPosition, copyAttribs, i, + ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + ion_lev ? ion_lev[i] : 0, + m, q, pusher_algo, do_crr, do_copy, +#ifdef WARPX_QED + do_sync, + t_chi_max, +#endif + dt); + }); } #ifdef WARPX_QED @@ -1771,41 +1682,6 @@ PhysicalParticleContainer::PushP ( } void -PhysicalParticleContainer::copy_attribs (WarpXParIter& pti) -{ - auto& attribs = pti.GetAttribs(); - ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); - ParticleReal* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); - ParticleReal* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - - const auto np = pti.numParticles(); - const auto lev = pti.GetLevel(); - const auto index = pti.GetPairIndex(); - ParticleReal* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); - ParticleReal* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); - ParticleReal* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); - ParticleReal* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); - ParticleReal* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); - ParticleReal* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); - - const auto GetPosition = GetParticlePosition(pti); - - ParallelFor( np, - [=] AMREX_GPU_DEVICE (long i) { - ParticleReal x, y, z; - GetPosition(i, x, y, z); - xpold[i]=x; - ypold[i]=y; - zpold[i]=z; - - uxpold[i]=uxp[i]; - uypold[i]=uyp[i]; - uzpold[i]=uzp[i]; - } - ); -} - -void PhysicalParticleContainer::GetParticleSlice ( const int direction, const Real z_old, const Real z_new, const Real t_boost, diff --git a/Source/Particles/Pusher/CopyParticleAttribs.H b/Source/Particles/Pusher/CopyParticleAttribs.H new file mode 100644 index 000000000..1f8a190ea --- /dev/null +++ b/Source/Particles/Pusher/CopyParticleAttribs.H @@ -0,0 +1,94 @@ +/* Copyright 2020 Andrew Myers + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PARTICLES_PUSHER_COPYPARTICLEATTRIBS_H_ +#define WARPX_PARTICLES_PUSHER_COPYPARTICLEATTRIBS_H_ + +#include "Particles/WarpXParticleContainer.H" + +#include <AMReX_REAL.H> + +#include <limits> + + +/** \brief Functor that creates copies of the current particle + * positions and momenta for later use. This is needed + * by the back-transformed diagnostics. + * + * \param a_pti iterator to the tile containing the macroparticles + * \param a_tmp holder for the temporary particle data +*/ +struct CopyParticleAttribs +{ + using TmpParticles = WarpXParticleContainer::TmpParticles; + + GetParticlePosition m_get_position; + + const amrex::ParticleReal* AMREX_RESTRICT uxp = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT uyp = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT uzp = nullptr; + + amrex::ParticleReal* AMREX_RESTRICT xpold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT ypold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT zpold = nullptr; + + amrex::ParticleReal* AMREX_RESTRICT uxpold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT uypold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT uzpold = nullptr; + + CopyParticleAttribs (const WarpXParIter& a_pti, TmpParticles& tmp_particle_data) noexcept + { + if (tmp_particle_data.size() == 0) return; + + auto& attribs = a_pti.GetAttribs(); + uxp = attribs[PIdx::ux].dataPtr(); + uyp = attribs[PIdx::uy].dataPtr(); + uzp = attribs[PIdx::uz].dataPtr(); + + const auto lev = a_pti.GetLevel(); + const auto index = a_pti.GetPairIndex(); + xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); + ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); + zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); + uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); + uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); + uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); + + m_get_position = GetParticlePosition(a_pti); + } + + /** \brief copy the position and momentum of particle i to the + * temporary data holder + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (const int i) const noexcept + { + AMREX_ASSERT(uxp != nullptr); + AMREX_ASSERT(uyp != nullptr); + AMREX_ASSERT(uzp != nullptr); + + AMREX_ASSERT(xpold != nullptr); + AMREX_ASSERT(ypold != nullptr); + AMREX_ASSERT(zpold != nullptr); + + AMREX_ASSERT(uxpold != nullptr); + AMREX_ASSERT(uypold != nullptr); + AMREX_ASSERT(uzpold != nullptr); + + amrex::ParticleReal x, y, z; + m_get_position(i, x, y, z); + + xpold[i] = x; + ypold[i] = y; + zpold[i] = z; + + uxpold[i] = uxp[i]; + uypold[i] = uyp[i]; + uzpold[i] = uzp[i]; + } +}; + +#endif // WARPX_PARTICLES_PUSHER_COPYPARTICLEATTRIBS_H_ diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H new file mode 100644 index 000000000..149342b74 --- /dev/null +++ b/Source/Particles/Pusher/PushSelector.H @@ -0,0 +1,144 @@ +/* Copyright 2020 Andrew Myers + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PARTICLES_PUSHER_PUSHSELECTOR_H_ +#define WARPX_PARTICLES_PUSHER_PUSHSELECTOR_H_ + +// Import low-level single-particle kernels +#include "Particles/Pusher/UpdatePosition.H" +#include "Particles/Pusher/UpdateMomentumBoris.H" +#include "Particles/Pusher/UpdateMomentumVay.H" +#include "Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H" +#include "Particles/Pusher/UpdateMomentumHigueraCary.H" +#include "Particles/WarpXParticleContainer.H" + +#include <AMReX_REAL.H> + +#include <limits> + +/** + * \brief Push position and momentum for a single particle + * + * /param GetPosition : A functor for returning the particle position. + * /param GetPosition : A functor for setting the particle position. + * /param copyAttribs : A functor for storing the old u and x + * /param i : The index of the particle to work on + * \param ux, uy, uz : Particle momentum + * \param Ex, Ey, Ez : Electric field on particles. + * \param Bx, By, Bz : Magnetic field on particles. + * \param ion_lev : Ionization level of this particle (0 if ioniziation not on) + * \param m : Mass of this species. + * \param q : Charge of this species. + * \param pusher_algo : 0: Boris, 1: Vay, 2: HigueraCary + * \param do_crr : Whether to do the classical radiation reaction + * \param do_crr : Whether to copy the old x and u for the BTD + * \param do_sync : Whether to include quantum synchrotron radiation (QSR) + * \param t_chi_max : Cutoff chi for QSR + * \param dt : Time step size + */ +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void doParticlePush(const GetParticlePosition& GetPosition, + const SetParticlePosition& SetPosition, + const CopyParticleAttribs& copyAttribs, + const long i, + amrex::ParticleReal& ux, + amrex::ParticleReal& uy, + amrex::ParticleReal& uz, + const amrex::ParticleReal Ex, + const amrex::ParticleReal Ey, + const amrex::ParticleReal Ez, + const amrex::ParticleReal Bx, + const amrex::ParticleReal By, + const amrex::ParticleReal Bz, + const int ion_lev, + const amrex::Real m, + const amrex::Real q, + const int pusher_algo, + const int do_crr, + const int do_copy, +#ifdef WARPX_QED + const int do_sync, + const amrex::Real t_chi_max, +#endif + const amrex::Real dt) +{ + if (do_copy) copyAttribs(i); + if (do_crr) { +#ifdef WARPX_QED + if (do_sync) { + auto chi = QedUtils::chi_lepton(m*ux, m*uy, m*uz, + Ex, Ey, Ez, + Bx, By, Bz); + if (chi < t_chi_max) { + UpdateMomentumBorisWithRadiationReaction(ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, q, m, dt); + } + else { + UpdateMomentumBoris( ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, q, m, dt); + } + amrex::ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux, uy, uz, dt ); + SetPosition(i, x, y, z); + } else { + UpdateMomentumBorisWithRadiationReaction(ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, q, m, dt); + amrex::ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux, uy, uz, dt ); + SetPosition(i, x, y, z); + } +#else + amrex::Real qp = q; + if (ion_lev) { qp *= ion_lev; } + UpdateMomentumBorisWithRadiationReaction(ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, qp, m, dt); + amrex::ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux, uy, uz, dt ); + SetPosition(i, x, y, z); +#endif + } else if (pusher_algo == ParticlePusherAlgo::Boris) { + amrex::Real qp = q; + if (ion_lev) { qp *= ion_lev; } + UpdateMomentumBoris( ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, qp, m, dt); + amrex::ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux, uy, uz, dt ); + SetPosition(i, x, y, z); + } else if (pusher_algo == ParticlePusherAlgo::Vay) { + amrex::Real qp = q; + if (ion_lev){ qp *= ion_lev; } + UpdateMomentumVay( ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, qp, m, dt); + amrex::ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux, uy, uz, dt ); + SetPosition(i, x, y, z); + } else if (pusher_algo == ParticlePusherAlgo::HigueraCary) { + amrex::Real qp = q; + if (ion_lev){ qp *= ion_lev; } + UpdateMomentumHigueraCary( ux, uy, uz, + Ex, Ey, Ez, Bx, + By, Bz, qp, m, dt); + amrex::ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux, uy, uz, dt ); + SetPosition(i, x, y, z); + } else { + amrex::Abort("Unknown particle pusher"); + } +} + +#endif // WARPX_PARTICLES_PUSHER_SELECTOR_H_ diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 477e70ce9..8bc4491d7 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -407,10 +407,14 @@ protected: amrex::Vector<amrex::FArrayBox> local_jy; amrex::Vector<amrex::FArrayBox> local_jz; +public: using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>; using PairIndex = std::pair<int, int>; + using TmpParticleTile = std::array<DataContainer, TmpIdx::nattribs>; + using TmpParticles = amrex::Vector<std::map<PairIndex, TmpParticleTile> >; - amrex::Vector<std::map<PairIndex, std::array<DataContainer, TmpIdx::nattribs> > > tmp_particle_data; +protected: + TmpParticles tmp_particle_data; /** * When using runtime components, AMReX requires to touch all tiles |