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.cpp93
1 files changed, 41 insertions, 52 deletions
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index 80f3882a0..1abbd747d 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -86,9 +86,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 +174,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 +230,15 @@ WarpXParticleContainer::AddNParticles (int lev,
#endif
p.pos(1) = z[i];
#endif
-
+
if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles)
{
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]);
+ particle_tile.push_back_real(particle_comps["zold"], z[i]);
}
-
+
particle_tile.push_back(p);
}
@@ -254,9 +254,9 @@ WarpXParticleContainer::AddNParticles (int lev,
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);
+ particle_tile.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 +327,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 +426,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 +477,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) {
@@ -648,7 +648,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp,
#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);
@@ -708,7 +708,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp,
#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);
@@ -1023,8 +1023,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;
@@ -1034,47 +1032,38 @@ 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();
+ // Loop over the particles and update the position
+ const long np = pti.numParticles();
+ const Real inv_c2 = 1./(PhysConst::c*PhysConst::c);
+ amrex::ParallelFor( np,
+ [=] AMREX_GPU_DEVICE (long i) {
+ // Compute inverse Lorentz factor
+ const Real inv_gamma = 1./std::sqrt(
+ 1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_c2);
+ // Update positions over one time step
+ pstructs[i].pos(0) += ux[i] * inv_gamma * dt;
+#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D
+ pstructs[i].pos(1) += uy[i] * inv_gamma * dt;
+#endif
+ pstructs[i].pos(2) += uz[i] * inv_gamma * dt;
+ }
+ );
if (cost) {
const Box& tbx = pti.tilebox();
@@ -1090,15 +1079,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 01:56:06 +0000