From 6f23f443585034a413059074830d7ffd3052520e Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Wed, 17 Jul 2019 09:28:50 -0700 Subject: Convert PushP and RigidPC to c++ --- .../Particles/RigidInjectedParticleContainer.cpp | 164 ++++++++++++++------- 1 file changed, 113 insertions(+), 51 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 9bd4cb4fc..a0aadb57d 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -10,6 +10,9 @@ #include #include #include +#include +#include +#include using namespace amrex; @@ -233,6 +236,20 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Cuda::ManagedDeviceVector xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; + Real* AMREX_RESTRICT x = xp.dataPtr(); + Real* AMREX_RESTRICT y = yp.dataPtr(); + Real* AMREX_RESTRICT z = zp.dataPtr(); + Real* AMREX_RESTRICT gi = giv.dataPtr(); + Real* AMREX_RESTRICT uxpp = uxp.dataPtr(); + Real* AMREX_RESTRICT uypp = uyp.dataPtr(); + Real* AMREX_RESTRICT uzpp = uzp.dataPtr(); + Real* AMREX_RESTRICT Expp = Exp.dataPtr(); + Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); + Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); + Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + if (!done_injecting_lev) { xp_save = xp; yp_save = yp; @@ -240,33 +257,34 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, uxp_save = uxp; uyp_save = uyp; uzp_save = uzp; + + Real* AMREX_RESTRICT xp_savep = xp_save.dataPtr(); + Real* AMREX_RESTRICT yp_savep = yp_save.dataPtr(); + Real* AMREX_RESTRICT zp_savep = zp_save.dataPtr(); + Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles // should by advanced a fraction of a time step instead. // Scaling the fields is much easier and may be good enough. - for (int i=0 ; i < zp.size() ; i++) { - const Real dtscale = dt - (zinject_plane_lev_previous - zp[i])/(vzbeam_ave_boosted + WarpX::beta_boost*PhysConst::c); + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + const Real dtscale = dt - (zinject_plane_lev_previous - z[i])/(vzbeam_ave_boosted + WarpX::beta_boost*PhysConst::c); if (0. < dtscale && dtscale < dt) { - Exp[i] *= dtscale; - Eyp[i] *= dtscale; - Ezp[i] *= dtscale; - Bxp[i] *= dtscale; - Byp[i] *= dtscale; - Bzp[i] *= dtscale; + Expp[i] *= dtscale; + Eypp[i] *= dtscale; + Ezpp[i] *= dtscale; + Bxpp[i] *= dtscale; + Bypp[i] *= dtscale; + Bzpp[i] *= dtscale; } } + ); } - warpx_particle_pusher(&np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - giv.dataPtr(), - Exp.dataPtr(), Eyp.dataPtr(), Ezp.dataPtr(), - Bxp.dataPtr(), Byp.dataPtr(), Bzp.dataPtr(), - &this->charge, &this->mass, &dt, - &WarpX::particle_pusher_algo); + PhysicalParticleContainer::PushPX(pti, xp, yp, zp, giv, dt); if (!done_injecting_lev) { #ifdef _OPENMP @@ -274,25 +292,35 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, #else const int tid = 0; #endif + + Real* AMREX_RESTRICT xp_savep = xp_save.dataPtr(); + Real* AMREX_RESTRICT yp_savep = yp_save.dataPtr(); + Real* AMREX_RESTRICT zp_savep = zp_save.dataPtr(); + Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. - for (int i=0 ; i < zp.size() ; i++) { - if (zp[i] <= zinject_plane_lev) { - uxp[i] = uxp_save[i]; - uyp[i] = uyp_save[i]; - uzp[i] = uzp_save[i]; - giv[i] = 1./std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/(PhysConst::c*PhysConst::c)); - xp[i] = xp_save[i]; - yp[i] = yp_save[i]; + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + if (z[i] <= zinject_plane_lev) { + uxpp[i] = uxp_savep[i]; + uypp[i] = uyp_savep[i]; + uzpp[i] = uzp_savep[i]; + gi[i] = 1./std::sqrt(1. + (uxpp[i]*uxpp[i] + uypp[i]*uypp[i] + uzpp[i]*uzpp[i])/(PhysConst::c*PhysConst::c)); + x[i] = xp_savep[i]; + y[i] = yp_savep[i]; if (rigid_advance) { - zp[i] = zp_save[i] + dt*vzbeam_ave_boosted; + z[i] = zp_savep[i] + dt*vzbeam_ave_boosted; } else { - zp[i] = zp_save[i] + dt*uzp[i]*giv[i]; + z[i] = zp_savep[i] + dt*uzpp[i]*gi[i]; } done_injecting_temp[tid] = 0; } } + ); } } @@ -343,6 +371,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) { + BL_PROFILE("RigidInjectedParticleContainer::PushP"); + if (do_not_push) return; const std::array& dx = WarpX::CellSize(lev); @@ -351,8 +381,11 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, #pragma omp parallel #endif { - Cuda::ManagedDeviceVector xp, yp, zp, giv; - +#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(); @@ -386,24 +419,24 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, Byp.assign(np,WarpX::B_external[1]); Bzp.assign(np,WarpX::B_external[2]); - giv.resize(np); + m_giv[thread_num].resize(np); // // copy data from particle container to temp arrays // - pti.GetPosition(xp, yp, zp); + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); const int* ixyzmin_grid = box.loVect(); const int ll4symtry = false; - const int l_lower_order_in_v = true; long lvect_fieldgathe = 64; + warpx_geteb_energy_conserving( &np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), + m_xp[thread_num].dataPtr(), + m_yp[thread_num].dataPtr(), + m_zp[thread_num].dataPtr(), Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), ixyzmin_grid, @@ -416,7 +449,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &l_lower_order_in_v, &WarpX::do_nodal, + &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); // Save the position and momenta, making copies @@ -424,27 +457,56 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, auto uyp_save = uyp; auto uzp_save = uzp; - warpx_particle_pusher_momenta(&np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - giv.dataPtr(), - Exp.dataPtr(), Eyp.dataPtr(), Ezp.dataPtr(), - Bxp.dataPtr(), Byp.dataPtr(), Bzp.dataPtr(), - &this->charge, &this->mass, &dt, - &WarpX::particle_pusher_algo); + // This wraps the momentum advance so that inheritors can modify the call. + // Extract pointers to the different particle quantities + Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + Real* AMREX_RESTRICT gi = m_giv[thread_num].dataPtr(); + Real* AMREX_RESTRICT uxpp = uxp.dataPtr(); + Real* AMREX_RESTRICT uypp = uyp.dataPtr(); + Real* AMREX_RESTRICT uzpp = uzp.dataPtr(); + Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + Real* AMREX_RESTRICT Expp = Exp.dataPtr(); + Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); + Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); + Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + + // Loop over the particles and update their momentum + const Real q = this->charge; + const Real m = this-> mass; + if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumBoris( uxpp[i], uypp[i], uzpp[i], gi[i], + Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumVay( uxpp[i], uypp[i], uzpp[i], gi[i], + Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + }; // Undo the push for particles not injected yet. // It is assumed that PushP will only be called on the first and last steps // and that no particles will cross zinject_plane. - for (int i=0 ; i < zp.size() ; i++) { + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { if (zp[i] <= zinject_plane_levels[lev]) { - uxp[i] = uxp_save[i]; - uyp[i] = uyp_save[i]; - uzp[i] = uzp_save[i]; + uxpp[i] = uxp_save[i]; + uypp[i] = uyp_save[i]; + uzpp[i] = uzp_save[i]; } } + ); } } -- cgit v1.2.3 From ea3d09e05401e2b5aa53668eb637f4e25d9fd15e Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Wed, 17 Jul 2019 09:40:27 -0700 Subject: Removed warpx_particle_pusher and other cleanup --- Source/FortranInterface/WarpX_picsar.F90 | 112 --------------------- .../Particles/RigidInjectedParticleContainer.cpp | 12 +-- 2 files changed, 6 insertions(+), 118 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 33f85c633..40b9f4fc9 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -417,118 +417,6 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n end subroutine warpx_current_deposition_rz_volume_scaling - ! _________________________________________________________________ - !> - !> @brief - !> Main subroutine for the particle pusher (velocity and position) - !> - !> @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] ex,ey,ez particle electric fields in each direction - !> @param[in] bx,by,bz particle magnetic fields in each direction - !> @param[in] q charge - !> @param[in] m masse - !> @param[in] dt time step - !> @param[in] particle_pusher_algo Particle pusher algorithm - subroutine warpx_particle_pusher(np,xp,yp,zp,uxp,uyp,uzp, & - gaminv,& - ex,ey,ez,bx,by,bz,q,m,dt, & - particle_pusher_algo) & - bind(C, name="warpx_particle_pusher") - - 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) :: ex(np),ey(np),ez(np) - REAL(amrex_real),INTENT(IN) :: bx(np),by(np),bz(np) - REAL(amrex_real),INTENT(IN) :: q,m,dt - INTEGER(c_long), INTENT(IN) :: particle_pusher_algo - - SELECT CASE (particle_pusher_algo) - - !! Vay pusher -- Full push - CASE (1_c_long) - CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) - - CALL pxr_ebcancelpush3d(np,uxp,uyp,uzp,gaminv, & - ex,ey,ez, & - bx,by,bz,q,m,dt,0_c_long) - CASE DEFAULT - - ! Momentum pusher in a single loop - CALL pxr_boris_push_u_3d(np,uxp,uyp,uzp,& - gaminv, & - ex,ey,ez, & - bx,by,bz, & - q,m,dt) - - END SELECT - - !!!! --- 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 - - - ! _________________________________________________________________ - !> - !> @brief - !> Main subroutine for the particle pusher (velocity) - !> - !> @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] ex,ey,ez particle electric fields in each direction - !> @param[in] bx,by,bz particle magnetic fields in each direction - !> @param[in] q charge - !> @param[in] m masse - !> @param[in] dt time step - !> @param[in] particle_pusher_algo Particle pusher algorithm - subroutine warpx_particle_pusher_momenta(np,xp,yp,zp,uxp,uyp,uzp, & - gaminv,& - ex,ey,ez,bx,by,bz,q,m,dt, & - particle_pusher_algo) & - bind(C, name="warpx_particle_pusher_momenta") - - INTEGER(c_long), INTENT(IN) :: np - REAL(amrex_real),INTENT(INOUT) :: gaminv(np) - REAL(amrex_real),INTENT(IN) :: xp(np),yp(np),zp(np) - REAL(amrex_real),INTENT(INOUT) :: uxp(np),uyp(np),uzp(np) - REAL(amrex_real),INTENT(IN) :: ex(np),ey(np),ez(np) - REAL(amrex_real),INTENT(IN) :: bx(np),by(np),bz(np) - REAL(amrex_real),INTENT(IN) :: q,m,dt - INTEGER(c_long), INTENT(IN) :: particle_pusher_algo - - SELECT CASE (particle_pusher_algo) - - !! Vay pusher -- Full push - CASE (1_c_long) - CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) - - CALL pxr_ebcancelpush3d(np,uxp,uyp,uzp,gaminv, & - ex,ey,ez, & - bx,by,bz,q,m,dt,0_c_long) - CASE DEFAULT - - ! Momentum pusher in a single loop - CALL pxr_boris_push_u_3d(np,uxp,uyp,uzp,& - gaminv, & - ex,ey,ez, & - bx,by,bz, & - q,m,dt) - - END SELECT - - end subroutine warpx_particle_pusher_momenta - ! _________________________________________________________________ !> !> @brief diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a0aadb57d..7e4434278 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -207,7 +207,7 @@ RigidInjectedParticleContainer::BoostandRemapParticles() void RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector& xp, + Cuda::ManagedDeviceVector& xp, Cuda::ManagedDeviceVector& yp, Cuda::ManagedDeviceVector& zp, Cuda::ManagedDeviceVector& giv, @@ -218,8 +218,8 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, { copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); } - - // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. + + // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; @@ -352,9 +352,9 @@ RigidInjectedParticleContainer::Evolve (int lev, done_injecting_lev = done_injecting[lev]; PhysicalParticleContainer::Evolve (lev, - Ex, Ey, Ez, - Bx, By, Bz, - jx, jy, jz, + Ex, Ey, Ez, + Bx, By, Bz, + jx, jy, jz, cjx, cjy, cjz, rho, crho, cEx, cEy, cEz, -- cgit v1.2.3 From 2b2d1c40ac995951661acfe03c8f13993ece741e Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 22 Jul 2019 16:08:51 -0700 Subject: Momentum push conversion clean up for clarity and optimization --- Source/Particles/PhysicalParticleContainer.cpp | 54 ++++---- .../Particles/RigidInjectedParticleContainer.cpp | 146 ++++++++++----------- 2 files changed, 97 insertions(+), 103 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 132224db2..5ee2ef5b6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1751,19 +1751,19 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, // 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 - Real* AMREX_RESTRICT x = xp.dataPtr(); - Real* AMREX_RESTRICT y = yp.dataPtr(); - Real* AMREX_RESTRICT z = zp.dataPtr(); - Real* AMREX_RESTRICT gi = giv.dataPtr(); - Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - Real* AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); - Real* AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); - Real* AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); - Real* AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); - Real* AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); - Real* AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); + Real* const AMREX_RESTRICT x = xp.dataPtr(); + Real* const AMREX_RESTRICT y = yp.dataPtr(); + Real* const AMREX_RESTRICT z = zp.dataPtr(); + Real* const AMREX_RESTRICT gi = giv.dataPtr(); + Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const Real* AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const Real* AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const Real* AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const Real* AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const Real* AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const Real* AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { @@ -1794,7 +1794,6 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, } else { amrex::Abort("Unknown particle pusher"); }; - } void @@ -1823,9 +1822,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt, auto& attribs = pti.GetAttribs(); - auto& uxp = attribs[PIdx::ux]; - auto& uyp = attribs[PIdx::uy]; - auto& uzp = attribs[PIdx::uz]; auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; auto& Ezp = attribs[PIdx::Ez]; @@ -1885,16 +1881,16 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - Real* AMREX_RESTRICT gi = m_giv[thread_num].dataPtr(); - Real* AMREX_RESTRICT uxpp = uxp.dataPtr(); - Real* AMREX_RESTRICT uypp = uyp.dataPtr(); - Real* AMREX_RESTRICT uzpp = uzp.dataPtr(); - Real* AMREX_RESTRICT Expp = Exp.dataPtr(); - Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); - Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); - Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + Real* const AMREX_RESTRICT gi = m_giv[thread_num].dataPtr(); + Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const Real* AMREX_RESTRICT Expp = Exp.dataPtr(); + const Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); + const Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; @@ -1902,14 +1898,14 @@ PhysicalParticleContainer::PushP (int lev, Real dt, if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBoris( uxpp[i], uypp[i], uzpp[i], gi[i], + UpdateMomentumBoris( ux[i], uy[i], uz[i], gi[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumVay( uxpp[i], uypp[i], uzpp[i], gi[i], + UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 7e4434278..09df32db4 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -214,56 +214,40 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real dt) { - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); - } - // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; - auto& Exp = attribs[PIdx::Ex]; - auto& Eyp = attribs[PIdx::Ey]; - auto& Ezp = attribs[PIdx::Ez]; - auto& Bxp = attribs[PIdx::Bx]; - auto& Byp = attribs[PIdx::By]; - auto& Bzp = attribs[PIdx::Bz]; - const long np = pti.numParticles(); // Save the position and momenta, making copies Cuda::ManagedDeviceVector xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; - Real* AMREX_RESTRICT x = xp.dataPtr(); - Real* AMREX_RESTRICT y = yp.dataPtr(); - Real* AMREX_RESTRICT z = zp.dataPtr(); - Real* AMREX_RESTRICT gi = giv.dataPtr(); - Real* AMREX_RESTRICT uxpp = uxp.dataPtr(); - Real* AMREX_RESTRICT uypp = uyp.dataPtr(); - Real* AMREX_RESTRICT uzpp = uzp.dataPtr(); - Real* AMREX_RESTRICT Expp = Exp.dataPtr(); - Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); - Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); - Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + Real* const AMREX_RESTRICT x = xp.dataPtr(); + Real* const AMREX_RESTRICT y = yp.dataPtr(); + Real* const AMREX_RESTRICT z = zp.dataPtr(); + Real* const AMREX_RESTRICT gi = giv.dataPtr(); + Real* const AMREX_RESTRICT ux = uxp.dataPtr(); + Real* const AMREX_RESTRICT uy = uyp.dataPtr(); + Real* const AMREX_RESTRICT uz = uzp.dataPtr(); + Real* const AMREX_RESTRICT Exp = attribs[PIdx::Ex].dataPtr(); + Real* const AMREX_RESTRICT Eyp = attribs[PIdx::Ey].dataPtr(); + Real* const AMREX_RESTRICT Ezp = attribs[PIdx::Ez].dataPtr(); + Real* const AMREX_RESTRICT Bxp = attribs[PIdx::Bx].dataPtr(); + Real* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr(); + Real* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; - - Real* AMREX_RESTRICT xp_savep = xp_save.dataPtr(); - Real* AMREX_RESTRICT yp_savep = yp_save.dataPtr(); - Real* AMREX_RESTRICT zp_savep = zp_save.dataPtr(); - Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); - Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); - Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { + // If the old values are not already saved, create copies here. + xp_save = xp; + yp_save = yp; + zp_save = zp; + uxp_save = uxp; + uyp_save = uyp; + uzp_save = uzp; + } // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -273,12 +257,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, [=] AMREX_GPU_DEVICE (long i) { const Real dtscale = dt - (zinject_plane_lev_previous - z[i])/(vzbeam_ave_boosted + WarpX::beta_boost*PhysConst::c); if (0. < dtscale && dtscale < dt) { - Expp[i] *= dtscale; - Eypp[i] *= dtscale; - Ezpp[i] *= dtscale; - Bxpp[i] *= dtscale; - Bypp[i] *= dtscale; - Bzpp[i] *= dtscale; + Exp[i] *= dtscale; + Eyp[i] *= dtscale; + Ezp[i] *= dtscale; + Bxp[i] *= dtscale; + Byp[i] *= dtscale; + Bzp[i] *= dtscale; } } ); @@ -293,36 +277,50 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, const int tid = 0; #endif - Real* AMREX_RESTRICT xp_savep = xp_save.dataPtr(); - Real* AMREX_RESTRICT yp_savep = yp_save.dataPtr(); - Real* AMREX_RESTRICT zp_savep = zp_save.dataPtr(); - Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); - Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); - Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + Real* AMREX_RESTRICT x_save; + Real* AMREX_RESTRICT y_save; + Real* AMREX_RESTRICT z_save; + Real* AMREX_RESTRICT ux_save; + Real* AMREX_RESTRICT uy_save; + Real* AMREX_RESTRICT uz_save; + if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { + x_save = xp_save.dataPtr(); + y_save = yp_save.dataPtr(); + z_save = zp_save.dataPtr(); + ux_save = uxp_save.dataPtr(); + uy_save = uyp_save.dataPtr(); + uz_save = uzp_save.dataPtr(); + } else { + x_save = pti.GetAttribs(particle_comps["xold"]).dataPtr(); + y_save = pti.GetAttribs(particle_comps["yold"]).dataPtr(); + z_save = pti.GetAttribs(particle_comps["zold"]).dataPtr(); + ux_save = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); + uy_save = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); + uz_save = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); + } // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { if (z[i] <= zinject_plane_lev) { - uxpp[i] = uxp_savep[i]; - uypp[i] = uyp_savep[i]; - uzpp[i] = uzp_savep[i]; - gi[i] = 1./std::sqrt(1. + (uxpp[i]*uxpp[i] + uypp[i]*uypp[i] + uzpp[i]*uzpp[i])/(PhysConst::c*PhysConst::c)); - x[i] = xp_savep[i]; - y[i] = yp_savep[i]; + ux[i] = ux_save[i]; + uy[i] = uy_save[i]; + uz[i] = uz_save[i]; + gi[i] = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])/(PhysConst::c*PhysConst::c)); + x[i] = x_save[i]; + y[i] = y_save[i]; if (rigid_advance) { - z[i] = zp_savep[i] + dt*vzbeam_ave_boosted; + z[i] = z_save[i] + dt*vzbeam_ave_boosted; } else { - z[i] = zp_savep[i] + dt*uzpp[i]*gi[i]; + z[i] = z_save[i] + dt*uz[i]*gi[i]; } done_injecting_temp[tid] = 0; } } ); } - } void @@ -459,24 +457,24 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); - Real* AMREX_RESTRICT gi = m_giv[thread_num].dataPtr(); - Real* AMREX_RESTRICT uxpp = uxp.dataPtr(); - Real* AMREX_RESTRICT uypp = uyp.dataPtr(); - Real* AMREX_RESTRICT uzpp = uzp.dataPtr(); - Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); - Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); - Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); - Real* AMREX_RESTRICT Expp = Exp.dataPtr(); - Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); - Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); - Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + const Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + Real* const AMREX_RESTRICT gi = m_giv[thread_num].dataPtr(); + Real* const AMREX_RESTRICT uxpp = uxp.dataPtr(); + Real* const AMREX_RESTRICT uypp = uyp.dataPtr(); + Real* const AMREX_RESTRICT uzpp = uzp.dataPtr(); + const Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); + const Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); + const Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + const Real* AMREX_RESTRICT Expp = Exp.dataPtr(); + const Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); + const Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; - const Real m = this-> mass; + const Real m = this->mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { -- cgit v1.2.3 From df262ec96f182446274d8d3ab6c627d94b910f26 Mon Sep 17 00:00:00 2001 From: grote Date: Wed, 24 Jul 2019 10:01:31 -0700 Subject: Bug fixes for rigid injection in conversion of momentum push --- .../Particles/RigidInjectedParticleContainer.cpp | 34 +++++++++++----------- 1 file changed, 17 insertions(+), 17 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 09df32db4..dc594f5a1 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -253,9 +253,10 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, // This only approximates what should be happening. The particles // should by advanced a fraction of a time step instead. // Scaling the fields is much easier and may be good enough. + const Real v_boost = WarpX::beta_boost*PhysConst::c; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const Real dtscale = dt - (zinject_plane_lev_previous - z[i])/(vzbeam_ave_boosted + WarpX::beta_boost*PhysConst::c); + const Real dtscale = dt - (zinject_plane_lev_previous - z[i])/(vzbeam_ave_boosted + v_boost); if (0. < dtscale && dtscale < dt) { Exp[i] *= dtscale; Eyp[i] *= dtscale; @@ -316,7 +317,6 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, else { z[i] = z_save[i] + dt*uz[i]*gi[i]; } - done_injecting_temp[tid] = 0; } } ); @@ -340,13 +340,13 @@ RigidInjectedParticleContainer::Evolve (int lev, zinject_plane_levels[lev] -= dt*WarpX::beta_boost*PhysConst::c; zinject_plane_lev = zinject_plane_levels[lev]; - // Setup check of whether more particles need to be injected -#ifdef _OPENMP - const int nthreads = omp_get_max_threads(); -#else - const int nthreads = 1; -#endif - done_injecting_temp.assign(nthreads, 1); // We do not use bool because vector is special. + // Set the done injecting flag whan the inject plane moves out of the + // simulation domain. + // It is much easier to do this check, rather than checking if all of the + // particles have crossed the inject plane. + const Real* plo = Geom(lev).ProbLo(); + const Real* phi = Geom(lev).ProbHi(); + done_injecting[lev] = (zinject_plane_levels[lev] < plo[2] || zinject_plane_levels[lev] > phi[2]); done_injecting_lev = done_injecting[lev]; PhysicalParticleContainer::Evolve (lev, @@ -358,10 +358,6 @@ RigidInjectedParticleContainer::Evolve (int lev, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); - - // Check if all done_injecting_temp are still true. - done_injecting[lev] = std::all_of(done_injecting_temp.begin(), done_injecting_temp.end(), - [](int i) -> bool { return i; }); } void @@ -496,12 +492,16 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // Undo the push for particles not injected yet. // It is assumed that PushP will only be called on the first and last steps // and that no particles will cross zinject_plane. + const Real* const AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + const Real* const AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + const Real* const AMREX_RESTRICT uz_save = uzp_save.dataPtr(); + const Real zz = zinject_plane_levels[lev]; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - if (zp[i] <= zinject_plane_levels[lev]) { - uxpp[i] = uxp_save[i]; - uypp[i] = uyp_save[i]; - uzpp[i] = uzp_save[i]; + if (zp[i] <= zz) { + uxpp[i] = ux_save[i]; + uypp[i] = uy_save[i]; + uzpp[i] = uz_save[i]; } } ); -- cgit v1.2.3 From 9105ca79ac7a115cc3348d937462044a7c55bf2e Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Wed, 24 Jul 2019 10:16:30 -0700 Subject: For momentum push conversion, added more consts --- Source/Particles/PhysicalParticleContainer.cpp | 24 +++++++++++----------- .../Particles/RigidInjectedParticleContainer.cpp | 20 +++++++++--------- 2 files changed, 22 insertions(+), 22 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 5ee2ef5b6..b913ed5cd 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1758,12 +1758,12 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - const Real* AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); - const Real* AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); - const Real* AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); - const Real* AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); - const Real* AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); - const Real* AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); + const Real* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const Real* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const Real* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const Real* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const Real* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const Real* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { @@ -1885,12 +1885,12 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - const Real* AMREX_RESTRICT Expp = Exp.dataPtr(); - const Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); - const Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - const Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - const Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); - const Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + const Real* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const Real* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const Real* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const Real* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const Real* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const Real* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index dc594f5a1..cad82d6d4 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -453,20 +453,20 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - const Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + const Real* const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); Real* const AMREX_RESTRICT gi = m_giv[thread_num].dataPtr(); Real* const AMREX_RESTRICT uxpp = uxp.dataPtr(); Real* const AMREX_RESTRICT uypp = uyp.dataPtr(); Real* const AMREX_RESTRICT uzpp = uzp.dataPtr(); - const Real* AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); - const Real* AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); - const Real* AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); - const Real* AMREX_RESTRICT Expp = Exp.dataPtr(); - const Real* AMREX_RESTRICT Eypp = Eyp.dataPtr(); - const Real* AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - const Real* AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - const Real* AMREX_RESTRICT Bypp = Byp.dataPtr(); - const Real* AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + const Real* const AMREX_RESTRICT uxp_savep = uxp_save.dataPtr(); + const Real* const AMREX_RESTRICT uyp_savep = uyp_save.dataPtr(); + const Real* const AMREX_RESTRICT uzp_savep = uzp_save.dataPtr(); + const Real* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const Real* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const Real* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const Real* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const Real* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const Real* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; -- cgit v1.2.3 From 30616270c7d7d0ca8a10912d624b445d6d249dc4 Mon Sep 17 00:00:00 2001 From: grote Date: Thu, 25 Jul 2019 09:51:15 -0700 Subject: For momentum push conversion, create copies of member variables for gpu --- Source/Particles/RigidInjectedParticleContainer.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index cad82d6d4..79396e6af 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -254,9 +254,11 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, // should by advanced a fraction of a time step instead. // Scaling the fields is much easier and may be good enough. const Real v_boost = WarpX::beta_boost*PhysConst::c; + const Real z_plane_previous = zinject_plane_lev_previous; + const Real vz_ave_boosted = vzbeam_ave_boosted; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const Real dtscale = dt - (zinject_plane_lev_previous - z[i])/(vzbeam_ave_boosted + v_boost); + const Real dtscale = dt - (z_plane_previous - z[i])/(vz_ave_boosted + v_boost); if (0. < dtscale && dtscale < dt) { Exp[i] *= dtscale; Eyp[i] *= dtscale; @@ -272,11 +274,6 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, PhysicalParticleContainer::PushPX(pti, xp, yp, zp, giv, dt); if (!done_injecting_lev) { -#ifdef _OPENMP - const int tid = omp_get_thread_num(); -#else - const int tid = 0; -#endif Real* AMREX_RESTRICT x_save; Real* AMREX_RESTRICT y_save; @@ -302,17 +299,21 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. + const Real z_plane_lev = zinject_plane_lev; + const Real vz_ave_boosted = vzbeam_ave_boosted; + const bool rigid = rigid_advance; + const Real inv_csq = 1./(PhysConst::c*PhysConst::c); amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - if (z[i] <= zinject_plane_lev) { + if (z[i] <= z_plane_lev) { ux[i] = ux_save[i]; uy[i] = uy_save[i]; uz[i] = uz_save[i]; - gi[i] = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])/(PhysConst::c*PhysConst::c)); + gi[i] = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_csq); x[i] = x_save[i]; y[i] = y_save[i]; - if (rigid_advance) { - z[i] = z_save[i] + dt*vzbeam_ave_boosted; + if (rigid) { + z[i] = z_save[i] + dt*vz_ave_boosted; } else { z[i] = z_save[i] + dt*uz[i]*gi[i]; -- cgit v1.2.3