diff options
Diffstat (limited to 'Source/Particles')
23 files changed, 545 insertions, 435 deletions
diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index b9210e67c..eec407a2b 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -18,10 +18,10 @@ * /param q : species charge. */ template <int depos_order> -void doChargeDepositionShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - const amrex::Real * const wp, +void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + const amrex::ParticleReal * const wp, const int * const ion_lev, const amrex::Array4<amrex::Real>& rho_arr, const long np_to_depose, diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 7a96dab9a..6da0f1155 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -24,13 +24,13 @@ * /param q : species charge. */ template <int depos_order> -void doDepositionShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - const amrex::Real * const wp, - const amrex::Real * const uxp, - const amrex::Real * const uyp, - const amrex::Real * const uzp, +void doDepositionShapeN(const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + const amrex::ParticleReal * const wp, + const amrex::ParticleReal * const uxp, + const amrex::ParticleReal * const uyp, + const amrex::ParticleReal * const uzp, const int * const ion_lev, const amrex::Array4<amrex::Real>& jx_arr, const amrex::Array4<amrex::Real>& jy_arr, @@ -189,13 +189,13 @@ void doDepositionShapeN(const amrex::Real * const xp, * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order> -void doEsirkepovDepositionShapeN (const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - const amrex::Real * const wp, - const amrex::Real * const uxp, - const amrex::Real * const uyp, - const amrex::Real * const uzp, +void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + const amrex::ParticleReal * const wp, + const amrex::ParticleReal * const uxp, + const amrex::ParticleReal * const uyp, + const amrex::ParticleReal * const uzp, const int * ion_lev, const amrex::Array4<amrex::Real>& Jx_arr, const amrex::Array4<amrex::Real>& Jy_arr, @@ -397,7 +397,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djr_cmplx = 2.*sdxi*xy_mid; + const Complex djr_cmplx = amrex::Real(2.)*sdxi*xy_mid; amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real()); amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag()); xy_mid = xy_mid*xy_mid0; @@ -418,7 +418,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -2.*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode* + const Complex djt_cmplx = -amrex::Real(2.)*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode* (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); @@ -438,7 +438,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djz_cmplx = 2.*sdzk*xy_mid; + const Complex djz_cmplx = amrex::Real(2.)*sdzk*xy_mid; amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real()); amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag()); xy_mid = xy_mid*xy_mid0; diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 6727b0aa9..c5ec6fb5b 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -19,12 +19,12 @@ * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order, int lower_in_v> -void doGatherShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - amrex::Real * const Exp, amrex::Real * const Eyp, - amrex::Real * const Ezp, amrex::Real * const Bxp, - amrex::Real * const Byp, amrex::Real * const Bzp, +void doGatherShapeN(const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, + amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, + amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, const amrex::Array4<const amrex::Real>& ex_arr, const amrex::Array4<const amrex::Real>& ey_arr, const amrex::Array4<const amrex::Real>& ez_arr, diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 95f36cf65..c9d520873 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -18,6 +18,7 @@ CEXE_headers += ShapeFactors.H include $(WARPX_HOME)/Source/Particles/Pusher/Make.package include $(WARPX_HOME)/Source/Particles/Deposition/Make.package include $(WARPX_HOME)/Source/Particles/Gather/Make.package +include $(WARPX_HOME)/Source/Particles/Sorting/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index bb795465e..715c97b99 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -545,12 +545,12 @@ namespace WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); // --- source SoA particle data auto& soa_source = ptile_source.GetStructOfArrays(); - GpuArray<Real*,PIdx::nattribs> attribs_source; + GpuArray<ParticleReal*,PIdx::nattribs> attribs_source; for (int ia = 0; ia < PIdx::nattribs; ++ia) { attribs_source[ia] = soa_source.GetRealData(ia).data(); } // --- source runtime attribs - GpuArray<Real*,3> runtime_uold_source; + GpuArray<ParticleReal*,3> runtime_uold_source; // Prepare arrays for boosted frame diagnostics. runtime_uold_source[0] = soa_source.GetRealData(PIdx::ux).data(); runtime_uold_source[1] = soa_source.GetRealData(PIdx::uy).data(); @@ -590,13 +590,13 @@ namespace WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; // --- product SoA particle data auto& soa_product = ptile_product.GetStructOfArrays(); - GpuArray<Real*,PIdx::nattribs> attribs_product; + GpuArray<ParticleReal*,PIdx::nattribs> attribs_product; for (int ia = 0; ia < PIdx::nattribs; ++ia) { // First element is the first newly-created product particle attribs_product[ia] = soa_product.GetRealData(ia).data() + np_product_old; } // --- product runtime attribs - GpuArray<Real*,6> runtime_attribs_product; + GpuArray<ParticleReal*,6> runtime_attribs_product; bool do_boosted_product = WarpX::do_boosted_frame_diagnostic && pc_product->DoBoostedFrameDiags(); if (do_boosted_product) { diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 3ac304bdc..a6ffd1d76 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -41,9 +41,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 9de441e5c..4a75ec9f3 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -34,27 +34,27 @@ void PhotonParticleContainer::InitData() void PhotonParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<Real>& xp, - Cuda::ManagedDeviceVector<Real>& yp, - Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<ParticleReal>& xp, + Cuda::ManagedDeviceVector<ParticleReal>& yp, + Cuda::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - 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 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* 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(); + ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); + ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); + ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ace1ec7f8..b558323a3 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -1,12 +1,13 @@ #ifndef WARPX_PhysicalParticleContainer_H_ #define WARPX_PhysicalParticleContainer_H_ -#include <map> +#include <PlasmaInjector.H> +#include <WarpXParticleContainer.H> +#include <NCIGodfreyFilter.H> #include <AMReX_IArrayBox.H> -#include <PlasmaInjector.H> -#include <WarpXParticleContainer.H> +#include <map> class PhysicalParticleContainer : public WarpXParticleContainer @@ -87,9 +88,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full); virtual void PushP (int lev, amrex::Real dt, @@ -100,8 +101,21 @@ public: const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, - const amrex::Real* yp, const amrex::Real* zp); + void PartitionParticlesInBuffers( + long& nfine_current, + long& nfine_gather, + long const np, + WarpXParIter& pti, + int const lev, + amrex::iMultiFab const* current_masks, + amrex::iMultiFab const* gather_masks, + RealVector& uxp, + RealVector& uyp, + RealVector& uzp, + RealVector& wp ); + + void copy_attribs(WarpXParIter& pti,const amrex::ParticleReal* xp, + const amrex::ParticleReal* yp, const amrex::ParticleReal* zp); virtual void PostRestart () final {} @@ -131,6 +145,33 @@ public: virtual void ConvertUnits (ConvertDirection convert_dir) override; +/** + * \brief Apply NCI Godfrey filter to all components of E and B before gather + * \param lev MR level + * \param box box onto which the filter is applied + * \param exeli safeguard to avoid destructing arrays between ParIter iterations on GPU + * \param filtered_Ex Array containing filtered value + * \param Ex Field array before filtering (not modified) + * \param ex_ptr pointer to the array used for field gather. + * + * The NCI Godfrey filter is applied on Ex, the result is stored in filtered_Ex + * and the pointer is modified (before this function is called, it points to Ex + * and after this function is called, it points to Ex_filtered) + */ + void applyNCIFilter ( + int lev, const amrex::Box& box, + amrex::Elixir& exeli, amrex::Elixir& eyeli, amrex::Elixir& ezeli, + amrex::Elixir& bxeli, amrex::Elixir& byeli, amrex::Elixir& bzeli, + amrex::FArrayBox& filtered_Ex, amrex::FArrayBox& filtered_Ey, + amrex::FArrayBox& filtered_Ez, amrex::FArrayBox& filtered_Bx, + amrex::FArrayBox& filtered_By, amrex::FArrayBox& filtered_Bz, + const amrex::FArrayBox& Ex, const amrex::FArrayBox& Ey, + const amrex::FArrayBox& Ez, const amrex::FArrayBox& Bx, + const amrex::FArrayBox& By, const amrex::FArrayBox& Bz, + amrex::FArrayBox const * & exfab, amrex::FArrayBox const * & eyfab, + amrex::FArrayBox const * & ezfab, amrex::FArrayBox const * & bxfab, + amrex::FArrayBox const * & byfab, amrex::FArrayBox const * & bzfab); + protected: std::string species_name; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b1e83d652..02dee1967 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -31,8 +31,11 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("boost_adjust_transverse_positions", boost_adjust_transverse_positions); pp.query("do_backward_propagation", do_backward_propagation); + + // Initialize splitting pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); + pp.query("do_continuous_injection", do_continuous_injection); // Whether to plot back-transformed (lab-frame) diagnostics // for this species. @@ -196,7 +199,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, std::array<Real, 3> u, Real weight) { - std::array<Real,PIdx::nattribs> attribs; + std::array<ParticleReal,PIdx::nattribs> attribs; attribs.fill(0.0); // update attribs with input arguments @@ -361,13 +364,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int dir=0; dir<AMREX_SPACEDIM; dir++) { if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0.) * dx[dir]); + overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]); } else { no_overlap = true; break; } if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0.) * dx[dir]); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]); } else { no_overlap = true; break; } @@ -437,7 +440,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size; auto& soa = particle_tile.GetStructOfArrays(); - GpuArray<Real*,PIdx::nattribs> pa; + GpuArray<ParticleReal*,PIdx::nattribs> pa; for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } @@ -948,7 +951,6 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy); BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); - BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); @@ -960,7 +962,6 @@ PhysicalParticleContainer::Evolve (int lev, BL_ASSERT(OnSameGrids(lev,jx)); MultiFab* cost = WarpX::getCosts(lev); - const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev); const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev); @@ -991,10 +992,6 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox filtered_Ex, filtered_Ey, filtered_Ez; FArrayBox filtered_Bx, filtered_By, filtered_Bz; - std::vector<bool> inexflag; - Vector<long> pid; - RealVector tmp; - ParticleVector particle_tmp; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -1026,56 +1023,18 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* bzfab = &(Bz[pti]); Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; + if (WarpX::use_fdtd_nci_corr) { -#if (AMREX_SPACEDIM == 2) - const Box& tbox = amrex::grow(pti.tilebox(),{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noz)}); -#else - const Box& tbox = amrex::grow(pti.tilebox(),{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noy), - static_cast<int>(WarpX::noz)}); -#endif - - // Filter Ex (Both 2D and 3D) - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - // Safeguard for GPU - exeli = filtered_Ex.elixir(); - // Apply filter on Ex, result stored in filtered_Ex - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); - // Update exfab reference - exfab = &filtered_Ex; - - // Filter Ez - filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); - ezfab = &filtered_Ez; - - // Filter By - filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); - byfab = &filtered_By; -#if (AMREX_SPACEDIM == 3) - // Filter Ey - filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); - eyfab = &filtered_Ey; - - // Filter Bx - filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); - bxfab = &filtered_Bx; - - // Filter Bz - filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); - bzfab = &filtered_Bz; -#endif + // Filter arrays Ex[pti], store the result in + // filtered_Ex and update pointer exfab so that it + // points to filtered_Ex (and do the same for all + // components of E and B). + applyNCIFilter(lev, pti.tilebox(), exeli, eyeli, ezeli, bxeli, byeli, bzeli, + filtered_Ex, filtered_Ey, filtered_Ez, + filtered_Bx, filtered_By, filtered_Bz, + Ex[pti], Ey[pti], Ez[pti], Bx[pti], By[pti], Bz[pti], + exfab, eyfab, ezfab, bxfab, byfab, bzfab); } Exp.assign(np,0.0); @@ -1085,99 +1044,21 @@ PhysicalParticleContainer::Evolve (int lev, Byp.assign(np,WarpX::B_external[1]); Bzp.assign(np,WarpX::B_external[2]); - long nfine_current = np; //! number of particles depositing to fine grid - long nfine_gather = np; //! number of particles gathering from fine grid - if (has_buffer && !do_not_push) - { - BL_PROFILE_VAR_START(blp_partition); - inexflag.resize(np); - auto& aos = pti.GetArrayOfStructs(); - // We need to partition the large buffer first - iMultiFab const* bmasks = (WarpX::n_field_gather_buffer >= WarpX::n_current_deposition_buffer) ? - gather_masks : current_masks; - int i = 0; - const auto& msk = (*bmasks)[pti]; - for (const auto& p : aos) { - const IntVect& iv = Index(p, lev); - inexflag[i++] = msk(iv); - } - - pid.resize(np); - std::iota(pid.begin(), pid.end(), 0L); - - auto sep = std::stable_partition(pid.begin(), pid.end(), - [&inexflag](long id) { return inexflag[id]; }); - - if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer) { - nfine_current = nfine_gather = std::distance(pid.begin(), sep); - } else if (sep != pid.end()) { - int n_buf; - if (bmasks == gather_masks) { - nfine_gather = std::distance(pid.begin(), sep); - bmasks = current_masks; - n_buf = WarpX::n_current_deposition_buffer; - } else { - nfine_current = std::distance(pid.begin(), sep); - bmasks = gather_masks; - n_buf = WarpX::n_field_gather_buffer; - } - if (n_buf > 0) - { - const auto& msk2 = (*bmasks)[pti]; - for (auto it = sep; it != pid.end(); ++it) { - const long id = *it; - const IntVect& iv = Index(aos[id], lev); - inexflag[id] = msk2(iv); - } - - auto sep2 = std::stable_partition(sep, pid.end(), - [&inexflag](long id) {return inexflag[id];}); - if (bmasks == gather_masks) { - nfine_gather = std::distance(pid.begin(), sep2); - } else { - nfine_current = std::distance(pid.begin(), sep2); - } - } - } - - // only deposit / gather to coarsest grid - if (m_deposit_on_main_grid && lev > 0) { - nfine_current = 0; - } - if (m_gather_from_main_grid && lev > 0) { - nfine_gather = 0; - } - - if (nfine_current != np || nfine_gather != np) - { - particle_tmp.resize(np); - for (long ip = 0; ip < np; ++ip) { - particle_tmp[ip] = aos[pid[ip]]; - } - std::swap(aos(), particle_tmp); - - tmp.resize(np); - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = wp[pid[ip]]; - } - std::swap(wp, tmp); - - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = uxp[pid[ip]]; - } - std::swap(uxp, tmp); - - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = uyp[pid[ip]]; - } - std::swap(uyp, tmp); - - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = uzp[pid[ip]]; - } - std::swap(uzp, tmp); - } - BL_PROFILE_VAR_STOP(blp_partition); + // Determine which particles deposit/gather in the buffer, and + // which particles deposit/gather in the fine patch + long nfine_current = np; + long nfine_gather = np; + if (has_buffer && !do_not_push) { + // - Modify `nfine_current` and `nfine_gather` (in place) + // so that they correspond to the number of particles + // that deposit/gather in the fine patch respectively. + // - Reorder the particle arrays, + // so that the `nfine_current`/`nfine_gather` first particles + // deposit/gather in the fine patch + // and (thus) the `np-nfine_current`/`np-nfine_gather` last particles + // deposit/gather in the buffer + PartitionParticlesInBuffers( nfine_current, nfine_gather, np, + pti, lev, current_masks, gather_masks, uxp, uyp, uzp, wp ); } const long np_current = (cjx) ? nfine_current : np; @@ -1235,54 +1116,16 @@ PhysicalParticleContainer::Evolve (int lev, if (WarpX::use_fdtd_nci_corr) { -#if (AMREX_SPACEDIM == 2) - const Box& tbox = amrex::grow(cbox,{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noz)}); -#else - const Box& tbox = amrex::grow(cbox,{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noy), - static_cast<int>(WarpX::noz)}); -#endif - - // Filter Ex (both 2D and 3D) - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - // Safeguard for GPU - exeli = filtered_Ex.elixir(); - // Apply filter on Ex, result stored in filtered_Ex - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); - // Update exfab reference - cexfab = &filtered_Ex; - - // Filter Ez - filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); - cezfab = &filtered_Ez; - - // Filter By - filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); - cbyfab = &filtered_By; -#if (AMREX_SPACEDIM == 3) - // Filter Ey - filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); - ceyfab = &filtered_Ey; - - // Filter Bx - filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); - cbxfab = &filtered_Bx; - - // Filter Bz - filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); - cbzfab = &filtered_Bz; -#endif + // Filter arrays (*cEx)[pti], store the result in + // filtered_Ex and update pointer cexfab so that it + // points to filtered_Ex (and do the same for all + // components of E and B) + applyNCIFilter(lev-1, cbox, exeli, eyeli, ezeli, bxeli, byeli, bzeli, + filtered_Ex, filtered_Ey, filtered_Ez, + filtered_Bx, filtered_By, filtered_Bz, + (*cEx)[pti], (*cEy)[pti], (*cEz)[pti], + (*cBx)[pti], (*cBy)[pti], (*cBz)[pti], + cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab); } // Field gather for particles in gather buffers @@ -1375,6 +1218,74 @@ PhysicalParticleContainer::Evolve (int lev, } } +void +PhysicalParticleContainer::applyNCIFilter ( + int lev, const Box& box, + Elixir& exeli, Elixir& eyeli, Elixir& ezeli, + Elixir& bxeli, Elixir& byeli, Elixir& bzeli, + FArrayBox& filtered_Ex, FArrayBox& filtered_Ey, FArrayBox& filtered_Ez, + FArrayBox& filtered_Bx, FArrayBox& filtered_By, FArrayBox& filtered_Bz, + const FArrayBox& Ex, const FArrayBox& Ey, const FArrayBox& Ez, + const FArrayBox& Bx, const FArrayBox& By, const FArrayBox& Bz, + FArrayBox const * & ex_ptr, FArrayBox const * & ey_ptr, + FArrayBox const * & ez_ptr, FArrayBox const * & bx_ptr, + FArrayBox const * & by_ptr, FArrayBox const * & bz_ptr) +{ + + // Get instances of NCI Godfrey filters + const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; + const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; + +#if (AMREX_SPACEDIM == 2) + const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox), + static_cast<int>(WarpX::noz)}); +#else + const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox), + static_cast<int>(WarpX::noy), + static_cast<int>(WarpX::noz)}); +#endif + + // Filter Ex (Both 2D and 3D) + filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex, filtered_Ex.box()); + // Update ex_ptr reference + ex_ptr = &filtered_Ex; + + // Filter Ez + filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez, filtered_Ez.box()); + ez_ptr = &filtered_Ez; + + // Filter By + filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By, filtered_By.box()); + by_ptr = &filtered_By; +#if (AMREX_SPACEDIM == 3) + // Filter Ey + filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey, filtered_Ey.box()); + ey_ptr = &filtered_Ey; + + // Filter Bx + filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx, filtered_Bx.box()); + bx_ptr = &filtered_Bx; + + // Filter Bz + filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz, filtered_Bz.box()); + bz_ptr = &filtered_Bz; +#endif +} + // Loop over all particles in the particle container and // split particles tagged with p.id()=DoSplitParticleID void @@ -1382,7 +1293,7 @@ PhysicalParticleContainer::SplitParticles(int lev) { auto& mypc = WarpX::GetInstance().GetPartContainer(); auto& pctmp_split = mypc.GetPCtmp(); - Cuda::ManagedDeviceVector<Real> xp, yp, zp; + Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; RealVector psplit_x, psplit_y, psplit_z, psplit_w; RealVector psplit_ux, psplit_uy, psplit_uz; long np_split_to_add = 0; @@ -1398,7 +1309,16 @@ PhysicalParticleContainer::SplitParticles(int lev) for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { pti.GetPosition(xp, yp, zp); + + // offset for split particles is computed as a function of cell size + // and number of particles per cell, so that a uniform distribution + // before splitting results in a uniform distribution after splitting + const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); + amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0], + dx[1]/2./ppc_nd[1], + dx[2]/2./ppc_nd[2]}; + // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); // particle Struct Of Arrays data @@ -1421,9 +1341,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int ishift = -1; ishift < 2; ishift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x and z - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + kshift*dx[2]/2 ); + psplit_z.push_back( zp[i] + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1435,7 +1355,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // 4 particles in 2d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); psplit_y.push_back( yp[i] ); psplit_z.push_back( zp[i] ); psplit_ux.push_back( uxp[i] ); @@ -1445,7 +1365,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // Add one particle with offset in z psplit_x.push_back( xp[i] ); psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*dx[2]/2 ); + psplit_z.push_back( zp[i] + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1460,9 +1380,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int jshift = -1; jshift < 2; jshift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x, y and z - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); - psplit_y.push_back( yp[i] + jshift*dx[1]/2 ); - psplit_z.push_back( zp[i] + kshift*dx[2]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); + psplit_y.push_back( yp[i] + jshift*split_offset[1] ); + psplit_z.push_back( zp[i] + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1475,7 +1395,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // 6 particles in 3d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); psplit_y.push_back( yp[i] ); psplit_z.push_back( zp[i] ); psplit_ux.push_back( uxp[i] ); @@ -1484,7 +1404,7 @@ PhysicalParticleContainer::SplitParticles(int lev) psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in y psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] + ishift*dx[1]/2 ); + psplit_y.push_back( yp[i] + ishift*split_offset[1] ); psplit_z.push_back( zp[i] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); @@ -1493,7 +1413,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // Add one particle with offset in z psplit_x.push_back( xp[i] ); psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*dx[2]/2 ); + psplit_z.push_back( zp[i] + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1531,27 +1451,27 @@ PhysicalParticleContainer::SplitParticles(int lev) void PhysicalParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<Real>& xp, - Cuda::ManagedDeviceVector<Real>& yp, - Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<ParticleReal>& xp, + Cuda::ManagedDeviceVector<ParticleReal>& yp, + Cuda::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - 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 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* 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(); + ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); + ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); + ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf)) { @@ -1660,15 +1580,15 @@ 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* 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* 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(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; @@ -1694,23 +1614,23 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, - const Real* yp, const Real* zp) +void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, + const ParticleReal* yp, const ParticleReal* zp) { auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); + ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); + ParticleReal* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); + ParticleReal* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); - Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); - Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); - Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); - Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); - Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); - Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); + ParticleReal* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); + ParticleReal* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); + ParticleReal* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); + ParticleReal* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); + ParticleReal* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); + ParticleReal* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { @@ -1929,9 +1849,9 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const Array4<const Real>& by_arr = byfab->array(); const Array4<const Real>& bz_arr = bzfab->array(); - const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + const ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + const ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev); @@ -2078,15 +1998,15 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Otherwise, resize ionization_mask, and get poiters to attribs arrays. ionization_mask.resize(np); int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data(); - const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); - const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); - const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); - const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); - const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); - const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); - const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); - const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); - const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); + const ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + const ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); + const ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); + const ParticleReal * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); + const ParticleReal * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); + const ParticleReal * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); + const ParticleReal * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); + const ParticleReal * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); + const ParticleReal * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); Real c = PhysConst::c; diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 3c74baeb2..f0dfa4c83 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -11,7 +11,7 @@ * and stores them in the variables `x`, `y`, `z`. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void GetPosition( - amrex::Real& x, amrex::Real& y, amrex::Real& z, + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, const WarpXParticleContainer::ParticleType& p) { #if (AMREX_SPACEDIM==3) @@ -20,7 +20,7 @@ void GetPosition( z = p.pos(2); #else x = p.pos(0); - y = std::numeric_limits<amrex::Real>::quiet_NaN(); + y = std::numeric_limits<amrex::ParticleReal>::quiet_NaN(); z = p.pos(1); #endif } @@ -30,7 +30,7 @@ void GetPosition( AMREX_GPU_HOST_DEVICE AMREX_INLINE void SetPosition( WarpXParticleContainer::ParticleType& p, - const amrex::Real x, const amrex::Real y, const amrex::Real z) + const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z) { #if (AMREX_SPACEDIM==3) p.pos(0) = x; @@ -49,10 +49,10 @@ void SetPosition( * and store them in the variables `x`, `y`, `z` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void GetCartesianPositionFromCylindrical( - amrex::Real& x, amrex::Real& y, amrex::Real& z, - const WarpXParticleContainer::ParticleType& p, const amrex::Real theta) + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const WarpXParticleContainer::ParticleType& p, const amrex::ParticleReal theta) { - const amrex::Real r = p.pos(0); + const amrex::ParticleReal r = p.pos(0); x = r*std::cos(theta); y = r*std::sin(theta); z = p.pos(1); @@ -63,8 +63,8 @@ void GetCartesianPositionFromCylindrical( * from the values of `x`, `y`, `z` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void SetCylindricalPositionFromCartesian( - WarpXParticleContainer::ParticleType& p, amrex::Real& theta, - const amrex::Real x, const amrex::Real y, const amrex::Real z ) + WarpXParticleContainer::ParticleType& p, amrex::ParticleReal& theta, + const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z ) { theta = std::atan2(y, x); p.pos(0) = std::sqrt(x*x + y*y); diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H index a33058347..205cc9a71 100644 --- a/Source/Particles/Pusher/UpdateMomentumBoris.H +++ b/Source/Particles/Pusher/UpdateMomentumBoris.H @@ -7,9 +7,9 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumBoris( - amrex::Real& ux, amrex::Real& uy, amrex::Real& uz, - const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez, - const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz, + amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz, + const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez, + const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz, const amrex::Real q, const amrex::Real m, const amrex::Real dt ) { const amrex::Real econst = 0.5*q*dt/m; diff --git a/Source/Particles/Pusher/UpdateMomentumVay.H b/Source/Particles/Pusher/UpdateMomentumVay.H index 1f0f19e63..433a891c5 100644 --- a/Source/Particles/Pusher/UpdateMomentumVay.H +++ b/Source/Particles/Pusher/UpdateMomentumVay.H @@ -9,9 +9,9 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumVay( - amrex::Real& ux, amrex::Real& uy, amrex::Real& uz, - const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez, - const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz, + amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz, + const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez, + const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz, const amrex::Real q, const amrex::Real m, const amrex::Real dt ) { // Constants diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index a9df63a30..da0e9cdf9 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -9,8 +9,8 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdatePosition( - amrex::Real& x, amrex::Real& y, amrex::Real& z, - const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::Real dt ) { diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H index bd6e6cd21..f95c2b09d 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -10,8 +10,8 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdatePositionPhoton( - amrex::Real& x, amrex::Real& y, amrex::Real& z, - const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::Real dt ) { // Compute speed of light over inverse of momentum modulus diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index be3dd21f9..3abbb4afe 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -44,9 +44,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; virtual void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 7d129fc01..891ade76d 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -76,7 +76,7 @@ RigidInjectedParticleContainer::RemapParticles() // Note that the particles are already in the boosted frame. // This value is saved to advance the particles not injected yet - Cuda::ManagedDeviceVector<Real> xp, yp, zp; + Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -136,7 +136,7 @@ RigidInjectedParticleContainer::BoostandRemapParticles() #pragma omp parallel #endif { - Cuda::ManagedDeviceVector<Real> xp, yp, zp; + Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; for (WarpXParIter pti(*this, 0); pti.isValid(); ++pti) { @@ -207,9 +207,9 @@ RigidInjectedParticleContainer::BoostandRemapParticles() void RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<Real>& xp, - Cuda::ManagedDeviceVector<Real>& yp, - Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<ParticleReal>& xp, + Cuda::ManagedDeviceVector<ParticleReal>& yp, + Cuda::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { @@ -220,21 +220,21 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& uzp = attribs[PIdx::uz]; // Save the position and momenta, making copies - Cuda::ManagedDeviceVector<Real> xp_save, yp_save, zp_save; + Cuda::ManagedDeviceVector<ParticleReal> xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; - 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 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(); + ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); + ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); + ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = uxp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = uyp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = uzp.dataPtr(); + ParticleReal* const AMREX_RESTRICT Exp = attribs[PIdx::Ex].dataPtr(); + ParticleReal* const AMREX_RESTRICT Eyp = attribs[PIdx::Ey].dataPtr(); + ParticleReal* const AMREX_RESTRICT Ezp = attribs[PIdx::Ez].dataPtr(); + ParticleReal* const AMREX_RESTRICT Bxp = attribs[PIdx::Bx].dataPtr(); + ParticleReal* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr(); + ParticleReal* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { // If the old values are not already saved, create copies here. @@ -271,12 +271,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (!done_injecting_lev) { - Real* AMREX_RESTRICT x_save = xp_save.dataPtr(); - Real* AMREX_RESTRICT y_save = yp_save.dataPtr(); - Real* AMREX_RESTRICT z_save = zp_save.dataPtr(); - Real* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); - Real* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); - Real* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT x_save = xp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT y_save = yp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT z_save = zp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. @@ -415,16 +415,16 @@ 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* const AMREX_RESTRICT zp = m_zp[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* 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(); + const ParticleReal* const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + ParticleReal* const AMREX_RESTRICT uxpp = uxp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uypp = uyp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uzpp = uzp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; @@ -450,10 +450,10 @@ 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]; + const ParticleReal* const AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + const ParticleReal* const AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + const ParticleReal* const AMREX_RESTRICT uz_save = uzp_save.dataPtr(); + const ParticleReal zz = zinject_plane_levels[lev]; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { if (zp[i] <= zz) { diff --git a/Source/Particles/Sorting/Make.package b/Source/Particles/Sorting/Make.package new file mode 100644 index 000000000..750d2607e --- /dev/null +++ b/Source/Particles/Sorting/Make.package @@ -0,0 +1,3 @@ +CEXE_sources += Partition.cpp +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Sorting +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Sorting diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp new file mode 100644 index 000000000..683dbbd04 --- /dev/null +++ b/Source/Particles/Sorting/Partition.cpp @@ -0,0 +1,141 @@ +#include <WarpX.H> +#include <PhysicalParticleContainer.H> +#include <AMReX_Particles.H> + +using namespace amrex; + +/* \brief Determine which particles deposit/gather in the buffer, and + * and reorder the particle arrays accordingly + * + * More specifically: + * - Modify `nfine_current` and `nfine_gather` (in place) + * so that they correspond to the number of particles + * that deposit/gather in the fine patch respectively. + * - Reorder the particle arrays, + * so that the `nfine_current`/`nfine_gather` first particles + * deposit/gather in the fine patch + * and (thus) the `np-nfine_current`/`np-nfine_gather` last particles + * deposit/gather in the buffer + * + * \param nfine_current number of particles that deposit to the fine patch + * (modified by this function) + * \param nfine_gather number of particles that gather into the fine patch + * (modified by this function) + * \param np total number of particles in this tile + * \param pti object that holds the particle information for this tile + * \param lev current refinement level + * \param current_masks indicates, for each cell, whether that cell is + * in the deposition buffers or in the interior of the fine patch + * \param gather_masks indicates, for each cell, whether that cell is + * in the gather buffers or in the interior of the fine patch + * \param uxp, uyp, uzp, wp references to the particle momenta and weight + * (modified by this function) + */ +void +PhysicalParticleContainer::PartitionParticlesInBuffers( + long& nfine_current, long& nfine_gather, long const np, + WarpXParIter& pti, int const lev, + iMultiFab const* current_masks, + iMultiFab const* gather_masks, + RealVector& uxp, RealVector& uyp, RealVector& uzp, RealVector& wp) +{ + BL_PROFILE("PPC::Evolve::partition"); + + std::vector<bool> inexflag; + inexflag.resize(np); + + auto& aos = pti.GetArrayOfStructs(); + + // We need to partition the large buffer first + iMultiFab const* bmasks = + (WarpX::n_field_gather_buffer >= WarpX::n_current_deposition_buffer) ? + gather_masks : current_masks; + int i = 0; + const auto& msk = (*bmasks)[pti]; + for (const auto& p : aos) { + const IntVect& iv = Index(p, lev); + inexflag[i++] = msk(iv); + } + + Vector<long> pid; + pid.resize(np); + std::iota(pid.begin(), pid.end(), 0L); + auto sep = std::stable_partition(pid.begin(), pid.end(), + [&inexflag](long id) { return inexflag[id]; }); + + if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer) { + nfine_current = nfine_gather = std::distance(pid.begin(), sep); + } else if (sep != pid.end()) { + int n_buf; + if (bmasks == gather_masks) { + nfine_gather = std::distance(pid.begin(), sep); + bmasks = current_masks; + n_buf = WarpX::n_current_deposition_buffer; + } else { + nfine_current = std::distance(pid.begin(), sep); + bmasks = gather_masks; + n_buf = WarpX::n_field_gather_buffer; + } + if (n_buf > 0) + { + const auto& msk2 = (*bmasks)[pti]; + for (auto it = sep; it != pid.end(); ++it) { + const long id = *it; + const IntVect& iv = Index(aos[id], lev); + inexflag[id] = msk2(iv); + } + + auto sep2 = std::stable_partition(sep, pid.end(), + [&inexflag](long id) {return inexflag[id];}); + if (bmasks == gather_masks) { + nfine_gather = std::distance(pid.begin(), sep2); + } else { + nfine_current = std::distance(pid.begin(), sep2); + } + } + } + + // only deposit / gather to coarsest grid + if (m_deposit_on_main_grid && lev > 0) { + nfine_current = 0; + } + if (m_gather_from_main_grid && lev > 0) { + nfine_gather = 0; + } + + if (nfine_current != np || nfine_gather != np) + { + + ParticleVector particle_tmp; + particle_tmp.resize(np); + + for (long ip = 0; ip < np; ++ip) { + particle_tmp[ip] = aos[pid[ip]]; + } + std::swap(aos(), particle_tmp); + + RealVector tmp; + tmp.resize(np); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = wp[pid[ip]]; + } + std::swap(wp, tmp); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = uxp[pid[ip]]; + } + std::swap(uxp, tmp); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = uyp[pid[ip]]; + } + std::swap(uyp, tmp); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = uzp[pid[ip]]; + } + std::swap(uzp, tmp); + } + +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ee75bc511..7b0d2d1d0 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -68,12 +68,12 @@ public: WarpXParIter (ContainerType& pc, int level); #if (AMREX_SPACEDIM == 2) - void GetPosition (amrex::Cuda::ManagedDeviceVector<amrex::Real>& x, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& y, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& z) const; - void SetPosition (const amrex::Cuda::ManagedDeviceVector<amrex::Real>& x, - const amrex::Cuda::ManagedDeviceVector<amrex::Real>& y, - const amrex::Cuda::ManagedDeviceVector<amrex::Real>& z); + void GetPosition (amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& x, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& y, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& z) const; + void SetPosition (const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& x, + const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& y, + const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& z); #endif const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { return GetStructOfArrays().GetRealData(); @@ -104,7 +104,7 @@ class WarpXParticleContainer public: friend MultiParticleContainer; - // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // amrex::StructOfArrays with DiagIdx::nattribs amrex::ParticleReal components // and 0 int components for the particle data. using DiagnosticParticleData = amrex::StructOfArrays<DiagIdx::nattribs, 0>; // DiagnosticParticles is a vector, with one element per MR level. @@ -232,17 +232,17 @@ public: amrex::Real maxParticleVelocity(bool local = false); void AddNParticles (int lev, - int n, const amrex::Real* x, const amrex::Real* y, const amrex::Real* z, - const amrex::Real* vx, const amrex::Real* vy, const amrex::Real* vz, - int nattr, const amrex::Real* attr, int uniqueparticles, int id=-1); + int n, const amrex::ParticleReal* x, const amrex::ParticleReal* y, const amrex::ParticleReal* z, + const amrex::ParticleReal* vx, const amrex::ParticleReal* vy, const amrex::ParticleReal* vz, + int nattr, const amrex::ParticleReal* attr, int uniqueparticles, int id=-1); void AddOneParticle (int lev, int grid, int tile, - amrex::Real x, amrex::Real y, amrex::Real z, - std::array<amrex::Real,PIdx::nattribs>& attribs); + amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, + std::array<amrex::ParticleReal,PIdx::nattribs>& attribs); void AddOneParticle (ParticleTileType& particle_tile, - amrex::Real x, amrex::Real y, amrex::Real z, - std::array<amrex::Real,PIdx::nattribs>& attribs); + amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, + std::array<amrex::ParticleReal,PIdx::nattribs>& attribs); virtual void ReadHeader (std::istream& is); @@ -326,7 +326,7 @@ protected: amrex::Vector<amrex::FArrayBox> local_jy; amrex::Vector<amrex::FArrayBox> local_jz; - using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::Real>; + using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>; using PairIndex = std::pair<int, int>; amrex::Vector<DataContainer> m_xp, m_yp, m_zp; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 176c147da..65a82f233 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -25,7 +25,9 @@ WarpXParIter::WarpXParIter (ContainerType& pc, int level) #if (AMREX_SPACEDIM == 2) void -WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDeviceVector<Real>& y, Cuda::ManagedDeviceVector<Real>& z) const +WarpXParIter::GetPosition (Gpu::ManagedDeviceVector<ParticleReal>& x, + Gpu::ManagedDeviceVector<ParticleReal>& y, + Gpu::ManagedDeviceVector<ParticleReal>& z) const { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); #ifdef WARPX_DIM_RZ @@ -38,17 +40,19 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDevi x[i] = x[i]*std::cos(theta[i]); } #else - y.resize(x.size(), std::numeric_limits<Real>::quiet_NaN()); + y.resize(x.size(), std::numeric_limits<ParticleReal>::quiet_NaN()); #endif } void -WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector<Real>& x, const Cuda::ManagedDeviceVector<Real>& y, const Cuda::ManagedDeviceVector<Real>& z) +WarpXParIter::SetPosition (const Gpu::ManagedDeviceVector<ParticleReal>& x, + const Gpu::ManagedDeviceVector<ParticleReal>& y, + const Gpu::ManagedDeviceVector<ParticleReal>& z) { #ifdef WARPX_DIM_RZ auto& attribs = GetAttribs(); auto& theta = attribs[PIdx::theta]; - Cuda::ManagedDeviceVector<Real> r(x.size()); + Gpu::ManagedDeviceVector<ParticleReal> r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { theta[i] = std::atan2(y[i], x[i]); r[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); @@ -132,8 +136,8 @@ WarpXParticleContainer::AllocData () void WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, - Real x, Real y, Real z, - std::array<Real,PIdx::nattribs>& attribs) + ParticleReal x, ParticleReal y, ParticleReal z, + std::array<ParticleReal,PIdx::nattribs>& attribs) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid, tile); AddOneParticle(particle_tile, x, y, z, attribs); @@ -141,8 +145,8 @@ WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, void WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, - Real x, Real y, Real z, - std::array<Real,PIdx::nattribs>& attribs) + ParticleReal x, ParticleReal y, ParticleReal z, + std::array<ParticleReal,PIdx::nattribs>& attribs) { ParticleType p; p.id() = ParticleType::NextID(); @@ -171,12 +175,12 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, void WarpXParticleContainer::AddNParticles (int lev, - int n, const Real* x, const Real* y, const Real* z, - const Real* vx, const Real* vy, const Real* vz, - int nattr, const Real* attr, int uniqueparticles, int id) + int n, const ParticleReal* x, const ParticleReal* y, const ParticleReal* z, + const ParticleReal* vx, const ParticleReal* vy, const ParticleReal* vz, + int nattr, const ParticleReal* attr, int uniqueparticles, int id) { BL_ASSERT(nattr == 1); - const Real* weight = attr; + const ParticleReal* weight = attr; int ibegin, iend; if (uniqueparticles) { @@ -204,7 +208,7 @@ WarpXParticleContainer::AddNParticles (int lev, std::size_t np = iend-ibegin; #ifdef WARPX_DIM_RZ - Vector<Real> theta(np); + Vector<ParticleReal> theta(np); #endif for (int i = ibegin; i < iend; ++i) @@ -253,7 +257,7 @@ WarpXParticleContainer::AddNParticles (int lev, { #ifdef WARPX_DIM_RZ if (comp == PIdx::theta) { - particle_tile.push_back_real(comp, theta.front(), theta.back()); + particle_tile.push_back_real(comp, theta.data(), theta.data() + np); } else { particle_tile.push_back_real(comp, np, 0.0); @@ -362,9 +366,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // CPU, tiling: deposit into local_jx // (same for jx and jz) - Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -503,9 +507,9 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // GPU, no tiling: deposit directly in rho // CPU, tiling: deposit into local_rho - Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -731,7 +735,7 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) { Real WarpXParticleContainer::maxParticleVelocity(bool local) { - amrex::Real max_v = 0.0; + amrex::ParticleReal max_v = 0.0; for (int lev = 0; lev <= finestLevel(); ++lev) { @@ -745,12 +749,12 @@ Real WarpXParticleContainer::maxParticleVelocity(bool local) { auto& uy = pti.GetAttribs(PIdx::uy); auto& uz = pti.GetAttribs(PIdx::uz); for (unsigned long i = 0; i < ux.size(); i++) { - max_v = std::max(max_v, sqrt(ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])); + max_v = std::max(max_v, std::sqrt(ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])); } } } - if (!local) ParallelDescriptor::ReduceRealMax(max_v); + if (!local) ParallelAllReduce::Max(max_v, ParallelDescriptor::Communicator()); return max_v; } @@ -819,17 +823,17 @@ WarpXParticleContainer::PushX (int lev, Real dt) ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); // - momenta are stored as a struct of array, in `attribs` auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + ParticleReal* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); #ifdef WARPX_DIM_RZ - Real* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); + ParticleReal* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); #endif // Loop over the particles and update their position amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { ParticleType& p = pstructs[i]; // Particle object that gets updated - Real x, y, z; // Temporary variables + ParticleReal x, y, z; // Temporary variables #ifndef WARPX_DIM_RZ GetPosition( x, y, z, p ); // Initialize x, y, z UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); diff --git a/Source/Particles/deposit_cic.F90 b/Source/Particles/deposit_cic.F90 index 1fe80016f..2831ce96b 100644 --- a/Source/Particles/deposit_cic.F90 +++ b/Source/Particles/deposit_cic.F90 @@ -1,7 +1,7 @@ module warpx_ES_deposit_cic use iso_c_binding - use amrex_fort_module, only : amrex_real + use amrex_fort_module, only : amrex_real, amrex_particle_real implicit none @@ -28,8 +28,8 @@ contains ng) & bind(c,name='warpx_deposit_cic_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(in) :: weights(np) + real(amrex_particle_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: weights(np) real(amrex_real), intent(in) :: charge integer, intent(in) :: lo(3) integer, intent(in) :: hi(3) @@ -86,8 +86,8 @@ contains ng) & bind(c,name='warpx_deposit_cic_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(in) :: weights(np) + real(amrex_particle_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: weights(np) real(amrex_real), intent(in) :: charge integer, intent(in) :: lo(2) integer, intent(in) :: hi(2) diff --git a/Source/Particles/interpolate_cic.F90 b/Source/Particles/interpolate_cic.F90 index 005ab35f4..3eb361d2f 100644 --- a/Source/Particles/interpolate_cic.F90 +++ b/Source/Particles/interpolate_cic.F90 @@ -1,7 +1,7 @@ module warpx_ES_interpolate_cic use iso_c_binding - use amrex_fort_module, only : amrex_real + use amrex_fort_module, only : amrex_real, amrex_particle_real implicit none @@ -31,7 +31,7 @@ contains lo, hi, plo, dx, ng) & bind(c,name='warpx_interpolate_cic_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) integer, intent(in) :: ng integer, intent(in) :: lo(3) @@ -103,7 +103,7 @@ contains lo, hi, plo, dx, ng) & bind(c,name='warpx_interpolate_cic_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) integer, intent(in) :: ng integer, intent(in) :: lo(2) @@ -157,7 +157,7 @@ contains plo, ng, lev) & bind(c,name='warpx_interpolate_cic_two_levels_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) integer, intent(in) :: ng, lev integer, intent(in) :: lo(3), hi(3) @@ -290,7 +290,7 @@ contains plo, ng, lev) & bind(c,name='warpx_interpolate_cic_two_levels_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) integer, intent(in) :: ng, lev integer, intent(in) :: lo(2), hi(2) diff --git a/Source/Particles/push_particles_ES.F90 b/Source/Particles/push_particles_ES.F90 index 53dd3c181..b84f48d5f 100644 --- a/Source/Particles/push_particles_ES.F90 +++ b/Source/Particles/push_particles_ES.F90 @@ -1,7 +1,7 @@ module warpx_ES_push_particles use iso_c_binding - use amrex_fort_module, only : amrex_real + use amrex_fort_module, only : amrex_real, amrex_particle_real implicit none @@ -36,8 +36,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np), Ez_p(np) real(amrex_real), intent(in) :: charge real(amrex_real), intent(in) :: mass @@ -100,8 +100,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np) real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np) real(amrex_real), intent(in) :: charge real(amrex_real), intent(in) :: mass @@ -167,8 +167,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_positions_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) real(amrex_real), intent(in) :: dt real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) @@ -219,8 +219,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_positions_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np) real(amrex_real), intent(in) :: dt real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) |