aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/FortranInterface/WarpX_f.H6
-rw-r--r--Source/FortranInterface/WarpX_picsar.F9032
-rw-r--r--Source/Particles/WarpXParticleContainer.cpp93
3 files changed, 41 insertions, 90 deletions
diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H
index 3d92b7651..4457e34d8 100644
--- a/Source/FortranInterface/WarpX_f.H
+++ b/Source/FortranInterface/WarpX_f.H
@@ -176,12 +176,6 @@ extern "C"
const amrex::Real* charge, const amrex::Real* mass, const amrex::Real* dt,
const long* particle_pusher_algo);
- // Particle pusher (position)
-
- void warpx_particle_pusher_positions(const long* np,
- amrex::Real* xp, amrex::Real* yp, amrex::Real* zp,
- amrex::Real* uxp, amrex::Real* uyp, amrex::Real* uzp, amrex::Real* gaminv,
- const amrex::Real* dt);
// Laser pusher
diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90
index 6d6a2e411..a4cc99926 100644
--- a/Source/FortranInterface/WarpX_picsar.F90
+++ b/Source/FortranInterface/WarpX_picsar.F90
@@ -522,38 +522,6 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n
! _________________________________________________________________
!>
!> @brief
- !> Main subroutine for the particle pusher of positions
- !>
- !> @param[in] np number of super-particles
- !> @param[in] xp,yp,zp particle position arrays
- !> @param[in] uxp,uyp,uzp normalized momentum in each direction
- !> @param[in] gaminv particle Lorentz factors
- !> @param[in] dt time step
- !> @param[in] particle_pusher_algo Particle pusher algorithm
- subroutine warpx_particle_pusher_positions(np,xp,yp,zp,uxp,uyp,uzp, &
- gaminv,dt) &
- bind(C, name="warpx_particle_pusher_positions")
-
- INTEGER(c_long), INTENT(IN) :: np
- REAL(amrex_real),INTENT(INOUT) :: gaminv(np)
- REAL(amrex_real),INTENT(INOUT) :: xp(np),yp(np),zp(np)
- REAL(amrex_real),INTENT(INOUT) :: uxp(np),uyp(np),uzp(np)
- REAL(amrex_real),INTENT(IN) :: dt
-
- CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv)
-
- !!!! --- push particle species positions a time step
-#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ)
- CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt)
-#elif (AMREX_SPACEDIM == 2)
- CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt)
-#endif
-
- end subroutine warpx_particle_pusher_positions
-
- ! _________________________________________________________________
- !>
- !> @brief
!> Main subroutine for the Maxwell solver (E field)
!>
!> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds
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-09-04 16:56:24 +0000