aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles/WarpXParticleContainer.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles/WarpXParticleContainer.cpp')
-rw-r--r--Source/Particles/WarpXParticleContainer.cpp167
1 files changed, 89 insertions, 78 deletions
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index 2edd3c636..47d57294d 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -7,6 +7,10 @@
#include <WarpX_f.H>
#include <WarpX.H>
+// Import low-level single-particle kernels
+#include <GetAndSetPosition.H>
+#include <UpdatePosition.H>
+
using namespace amrex;
int WarpXParticleContainer::do_not_push = 0;
@@ -78,7 +82,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies)
particle_comps["theta"] = PIdx::theta;
#endif
- if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles)
+ if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
{
particle_comps["xold"] = PIdx::nattribs;
particle_comps["yold"] = PIdx::nattribs+1;
@@ -86,9 +90,9 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies)
particle_comps["uxold"] = PIdx::nattribs+3;
particle_comps["uyold"] = PIdx::nattribs+4;
particle_comps["uzold"] = PIdx::nattribs+5;
-
+
}
-
+
// Initialize temporary local arrays for charge/current deposition
int num_threads = 1;
#ifdef _OPENMP
@@ -174,7 +178,7 @@ WarpXParticleContainer::AddNParticles (int lev,
int n, const Real* x, const Real* y, const Real* z,
const Real* vx, const Real* vy, const Real* vz,
int nattr, const Real* attr, int uniqueparticles, int id)
-{
+{
BL_ASSERT(nattr == 1);
const Real* weight = attr;
@@ -230,15 +234,15 @@ WarpXParticleContainer::AddNParticles (int lev,
#endif
p.pos(1) = z[i];
#endif
-
- if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles)
+
+ if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
{
- auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0);
- particle_tile.push_back_real(particle_comps["xold"], x[i]);
- particle_tile.push_back_real(particle_comps["yold"], y[i]);
- particle_tile.push_back_real(particle_comps["zold"], z[i]);
+ auto& ptile = DefineAndReturnParticleTile(0, 0, 0);
+ ptile.push_back_real(particle_comps["xold"], x[i]);
+ ptile.push_back_real(particle_comps["yold"], y[i]);
+ ptile.push_back_real(particle_comps["zold"], z[i]);
}
-
+
particle_tile.push_back(p);
}
@@ -249,14 +253,14 @@ WarpXParticleContainer::AddNParticles (int lev,
particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend);
particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend);
- if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles)
+ if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
{
- auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0);
- particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend);
- particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend);
- particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend);
+ auto& ptile = DefineAndReturnParticleTile(0, 0, 0);
+ ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend);
+ ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend);
+ ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend);
}
-
+
for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp)
{
#ifdef WARPX_RZ
@@ -327,7 +331,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti,
jx_ptr = local_jx[thread_num].dataPtr();
jy_ptr = local_jy[thread_num].dataPtr();
jz_ptr = local_jz[thread_num].dataPtr();
-
+
local_jx[thread_num].setVal(0.0);
local_jy[thread_num].setVal(0.0);
local_jz[thread_num].setVal(0.0);
@@ -426,7 +430,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti,
BL_PROFILE_VAR_STOP(blp_pxr_cd);
-#ifndef AMREX_USE_GPU
+#ifndef AMREX_USE_GPU
BL_PROFILE_VAR_START(blp_accumulate);
jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1);
@@ -477,7 +481,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti,
auto jyntot = local_jy[thread_num].length();
auto jzntot = local_jz[thread_num].length();
#endif
-
+
long ncrse = np - np_current;
BL_PROFILE_VAR_START(blp_pxr_cd);
if (j_is_nodal) {
@@ -608,7 +612,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp,
const std::array<Real, 3>& xyzmin = xyzmin_tile;
#ifdef AMREX_USE_GPU
- data_ptr = (*rhomf)[pti].dataPtr();
+ data_ptr = (*rhomf)[pti].dataPtr(icomp);
auto rholen = (*rhomf)[pti].length();
#else
tile_box.grow(ngRho);
@@ -641,9 +645,14 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp,
&ngRho, &ngRho, &ngRho,
&WarpX::nox,&WarpX::noy,&WarpX::noz,
&lvect, &WarpX::charge_deposition_algo);
+#ifdef WARPX_RZ
+ warpx_charge_deposition_rz_volume_scaling(
+ data_ptr, &ngRho, rholen.getVect(),
+ &xyzmin[0], &dx[0]);
+#endif
BL_PROFILE_VAR_STOP(blp_pxr_chd);
-#ifndef AMREX_USE_GPU
+#ifndef AMREX_USE_GPU
BL_PROFILE_VAR_START(blp_accumulate);
(*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1);
@@ -660,7 +669,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp,
const std::array<Real,3>& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1);
#ifdef AMREX_USE_GPU
- data_ptr = (*crhomf)[pti].dataPtr();
+ data_ptr = (*crhomf)[pti].dataPtr(icomp);
auto rholen = (*crhomf)[pti].length();
#else
tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector());
@@ -696,9 +705,14 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp,
&ngRho, &ngRho, &ngRho,
&WarpX::nox,&WarpX::noy,&WarpX::noz,
&lvect, &WarpX::charge_deposition_algo);
+#ifdef WARPX_RZ
+ warpx_charge_deposition_rz_volume_scaling(
+ data_ptr, &ngRho, rholen.getVect(),
+ &cxyzmin_tile[0], &cdx[0]);
+#endif
BL_PROFILE_VAR_STOP(blp_pxr_chd);
-#ifndef AMREX_USE_GPU
+#ifndef AMREX_USE_GPU
BL_PROFILE_VAR_START(blp_accumulate);
(*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1);
@@ -723,7 +737,6 @@ WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho,
const auto& gm = m_gdb->Geom(lev);
const auto& ba = m_gdb->ParticleBoxArray(lev);
- const auto& dm = m_gdb->DistributionMap(lev);
const Real* dx = gm.CellSize();
const Real* plo = gm.ProbLo();
@@ -793,36 +806,36 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local)
#ifdef _OPENMP
#pragma omp parallel
-#endif
{
+#endif
Cuda::ManagedDeviceVector<Real> xp, yp, zp;
- FArrayBox local_rho;
+#ifdef _OPENMP
+ FArrayBox rho_loc;
+#endif
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
- const Box& box = pti.validbox();
-
auto& wp = pti.GetAttribs(PIdx::w);
const long np = pti.numParticles();
pti.GetPosition(xp, yp, zp);
- const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev);
- const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev);
-
// Data on the grid
Real* data_ptr;
FArrayBox& rhofab = (*rho)[pti];
#ifdef _OPENMP
+ const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev);
Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector());
const std::array<Real, 3>& xyzmin = xyzmin_tile;
tile_box.grow(ng);
- local_rho.resize(tile_box);
- local_rho = 0.0;
- data_ptr = local_rho.dataPtr();
- auto rholen = local_rho.length();
+ rho_loc.resize(tile_box);
+ rho_loc = 0.0;
+ data_ptr = rho_loc.dataPtr();
+ auto rholen = rho_loc.length();
#else
+ const Box& box = pti.validbox();
+ const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev);
const std::array<Real, 3>& xyzmin = xyzmin_grid;
data_ptr = rhofab.dataPtr();
auto rholen = rhofab.length();
@@ -852,12 +865,17 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local)
&dx[0], &dx[1], &dx[2], &nx, &ny, &nz,
&nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz,
&lvect, &WarpX::charge_deposition_algo);
+#ifdef WARPX_RZ
+ long ngRho = WarpX::nox;
+ warpx_charge_deposition_rz_volume_scaling(
+ data_ptr, &ngRho, rholen.getVect(),
+ &xyzmin[0], &dx[0]);
+#endif
#ifdef _OPENMP
- rhofab.atomicAdd(local_rho);
-#endif
+ rhofab.atomicAdd(rho_loc);
}
-
+#endif
}
if (!local) rho->SumBoundary(gm.periodicity());
@@ -1007,8 +1025,6 @@ void
WarpXParticleContainer::PushX (int lev, Real dt)
{
BL_PROFILE("WPC::PushX()");
- BL_PROFILE_VAR_NS("WPC::PushX::Copy", blp_copy);
- BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pxr_pp);
if (do_not_push) return;
@@ -1018,47 +1034,42 @@ WarpXParticleContainer::PushX (int lev, Real dt)
#pragma omp parallel
#endif
{
- Cuda::ManagedDeviceVector<Real> xp, yp, zp, giv;
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
Real wt = amrex::second();
- auto& attribs = pti.GetAttribs();
-
- auto& uxp = attribs[PIdx::ux];
- auto& uyp = attribs[PIdx::uy];
- auto& uzp = attribs[PIdx::uz];
-
- const long np = pti.numParticles();
-
- giv.resize(np);
-
- //
- // copy data from particle container to temp arrays
- //
- BL_PROFILE_VAR_START(blp_copy);
- pti.GetPosition(xp, yp, zp);
- BL_PROFILE_VAR_STOP(blp_copy);
-
//
// Particle Push
//
- BL_PROFILE_VAR_START(blp_pxr_pp);
- warpx_particle_pusher_positions(&np,
- xp.dataPtr(),
- yp.dataPtr(),
- zp.dataPtr(),
- uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(),
- giv.dataPtr(), &dt);
- BL_PROFILE_VAR_STOP(blp_pxr_pp);
-
- //
- // copy particle data back
- //
- BL_PROFILE_VAR_START(blp_copy);
- pti.SetPosition(xp, yp, zp);
- BL_PROFILE_VAR_STOP(blp_copy);
+ // 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]);
+ // - momenta are stored as a struct of array, in `attribs`
+ auto& attribs = pti.GetAttribs();
+ Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+ Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+ Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+#ifdef WARPX_RZ
+ Real* 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; // Temporary variables
+#ifndef WARPX_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_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
+ }
+ );
if (cost) {
const Box& tbx = pti.tilebox();
@@ -1074,15 +1085,15 @@ WarpXParticleContainer::PushX (int lev, Real dt)
}
// This function is called in Redistribute, just after locate
-void
-WarpXParticleContainer::particlePostLocate(ParticleType& p,
+void
+WarpXParticleContainer::particlePostLocate(ParticleType& p,
const ParticleLocData& pld,
const int lev)
{
// Tag particle if goes to higher level.
// It will be split later in the loop
- if (pld.m_lev == lev+1
- and p.m_idata.id != NoSplitParticleID
+ if (pld.m_lev == lev+1
+ and p.m_idata.id != NoSplitParticleID
and p.m_idata.id >= 0)
{
p.m_idata.id = DoSplitParticleID;
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-15 20:43:35 +0000