Merge branch 'dev' into elementary_process

1 files changed, 111 insertions, 152 deletions

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,