From eee063ceca2707824b63ac0966acf265c2b67457 Mon Sep 17 00:00:00 2001 From: Diana Amorim Date: Tue, 9 Jul 2019 16:22:09 -0700 Subject: Changed warpx_copy_attribs function from F90 to CPP inside PhysicalParticleContainer class. --- Source/Particles/RigidInjectedParticleContainer.cpp | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 2a3e8dd0d..ca02d1458 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -211,6 +211,11 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real dt) { + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + } + // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); auto& uxp = attribs[PIdx::ux]; @@ -224,21 +229,6 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& xpold = pti.GetAttribs(particle_comps["xold"]); - auto& ypold = pti.GetAttribs(particle_comps["yold"]); - auto& zpold = pti.GetAttribs(particle_comps["zold"]); - auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); - auto& uypold = pti.GetAttribs(particle_comps["uyold"]); - auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - } - // Save the position and momenta, making copies Cuda::ManagedDeviceVector xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; -- cgit v1.2.3 From fb25dde6d6770f01bdfee71e10137c5fe8a8e035 Mon Sep 17 00:00:00 2001 From: Diana Amorim Date: Wed, 10 Jul 2019 16:56:58 -0700 Subject: Changed name of function to copy_attribs() --- Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- Source/Particles/RigidInjectedParticleContainer.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 99fc0f8da..d55764682 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -78,7 +78,7 @@ public: const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - void warpx_copy_attribs(WarpXParIter& pti,const amrex::Real* xp, + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c45f73f0f..addde5e37 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1731,7 +1731,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); } // The following attributes should be included in CPP version of warpx_particle_pusher @@ -1861,7 +1861,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::warpx_copy_attribs(WarpXParIter& pti,const Real* xp, +void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const Real* yp, const Real* zp) { diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index ca02d1458..9bd4cb4fc 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -213,7 +213,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); } // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. -- cgit v1.2.3 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++ --- Source/Particles/PhysicalParticleContainer.cpp | 45 ++++-- .../Particles/RigidInjectedParticleContainer.cpp | 164 ++++++++++++++------- 2 files changed, 147 insertions(+), 62 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3a512d9dc..4b4747d0e 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1748,7 +1748,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, Real dt) { - // 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(); // Extract pointers to the different particle quantities Real* AMREX_RESTRICT x = xp.dataPtr(); @@ -1883,16 +1883,39 @@ PhysicalParticleContainer::PushP (int lev, Real dt, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); - warpx_particle_pusher_momenta(&np, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - m_giv[thread_num].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 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(); + + // 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"); + }; } } } 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 From 2a45d27a5aed1de8531f07e0756c14533e09bc1a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 1 Aug 2019 15:32:28 -0700 Subject: remove unused variables --- Source/Particles/RigidInjectedParticleContainer.cpp | 3 --- 1 file changed, 3 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 79396e6af..84c14dd7c 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -459,9 +459,6 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, Real* const AMREX_RESTRICT uxpp = uxp.dataPtr(); Real* const AMREX_RESTRICT uypp = uyp.dataPtr(); Real* const AMREX_RESTRICT uzpp = uzp.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(); -- cgit v1.2.3 From 43ef6e29a1eeba61b9fdb95efb1bb1e301ba1b6d Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 2 Aug 2019 17:25:35 -0700 Subject: PushP now uses FIeldGather --- Source/Particles/PhysicalParticleContainer.cpp | 7 +++++++ Source/Particles/RigidInjectedParticleContainer.cpp | 7 +++++++ 2 files changed, 14 insertions(+) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a34fb69e2..89ea3c9f3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1647,6 +1647,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const int ll4symtry = false; long lvect_fieldgathe = 64; +#ifdef WARPX_RZ warpx_geteb_energy_conserving( &np, m_xp[thread_num].dataPtr(), @@ -1666,6 +1667,12 @@ PhysicalParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bzfab), &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); +#else + 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); +#endif // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 84c14dd7c..df809a5f0 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -427,6 +427,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const int ll4symtry = false; long lvect_fieldgathe = 64; +#ifdef WARPX_RZ warpx_geteb_energy_conserving( &np, m_xp[thread_num].dataPtr(), @@ -446,6 +447,12 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bzfab), &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); +#else + 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); +#endif // Save the position and momenta, making copies auto uxp_save = uxp; -- cgit v1.2.3 From 69a94bfb22c09e96ffb31f3bea95a69f7f949319 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 2 Aug 2019 17:44:26 -0700 Subject: Minor fix in RigidParticleContainer::PushP --- Source/Particles/RigidInjectedParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index df809a5f0..83556e69f 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -421,13 +421,13 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); +#ifdef WARPX_RZ const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); const int* ixyzmin_grid = box.loVect(); const int ll4symtry = false; long lvect_fieldgathe = 64; -#ifdef WARPX_RZ warpx_geteb_energy_conserving( &np, m_xp[thread_num].dataPtr(), -- cgit v1.2.3 From 411d953482a42070484fe780deeff4f7593c4b9e Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 5 Aug 2019 09:23:16 -0700 Subject: Removed warpx_geteb_energy_conserving --- Source/FortranInterface/WarpX_f.H | 28 ----- Source/FortranInterface/WarpX_picsar.F90 | 81 --------------- Source/Particles/PhysicalParticleContainer.cpp | 113 +-------------------- .../Particles/RigidInjectedParticleContainer.cpp | 28 ----- 4 files changed, 4 insertions(+), 246 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 0440148eb..689029059 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -106,34 +106,6 @@ extern "C" const long* nox, const long* noy,const long* noz, const int* l_nodal, const long* lvect, const long* current_depo_algo); - // Current deposition finalize for RZ - void warpx_current_deposition_rz_volume_scaling( - amrex::Real* jx, const long* jx_ng, const int* jx_ntot, - amrex::Real* jy, const long* jy_ng, const int* jy_ntot, - amrex::Real* jz, const long* jz_ng, const int* jz_ntot, - const amrex::Real* rmin, - const amrex::Real* dr); - - // Field gathering - - void warpx_geteb_energy_conserving(const long* np, - const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, - amrex::Real* exp, amrex::Real* eyp, amrex::Real* ezp, - amrex::Real* bxp, amrex::Real* byp, amrex::Real* bzp, - const int* ixyzmin, - const amrex::Real* xmin, const amrex::Real* ymin, const amrex::Real* zmin, - const amrex::Real* dx, const amrex::Real* dy, const amrex::Real* dz, - const long* nox, const long* noy, const long* noz, - const amrex::Real* exg, const int* exg_lo, const int* exg_hi, - const amrex::Real* eyg, const int* eyg_lo, const int* eyg_hi, - const amrex::Real* ezg, const int* ezg_lo, const int* ezg_hi, - const amrex::Real* bxg, const int* bxg_lo, const int* bxg_hi, - const amrex::Real* byg, const int* byg_lo, const int* byg_hi, - const amrex::Real* bzg, const int* bzg_lo, const int* bzg_hi, - const int* ll4symtry, const int* l_lower_order_in_v, - const int* l_nodal, const long* lvect, - const long* field_gathe_algo); - // Particle pusher (velocity and position) void warpx_particle_pusher(const long* np, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 14bd79ad4..e65c30e42 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -1,19 +1,16 @@ #if (AMREX_SPACEDIM == 3) -#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb3d_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic #elif (AMREX_SPACEDIM == 2) #ifdef WARPX_RZ -#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2drz_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz #define WRPX_PXR_RZ_VOLUME_SCALING_RHO apply_rz_volume_scaling_rho #else -#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2dxz_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic_2d #endif @@ -52,84 +49,6 @@ module warpx_to_pxr_module contains - ! _________________________________________________________________ - !> - !> @brief - !> Main subroutine for the field gathering process - !> - !> @param[in] np number of particles - !> @param[in] xp,yp,zp particle position arrays - !> @param[in] ex,ey,ez particle electric fields in each direction - !> @param[in] bx,by,bz particle magnetic fields in each direction - !> @param[in] ixyzmin tile grid minimum index - !> @param[in] xmin,ymin,zmin tile grid minimum position - !> @param[in] dx,dy,dz space discretization steps - !> @param[in] xyzmin grid minimum position - !> @param[in] dxyz space discretization steps - !> @param[in] nox,noy,noz interpolation order - !> @param[in] exg,eyg,ezg electric field grid arrays - !> @param[in] bxg,byg,bzg electric field grid arrays - !> @param[in] lvect vector length - !> - subroutine warpx_geteb_energy_conserving(np,xp,yp,zp, & - ex,ey,ez,bx,by,bz,ixyzmin,xmin,ymin,zmin,dx,dy,dz,nox,noy,noz, & - exg,exg_lo,exg_hi,eyg,eyg_lo,eyg_hi,ezg,ezg_lo,ezg_hi, & - bxg,bxg_lo,bxg_hi,byg,byg_lo,byg_hi,bzg,bzg_lo,bzg_hi, & - ll4symtry,l_lower_order_in_v, l_nodal,& - lvect,field_gathe_algo) & - bind(C, name="warpx_geteb_energy_conserving") - - integer, intent(in) :: exg_lo(AMREX_SPACEDIM), eyg_lo(AMREX_SPACEDIM), ezg_lo(AMREX_SPACEDIM), & - bxg_lo(AMREX_SPACEDIM), byg_lo(AMREX_SPACEDIM), bzg_lo(AMREX_SPACEDIM) - integer, intent(in) :: exg_hi(AMREX_SPACEDIM), eyg_hi(AMREX_SPACEDIM), ezg_hi(AMREX_SPACEDIM), & - bxg_hi(AMREX_SPACEDIM), byg_hi(AMREX_SPACEDIM), bzg_hi(AMREX_SPACEDIM) - integer, intent(in) :: ixyzmin(AMREX_SPACEDIM) - real(amrex_real), intent(in) :: xmin,ymin,zmin,dx,dy,dz - integer(c_long), intent(in) :: field_gathe_algo - integer(c_long), intent(in) :: np,nox,noy,noz - integer(c_int), intent(in) :: ll4symtry,l_lower_order_in_v, l_nodal - integer(c_long),intent(in) :: lvect - real(amrex_real), intent(in), dimension(np) :: xp,yp,zp - real(amrex_real), intent(out), dimension(np) :: ex,ey,ez,bx,by,bz - real(amrex_real),intent(in):: exg(*), eyg(*), ezg(*), bxg(*), byg(*), bzg(*) - logical(pxr_logical) :: pxr_ll4symtry, pxr_l_lower_order_in_v, pxr_l_nodal - - ! Compute the number of valid cells and guard cells - integer(c_long) :: exg_nvalid(AMREX_SPACEDIM), eyg_nvalid(AMREX_SPACEDIM), ezg_nvalid(AMREX_SPACEDIM), & - bxg_nvalid(AMREX_SPACEDIM), byg_nvalid(AMREX_SPACEDIM), bzg_nvalid(AMREX_SPACEDIM), & - exg_nguards(AMREX_SPACEDIM), eyg_nguards(AMREX_SPACEDIM), ezg_nguards(AMREX_SPACEDIM), & - bxg_nguards(AMREX_SPACEDIM), byg_nguards(AMREX_SPACEDIM), bzg_nguards(AMREX_SPACEDIM) - - pxr_ll4symtry = ll4symtry .eq. 1 - pxr_l_lower_order_in_v = l_lower_order_in_v .eq. 1 - pxr_l_nodal = l_nodal .eq. 1 - - exg_nguards = ixyzmin - exg_lo - eyg_nguards = ixyzmin - eyg_lo - ezg_nguards = ixyzmin - ezg_lo - bxg_nguards = ixyzmin - bxg_lo - byg_nguards = ixyzmin - byg_lo - bzg_nguards = ixyzmin - bzg_lo - exg_nvalid = exg_lo + exg_hi - 2_c_long*ixyzmin + 1_c_long - eyg_nvalid = eyg_lo + eyg_hi - 2_c_long*ixyzmin + 1_c_long - ezg_nvalid = ezg_lo + ezg_hi - 2_c_long*ixyzmin + 1_c_long - bxg_nvalid = bxg_lo + bxg_hi - 2_c_long*ixyzmin + 1_c_long - byg_nvalid = byg_lo + byg_hi - 2_c_long*ixyzmin + 1_c_long - bzg_nvalid = bzg_lo + bzg_hi - 2_c_long*ixyzmin + 1_c_long - - CALL WRPX_PXR_GETEB_ENERGY_CONSERVING(np,xp,yp,zp, & - ex,ey,ez,bx,by,bz,xmin,ymin,zmin,dx,dy,dz,nox,noy,noz, & - exg,exg_nguards,exg_nvalid,& - eyg,eyg_nguards,eyg_nvalid,& - ezg,ezg_nguards,ezg_nvalid,& - bxg,bxg_nguards,bxg_nvalid,& - byg,byg_nguards,byg_nvalid,& - bzg,bzg_nguards,bzg_nvalid,& - pxr_ll4symtry, pxr_l_lower_order_in_v, pxr_l_nodal, & - lvect, field_gathe_algo ) - - end subroutine warpx_geteb_energy_conserving - ! _________________________________________________________________ !> !> @brief diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7d63bd8e7..08f8a77b4 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -893,40 +893,13 @@ PhysicalParticleContainer::FieldGather (int lev, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - const std::array& xyzmin = WarpX::LowerCorner(box, lev); - const int* ixyzmin = box.loVect(); - // // Field Gather // -#ifdef WARPX_RZ - const int ll4symtry = false; - long lvect_fieldgathe = 64; - warpx_geteb_energy_conserving( - &np, - 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, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else 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); -#endif if (cost) { const Box& tbx = pti.tilebox(); @@ -1195,35 +1168,9 @@ PhysicalParticleContainer::Evolve (int lev, // Field Gather of Aux Data (i.e., the full solution) // BL_PROFILE_VAR_START(blp_pxr_fg); -#ifdef WARPX_RZ - const int ll4symtry = false; - long lvect_fieldgathe = 64; - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - const int* ixyzmin_grid = box.loVect(); - warpx_geteb_energy_conserving( - &np_gather, - 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, - &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(*exfab), - BL_TO_FORTRAN_ANYD(*eyfab), - BL_TO_FORTRAN_ANYD(*ezfab), - BL_TO_FORTRAN_ANYD(*bxfab), - BL_TO_FORTRAN_ANYD(*byfab), - BL_TO_FORTRAN_ANYD(*bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else 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); -#endif if (np_gather < np) { @@ -1293,29 +1240,6 @@ PhysicalParticleContainer::Evolve (int lev, } // Field gather for particles in gather buffers -#ifdef WARPX_RZ - - long ncrse = np - nfine_gather; - warpx_geteb_energy_conserving( - &ncrse, - m_xp[thread_num].dataPtr()+nfine_gather, - m_yp[thread_num].dataPtr()+nfine_gather, - m_zp[thread_num].dataPtr()+nfine_gather, - Exp.dataPtr()+nfine_gather, Eyp.dataPtr()+nfine_gather, Ezp.dataPtr()+nfine_gather, - Bxp.dataPtr()+nfine_gather, Byp.dataPtr()+nfine_gather, Bzp.dataPtr()+nfine_gather, - cixyzmin_grid, - &cxyzmin_grid[0], &cxyzmin_grid[1], &cxyzmin_grid[2], - &cdx[0], &cdx[1], &cdx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(*cexfab), - BL_TO_FORTRAN_ANYD(*ceyfab), - BL_TO_FORTRAN_ANYD(*cezfab), - BL_TO_FORTRAN_ANYD(*cbxfab), - BL_TO_FORTRAN_ANYD(*cbyfab), - BL_TO_FORTRAN_ANYD(*cbzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, @@ -1323,7 +1247,6 @@ PhysicalParticleContainer::Evolve (int lev, cEx->nGrow(), e_is_nodal, nfine_gather, np-nfine_gather, thread_num, lev, lev-1); -#endif } BL_PROFILE_VAR_STOP(blp_pxr_fg); @@ -1657,38 +1580,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); -#ifdef WARPX_RZ - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - const int* ixyzmin_grid = box.loVect(); - - const int ll4symtry = false; - long lvect_fieldgathe = 64; - - warpx_geteb_energy_conserving( - &np, - 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, - &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else - 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); -#endif + 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); // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 83556e69f..f049fdb7c 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -421,38 +421,10 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); -#ifdef WARPX_RZ - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - const int* ixyzmin_grid = box.loVect(); - - const int ll4symtry = false; - long lvect_fieldgathe = 64; - - warpx_geteb_energy_conserving( - &np, - 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, - &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else 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); -#endif // Save the position and momenta, making copies auto uxp_save = uxp; -- cgit v1.2.3