From 4b4f13857bd4fd623096a367b784e30fe15a8810 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 27 Jan 2020 13:42:39 -0800 Subject: remove the copies between soa and aos for the particle positions --- Source/Particles/PhysicalParticleContainer.H | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.H') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 31d3cbbf3..99665f98a 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -62,10 +62,7 @@ public: void AssignExternalFieldOnParticles ( WarpXParIter& pti, RealVector& Exp, RealVector& Eyp, RealVector& Ezp, RealVector& Bxp, - RealVector& Byp, RealVector& Bzp, - const amrex::Gpu::ManagedDeviceVector& xp, - const amrex::Gpu::ManagedDeviceVector& yp, - const amrex::Gpu::ManagedDeviceVector& zp, int lev); + RealVector& Byp, RealVector& Bzp, int lev); virtual void FieldGather (int lev, const amrex::MultiFab& Ex, @@ -153,11 +150,7 @@ public: amrex::Real dt, DtType a_dt_type=DtType::Full) override; - virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector& xp, - amrex::Gpu::ManagedDeviceVector& yp, - amrex::Gpu::ManagedDeviceVector& 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, @@ -180,8 +173,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, const ParticleType* pstruct); virtual void PostRestart () final {} -- cgit v1.2.3 From f5328965e6bcb338e50be09bb9fa6af8ccc6a6e6 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 28 Jan 2020 11:28:25 -0800 Subject: remove the thread_num variable from places where it is no longer needed --- Source/Particles/PhysicalParticleContainer.H | 1 - Source/Particles/PhysicalParticleContainer.cpp | 25 ++++--------------------- 2 files changed, 4 insertions(+), 22 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.H') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 99665f98a..4cb71a5d2 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -88,7 +88,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); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4a35e5e12..a660adadd 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1056,11 +1056,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(); @@ -1093,7 +1088,7 @@ PhysicalParticleContainer::FieldGather (int lev, 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(); @@ -1259,7 +1254,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) { @@ -1295,7 +1290,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); @@ -1802,11 +1797,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(); @@ -1834,7 +1824,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, 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 @@ -1990,11 +1980,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(); @@ -2157,7 +2142,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. @@ -2179,7 +2163,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) { -- cgit v1.2.3 From 80074a19c67a17a9b8169dc64004322bbd0c31b5 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 28 Jan 2020 13:24:49 -0800 Subject: switch deposition, gather, and pushers to use the get / set position functor --- Source/Particles/Deposition/ChargeDeposition.H | 17 +++--- Source/Particles/Deposition/CurrentDeposition.H | 49 ++++++++------- Source/Particles/Gather/FieldGather.H | 22 ++++--- Source/Particles/PhotonParticleContainer.cpp | 13 ++-- Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 69 ++++++++++++++-------- Source/Particles/Pusher/GetAndSetPosition.H | 12 ++-- Source/Particles/Pusher/UpdatePosition.H | 8 +-- Source/Particles/Pusher/UpdatePositionPhoton.H | 9 ++- .../Particles/RigidInjectedParticleContainer.cpp | 7 +-- Source/Particles/WarpXParticleContainer.cpp | 45 +++++++------- 11 files changed, 144 insertions(+), 109 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.H') diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 53861984e..9e3df600b 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -1,10 +1,11 @@ #ifndef CHARGEDEPOSITION_H_ #define CHARGEDEPOSITION_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" /* \brief Charge Deposition for thread thread_num - * /param pstruct : Pointer to array of particle structs + * /param get_position : 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 @@ -18,7 +19,7 @@ * /param q : species charge. */ template -void doChargeDepositionShapeN(const WarpXParticleContainer::ParticleType* const pstruct, +void doChargeDepositionShapeN(const GetPosition& get_position, const amrex::ParticleReal * const wp, const int * const ion_lev, const amrex::Array4& rho_arr, @@ -54,15 +55,17 @@ void doChargeDepositionShapeN(const WarpXParticleContainer::ParticleType* const wq *= ion_lev[ip]; } + amrex::Real xp, yp, zp; + get_position(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(pstruct[ip].pos(0)*pstruct[ip].pos(0) - + pstruct[ip].pos(1)*pstruct[ip].pos(1)); + const amrex::Real r = std::sqrt(xp*xp + yp*yp); const amrex::Real x = (r - xmin)*dxi; #else - const amrex::Real x = (pstruct[ip].pos(0) - xmin)*dxi; + const amrex::Real x = (xp - xmin)*dxi; #endif // Compute shape factors for node-centered quantities amrex::Real sx[depos_order + 1]; @@ -71,12 +74,12 @@ void doChargeDepositionShapeN(const WarpXParticleContainer::ParticleType* const #if (defined WARPX_DIM_3D) // y direction - const amrex::Real y = (pstruct[ip].pos(1) - ymin)*dyi; + const amrex::Real y = (yp - ymin)*dyi; amrex::Real sy[depos_order + 1]; const int j = compute_shape_factor(sy, y); #endif // z direction - const amrex::Real z = (pstruct[ip].pos(2) - zmin)*dzi; + const amrex::Real z = (zp - zmin)*dzi; amrex::Real sz[depos_order + 1]; const int k = compute_shape_factor(sz, z); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index b5e5c13a6..8e7f1a552 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -1,6 +1,7 @@ #ifndef CURRENTDEPOSITION_H_ #define CURRENTDEPOSITION_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" #include @@ -9,7 +10,7 @@ /** * \brief Current Deposition for thread thread_num - * /param pstruct : Pointer to array of particle structs. + * /param get_position : 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 @@ -27,7 +28,7 @@ * /param q : species charge. */ template -void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstruct, +void doDepositionShapeN(const GetPosition& get_position, const amrex::ParticleReal * const wp, const amrex::ParticleReal * const uxp, const amrex::ParticleReal * const uyp, @@ -86,6 +87,10 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru if (do_ionization){ wq *= ion_lev[ip]; } + + amrex::Real xp, yp, zp; + get_position(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; @@ -93,8 +98,8 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru #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 = pstruct[ip].pos(0) - 0.5*dt*vx; - const amrex::Real ypmid = pstruct[ip].pos(1) - 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; @@ -119,7 +124,7 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru #if (defined WARPX_DIM_RZ) const amrex::Real xmid = (rpmid - xmin)*dxi; #else - const amrex::Real xmid = (pstruct[ip].pos(0) - 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 @@ -144,7 +149,7 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru #if (defined WARPX_DIM_3D) // y direction - const amrex::Real ymid = (pstruct[ip].pos(1) - 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; @@ -164,7 +169,7 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru #endif // z direction - const amrex::Real zmid = (pstruct[ip].pos(2) - 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; @@ -221,7 +226,7 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru /** * \brief Esirkepov Current Deposition for thread thread_num * - * \param pstruct : Pointer to array of particle structs. + * /param get_position : 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 @@ -240,7 +245,7 @@ void doDepositionShapeN(const WarpXParticleContainer::ParticleType * const pstru * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template -void doEsirkepovDepositionShapeN (const WarpXParticleContainer::ParticleType * pstruct, +void doEsirkepovDepositionShapeN (const GetPosition& get_position, const amrex::ParticleReal * const wp, const amrex::ParticleReal * const uxp, const amrex::ParticleReal * const uyp, @@ -305,6 +310,10 @@ void doEsirkepovDepositionShapeN (const WarpXParticleContainer::ParticleType * p if (do_ionization){ wq *= ion_lev[ip]; } + + Real xp, yp, zp; + get_position(ip, xp, yp, zp); + Real const wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) Real const wqy = wq*invdtdy; @@ -313,18 +322,18 @@ void doEsirkepovDepositionShapeN (const WarpXParticleContainer::ParticleType * p // computes current and old position in grid units #if (defined WARPX_DIM_RZ) - Real const xp_mid = pstruct[ip].pos(0) - 0.5_rt * dt*uxp[ip]*gaminv; - Real const yp_mid = pstruct[ip].pos(1) - 0.5_rt * dt*uyp[ip]*gaminv; - Real const xp_old = pstruct[ip].pos(0) - dt*uxp[ip]*gaminv; - Real const yp_old = pstruct[ip].pos(1) - dt*uyp[ip]*gaminv; - Real const rp_new = std::sqrt(pstruct[ip].pos(0)*pstruct[ip].pos(0) - + pstruct[ip].pos(1)*pstruct[ip].pos(1)); + 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 = pstruct[ip].pos(0)/rp_new; - sintheta_new = pstruct[ip].pos(1)/rp_new; + costheta_new = xp/rp_new; + sintheta_new = yp/rp_new; } else { costheta_new = 1.; sintheta_new = 0.; @@ -351,14 +360,14 @@ void doEsirkepovDepositionShapeN (const WarpXParticleContainer::ParticleType * p Real const x_new = (rp_new - xmin)*dxi; Real const x_old = (rp_old - xmin)*dxi; #else - Real const x_new = (pstruct[ip].pos(0) - 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 = (pstruct[ip].pos(1) - ymin)*dyi; + Real const y_new = (yp - ymin)*dyi; Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv; #endif - Real const z_new = (pstruct[ip].pos(2) - 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 0af1e362a..369b648c1 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -1,12 +1,13 @@ #ifndef FIELDGATHER_H_ #define FIELDGATHER_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" #include /** * \brief Field gather for particles handled by thread thread_num - * \param pstruct : Pointer to array of particle structs. + * /param get_position : 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. @@ -20,7 +21,7 @@ * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template -void doGatherShapeN(const WarpXParticleContainer::ParticleType * const pstruct, +void doGatherShapeN(const GetPosition& get_position, amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, @@ -71,15 +72,18 @@ void doGatherShapeN(const WarpXParticleContainer::ParticleType * const pstruct, amrex::ParallelFor( np_to_gather, [=] AMREX_GPU_DEVICE (long ip) { + + amrex::Real xp, yp, zp; + get_position(ip, xp, yp, zp); + // --- Compute shape factors // x direction // Get particle position #ifdef WARPX_DIM_RZ - const amrex::Real rp = std::sqrt(pstruct[ip].pos(0)*pstruct[ip].pos(0) - + pstruct[ip].pos(1)*pstruct[ip].pos(1)); + const amrex::Real rp = std::sqrt(xp*xp + yp*yp); const amrex::Real x = (rp - xmin)*dxi; #else - const amrex::Real x = (pstruct[ip].pos(0)-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 @@ -121,7 +125,7 @@ void doGatherShapeN(const WarpXParticleContainer::ParticleType * const pstruct, #if (AMREX_SPACEDIM == 3) // y direction - const amrex::Real y = (pstruct[ip].pos(1)-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]; @@ -157,7 +161,7 @@ void doGatherShapeN(const WarpXParticleContainer::ParticleType * const pstruct, #endif // z direction - const amrex::Real z = (pstruct[ip].pos(2)-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]; @@ -237,8 +241,8 @@ void doGatherShapeN(const WarpXParticleContainer::ParticleType * const pstruct, amrex::Real costheta; amrex::Real sintheta; if (rp > 0.) { - costheta = pstruct[ip].pos(0)/rp; - sintheta = pstruct[ip].pos(1)/rp; + costheta = xp/rp; + sintheta = yp/rp; } else { costheta = 1.; sintheta = 0.; diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index a4d738f93..86cc206a7 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -14,7 +14,7 @@ // Import low-level single-particle kernels #include - +#include using namespace amrex; @@ -72,14 +72,19 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type) if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { - copy_attribs(pti, pstruct); + copy_attribs(pti); } + const auto get_position = GetPosition(pti); + auto set_position = SetPosition(pti); + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - - UpdatePositionPhoton( pstruct[i], ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePositionPhoton( x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 4cb71a5d2..392d13244 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -172,7 +172,7 @@ public: RealVector& uzp, RealVector& wp ); - void copy_attribs (WarpXParIter& pti, const ParticleType* pstruct); + void copy_attribs (WarpXParIter& pti); virtual void PostRestart () final {} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a660adadd..110d4294c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1618,8 +1618,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - auto& aos = pti.GetArrayOfStructs(); - ParticleType* AMREX_RESTRICT const pstruct = aos().dataPtr(); + const auto get_position = GetPosition(pti); + auto set_position = SetPosition(pti); ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); @@ -1633,7 +1633,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) { - copy_attribs(pti, pstruct); + copy_attribs(pti); } int* AMREX_RESTRICT ion_lev = nullptr; @@ -1666,7 +1666,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); } - UpdatePosition( pstruct[i], ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); }else{ @@ -1676,7 +1679,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - UpdatePosition( pstruct[i], ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); } @@ -1690,8 +1696,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( pstruct[i], - ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); #endif @@ -1704,8 +1712,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) UpdateMomentumBoris( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( pstruct[i], - ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { @@ -1717,8 +1727,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) UpdateMomentumVay( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( pstruct[i], - ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { @@ -1730,8 +1742,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( pstruct[i], - ux[i], uy[i], uz[i], dt ); + Real x, y, z; + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + set_position(i, x, y, z); } ); } else { @@ -1901,7 +1915,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::copy_attribs (WarpXParIter& pti, const ParticleType* pstruct) +void PhysicalParticleContainer::copy_attribs (WarpXParIter& pti) { auto& attribs = pti.GetAttribs(); ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); @@ -1918,11 +1932,15 @@ void PhysicalParticleContainer::copy_attribs (WarpXParIter& pti, const ParticleT 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 get_position = GetPosition(pti); + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { - xpold[i]=pstruct[i].pos(0); - ypold[i]=pstruct[i].pos(1); - zpold[i]=pstruct[i].pos(2); + Real x, y, z; + get_position(i, x, y, z); + xpold[i]=x; + ypold[i]=y; + zpold[i]=z; uxpold[i]=uxp[i]; uypold[i]=uyp[i]; @@ -2193,9 +2211,8 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // Add guard cells to the box. box.grow(ngE); - const auto& aos = pti.GetArrayOfStructs(); - const ParticleType* AMREX_RESTRICT pstruct = aos().dataPtr() + offset; - + const auto get_position = GetPosition(pti, offset); + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); @@ -2205,7 +2222,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>(pstruct, + doGatherShapeN<1,1>(get_position, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2213,7 +2230,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doGatherShapeN<2,1>(pstruct, + doGatherShapeN<2,1>(get_position, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2221,7 +2238,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doGatherShapeN<3,1>(pstruct, + doGatherShapeN<3,1>(get_position, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2231,7 +2248,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } else { if (WarpX::nox == 1){ - doGatherShapeN<1,0>(pstruct, + doGatherShapeN<1,0>(get_position, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2239,7 +2256,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doGatherShapeN<2,0>(pstruct, + doGatherShapeN<2,0>(get_position, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2247,7 +2264,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doGatherShapeN<3,0>(pstruct, + doGatherShapeN<3,0>(get_position, 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 e3c017aa6..a429a98fc 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -16,13 +16,13 @@ struct GetPosition #elif (AMREX_SPACEDIM == 2) static constexpr RType m_snan = std::numeric_limits::quiet_NaN(); #endif - GetPosition (const WarpXParIter& a_pti) noexcept + GetPosition (const WarpXParIter& a_pti, int a_offset = 0) noexcept { const auto& aos = a_pti.GetArrayOfStructs(); - m_structs = aos().dataPtr(); + m_structs = aos().dataPtr() + a_offset; #if (defined WARPX_DIM_RZ) const auto& soa = a_pti.GetStructOfArrays(); - m_theta = soa.GetRealData(PIdx::theta).dataPtr(); + m_theta = soa.GetRealData(PIdx::theta).dataPtr() + a_offset; #elif (AMREX_SPACEDIM == 2) static constexpr RType m_snan = std::numeric_limits::quiet_NaN(); #endif @@ -57,13 +57,13 @@ struct SetPosition #if (defined WARPX_DIM_RZ) RType* AMREX_RESTRICT m_theta; #endif - SetPosition (WarpXParIter& a_pti) noexcept + SetPosition (WarpXParIter& a_pti, int a_offset = 0) noexcept { auto& aos = a_pti.GetArrayOfStructs(); - m_structs = aos().dataPtr(); + m_structs = aos().dataPtr() + a_offset; #if (defined WARPX_DIM_RZ) auto& soa = a_pti.GetStructOfArrays(); - m_theta = soa.GetRealData(PIdx::theta).dataPtr(); + m_theta = soa.GetRealData(PIdx::theta).dataPtr() + a_offset; #endif } diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 564355887..ad9f0f07e 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -8,7 +8,7 @@ /** \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(WarpXParticleContainer::ParticleType& p, +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 ) { @@ -18,11 +18,11 @@ void UpdatePosition(WarpXParticleContainer::ParticleType& p, // Compute inverse Lorentz factor const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2); // Update positions over one time step - p.pos(0) += ux * inv_gamma * dt; + x += ux * inv_gamma * dt; #if (AMREX_SPACEDIM == 3) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D - p.pos(1) += uy * inv_gamma * dt; + y += uy * inv_gamma * dt; #endif - p.pos(2) += uz * inv_gamma * dt; + 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 fd9597b3f..325e4b748 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -12,7 +12,7 @@ */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdatePositionPhoton( - WarpXParticleContainer::ParticleType& p, + 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 ) { @@ -20,12 +20,11 @@ void UpdatePositionPhoton( const amrex::Real c_over_umod = PhysConst::c/std::sqrt(ux*ux + uy*uy + uz*uz); // Update positions over one time step - p.pos(0) += ux * c_over_umod * dt; + x += ux * c_over_umod * dt; #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D - p.pos(1) += uy * c_over_umod * dt; + y += uy * c_over_umod * dt; #endif - p.pos(2) += uz * c_over_umod * dt; - + z += uz * c_over_umod * dt; } #endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 4d4aedb0e..fc41fd3d2 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -366,11 +366,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,7 +396,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, 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; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 170770c58..779771360 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -286,9 +286,8 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // CPU, tiling: deposit into local_jx // (same for jx and jz) - const auto& aos = pti.GetArrayOfStructs(); - const ParticleType* AMREX_RESTRICT pstruct = aos().dataPtr() + offset; - + const auto get_position = GetPosition(pti, offset); + // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow const Dim3 lo = lbound(tilebox); @@ -298,19 +297,19 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>( - pstruct, wp.dataPtr() + offset, uxp.dataPtr() + offset, + get_position, 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>( - pstruct, wp.dataPtr() + offset, uxp.dataPtr() + offset, + get_position, 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>( - pstruct, wp.dataPtr() + offset, uxp.dataPtr() + offset, + get_position, 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); @@ -318,19 +317,19 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } else { if (WarpX::nox == 1){ doDepositionShapeN<1>( - pstruct, wp.dataPtr() + offset, uxp.dataPtr() + offset, + get_position, 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>( - pstruct, wp.dataPtr() + offset, uxp.dataPtr() + offset, + get_position, 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>( - pstruct, wp.dataPtr() + offset, uxp.dataPtr() + offset, + get_position, 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); @@ -424,8 +423,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // GPU, no tiling: deposit directly in rho // CPU, tiling: deposit into local_rho - const auto& aos = pti.GetArrayOfStructs(); - const ParticleType* AMREX_RESTRICT pstruct = aos().dataPtr() + offset; + const auto get_position = GetPosition(pti, offset); // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -435,13 +433,13 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(pstruct, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<1>(get_position, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(pstruct, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<2>(get_position, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(pstruct, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<3>(get_position, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -734,28 +732,33 @@ 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 get_position = GetPosition(pti); + auto set_position = SetPosition(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 + Real x, y, z; #ifndef WARPX_DIM_RZ - UpdatePosition( p, ux[i], uy[i], uz[i], dt); + get_position(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt); + set_position(i, x, y, z); #else // For WARPX_DIM_RZ, the particles are still pushed in 3D Cartesian - ParticleReal x, y, z; // Temporary variables + ParticleType& p = pstruct[i]; // Particle object that gets updated GetCartesianPositionFromCylindrical( x, y, z, p, theta[i] ); - UpdatePosition( p, ux[i], uy[i], uz[i], dt); + UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); SetCylindricalPositionFromCartesian( p, theta[i], x, y, z ); #endif } -- cgit v1.2.3