aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles')
-rwxr-xr-xSource/Particles/Deposition/ChargeDeposition.H18
-rw-r--r--Source/Particles/Deposition/CurrentDeposition.H56
-rw-r--r--Source/Particles/Gather/FieldGather.H23
-rw-r--r--Source/Particles/PhotonParticleContainer.H3
-rw-r--r--Source/Particles/PhotonParticleContainer.cpp24
-rw-r--r--Source/Particles/PhysicalParticleContainer.H15
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp263
-rw-r--r--Source/Particles/Pusher/GetAndSetPosition.H143
-rw-r--r--Source/Particles/Pusher/UpdatePosition.H8
-rw-r--r--Source/Particles/Pusher/UpdatePositionPhoton.H1
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.H6
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.cpp197
-rw-r--r--Source/Particles/WarpXParticleContainer.H10
-rw-r--r--Source/Particles/WarpXParticleContainer.cpp109
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);
}
);