diff options
Diffstat (limited to 'Source/Particles')
-rwxr-xr-x | Source/Particles/Deposition/ChargeDeposition.H | 18 | ||||
-rw-r--r-- | Source/Particles/Deposition/CurrentDeposition.H | 56 | ||||
-rw-r--r-- | Source/Particles/Gather/FieldGather.H | 23 | ||||
-rw-r--r-- | Source/Particles/PhotonParticleContainer.H | 3 | ||||
-rw-r--r-- | Source/Particles/PhotonParticleContainer.cpp | 24 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 15 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 263 | ||||
-rw-r--r-- | Source/Particles/Pusher/GetAndSetPosition.H | 143 | ||||
-rw-r--r-- | Source/Particles/Pusher/UpdatePosition.H | 8 | ||||
-rw-r--r-- | Source/Particles/Pusher/UpdatePositionPhoton.H | 1 | ||||
-rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.H | 6 | ||||
-rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 197 | ||||
-rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 10 | ||||
-rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 109 |
14 files changed, 394 insertions, 482 deletions
diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 669fc4c57..b03e4224f 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -8,10 +8,11 @@ #ifndef CHARGEDEPOSITION_H_ #define CHARGEDEPOSITION_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" /* \brief Charge Deposition for thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param wp : Pointer to array of particle weights. * \param ion_lev : Pointer to array of particle ionization level. This is required to have the charge of each macroparticle @@ -25,9 +26,7 @@ * /param q : species charge. */ template <int depos_order> -void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doChargeDepositionShapeN(const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, const int * const ion_lev, const amrex::Array4<amrex::Real>& rho_arr, @@ -63,14 +62,17 @@ void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, wq *= ion_lev[ip]; } + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + // --- Compute shape factors // x direction // Get particle position in grid coordinates #if (defined WARPX_DIM_RZ) - const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real r = std::sqrt(xp*xp + yp*yp); const amrex::Real x = (r - xmin)*dxi; #else - const amrex::Real x = (xp[ip] - xmin)*dxi; + const amrex::Real x = (xp - xmin)*dxi; #endif // Compute shape factors for node-centered quantities amrex::Real sx[depos_order + 1]; @@ -79,12 +81,12 @@ void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_3D) // y direction - const amrex::Real y = (yp[ip] - ymin)*dyi; + const amrex::Real y = (yp - ymin)*dyi; amrex::Real sy[depos_order + 1]; const int j = compute_shape_factor<depos_order>(sy, y); #endif // z direction - const amrex::Real z = (zp[ip] - zmin)*dzi; + const amrex::Real z = (zp - zmin)*dzi; amrex::Real sz[depos_order + 1]; const int k = compute_shape_factor<depos_order>(sz, z); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 90039dea2..af3b0006b 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -8,6 +8,7 @@ #ifndef CURRENTDEPOSITION_H_ #define CURRENTDEPOSITION_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" #include <WarpX_Complex.H> @@ -16,7 +17,7 @@ /** * \brief Current Deposition for thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. * \param ion_lev : Pointer to array of particle ionization level. This is @@ -34,9 +35,7 @@ * /param q : species charge. */ template <int depos_order> -void doDepositionShapeN(const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doDepositionShapeN(const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, const amrex::ParticleReal * const uxp, const amrex::ParticleReal * const uyp, @@ -95,6 +94,10 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, if (do_ionization){ wq *= ion_lev[ip]; } + + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + const amrex::Real vx = uxp[ip]*gaminv; const amrex::Real vy = uyp[ip]*gaminv; const amrex::Real vz = uzp[ip]*gaminv; @@ -102,8 +105,8 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_RZ) // In RZ, wqx is actually wqr, and wqy is wqtheta // Convert to cylinderical at the mid point - const amrex::Real xpmid = xp[ip] - 0.5*dt*vx; - const amrex::Real ypmid = yp[ip] - 0.5*dt*vy; + const amrex::Real xpmid = xp - 0.5*dt*vx; + const amrex::Real ypmid = yp - 0.5*dt*vy; const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid); amrex::Real costheta; amrex::Real sintheta; @@ -128,7 +131,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_RZ) const amrex::Real xmid = (rpmid - xmin)*dxi; #else - const amrex::Real xmid = (xp[ip] - xmin)*dxi - dts2dx*vx; + const amrex::Real xmid = (xp - xmin)*dxi - dts2dx*vx; #endif // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current // sx_j[xyz] shape factor along x for the centering of each current @@ -153,7 +156,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_3D) // y direction - const amrex::Real ymid = (yp[ip] - ymin)*dyi - dts2dy*vy; + const amrex::Real ymid = (yp - ymin)*dyi - dts2dy*vy; amrex::Real sy_node[depos_order + 1]; amrex::Real sy_cell[depos_order + 1]; int k_node; @@ -173,7 +176,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #endif // z direction - const amrex::Real zmid = (zp[ip] - zmin)*dzi - dts2dz*vz; + const amrex::Real zmid = (zp - zmin)*dzi - dts2dz*vz; amrex::Real sz_node[depos_order + 1]; amrex::Real sz_cell[depos_order + 1]; int l_node; @@ -230,7 +233,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, /** * \brief Esirkepov Current Deposition for thread thread_num * - * \param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. * \param ion_lev : Pointer to array of particle ionization level. This is @@ -249,9 +252,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order> -void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, const amrex::ParticleReal * const uxp, const amrex::ParticleReal * const uyp, @@ -308,14 +309,18 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, // --- Get particle quantities Real const gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq - + uyp[ip]*uyp[ip]*clightsq - + uzp[ip]*uzp[ip]*clightsq); + + uyp[ip]*uyp[ip]*clightsq + + uzp[ip]*uzp[ip]*clightsq); // wqx, wqy wqz are particle current in each direction Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; } + + ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + Real const wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) Real const wqy = wq*invdtdy; @@ -324,17 +329,18 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, // computes current and old position in grid units #if (defined WARPX_DIM_RZ) - Real const xp_mid = xp[ip] - 0.5_rt * dt*uxp[ip]*gaminv; - Real const yp_mid = yp[ip] - 0.5_rt * dt*uyp[ip]*gaminv; - Real const xp_old = xp[ip] - dt*uxp[ip]*gaminv; - Real const yp_old = yp[ip] - dt*uyp[ip]*gaminv; - Real const rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + Real const xp_mid = xp - 0.5_rt * dt*uxp[ip]*gaminv; + Real const yp_mid = yp - 0.5_rt * dt*uyp[ip]*gaminv; + Real const xp_old = xp - dt*uxp[ip]*gaminv; + Real const yp_old = yp - dt*uyp[ip]*gaminv; + Real const rp_new = std::sqrt(xp*xp + + yp*yp); Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid); Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); Real costheta_new, sintheta_new; if (rp_new > 0._rt) { - costheta_new = xp[ip]/rp_new; - sintheta_new = yp[ip]/rp_new; + costheta_new = xp/rp_new; + sintheta_new = yp/rp_new; } else { costheta_new = 1.; sintheta_new = 0.; @@ -361,14 +367,14 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, Real const x_new = (rp_new - xmin)*dxi; Real const x_old = (rp_old - xmin)*dxi; #else - Real const x_new = (xp[ip] - xmin)*dxi; + Real const x_new = (xp - xmin)*dxi; Real const x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif #if (defined WARPX_DIM_3D) - Real const y_new = (yp[ip] - ymin)*dyi; + Real const y_new = (yp - ymin)*dyi; Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv; #endif - Real const z_new = (zp[ip] - zmin)*dzi; + Real const z_new = (zp - zmin)*dzi; Real const z_old = z_new - dtsdz0*uzp[ip]*gaminv; #if (defined WARPX_DIM_RZ) diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 0a58f8425..12d9b6291 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -8,12 +8,13 @@ #ifndef FIELDGATHER_H_ #define FIELDGATHER_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" #include <WarpX_Complex.H> /** * \brief Field gather for particles handled by thread thread_num - * \param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param Exp, Eyp, Ezp: Pointer to array of electric field on particles. * \param Bxp, Byp, Bzp: Pointer to array of magnetic field on particles. * \param ex_arr ey_arr: Array4 of current density, either full array or tile. @@ -27,9 +28,7 @@ * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order, int lower_in_v> -void doGatherShapeN(const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doGatherShapeN(const GetParticlePosition& GetPosition, amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, @@ -80,14 +79,18 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, amrex::ParallelFor( np_to_gather, [=] AMREX_GPU_DEVICE (long ip) { + + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + // --- Compute shape factors // x direction // Get particle position #ifdef WARPX_DIM_RZ - const amrex::Real rp = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real rp = std::sqrt(xp*xp + yp*yp); const amrex::Real x = (rp - xmin)*dxi; #else - const amrex::Real x = (xp[ip]-xmin)*dxi; + const amrex::Real x = (xp-xmin)*dxi; #endif // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current @@ -129,7 +132,7 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, #if (AMREX_SPACEDIM == 3) // y direction - const amrex::Real y = (yp[ip]-ymin)*dyi; + const amrex::Real y = (yp-ymin)*dyi; amrex::Real sy_node[depos_order + 1]; amrex::Real sy_cell[depos_order + 1]; amrex::Real sy_node_v[depos_order + 1 - lower_in_v]; @@ -165,7 +168,7 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, #endif // z direction - const amrex::Real z = (zp[ip]-zmin)*dzi; + const amrex::Real z = (zp-zmin)*dzi; amrex::Real sz_node[depos_order + 1]; amrex::Real sz_cell[depos_order + 1]; amrex::Real sz_node_v[depos_order + 1 - lower_in_v]; @@ -245,8 +248,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, amrex::Real costheta; amrex::Real sintheta; if (rp > 0.) { - costheta = xp[ip]/rp; - sintheta = yp[ip]/rp; + costheta = xp/rp; + sintheta = yp/rp; } else { costheta = 1.; sintheta = 0.; diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 7750d5ce8..ef32ebcb0 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -55,9 +55,6 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; // Do nothing diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index ab85170ac..d007d3c1c 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -21,7 +21,7 @@ // Import low-level single-particle kernels #include <UpdatePositionPhoton.H> - +#include <GetAndSetPosition.H> using namespace amrex; @@ -57,19 +57,13 @@ void PhotonParticleContainer::InitData() } void -PhotonParticleContainer::PushPX(WarpXParIter& pti, - Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp, - Real dt, DtType a_dt_type) +PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); + // Extract pointers to the different particle quantities - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); @@ -82,15 +76,19 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { - copy_attribs(pti, x, y, z); + copy_attribs(pti); } + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - - UpdatePositionPhoton( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePositionPhoton( x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 0500d091a..421dcd842 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -72,10 +72,7 @@ public: void AssignExternalFieldOnParticles ( WarpXParIter& pti, RealVector& Exp, RealVector& Eyp, RealVector& Ezp, RealVector& Bxp, - RealVector& Byp, RealVector& Bzp, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, int lev); + RealVector& Byp, RealVector& Bzp, int lev); virtual void FieldGather (int lev, const amrex::MultiFab& Ex, @@ -101,7 +98,6 @@ public: const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, - int thread_num, int lev, int depos_lev); @@ -163,11 +159,7 @@ public: amrex::Real dt, DtType a_dt_type=DtType::Full) override; - virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, - amrex::Real dt, DtType a_dt_type=DtType::Full); + virtual void PushPX (WarpXParIter& pti, amrex::Real dt, DtType a_dt_type=DtType::Full); virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, @@ -190,8 +182,7 @@ public: RealVector& uzp, RealVector& wp ); - void copy_attribs(WarpXParIter& pti,const amrex::ParticleReal* xp, - const amrex::ParticleReal* yp, const amrex::ParticleReal* zp); + void copy_attribs (WarpXParIter& pti); virtual void PostRestart () final {} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4c98f436e..6572657ff 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -20,6 +20,7 @@ #include <WarpXWrappers.h> #include <IonizationEnergiesTable.H> #include <FieldGather.H> +#include <GetAndSetPosition.H> #include <WarpXAlgorithmSelection.H> @@ -967,10 +968,7 @@ PhysicalParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<Mul void PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, RealVector& Exp, RealVector& Eyp, RealVector& Ezp, - RealVector& Bxp, RealVector& Byp, RealVector& Bzp, - const Gpu::ManagedDeviceVector<ParticleReal>& xp, - const Gpu::ManagedDeviceVector<ParticleReal>& yp, - const Gpu::ManagedDeviceVector<ParticleReal>& zp, int lev) + RealVector& Bxp, RealVector& Byp, RealVector& Bzp, int lev) { const long np = pti.numParticles(); /// get WarpX class object @@ -990,9 +988,7 @@ PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, Bzp.assign(np,mypc.m_B_external_particle[2]); } if (mypc.m_E_ext_particle_s=="parse_e_ext_particle_function") { - const Real* const AMREX_RESTRICT xp_data = xp.dataPtr(); - const Real* const AMREX_RESTRICT yp_data = yp.dataPtr(); - const Real* const AMREX_RESTRICT zp_data = zp.dataPtr(); + const auto GetPosition = GetParticlePosition(pti); Real* const AMREX_RESTRICT Exp_data = Exp.dataPtr(); Real* const AMREX_RESTRICT Eyp_data = Eyp.dataPtr(); Real* const AMREX_RESTRICT Ezp_data = Ezp.dataPtr(); @@ -1001,20 +997,20 @@ PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, ParserWrapper *zfield_partparser = mypc.m_Ez_particle_parser.get(); Real time = warpx.gett_new(lev); amrex::ParallelFor(pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Exp_data[i] = xfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Eyp_data[i] = yfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Ezp_data[i] = zfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - }, + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal x, y, z; + GetPosition(i, x, y, z); + Exp_data[i] = xfield_partparser->getField(x, y, z, time); + Eyp_data[i] = yfield_partparser->getField(x, y, z, time); + Ezp_data[i] = zfield_partparser->getField(x, y, z, time); + }, /* To allocate shared memory for the GPU threads. */ /* But, for now only 4 doubles (x,y,z,t) are allocated. */ amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 ); } if (mypc.m_B_ext_particle_s=="parse_b_ext_particle_function") { - const Real* const AMREX_RESTRICT xp_data = xp.dataPtr(); - const Real* const AMREX_RESTRICT yp_data = yp.dataPtr(); - const Real* const AMREX_RESTRICT zp_data = zp.dataPtr(); + const auto GetPosition = GetParticlePosition(pti); Real* const AMREX_RESTRICT Bxp_data = Bxp.dataPtr(); Real* const AMREX_RESTRICT Byp_data = Byp.dataPtr(); Real* const AMREX_RESTRICT Bzp_data = Bzp.dataPtr(); @@ -1024,10 +1020,12 @@ PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, Real time = warpx.gett_new(lev); amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - Bxp_data[i] = xfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Byp_data[i] = yfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Bzp_data[i] = zfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - }, + ParticleReal x, y, z; + GetPosition(i, x, y, z); + Bxp_data[i] = xfield_partparser->getField(x, y, z, time); + Byp_data[i] = yfield_partparser->getField(x, y, z, time); + Bzp_data[i] = zfield_partparser->getField(x, y, z, time); + }, /* To allocate shared memory for the GPU threads. */ /* But, for now only 4 doubles (x,y,z,t) are allocated. */ amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 @@ -1056,11 +1054,6 @@ PhysicalParticleContainer::FieldGather (int lev, #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { Real wt = amrex::second(); @@ -1087,18 +1080,13 @@ PhysicalParticleContainer::FieldGather (int lev, const FArrayBox& bzfab = Bz[pti]; // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - - // // Field Gather // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, - 0, np, thread_num, lev, lev); + 0, np, lev, lev); if (cost) { const Box& tbx = pti.tilebox(); @@ -1110,9 +1098,6 @@ PhysicalParticleContainer::FieldGather (int lev, costarr(i,j,k) += wt; }); } - // synchronize avoids cudaStreams from over-writing the temporary arrays used to - // store positions - Gpu::synchronize(); } } } @@ -1238,13 +1223,6 @@ PhysicalParticleContainer::Evolve (int lev, const long np_current = (cjx) ? nfine_current : np; - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); - if (rho) { // Deposit charge before particle push, in component 0 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1274,7 +1252,7 @@ PhysicalParticleContainer::Evolve (int lev, FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, - 0, np_gather, thread_num, lev, lev); + 0, np_gather, lev, lev); if (np_gather < np) { @@ -1310,7 +1288,7 @@ PhysicalParticleContainer::Evolve (int lev, cbxfab, cbyfab, cbzfab, cEx->nGrow(), e_is_nodal, nfine_gather, np-nfine_gather, - thread_num, lev, lev-1); + lev, lev-1); } BL_PROFILE_VAR_STOP(blp_fg); @@ -1328,7 +1306,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], dt, a_dt_type); + PushPX(pti, dt, a_dt_type); BL_PROFILE_VAR_STOP(blp_ppc_pp); // @@ -1352,14 +1330,6 @@ PhysicalParticleContainer::Evolve (int lev, np_current, np-np_current, thread_num, lev, lev-1, dt); } - - - // - // copy particle data back - // - BL_PROFILE_VAR_START(blp_copy); - pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); } if (rho) { @@ -1478,7 +1448,6 @@ PhysicalParticleContainer::SplitParticles(int lev) { auto& mypc = WarpX::GetInstance().GetPartContainer(); auto& pctmp_split = mypc.GetPCtmp(); - Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; RealVector psplit_x, psplit_y, psplit_z, psplit_w; RealVector psplit_ux, psplit_uy, psplit_uz; long np_split_to_add = 0; @@ -1493,7 +1462,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // Loop over particle interator for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - pti.GetPosition(xp, yp, zp); + const auto GetPosition = GetParticlePosition(pti); const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); @@ -1518,6 +1487,8 @@ PhysicalParticleContainer::SplitParticles(int lev) auto& uzp = attribs[PIdx::uz]; const long np = pti.numParticles(); for(int i=0; i<np; i++){ + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); auto& p = particles[i]; if (p.id() == DoSplitParticleID){ // If particle is tagged, split it and put the @@ -1530,9 +1501,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int ishift = -1; ishift < 2; ishift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x and z - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + kshift*split_offset[2] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1544,17 +1515,17 @@ PhysicalParticleContainer::SplitParticles(int lev) // 4 particles in 2d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in z - psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*split_offset[2] ); + psplit_x.push_back( xp ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1569,9 +1540,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int jshift = -1; jshift < 2; jshift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x, y and z - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] + jshift*split_offset[1] ); - psplit_z.push_back( zp[i] + kshift*split_offset[2] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp + jshift*split_offset[1] ); + psplit_z.push_back( zp + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1584,25 +1555,25 @@ PhysicalParticleContainer::SplitParticles(int lev) // 6 particles in 3d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in y - psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] + ishift*split_offset[1] ); - psplit_z.push_back( zp[i] ); + psplit_x.push_back( xp ); + psplit_y.push_back( yp + ishift*split_offset[1] ); + psplit_z.push_back( zp ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in z - psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*split_offset[2] ); + psplit_x.push_back( xp ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1639,19 +1610,16 @@ PhysicalParticleContainer::SplitParticles(int lev) } void -PhysicalParticleContainer::PushPX(WarpXParIter& pti, - Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp, - Real dt, DtType a_dt_type) +PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); @@ -1664,7 +1632,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) { - copy_attribs(pti, x, y, z); + copy_attribs(pti); } int* AMREX_RESTRICT ion_lev = nullptr; @@ -1697,8 +1665,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); } - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], 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{ @@ -1708,8 +1678,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], 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); } ); } @@ -1723,8 +1695,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], 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); } ); #endif @@ -1737,8 +1711,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumBoris( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], 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) { @@ -1750,8 +1726,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumVay( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], 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) { @@ -1763,8 +1741,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], 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 { @@ -1830,11 +1810,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt, #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1858,16 +1833,11 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, - 0, np, thread_num, lev, lev); + 0, np, lev, lev); // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities @@ -1944,8 +1914,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, - const ParticleReal* yp, const ParticleReal* zp) +void PhysicalParticleContainer::copy_attribs (WarpXParIter& pti) { auto& attribs = pti.GetAttribs(); ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); @@ -1962,11 +1931,15 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleRea 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) { - xpold[i]=xp[i]; - ypold[i]=yp[i]; - zpold[i]=zp[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]; @@ -2024,11 +1997,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -2037,10 +2005,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real if ( !slice_box.intersects(tile_real_box) ) continue; - pti.GetPosition(m_xp[thread_num],m_yp[thread_num],m_zp[thread_num]); - Real *const AMREX_RESTRICT xpnew = m_xp[thread_num].dataPtr(); - Real *const AMREX_RESTRICT ypnew = m_yp[thread_num].dataPtr(); - Real *const AMREX_RESTRICT zpnew = m_zp[thread_num].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); auto& attribs = pti.GetAttribs(); Real* const AMREX_RESTRICT wpnew = attribs[PIdx::w].dataPtr(); @@ -2078,12 +2043,14 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) { - Flag[i] = 0; - if ( (((zpnew[i] >= z_new) && (zpold[i] <= z_old)) || - ((zpnew[i] <= z_new) && (zpold[i] >= z_old))) ) - { - Flag[i] = 1; - } + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + Flag[i] = 0; + if ( (((zp >= z_new) && (zpold[i] <= z_old)) || + ((zp <= z_new) && (zpold[i] >= z_old))) ) + { + Flag[i] = 1; + } }); // exclusive scan to obtain location indices using flag values @@ -2120,6 +2087,8 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) { + ParticleReal xp_new, yp_new, zp_new; + GetPosition(i, xp_new, yp_new, zp_new); if (Flag[i] == 1) { // Lorentz Transform particles to lab-frame @@ -2127,12 +2096,10 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real (uxpnew[i]*uxpnew[i] + uypnew[i]*uypnew[i] + uzpnew[i]*uzpnew[i])); - const Real t_new_p = gammaboost*t_boost - - uzfrm*zpnew[i]*inv_c2; - const Real z_new_p = gammaboost*(zpnew[i] - + betaboost*Phys_c*t_boost); - const Real uz_new_p = gammaboost*uzpnew[i] - - gamma_new_p*uzfrm; + const Real t_new_p = gammaboost*t_boost - uzfrm*zp_new*inv_c2; + const Real z_new_p = gammaboost*(zp_new + betaboost*Phys_c*t_boost); + const Real uz_new_p = gammaboost*uzpnew[i] - gamma_new_p*uzfrm; + const Real gamma_old_p = std::sqrt(1.0 + inv_c2* (uxpold[i]*uxpold[i] + uypold[i]*uypold[i] @@ -2150,12 +2117,10 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real const Real weight_new = (t_lab - t_old_p) / (t_new_p - t_old_p); - const Real xp = xpold[i]*weight_old - + xpnew[i]*weight_new; - const Real yp = ypold[i]*weight_old - + ypnew[i]*weight_new; - const Real zp = z_old_p*weight_old - + z_new_p*weight_new; + const Real xp = xpold[i]*weight_old + xp_new*weight_new; + const Real yp = ypold[i]*weight_old + yp_new*weight_new; + const Real zp = z_old_p*weight_old + z_new_p*weight_new; + const Real uxp = uxpold[i]*weight_old + uxpnew[i]*weight_new; const Real uyp = uypold[i]*weight_old @@ -2197,7 +2162,6 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param e_is_nodal: 0 if E is staggered, 1 if E is nodal * \param offset: index of first particle for which fields are gathered * \param np_to_gather: number of particles onto which fields are gathered - * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. @@ -2219,7 +2183,6 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, - int thread_num, int lev, int gather_lev) { @@ -2231,9 +2194,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // initializing the field value to the externally applied field before // gathering fields from the grid to the particles. - AssignExternalFieldOnParticles(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - m_xp[thread_num], m_yp[thread_num], - m_zp[thread_num], lev); + AssignExternalFieldOnParticles(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, lev); // Get cell size on gather_lev @@ -2252,9 +2213,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // Add guard cells to the box. box.grow(ngE); - const ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - const ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - const ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const auto GetPosition = GetParticlePosition(pti, offset); // Lower corner of tile box physical domain const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev); @@ -2265,7 +2224,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ if (WarpX::nox == 1){ - doGatherShapeN<1,1>(xp, yp, zp, + doGatherShapeN<1,1>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2273,7 +2232,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doGatherShapeN<2,1>(xp, yp, zp, + doGatherShapeN<2,1>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2281,7 +2240,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doGatherShapeN<3,1>(xp, yp, zp, + doGatherShapeN<3,1>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2291,7 +2250,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } else { if (WarpX::nox == 1){ - doGatherShapeN<1,0>(xp, yp, zp, + doGatherShapeN<1,0>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2299,7 +2258,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doGatherShapeN<2,0>(xp, yp, zp, + doGatherShapeN<2,0>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2307,7 +2266,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doGatherShapeN<3,0>(xp, yp, zp, + doGatherShapeN<3,0>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 7180f3c55..594260703 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -12,72 +12,99 @@ #include <WarpXParticleContainer.H> #include <AMReX_REAL.H> -#ifndef WARPX_DIM_RZ - -/** \brief Extract the particle's coordinates from the ParticleType struct `p`, - * and stores them in the variables `x`, `y`, `z`. */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void GetPosition( - amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const WarpXParticleContainer::ParticleType& p) +/** \brief Functor that can be used to extract the positions of the macroparticles + * inside a ParallelFor kernel + * + * \param a_pti iterator to the tile containing the macroparticles + * \param a_offset offset to apply to the particle indices +*/ +struct GetParticlePosition { -#if (AMREX_SPACEDIM==3) - x = p.pos(0); - y = p.pos(1); - z = p.pos(2); -#else - x = p.pos(0); - y = std::numeric_limits<amrex::ParticleReal>::quiet_NaN(); - z = p.pos(1); + using PType = WarpXParticleContainer::ParticleType; + using RType = amrex::ParticleReal; + + const PType* AMREX_RESTRICT m_structs; +#if (defined WARPX_DIM_RZ) + const RType* m_theta; +#elif (AMREX_SPACEDIM == 2) + static constexpr RType m_snan = std::numeric_limits<RType>::quiet_NaN(); #endif -} + GetParticlePosition (const WarpXParIter& a_pti, int a_offset = 0) noexcept + { + const auto& aos = a_pti.GetArrayOfStructs(); + m_structs = aos().dataPtr() + a_offset; +#if (defined WARPX_DIM_RZ) + const auto& soa = a_pti.GetStructOfArrays(); + m_theta = soa.GetRealData(PIdx::theta).dataPtr() + a_offset; +#endif + } -/** \brief Set the particle's coordinates in the ParticleType struct `p`, - * from their values in the variables `x`, `y`, `z`. */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void SetPosition( - WarpXParticleContainer::ParticleType& p, - const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z) -{ -#if (AMREX_SPACEDIM==3) - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; + /** \brief Extract the cartesian position coordinates of the particle + * located at index `i + a_offset` and store them in the variables + * `x`, `y`, `z` */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (const int i, RType& x, RType& y, RType& z) const noexcept + { +#ifdef WARPX_DIM_RZ + RType r = m_structs[i].pos(0); + x = r*std::cos(m_theta[i]); + y = r*std::sin(m_theta[i]); + z = m_structs[i].pos(1); +#elif WARPX_DIM_3D + x = m_structs[i].pos(0); + y = m_structs[i].pos(1); + z = m_structs[i].pos(2); #else - p.pos(0) = x; - p.pos(1) = z; + x = m_structs[i].pos(0); + y = m_snan; + z = m_structs[i].pos(1); #endif -} - -# elif defined WARPX_DIM_RZ + } +}; -/** \brief Extract the particle's coordinates from `theta` and the attributes - * of the ParticleType struct `p` (which contains the radius), - * and store them in the variables `x`, `y`, `z` */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void GetCartesianPositionFromCylindrical( - amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const WarpXParticleContainer::ParticleType& p, const amrex::ParticleReal theta) +/** \brief Functor that can be used to modify the positions of the macroparticles, + * inside a ParallelFor kernel. + * + * \param a_pti iterator to the tile being modified + * \param a_offset offset to apply to the particle indices +*/ +struct SetParticlePosition { - const amrex::ParticleReal r = p.pos(0); - x = r*std::cos(theta); - y = r*std::sin(theta); - z = p.pos(1); -} + using PType = WarpXParticleContainer::ParticleType; + using RType = amrex::ParticleReal; -/** \brief Set the particle's cylindrical coordinates by setting `theta` - * and the attributes of the ParticleType struct `p` (which stores the radius), - * from the values of `x`, `y`, `z` */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void SetCylindricalPositionFromCartesian( - WarpXParticleContainer::ParticleType& p, amrex::ParticleReal& theta, - const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z ) -{ - theta = std::atan2(y, x); - p.pos(0) = std::sqrt(x*x + y*y); - p.pos(1) = z; -} + PType* AMREX_RESTRICT m_structs; +#if (defined WARPX_DIM_RZ) + RType* AMREX_RESTRICT m_theta; +#endif + SetParticlePosition (WarpXParIter& a_pti, int a_offset = 0) noexcept + { + auto& aos = a_pti.GetArrayOfStructs(); + m_structs = aos().dataPtr() + a_offset; +#if (defined WARPX_DIM_RZ) + auto& soa = a_pti.GetStructOfArrays(); + m_theta = soa.GetRealData(PIdx::theta).dataPtr() + a_offset; +#endif + } -#endif // WARPX_DIM_RZ + /** \brief Set the position of the particle at index `i + a_offset` + * to `x`, `y`, `z` */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (const int i, RType x, RType y, RType z) const noexcept + { +#ifdef WARPX_DIM_RZ + m_theta[i] = std::atan2(y, x); + m_structs[i].pos(0) = std::sqrt(x*x + y*y); + m_structs[i].pos(1) = z; +#elif WARPX_DIM_3D + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = y; + m_structs[i].pos(2) = z; +#else + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = z; +#endif + } +}; #endif // WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 7f86a106d..1968a1439 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -15,10 +15,9 @@ /** \brief Push the particle's positions over one timestep, * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE -void UpdatePosition( - amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, - const amrex::Real dt ) +void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, + const amrex::Real dt ) { constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); @@ -31,7 +30,6 @@ void UpdatePosition( y += uy * inv_gamma * dt; #endif z += uz * inv_gamma * dt; - } #endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H index 3a28e87a1..44c0afda9 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -32,7 +32,6 @@ void UpdatePositionPhoton( y += uy * c_over_umod * dt; #endif z += uz * c_over_umod * dt; - } #endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index 5e5749093..7ad7fd82e 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -68,11 +68,7 @@ public: amrex::Real dt, DtType a_dt_type=DtType::Full) override; - virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, - amrex::Real dt, DtType a_dt_type=DtType::Full) override; + virtual void PushPX (WarpXParIter& pti, amrex::Real dt, DtType a_dt_type=DtType::Full) override; virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index f9db682d7..953874c0f 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -24,6 +24,7 @@ #include <UpdateMomentumVay.H> #include <UpdateMomentumBorisWithRadiationReaction.H> #include <UpdateMomentumHigueraCary.H> +#include <GetAndSetPosition.H> using namespace amrex; @@ -87,8 +88,6 @@ RigidInjectedParticleContainer::RemapParticles() // Note that the particles are already in the boosted frame. // This value is saved to advance the particles not injected yet - Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { auto& attribs = pti.GetAttribs(); @@ -96,31 +95,32 @@ RigidInjectedParticleContainer::RemapParticles() auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; - // Copy data from particle container to temp arrays - pti.GetPosition(xp, yp, zp); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); // Loop over particles const long np = pti.numParticles(); for (int i=0 ; i < np ; i++) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + const Real gammapr = std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/csq); const Real vzpr = uzp[i]/gammapr; // Back out the value of z_lab - const Real z_lab = (zp[i] + uz_boost*t_lab + WarpX::gamma_boost*t_lab*vzpr)/(WarpX::gamma_boost + uz_boost*vzpr/csq); + const Real z_lab = (zp + uz_boost*t_lab + WarpX::gamma_boost*t_lab*vzpr)/(WarpX::gamma_boost + uz_boost*vzpr/csq); // Time of the particle in the boosted frame given its position in the lab frame at t=0. const Real tpr = WarpX::gamma_boost*t_lab - uz_boost*z_lab/csq; // Adjust the position, taking away its motion from its own velocity and adding // the motion from the average velocity - zp[i] = zp[i] + tpr*vzpr - tpr*vzbeam_ave_boosted; - - } + zp += tpr*vzpr - tpr*vzbeam_ave_boosted; - // Copy the data back to the particle container - pti.SetPosition(xp, yp, zp); + SetPosition(i, xp, yp, zp); + } } } } @@ -147,8 +147,6 @@ RigidInjectedParticleContainer::BoostandRemapParticles() #pragma omp parallel #endif { - Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; - for (WarpXParIter pti(*this, 0); pti.isValid(); ++pti) { @@ -157,13 +155,16 @@ RigidInjectedParticleContainer::BoostandRemapParticles() auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; - // Copy data from particle container to temp arrays - pti.GetPosition(xp, yp, zp); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); // Loop over particles const long np = pti.numParticles(); for (int i=0 ; i < np ; i++) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + const Real gamma_lab = std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/(PhysConst::c*PhysConst::c)); const Real vx_lab = uxp[i]/gamma_lab; @@ -173,24 +174,24 @@ RigidInjectedParticleContainer::BoostandRemapParticles() // t0_lab is the time in the lab frame that the particles reaches z=0 // The location and time (z=0, t=0) is a synchronization point between the // lab and boosted frames. - const Real t0_lab = -zp[i]/vz_lab; + const Real t0_lab = -zp/vz_lab; if (!projected) { - xp[i] += t0_lab*vx_lab; - yp[i] += t0_lab*vy_lab; + xp += t0_lab*vx_lab; + yp += t0_lab*vy_lab; } if (focused) { // Correct for focusing effect from shift from z=0 to zinject const Real tfocus = -zinject_plane*WarpX::gamma_boost/vz_lab; - xp[i] -= tfocus*vx_lab; - yp[i] -= tfocus*vy_lab; + xp -= tfocus*vx_lab; + yp -= tfocus*vy_lab; } // Time of the particle in the boosted frame given its position in the lab frame at t=0. - const Real tpr = -WarpX::gamma_boost*WarpX::beta_boost*zp[i]/PhysConst::c; + const Real tpr = -WarpX::gamma_boost*WarpX::beta_boost*zp/PhysConst::c; // Position of the particle in the boosted frame given its position in the lab frame at t=0. - const Real zpr = WarpX::gamma_boost*zp[i]; + const Real zpr = WarpX::gamma_boost*zp; // Momentum of the particle in the boosted frame (assuming that it is fixed). uzp[i] = WarpX::gamma_boost*(uzp[i] - WarpX::beta_boost*PhysConst::c*gamma_lab); @@ -198,30 +199,23 @@ RigidInjectedParticleContainer::BoostandRemapParticles() // Put the particle at the location in the boosted frame at boost frame t=0, if (rigid_advance) { // with the particle moving at the average velocity - zp[i] = zpr - vzbeam_ave_boosted*tpr; + zp = zpr - vzbeam_ave_boosted*tpr; } else { // with the particle moving with its own velocity const Real gammapr = std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/(PhysConst::c*PhysConst::c)); const Real vzpr = uzp[i]/gammapr; - zp[i] = zpr - vzpr*tpr; + zp = zpr - vzpr*tpr; } + SetPosition(i, xp, yp, zp); } - - // Copy the data back to the particle container - pti.SetPosition(xp, yp, zp); - } } } void -RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, - Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp, - Real dt, DtType a_dt_type) +RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. @@ -234,9 +228,9 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Gpu::ManagedDeviceVector<ParticleReal> xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + ParticleReal* const AMREX_RESTRICT ux = uxp.dataPtr(); ParticleReal* const AMREX_RESTRICT uy = uyp.dataPtr(); ParticleReal* const AMREX_RESTRICT uz = uzp.dataPtr(); @@ -247,14 +241,38 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, ParticleReal* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr(); ParticleReal* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); - if (!done_injecting_lev) { + if (!done_injecting_lev) + { // If the old values are not already saved, create copies here. - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; + const auto np = pti.numParticles(); + + xp_save.resize(np); + yp_save.resize(np); + zp_save.resize(np); + + uxp_save.resize(np); + uyp_save.resize(np); + uzp_save.resize(np); + + amrex::Real* const AMREX_RESTRICT xp_save_ptr = xp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT yp_save_ptr = yp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT zp_save_ptr = zp_save.dataPtr(); + + amrex::Real* const AMREX_RESTRICT uxp_save_ptr = uxp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT uyp_save_ptr = uyp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT uzp_save_ptr = uzp_save.dataPtr(); + + amrex::ParallelFor( np, + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + xp_save_ptr[i] = xp; + yp_save_ptr[i] = yp; + zp_save_ptr[i] = zp; + uxp_save_ptr[i] = ux[i]; + uyp_save_ptr[i] = uy[i]; + uzp_save_ptr[i] = uz[i]; + }); // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -265,20 +283,22 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, const Real vz_ave_boosted = vzbeam_ave_boosted; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const Real dtscale = dt - (z_plane_previous - z[i])/(vz_ave_boosted + v_boost); - if (0. < dtscale && dtscale < dt) { - Exp[i] *= dtscale; - Eyp[i] *= dtscale; - Ezp[i] *= dtscale; - Bxp[i] *= dtscale; - Byp[i] *= dtscale; - Bzp[i] *= dtscale; - } - } - ); + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + const Real dtscale = dt - (z_plane_previous - zp)/(vz_ave_boosted + v_boost); + if (0. < dtscale && dtscale < dt) { + Exp[i] *= dtscale; + Eyp[i] *= dtscale; + Ezp[i] *= dtscale; + Bxp[i] *= dtscale; + Byp[i] *= dtscale; + Bzp[i] *= dtscale; + } + } + ); } - PhysicalParticleContainer::PushPX(pti, xp, yp, zp, dt, a_dt_type); + PhysicalParticleContainer::PushPX(pti, dt, a_dt_type); if (!done_injecting_lev) { @@ -296,23 +316,25 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, const bool rigid = rigid_advance; const Real inv_csq = 1./(PhysConst::c*PhysConst::c); amrex::ParallelFor( pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - if (z[i] <= z_plane_lev) { - ux[i] = ux_save[i]; - uy[i] = uy_save[i]; - uz[i] = uz_save[i]; - x[i] = x_save[i]; - y[i] = y_save[i]; - if (rigid) { - z[i] = z_save[i] + dt*vz_ave_boosted; - } - else { - const Real gi = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_csq); - z[i] = z_save[i] + dt*uz[i]*gi; - } - } - } - ); + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + if (zp <= z_plane_lev) { + ux[i] = ux_save[i]; + uy[i] = uy_save[i]; + uz[i] = uz_save[i]; + xp = x_save[i]; + yp = y_save[i]; + if (rigid) { + zp = z_save[i] + dt*vz_ave_boosted; + } + else { + const Real gi = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_csq); + zp = z_save[i] + dt*uz[i]*gi; + } + SetPosition(i, xp, yp, zp); + } + }); } } @@ -370,11 +392,6 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -401,16 +418,11 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, - 0, np, thread_num, lev, lev); + 0, np, lev, lev); // Save the position and momenta, making copies auto uxp_save = uxp; @@ -419,7 +431,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - const ParticleReal* const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); ParticleReal* const AMREX_RESTRICT uxpp = uxp.dataPtr(); ParticleReal* const AMREX_RESTRICT uypp = uyp.dataPtr(); ParticleReal* const AMREX_RESTRICT uzpp = uzp.dataPtr(); @@ -485,15 +497,16 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const ParticleReal* const AMREX_RESTRICT uz_save = uzp_save.dataPtr(); const ParticleReal zz = zinject_plane_levels[lev]; amrex::ParallelFor( pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - if (zp[i] <= zz) { - uxpp[i] = ux_save[i]; - uypp[i] = uy_save[i]; - uzpp[i] = uz_save[i]; - } - } - ); - + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + if (zp <= zz) { + uxpp[i] = ux_save[i]; + uypp[i] = uy_save[i]; + uzpp[i] = uz_save[i]; + } + } + ); } } } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index cac1d9883..b1def39d3 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -81,14 +81,6 @@ public: WarpXParIter (ContainerType& pc, int level); -#if (AMREX_SPACEDIM == 2) - void GetPosition (amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& x, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& y, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& z) const; - void SetPosition (const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& x, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& y, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& z); -#endif const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { return GetStructOfArrays().GetRealData(); } @@ -379,8 +371,6 @@ protected: using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>; using PairIndex = std::pair<int, int>; - amrex::Vector<DataContainer> m_xp, m_yp, m_zp; - // Whether to dump particle quantities. // If true, particle position is always dumped. int plot_species = 1; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 85cac0dbe..cd5bcaaf9 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -30,56 +30,6 @@ WarpXParIter::WarpXParIter (ContainerType& pc, int level) { } -#if (AMREX_SPACEDIM == 2) -void -WarpXParIter::GetPosition (Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp) const -{ - amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(xp, zp); -#ifdef WARPX_DIM_RZ - const auto& attribs = GetAttribs(); - const auto& thetap = attribs[PIdx::theta]; - yp.resize(xp.size()); - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT theta = thetap.dataPtr(); - amrex::ParallelFor( xp.size(), - [=] AMREX_GPU_DEVICE (long i) { - // The x stored in the particles is actually the radius - y[i] = x[i]*std::sin(theta[i]); - x[i] = x[i]*std::cos(theta[i]); - }); -#else - yp.resize(xp.size(), std::numeric_limits<ParticleReal>::quiet_NaN()); -#endif -} - -void -WarpXParIter::SetPosition (const Gpu::ManagedDeviceVector<ParticleReal>& xp, - const Gpu::ManagedDeviceVector<ParticleReal>& yp, - const Gpu::ManagedDeviceVector<ParticleReal>& zp) -{ -#ifdef WARPX_DIM_RZ - auto& attribs = GetAttribs(); - auto& thetap = attribs[PIdx::theta]; - Gpu::ManagedDeviceVector<ParticleReal> rp(xp.size()); - const ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT r = rp.dataPtr(); - ParticleReal* const AMREX_RESTRICT theta = thetap.dataPtr(); - amrex::ParallelFor( xp.size(), - [=] AMREX_GPU_DEVICE (long i) { - theta[i] = std::atan2(y[i], x[i]); - r[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); - }); - amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(rp, zp); -#else - amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(xp, zp); -#endif -} -#endif - WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) : ParticleContainer<0,0,PIdx::nattribs>(amr_core->GetParGDB()) , species_id(ispecies) @@ -116,10 +66,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) local_jx.resize(num_threads); local_jy.resize(num_threads); local_jz.resize(num_threads); - m_xp.resize(num_threads); - m_yp.resize(num_threads); - m_zp.resize(num_threads); - } void @@ -349,9 +295,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // CPU, tiling: deposit into local_jx // (same for jx and jz) - ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const auto GetPosition = GetParticlePosition(pti, offset); // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -362,19 +306,19 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doEsirkepovDepositionShapeN<2>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doEsirkepovDepositionShapeN<3>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); @@ -382,19 +326,19 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } else { if (WarpX::nox == 1){ doDepositionShapeN<1>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ doDepositionShapeN<2>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ doDepositionShapeN<3>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q); @@ -488,9 +432,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // GPU, no tiling: deposit directly in rho // CPU, tiling: deposit into local_rho - ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const auto GetPosition = GetParticlePosition(pti, offset); // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -500,13 +442,13 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<1>(GetPosition, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<2>(GetPosition, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<3>(GetPosition, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -545,8 +487,6 @@ WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::Mult const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -615,8 +555,6 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -772,32 +710,27 @@ WarpXParticleContainer::PushX (int lev, amrex::Real dt) // // Particle Push // - // Extract pointers to particle position and momenta, for this particle tile - // - positions are stored as an array of struct, in `ParticleType` - ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); + + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + // - momenta are stored as a struct of array, in `attribs` auto& attribs = pti.GetAttribs(); ParticleReal* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); #ifdef WARPX_DIM_RZ + auto& aos = pti.GetArrayOfStructs(); + ParticleType* AMREX_RESTRICT const pstruct = aos().dataPtr(); ParticleReal* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); #endif // Loop over the particles and update their position amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - ParticleType& p = pstructs[i]; // Particle object that gets updated - ParticleReal x, y, z; // Temporary variables -#ifndef WARPX_DIM_RZ - GetPosition( x, y, z, p ); // Initialize x, y, z - UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); - SetPosition( p, x, y, z ); // Update the object p -#else - // For WARPX_DIM_RZ, the particles are still pushed in 3D Cartesian - GetCartesianPositionFromCylindrical( x, y, z, p, theta[i] ); - UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); - SetCylindricalPositionFromCartesian( p, theta[i], x, y, z ); -#endif + 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); } ); |