diff options
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 263 |
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, |