diff options
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 155 |
1 files changed, 89 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); |