diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 155 | ||||
| -rw-r--r-- | Source/Particles/Pusher/Make.package | 2 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdateMomentumBoris.H | 47 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdateMomentumVay.H | 54 |
4 files changed, 192 insertions, 66 deletions
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d47a7b220..3a512d9dc 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -7,6 +7,12 @@ #include <WarpXConst.H> #include <WarpXWrappers.h> +#include <WarpXAlgorithmSelection.H> + +// Import low-level single-particle kernels +#include <UpdatePosition.H> +#include <UpdateMomentumBoris.H> +#include <UpdateMomentumVay.H> using namespace amrex; @@ -62,7 +68,7 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, ++np; } } - + return np; } @@ -82,7 +88,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp 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 + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -91,7 +97,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); @@ -108,9 +114,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -181,7 +187,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -195,7 +201,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, if (ParallelDescriptor::IOProcessor()) { std::array<Real, 3> u; Real weight; - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -226,7 +232,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, u_tmp[0] *= std::pow(-1,ix); y_tmp = y*std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); - CheckAndAddParticle(x_tmp, y_tmp, z, + CheckAndAddParticle(x_tmp, y_tmp, z, u_tmp, weight/4); } } @@ -263,7 +269,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); particle_tile.push_back_real(particle_comps["zold"], z); - + particle_tile.push_back_real(particle_comps["uxold"], u[0]); particle_tile.push_back_real(particle_comps["uyold"], u[1]); particle_tile.push_back_real(particle_comps["uzold"], u[2]); @@ -537,7 +543,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::ux] = u[0]; attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); @@ -838,7 +844,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) costarr(i,j,k) += wt; }); } - } + } } } #endif @@ -1166,11 +1172,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PICSAR::FieldGather", blp_pxr_fg); BL_PROFILE_VAR_NS("PICSAR::ParticlePush", blp_pxr_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)); - // Get instances of NCI Godfrey filters + // 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; @@ -1184,7 +1190,7 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1392,7 +1398,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); - + if (! do_not_push) { // @@ -1403,7 +1409,7 @@ PhysicalParticleContainer::Evolve (int lev, const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); const int* ixyzmin_grid = box.loVect(); - + const long np_gather = (cEx) ? nfine_gather : np; BL_PROFILE_VAR_START(blp_pxr_fg); @@ -1479,13 +1485,13 @@ PhysicalParticleContainer::Evolve (int lev, 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(); @@ -1493,7 +1499,7 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + long ncrse = np - nfine_gather; warpx_geteb_energy_conserving( &ncrse, @@ -1522,7 +1528,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_pxr_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_pxr_pp); @@ -1560,7 +1566,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev); if (cost) { @@ -1712,21 +1718,21 @@ PhysicalParticleContainer::SplitParticles(int lev) } // Add local arrays psplit_x etc. to the temporary // particle container pctmp_split. Split particles - // are tagged with p.id()=NoSplitParticleID so that + // are tagged with p.id()=NoSplitParticleID so that // they are not re-split when entering a higher level // AddNParticles calls Redistribute, so that particles // in pctmp_split are in the proper grids and tiles - pctmp_split.AddNParticles(lev, - np_split_to_add, - psplit_x.dataPtr(), - psplit_y.dataPtr(), - psplit_z.dataPtr(), - psplit_ux.dataPtr(), - psplit_uy.dataPtr(), - psplit_uz.dataPtr(), - 1, - psplit_w.dataPtr(), - 1, NoSplitParticleID); + pctmp_split.AddNParticles(lev, + np_split_to_add, + psplit_x.dataPtr(), + psplit_y.dataPtr(), + psplit_z.dataPtr(), + psplit_ux.dataPtr(), + psplit_uy.dataPtr(), + psplit_uz.dataPtr(), + 1, + psplit_w.dataPtr(), + 1, NoSplitParticleID); // Copy particles from tmp to current particle container addParticles(pctmp_split,1); // Clear tmp container @@ -1742,35 +1748,52 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, Real dt) { + // This wraps the call to warpx_particle_pusher 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(); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + copy_attribs(pti, x, y, z); } - // The following attributes should be included in CPP version of warpx_particle_pusher - // This wraps the call to warpx_particle_pusher 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(); - - 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); + // 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( ux[i], uy[i], uz[i], gi[i], + Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i], + Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + }; } @@ -1793,7 +1816,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1879,20 +1902,20 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, { 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(); - + Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["yold"]).dataPtr(); Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["zold"]).dataPtr(); Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - + const long np = pti.numParticles(); - + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { xpold[i]=xp[i]; @@ -1938,7 +1961,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -2020,7 +2043,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package index 8c8e77905..95a38fa2d 100644 --- a/Source/Particles/Pusher/Make.package +++ b/Source/Particles/Pusher/Make.package @@ -1,4 +1,6 @@ CEXE_headers += GetAndSetPosition.H CEXE_headers += UpdatePosition.H +CEXE_headers += UpdateMomentumBoris.H +CEXE_headers += UpdateMomentumVay.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H new file mode 100644 index 000000000..71e9a8ed1 --- /dev/null +++ b/Source/Particles/Pusher/UpdateMomentumBoris.H @@ -0,0 +1,47 @@ +#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_H_ +#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_H_ + +#include <AMReX_REAL.H> + +/* \brief Push the particle's positions over one timestep, + * 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, amrex::Real& gaminv, + const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez, + const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz, + const amrex::Real q, const amrex::Real m, const amrex::Real dt ) +{ + const amrex::Real econst = 0.5*q*dt/m; + + // First half-push for E + ux += econst*Ex; + uy += econst*Ey; + uz += econst*Ez; + // Compute temporary gamma factor + constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2); + // Magnetic rotation + // - Compute temporary variables + const amrex::Real tx = econst*inv_gamma*Bx; + const amrex::Real ty = econst*inv_gamma*By; + const amrex::Real tz = econst*inv_gamma*Bz; + const amrex::Real tsqi = 2./(1. + tx*tx + ty*ty + tz*tz); + const amrex::Real sx = tx*tsqi; + const amrex::Real sy = ty*tsqi; + const amrex::Real sz = tz*tsqi; + const amrex::Real ux_p = ux + uy*tz - uz*ty; + const amrex::Real uy_p = uy + uz*tx - ux*tz; + const amrex::Real uz_p = uz + ux*ty - uy*tx; + // - Update momentum + ux += uy_p*sz - uz_p*sy; + uy += uz_p*sx - ux_p*sz; + uz += ux_p*sy - uy_p*sx; + // Second half-push for E + ux += econst*Ex; + uy += econst*Ey; + uz += econst*Ez; + gaminv = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2); +} + +#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_H_ diff --git a/Source/Particles/Pusher/UpdateMomentumVay.H b/Source/Particles/Pusher/UpdateMomentumVay.H new file mode 100644 index 000000000..044297e22 --- /dev/null +++ b/Source/Particles/Pusher/UpdateMomentumVay.H @@ -0,0 +1,54 @@ +#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_ +#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_ + +#include <AMReX_FArrayBox.H> +#include <WarpXConst.H> +#include <AMReX_REAL.H> + +/* \brief Push the particle's positions over one timestep, + * 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, amrex::Real& gaminv, + const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez, + const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz, + const amrex::Real q, const amrex::Real m, const amrex::Real dt ) +{ + // Constants + const amrex::Real econst = q*dt/m; + const amrex::Real bconst = 0.5*q*dt/m; + constexpr amrex::Real invclight = 1./PhysConst::c; + constexpr amrex::Real invclightsq = 1./(PhysConst::c*PhysConst::c); + // Compute initial gamma + const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*invclightsq); + // Get tau + const amrex::Real taux = bconst*Bx; + const amrex::Real tauy = bconst*By; + const amrex::Real tauz = bconst*Bz; + const amrex::Real tausq = taux*taux+tauy*tauy+tauz*tauz; + // Get U', gamma'^2 + const amrex::Real uxpr = ux + econst*Ex + (uy*tauz-uz*tauy)*inv_gamma; + const amrex::Real uypr = uy + econst*Ey + (uz*taux-ux*tauz)*inv_gamma; + const amrex::Real uzpr = uz + econst*Ez + (ux*tauy-uy*taux)*inv_gamma; + const amrex::Real gprsq = (1. + (uxpr*uxpr + uypr*uypr + uzpr*uzpr)*invclightsq); + // Get u* + const amrex::Real ust = (uxpr*taux + uypr*tauy + uzpr*tauz)*invclight; + // Get new gamma + const amrex::Real sigma = gprsq-tausq; + const amrex::Real gisq = 2./(sigma + std::sqrt(sigma*sigma + 4.*(tausq + ust*ust)) ); + // Get t, s + const amrex::Real bg = bconst*std::sqrt(gisq); + const amrex::Real tx = bg*Bx; + const amrex::Real ty = bg*By; + const amrex::Real tz = bg*Bz; + const amrex::Real s = 1./(1.+tausq*gisq); + // Get t.u' + const amrex::Real tu = tx*uxpr + ty*uypr + tz*uzpr; + // Get new U + ux = s*(uxpr+tx*tu+uypr*tz-uzpr*ty); + uy = s*(uypr+ty*tu+uzpr*tx-uxpr*tz); + uz = s*(uzpr+tz*tu+uxpr*ty-uypr*tx); + gaminv = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*invclightsq); +} + +#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_ |