aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles/PhysicalParticleContainer.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp1323
1 files changed, 571 insertions, 752 deletions
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index d47a7b220..d10390204 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -6,65 +6,16 @@
#include <WarpX.H>
#include <WarpXConst.H>
#include <WarpXWrappers.h>
+#include <FieldGather.H>
+#include <WarpXAlgorithmSelection.H>
-using namespace amrex;
-
-long PhysicalParticleContainer::
-NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox,
- const RealBox& tile_realbox, const RealBox& particle_real_box)
-{
- const int lev = 0;
- const Geometry& geom = Geom(lev);
- int num_ppc = plasma_injector->num_particles_per_cell;
- const Real* dx = geom.CellSize();
+// Import low-level single-particle kernels
+#include <UpdatePosition.H>
+#include <UpdateMomentumBoris.H>
+#include <UpdateMomentumVay.H>
- long np = 0;
- const auto& overlap_corner = overlap_realbox.lo();
- for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv))
- {
- int fac;
- if (do_continuous_injection) {
-#if ( AMREX_SPACEDIM == 3 )
- Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0];
- Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1];
- Real z = overlap_corner[2] + (iv[2] + 0.5)*dx[2];
-#elif ( AMREX_SPACEDIM == 2 )
- Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0];
- Real y = 0;
- Real z = overlap_corner[1] + (iv[1] + 0.5)*dx[1];
-#endif
- fac = GetRefineFac(x, y, z);
- } else {
- fac = 1.0;
- }
-
- int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac);
- for (int i_part=0; i_part<ref_num_ppc;i_part++) {
- std::array<Real, 3> r;
- plasma_injector->getPositionUnitBox(r, i_part, fac);
-#if ( AMREX_SPACEDIM == 3 )
- Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0];
- Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1];
- Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2];
-#elif ( AMREX_SPACEDIM == 2 )
- Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0];
- Real y = 0;
- Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1];
-#endif
- // If the new particle is not inside the tile box,
- // go to the next generated particle.
-#if ( AMREX_SPACEDIM == 3 )
- if(!tile_realbox.contains( RealVect{x, y, z} )) continue;
-#elif ( AMREX_SPACEDIM == 2 )
- if(!tile_realbox.contains( RealVect{x, z} )) continue;
-#endif
- ++np;
- }
- }
-
- return np;
-}
+using namespace amrex;
PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies,
const std::string& name)
@@ -127,9 +78,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core)
void PhysicalParticleContainer::InitData()
{
AddParticles(0); // Note - add on level 0
- if (maxLevel() > 0) {
- Redistribute(); // We then redistribute
- }
+ Redistribute(); // We then redistribute
}
void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real& z, std::array<Real, 3>& u)
@@ -193,45 +142,36 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
std::normal_distribution<double> distz(z_m, z_rms);
if (ParallelDescriptor::IOProcessor()) {
- std::array<Real, 3> u;
- Real weight;
// 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;
}
for (long i = 0; i < npart; ++i) {
-#if ( AMREX_SPACEDIM == 3 | WARPX_RZ)
- weight = q_tot/npart/charge;
+#if ( AMREX_SPACEDIM == 3 | WARPX_DIM_RZ)
+ Real weight = q_tot/npart/charge;
Real x = distx(mt);
Real y = disty(mt);
Real z = distz(mt);
#elif ( AMREX_SPACEDIM == 2 )
- weight = q_tot/npart/charge/y_rms;
+ Real weight = q_tot/npart/charge/y_rms;
Real x = distx(mt);
Real y = 0.;
Real z = distz(mt);
#endif
if (plasma_injector->insideBounds(x, y, z)) {
- plasma_injector->getMomentum(u, x, y, z);
+ XDim3 u = plasma_injector->getMomentum(x, y, z);
+ u.x *= PhysConst::c;
+ u.y *= PhysConst::c;
+ u.z *= PhysConst::c;
if (do_symmetrize){
- std::array<Real, 3> u_tmp;
- Real x_tmp, y_tmp;
// Add four particles to the beam:
- // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy)
- for (int ix=0; ix<2; ix++){
- for (int iy=0; iy<2; iy++){
- u_tmp = u;
- x_tmp = x*std::pow(-1,ix);
- 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,
- u_tmp, weight/4);
- }
- }
+ CheckAndAddParticle( x, y, z, { u.x, u.y, u.z}, weight/4. );
+ CheckAndAddParticle( x,-y, z, { u.x,-u.y, u.z}, weight/4. );
+ CheckAndAddParticle(-x, y, z, {-u.x, u.y, u.z}, weight/4. );
+ CheckAndAddParticle(-x,-y, z, {-u.x,-u.y, u.z}, weight/4. );
} else {
- CheckAndAddParticle(x, y, z, u, weight);
+ CheckAndAddParticle(x, y, z, {u.x,u.y,u.z}, weight);
}
}
}
@@ -322,28 +262,19 @@ PhysicalParticleContainer::AddParticles (int lev)
void
PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
{
-#ifdef AMREX_USE_GPU
- AddPlasmaGPU(lev, part_realbox);
-#else
- AddPlasmaCPU(lev, part_realbox);
-#endif
-}
-
-void
-PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox)
-{
- BL_PROFILE("PhysicalParticleContainer::AddPlasmaCPU");
+ BL_PROFILE("PhysicalParticleContainer::AddPlasma");
// If no part_realbox is provided, initialize particles in the whole domain
const Geometry& geom = Geom(lev);
if (!part_realbox.ok()) part_realbox = geom.ProbDomain();
int num_ppc = plasma_injector->num_particles_per_cell;
-#ifdef WARPX_RZ
+#ifdef WARPX_DIM_RZ
Real rmax = std::min(plasma_injector->xmax, part_realbox.hi(0));
#endif
- const Real* dx = geom.CellSize();
+ const auto dx = geom.CellSizeArray();
+ const auto problo = geom.ProbLoArray();
Real scale_fac;
#if AMREX_SPACEDIM==3
@@ -358,490 +289,341 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox)
const int grid_id = mfi.index();
const int tile_id = mfi.LocalTileIndex();
GetParticles(lev)[std::make_pair(grid_id, tile_id)];
+ if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) {
+ DefineAndReturnParticleTile(lev, grid_id, tile_id);
+ }
}
#endif
MultiFab* cost = WarpX::getCosts(lev);
- if ( (not m_refined_injection_mask) and WarpX::do_moving_window)
+ const int nlevs = numLevels();
+ static bool refine_injection = false;
+ static Box fine_injection_box;
+ static int rrfac = 1;
+ // This does not work if the mesh is dynamic. But in that case, we should
+ // not use refined injected either. We also assume there is only one fine level.
+ if (WarpX::do_moving_window and WarpX::refine_plasma
+ and do_continuous_injection and nlevs == 2)
{
- Box mask_box = geom.Domain();
- mask_box.setSmall(WarpX::moving_window_dir, 0);
- mask_box.setBig(WarpX::moving_window_dir, 0);
- m_refined_injection_mask.reset( new IArrayBox(mask_box));
- m_refined_injection_mask->setVal(-1);
+ refine_injection = true;
+ fine_injection_box = ParticleBoxArray(1).minimalBox();
+ fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits<int>::lowest());
+ fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits<int>::max());
+ rrfac = m_gdb->refRatio(0)[0];
+ fine_injection_box.coarsen(rrfac);
}
+ InjectorPosition* inj_pos = plasma_injector->getInjectorPosition();
+ InjectorDensity* inj_rho = plasma_injector->getInjectorDensity();
+ InjectorMomentum* inj_mom = plasma_injector->getInjectorMomentum();
+ Real gamma_boost = WarpX::gamma_boost;
+ Real beta_boost = WarpX::beta_boost;
+ Real t = WarpX::GetInstance().gett_new(lev);
+ Real density_min = plasma_injector->density_min;
+ Real density_max = plasma_injector->density_max;
+
+#ifdef WARPX_DIM_RZ
+ bool radially_weighted = plasma_injector->radially_weighted;
+#endif
+
MFItInfo info;
- if (do_tiling) {
+ if (do_tiling && Gpu::notInLaunchRegion()) {
info.EnableTiling(tile_size);
}
- info.SetDynamic(true);
-
#ifdef _OPENMP
+ info.SetDynamic(true);
#pragma omp parallel if (not WarpX::serialize_ics)
#endif
+ for (MFIter mfi = MakeMFIter(lev, info); mfi.isValid(); ++mfi)
{
- std::array<Real,PIdx::nattribs> attribs;
- attribs.fill(0.0);
-
- // Loop through the tiles
- for (MFIter mfi = MakeMFIter(lev, info); mfi.isValid(); ++mfi) {
-
- Real wt = amrex::second();
-
- const Box& tile_box = mfi.tilebox();
- const RealBox tile_realbox = WarpX::getRealBox(tile_box, lev);
-
- // Find the cells of part_box that overlap with tile_realbox
- // If there is no overlap, just go to the next tile in the loop
- RealBox overlap_realbox;
- Box overlap_box;
- Real ncells_adjust;
- bool no_overlap = 0;
-
- for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
- if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) {
- 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]);
- } else {
- no_overlap = 1; break;
- }
- if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) {
- 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]);
- } else {
- no_overlap = 1; break;
- }
- // Count the number of cells in this direction in overlap_realbox
- overlap_box.setSmall( dir, 0 );
- overlap_box.setBig( dir,
- int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1);
+ Real wt = amrex::second();
+
+ const Box& tile_box = mfi.tilebox();
+ const RealBox tile_realbox = WarpX::getRealBox(tile_box, lev);
+
+ // Find the cells of part_box that overlap with tile_realbox
+ // If there is no overlap, just go to the next tile in the loop
+ RealBox overlap_realbox;
+ Box overlap_box;
+ IntVect shifted;
+ bool no_overlap = false;
+
+ 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]);
+ } else {
+ no_overlap = true; break;
}
- if (no_overlap == 1) {
- continue; // Go to the next tile
+ 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]);
+ } else {
+ no_overlap = true; break;
}
+ // Count the number of cells in this direction in overlap_realbox
+ overlap_box.setSmall( dir, 0 );
+ overlap_box.setBig( dir,
+ int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))
+ /dx[dir] )) - 1);
+ shifted[dir] = std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir]);
+ // shifted is exact in non-moving-window direction. That's all we care.
+ }
+ if (no_overlap == 1) {
+ continue; // Go to the next tile
+ }
- const int grid_id = mfi.index();
- const int tile_id = mfi.LocalTileIndex();
-
- // Loop through the cells of overlap_box and inject
- // the corresponding particles
- const auto& overlap_corner = overlap_realbox.lo();
- for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv))
- {
- int fac;
- if (do_continuous_injection) {
-#if ( AMREX_SPACEDIM == 3 )
- Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0];
- Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1];
- Real z = overlap_corner[2] + (iv[2] + 0.5)*dx[2];
-#elif ( AMREX_SPACEDIM == 2 )
- Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0];
- Real y = 0;
- Real z = overlap_corner[1] + (iv[1] + 0.5)*dx[1];
-#endif
- fac = GetRefineFac(x, y, z);
- } else {
- fac = 1.0;
- }
-
- int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac);
- for (int i_part=0; i_part<ref_num_ppc;i_part++) {
- std::array<Real, 3> r;
- plasma_injector->getPositionUnitBox(r, i_part, fac);
-#if ( AMREX_SPACEDIM == 3 )
- Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0];
- Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1];
- Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2];
-#elif ( AMREX_SPACEDIM == 2 )
- Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0];
- Real y = 0;
- Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1];
-#endif
- // If the new particle is not inside the tile box,
- // go to the next generated particle.
-#if ( AMREX_SPACEDIM == 3 )
- if(!tile_realbox.contains( RealVect{x, y, z} )) continue;
-#elif ( AMREX_SPACEDIM == 2 )
- if(!tile_realbox.contains( RealVect{x, z} )) continue;
-#endif
+ const int grid_id = mfi.index();
+ const int tile_id = mfi.LocalTileIndex();
- // Save the x and y values to use in the insideBounds checks.
- // This is needed with WARPX_RZ since x and y are modified.
- Real xb = x;
- Real yb = y;
-
-#ifdef WARPX_RZ
- // Replace the x and y, choosing the angle randomly.
- // These x and y are used to get the momentum and density
- Real theta = 2.*MathConst::pi*amrex::Random();
- y = x*std::sin(theta);
- x = x*std::cos(theta);
-#endif
+ // Max number of new particles, if particles are created in the whole
+ // overlap_box. All of them are created, and invalid ones are then
+ // discaded
+ int max_new_particles = overlap_box.numPts() * num_ppc;
- Real dens;
- std::array<Real, 3> u;
- if (WarpX::gamma_boost == 1.){
- // Lab-frame simulation
- // If the particle is not within the species's
- // xmin, xmax, ymin, ymax, zmin, zmax, go to
- // the next generated particle.
- if (!plasma_injector->insideBounds(xb, yb, z)) continue;
- plasma_injector->getMomentum(u, x, y, z);
- dens = plasma_injector->getDensity(x, y, z);
- } else {
- // Boosted-frame simulation
- Real c = PhysConst::c;
- Real gamma_boost = WarpX::gamma_boost;
- Real beta_boost = WarpX::beta_boost;
- // Since the user provides the density distribution
- // at t_lab=0 and in the lab-frame coordinates,
- // we need to find the lab-frame position of this
- // particle at t_lab=0, from its boosted-frame coordinates
- // Assuming ballistic motion, this is given by:
- // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) )
- // where betaz_lab is the speed of the particle in the lab frame
- //
- // In order for this equation to be solvable, betaz_lab
- // is explicitly assumed to have no dependency on z0_lab
- plasma_injector->getMomentum(u, x, y, 0.); // No z0_lab dependency
- // At this point u is the lab-frame momentum
- // => Apply the above formula for z0_lab
- Real gamma_lab = std::sqrt( 1 + (u[0]*u[0] + u[1]*u[1] + u[2]*u[2])/(c*c) );
- Real betaz_lab = u[2]/gamma_lab/c;
- Real t = WarpX::GetInstance().gett_new(lev);
- Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) );
- // If the particle is not within the lab-frame zmin, zmax, etc.
- // go to the next generated particle.
- if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue;
- // call `getDensity` with lab-frame parameters
- dens = plasma_injector->getDensity(x, y, z0_lab);
- // At this point u and dens are the lab-frame quantities
- // => Perform Lorentz transform
- dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab );
- u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab );
- }
- Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac));
-#ifdef WARPX_RZ
- if (plasma_injector->radially_weighted) {
- weight *= 2*MathConst::pi*xb;
- } else {
- // This is not correct since it might shift the particle
- // out of the local grid
- x = std::sqrt(xb*rmax);
- weight *= dx[0];
- }
-#endif
- attribs[PIdx::w ] = weight;
- 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);
- 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]);
- }
-
- AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs);
+ // If refine injection, build pointer dp_cellid that holds pointer to
+ // array of refined cell IDs.
+ Vector<int> cellid_v;
+ if (refine_injection and lev == 0)
+ {
+ // then how many new particles will be injected is not that simple
+ // We have to shift fine_injection_box because overlap_box has been shifted.
+ Box fine_overlap_box = overlap_box & amrex::shift(fine_injection_box,shifted);
+ max_new_particles += fine_overlap_box.numPts() * num_ppc
+ * (AMREX_D_TERM(rrfac,*rrfac,*rrfac)-1);
+ for (int icell = 0, ncells = overlap_box.numPts(); icell < ncells; ++icell) {
+ IntVect iv = overlap_box.atOffset(icell);
+ int r = (fine_overlap_box.contains(iv)) ? AMREX_D_TERM(rrfac,*rrfac,*rrfac) : 1;
+ for (int ipart = 0; ipart < r; ++ipart) {
+ cellid_v.push_back(icell);
+ cellid_v.push_back(ipart);
}
}
+ }
+ int const* hp_cellid = (cellid_v.empty()) ? nullptr : cellid_v.data();
+ amrex::AsyncArray<int> cellid_aa(hp_cellid, cellid_v.size());
+ int const* dp_cellid = cellid_aa.data();
- if (cost) {
- wt = (amrex::second() - wt) / tile_box.d_numPts();
- Array4<Real> const& costarr = cost->array(mfi);
- amrex::ParallelFor(tile_box,
- [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
- {
- costarr(i,j,k) += wt;
- });
- }
+ // Update NextID to include particles created in this function
+ int pid;
+#pragma omp critical (add_plasma_nextid)
+ {
+ pid = ParticleType::NextID();
+ ParticleType::NextID(pid+max_new_particles);
}
- }
-}
+ const int cpuid = ParallelDescriptor::MyProc();
-#ifdef AMREX_USE_GPU
-void
-PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox)
-{
- BL_PROFILE("PhysicalParticleContainer::AddPlasmaGPU");
+ auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)];
+ bool do_boosted = false;
+ if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) {
+ do_boosted = true;
+ DefineAndReturnParticleTile(lev, grid_id, tile_id);
+ }
+ auto old_size = particle_tile.GetArrayOfStructs().size();
+ auto new_size = old_size + max_new_particles;
+ particle_tile.resize(new_size);
+
+ ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size;
+ auto& soa = particle_tile.GetStructOfArrays();
+ GpuArray<Real*,PIdx::nattribs> pa;
+ for (int ia = 0; ia < PIdx::nattribs; ++ia) {
+ pa[ia] = soa.GetRealData(ia).data() + old_size;
+ }
+ GpuArray<Real*,6> pb;
+ if (do_boosted) {
+ pb[0] = soa.GetRealData(particle_comps[ "xold"]).data() + old_size;
+ pb[1] = soa.GetRealData(particle_comps[ "yold"]).data() + old_size;
+ pb[2] = soa.GetRealData(particle_comps[ "zold"]).data() + old_size;
+ pb[3] = soa.GetRealData(particle_comps["uxold"]).data() + old_size;
+ pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size;
+ pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size;
+ }
- // If no part_realbox is provided, initialize particles in the whole domain
- const Geometry& geom = Geom(lev);
- if (!part_realbox.ok()) part_realbox = geom.ProbDomain();
+ const GpuArray<Real,AMREX_SPACEDIM> overlap_corner
+ {AMREX_D_DECL(overlap_realbox.lo(0),
+ overlap_realbox.lo(1),
+ overlap_realbox.lo(2))};
- int num_ppc = plasma_injector->num_particles_per_cell;
-#ifdef WARPX_RZ
- Real rmax = std::min(plasma_injector->xmax, part_realbox.hi(0));
-#endif
+ std::size_t shared_mem_bytes = plasma_injector->sharedMemoryNeeded();
+ int lrrfac = rrfac;
- const Real* dx = geom.CellSize();
+ // Loop over all new particles and inject them (creates too many
+ // particles, in particular does not consider xmin, xmax etc.).
+ // The invalid ones are given negative ID and are deleted during the
+ // next redistribute.
+ amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept
+ {
+ ParticleType& p = pp[ip];
+ p.id() = pid+ip;
+ p.cpu() = cpuid;
+
+ int cellid, i_part;
+ Real fac;
+ if (dp_cellid == nullptr) {
+ cellid = ip/num_ppc;
+ i_part = ip - cellid*num_ppc;
+ fac = 1.0;
+ } else {
+ cellid = dp_cellid[2*ip];
+ i_part = dp_cellid[2*ip+1];
+ fac = lrrfac;
+ }
- Real scale_fac;
-#if AMREX_SPACEDIM==3
- scale_fac = dx[0]*dx[1]*dx[2]/num_ppc;
-#elif AMREX_SPACEDIM==2
- scale_fac = dx[0]*dx[1]/num_ppc;
-#endif
+ IntVect iv = overlap_box.atOffset(cellid);
-#ifdef _OPENMP
- // First touch all tiles in the map in serial
- for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
- const int grid_id = mfi.index();
- const int tile_id = mfi.LocalTileIndex();
- GetParticles(lev)[std::make_pair(grid_id, tile_id)];
- }
+ const XDim3 r = inj_pos->getPositionUnitBox(i_part, fac);
+#if (AMREX_SPACEDIM == 3)
+ Real x = overlap_corner[0] + (iv[0]+r.x)*dx[0];
+ Real y = overlap_corner[1] + (iv[1]+r.y)*dx[1];
+ Real z = overlap_corner[2] + (iv[2]+r.z)*dx[2];
+#else
+ Real x = overlap_corner[0] + (iv[0]+r.x)*dx[0];
+ Real y = 0.0;
+ Real z = overlap_corner[1] + (iv[1]+r.y)*dx[1];
#endif
- MultiFab* cost = WarpX::getCosts(lev);
-
- if ( (not m_refined_injection_mask) and WarpX::do_moving_window)
- {
- Box mask_box = geom.Domain();
- mask_box.setSmall(WarpX::moving_window_dir, 0);
- mask_box.setBig(WarpX::moving_window_dir, 0);
- m_refined_injection_mask.reset( new IArrayBox(mask_box));
- m_refined_injection_mask->setVal(-1);
- }
-
- MFItInfo info;
- if (do_tiling) {
- info.EnableTiling(tile_size);
- }
- info.SetDynamic(true);
-
-#ifdef _OPENMP
-#pragma omp parallel if (not WarpX::serialize_ics)
+#if (AMREX_SPACEDIM == 3)
+ if (!tile_realbox.contains(XDim3{x,y,z})) {
+ p.id() = -1;
+ return;
+ }
+#else
+ if (!tile_realbox.contains(XDim3{x,z,0.0})) {
+ p.id() = -1;
+ return;
+ }
#endif
- {
- std::array<Real,PIdx::nattribs> attribs;
- attribs.fill(0.0);
-
- // Loop through the tiles
- for (MFIter mfi = MakeMFIter(lev, info); mfi.isValid(); ++mfi) {
- Real wt = amrex::second();
-
- const Box& tile_box = mfi.tilebox();
- const RealBox tile_realbox = WarpX::getRealBox(tile_box, lev);
-
- // Find the cells of part_box that overlap with tile_realbox
- // If there is no overlap, just go to the next tile in the loop
- RealBox overlap_realbox;
- Box overlap_box;
- Real ncells_adjust;
- bool no_overlap = 0;
+ // Save the x and y values to use in the insideBounds checks.
+ // This is needed with WARPX_DIM_RZ since x and y are modified.
+ Real xb = x;
+ Real yb = y;
+
+#ifdef WARPX_DIM_RZ
+ // Replace the x and y, choosing the angle randomly.
+ // These x and y are used to get the momentum and density
+ Real theta = 2.*MathConst::pi*amrex::Random();
+ x = xb*std::cos(theta);
+ y = xb*std::sin(theta);
+#endif
- for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
- if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) {
- 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]);
- } else {
- no_overlap = 1; break;
+ Real dens;
+ XDim3 u;
+ if (gamma_boost == 1.) {
+ // Lab-frame simulation
+ // If the particle is not within the species's
+ // xmin, xmax, ymin, ymax, zmin, zmax, go to
+ // the next generated particle.
+ if (!inj_pos->insideBounds(xb, yb, z)) {
+ p.id() = -1;
+ return;
}
- if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) {
- 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]);
- } else {
- no_overlap = 1; break;
+ u = inj_mom->getMomentum(x, y, z);
+ dens = inj_rho->getDensity(x, y, z);
+ // Remove particle if density below threshold
+ if ( dens < density_min ){
+ p.id() = -1;
+ return;
}
- // Count the number of cells in this direction in overlap_realbox
- overlap_box.setSmall( dir, 0 );
- overlap_box.setBig( dir,
- int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1);
- }
- if (no_overlap == 1) {
- continue; // Go to the next tile
- }
-
- const int grid_id = mfi.index();
- const int tile_id = mfi.LocalTileIndex();
-
- Cuda::HostVector<ParticleType> host_particles;
- std::array<Cuda::HostVector<Real>, PIdx::nattribs> host_attribs;
-
- // Loop through the cells of overlap_box and inject
- // the corresponding particles
- const auto& overlap_corner = overlap_realbox.lo();
- for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv))
- {
- int fac;
- if (do_continuous_injection) {
-#if ( AMREX_SPACEDIM == 3 )
- Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0];
- Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1];
- Real z = overlap_corner[2] + (iv[2] + 0.5)*dx[2];
-#elif ( AMREX_SPACEDIM == 2 )
- Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0];
- Real y = 0;
- Real z = overlap_corner[1] + (iv[1] + 0.5)*dx[1];
-#endif
- fac = GetRefineFac(x, y, z);
- } else {
- fac = 1.0;
+ // Cut density if above threshold
+ dens = amrex::min(dens, density_max);
+ } else {
+ // Boosted-frame simulation
+ // Since the user provides the density distribution
+ // at t_lab=0 and in the lab-frame coordinates,
+ // we need to find the lab-frame position of this
+ // particle at t_lab=0, from its boosted-frame coordinates
+ // Assuming ballistic motion, this is given by:
+ // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) )
+ // where betaz_lab is the speed of the particle in the lab frame
+ //
+ // In order for this equation to be solvable, betaz_lab
+ // is explicitly assumed to have no dependency on z0_lab
+ u = inj_mom->getMomentum(x, y, 0.); // No z0_lab dependency
+ // At this point u is the lab-frame momentum
+ // => Apply the above formula for z0_lab
+ Real gamma_lab = std::sqrt( 1.+(u.x*u.x+u.y*u.y+u.z*u.z) );
+ Real betaz_lab = u.z/(gamma_lab);
+ Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab)
+ - PhysConst::c*t*(betaz_lab-beta_boost) );
+ // If the particle is not within the lab-frame zmin, zmax, etc.
+ // go to the next generated particle.
+ if (!inj_pos->insideBounds(xb, yb, z0_lab)) {
+ p.id() = -1;
+ return;
}
+ // call `getDensity` with lab-frame parameters
+ dens = inj_rho->getDensity(x, y, z0_lab);
+ // Remove particle if density below threshold
+ if ( dens < density_min ){
+ p.id() = -1;
+ return;
+ }
+ // Cut density if above threshold
+ dens = amrex::min(dens, density_max);
+ // At this point u and dens are the lab-frame quantities
+ // => Perform Lorentz transform
+ dens = gamma_boost * dens * ( 1.0 - beta_boost*betaz_lab );
+ u.z = gamma_boost * ( u.z -beta_boost*gamma_lab );
+ }
- int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac);
- for (int i_part=0; i_part<ref_num_ppc;i_part++) {
- std::array<Real, 3> r;
- plasma_injector->getPositionUnitBox(r, i_part, fac);
-#if ( AMREX_SPACEDIM == 3 )
- Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0];
- Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1];
- Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2];
-#elif ( AMREX_SPACEDIM == 2 )
- Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0];
- Real y = 0;
- Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1];
-#endif
- // If the new particle is not inside the tile box,
- // go to the next generated particle.
-#if ( AMREX_SPACEDIM == 3 )
- if(!tile_realbox.contains( RealVect{x, y, z} )) continue;
-#elif ( AMREX_SPACEDIM == 2 )
- if(!tile_realbox.contains( RealVect{x, z} )) continue;
-#endif
-
- // Save the x and y values to use in the insideBounds checks.
- // This is needed with WARPX_RZ since x and y are modified.
- Real xb = x;
- Real yb = y;
-
-#ifdef WARPX_RZ
- // Replace the x and y, choosing the angle randomly.
- // These x and y are used to get the momentum and density
- Real theta = 2.*MathConst::pi*amrex::Random();
- x = xb*std::cos(theta);
- y = xb*std::sin(theta);
-#endif
+ u.x *= PhysConst::c;
+ u.y *= PhysConst::c;
+ u.z *= PhysConst::c;
- Real dens;
- std::array<Real, 3> u;
- if (WarpX::gamma_boost == 1.){
- // Lab-frame simulation
- // If the particle is not within the species's
- // xmin, xmax, ymin, ymax, zmin, zmax, go to
- // the next generated particle.
- if (!plasma_injector->insideBounds(xb, yb, z)) continue;
- plasma_injector->getMomentum(u, x, y, z);
- dens = plasma_injector->getDensity(x, y, z);
- } else {
- // Boosted-frame simulation
- Real c = PhysConst::c;
- Real gamma_boost = WarpX::gamma_boost;
- Real beta_boost = WarpX::beta_boost;
- // Since the user provides the density distribution
- // at t_lab=0 and in the lab-frame coordinates,
- // we need to find the lab-frame position of this
- // particle at t_lab=0, from its boosted-frame coordinates
- // Assuming ballistic motion, this is given by:
- // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) )
- // where betaz_lab is the speed of the particle in the lab frame
- //
- // In order for this equation to be solvable, betaz_lab
- // is explicitly assumed to have no dependency on z0_lab
- plasma_injector->getMomentum(u, x, y, 0.); // No z0_lab dependency
- // At this point u is the lab-frame momentum
- // => Apply the above formula for z0_lab
- Real gamma_lab = std::sqrt( 1 + (u[0]*u[0] + u[1]*u[1] + u[2]*u[2])/(c*c) );
- Real betaz_lab = u[2]/gamma_lab/c;
- Real t = WarpX::GetInstance().gett_new(lev);
- Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) );
- // If the particle is not within the lab-frame zmin, zmax, etc.
- // go to the next generated particle.
- if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue;
- // call `getDensity` with lab-frame parameters
- dens = plasma_injector->getDensity(x, y, z0_lab);
- // At this point u and dens are the lab-frame quantities
- // => Perform Lorentz transform
- dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab );
- u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab );
- }
- Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac));
-#ifdef WARPX_RZ
- if (plasma_injector->radially_weighted) {
- weight *= 2*MathConst::pi*xb;
- } else {
- // This is not correct since it might shift the particle
- // out of the local grid
- x = std::sqrt(xb*rmax);
- weight *= dx[0];
- }
+ // Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac));
+ Real weight = dens * scale_fac;
+#ifdef WARPX_DIM_RZ
+ if (radially_weighted) {
+ weight *= 2.*MathConst::pi*xb;
+ } else {
+ // This is not correct since it might shift the particle
+ // out of the local grid
+ x = std::sqrt(xb*rmax);
+ weight *= dx[0];
+ }
#endif
- attribs[PIdx::w ] = weight;
- attribs[PIdx::ux] = u[0];
- attribs[PIdx::uy] = u[1];
- attribs[PIdx::uz] = u[2];
-
- // note - this will be slow on the GPU, need to revisit
- if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
- {
- auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id);
- 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]);
- }
+ pa[PIdx::w ][ip] = weight;
+ pa[PIdx::ux][ip] = u.x;
+ pa[PIdx::uy][ip] = u.y;
+ pa[PIdx::uz][ip] = u.z;
+
+ if (do_boosted) {
+ pb[0][ip] = x;
+ pb[1][ip] = y;
+ pb[2][ip] = z;
+ pb[3][ip] = u.x;
+ pb[4][ip] = u.y;
+ pb[5][ip] = u.z;
+ }
- ParticleType p;
- p.id() = ParticleType::NextID();
- p.cpu() = ParallelDescriptor::MyProc();
#if (AMREX_SPACEDIM == 3)
- p.pos(0) = x;
- p.pos(1) = y;
- p.pos(2) = z;
+ p.pos(0) = x;
+ p.pos(1) = y;
+ p.pos(2) = z;
#elif (AMREX_SPACEDIM == 2)
-#ifdef WARPX_RZ
- attribs[PIdx::theta] = theta;
+#ifdef WARPX_DIM_RZ
+ pa[PIdx::theta][ip] = theta;
#endif
- p.pos(0) = xb;
- p.pos(1) = z;
+ p.pos(0) = xb;
+ p.pos(1) = z;
#endif
-
- host_particles.push_back(p);
- for (int kk = 0; kk < PIdx::nattribs; ++kk)
- host_attribs[kk].push_back(attribs[kk]);
- }
- }
-
- auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)];
- auto old_size = particle_tile.GetArrayOfStructs().size();
- auto new_size = old_size + host_particles.size();
- particle_tile.resize(new_size);
-
- Cuda::thrust_copy(host_particles.begin(),
- host_particles.end(),
- particle_tile.GetArrayOfStructs().begin() + old_size);
-
- for (int kk = 0; kk < PIdx::nattribs; ++kk) {
- Cuda::thrust_copy(host_attribs[kk].begin(),
- host_attribs[kk].end(),
- particle_tile.GetStructOfArrays().GetRealData(kk).begin() + old_size);
- }
-
- if (cost) {
- wt = (amrex::second() - wt) / tile_box.d_numPts();
- Array4<Real> const& costarr = cost->array(mfi);
- amrex::ParallelFor(tile_box,
- [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
- {
- costarr(i,j,k) += wt;
- });
- }
- }
+ }, shared_mem_bytes);
+
+ if (cost) {
+ wt = (amrex::second() - wt) / tile_box.d_numPts();
+ Array4<Real> const& costarr = cost->array(mfi);
+ amrex::ParallelFor(tile_box,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ {
+ costarr(i,j,k) += wt;
+ });
+ }
}
+
+ // The function that calls this is responsible for redistributing particles.
}
-#endif
#ifdef WARPX_DO_ELECTROSTATIC
void
@@ -1066,11 +848,14 @@ PhysicalParticleContainer::FieldGather (int lev,
MultiFab* cost = WarpX::getCosts(lev);
#ifdef _OPENMP
-#pragma omp parallel
+#pragma omp parallel
#endif
{
- Cuda::ManagedDeviceVector<Real> xp, yp, zp;
-
+#ifdef _OPENMP
+ int thread_num = omp_get_thread_num();
+#else
+ int thread_num = 0;
+#endif
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
Real wt = amrex::second();
@@ -1106,35 +891,15 @@ PhysicalParticleContainer::FieldGather (int lev,
//
// copy data from particle container to temp arrays
//
- pti.GetPosition(xp, yp, zp);
-
- const std::array<Real,3>& xyzmin = WarpX::LowerCorner(box, lev);
- const int* ixyzmin = box.loVect();
+ pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]);
//
// Field Gather
//
- const int ll4symtry = false;
- long lvect_fieldgathe = 64;
- warpx_geteb_energy_conserving(
- &np,
- xp.dataPtr(),
- yp.dataPtr(),
- zp.dataPtr(),
- Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(),
- Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(),
- ixyzmin,
- &xyzmin[0], &xyzmin[1], &xyzmin[2],
- &dx[0], &dx[1], &dx[2],
- &WarpX::nox, &WarpX::noy, &WarpX::noz,
- BL_TO_FORTRAN_ANYD(exfab),
- BL_TO_FORTRAN_ANYD(eyfab),
- BL_TO_FORTRAN_ANYD(ezfab),
- BL_TO_FORTRAN_ANYD(bxfab),
- BL_TO_FORTRAN_ANYD(byfab),
- BL_TO_FORTRAN_ANYD(bzfab),
- &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal,
- &lvect_fieldgathe, &WarpX::field_gathering_algo);
+ int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
+ FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+ &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
+ Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev);
if (cost) {
const Box& tbx = pti.tilebox();
@@ -1164,7 +929,7 @@ PhysicalParticleContainer::Evolve (int lev,
BL_PROFILE("PPC::Evolve()");
BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy);
BL_PROFILE_VAR_NS("PICSAR::FieldGather", blp_pxr_fg);
- BL_PROFILE_VAR_NS("PICSAR::ParticlePush", blp_pxr_pp);
+ 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);
@@ -1391,57 +1156,40 @@ PhysicalParticleContainer::Evolve (int lev,
pti.GetPosition(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, 0, np_current, np, thread_num, lev);
+ if (rho) {
+ DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev);
+ if (has_buffer){
+ DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1);
+ }
+ }
if (! do_not_push)
{
+ const long np_gather = (cEx) ? nfine_gather : np;
+
+ int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
+
//
// Field Gather of Aux Data (i.e., the full solution)
//
- const int ll4symtry = false;
- long lvect_fieldgathe = 64;
-
- 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);
-
- warpx_geteb_energy_conserving(
- &np_gather,
- m_xp[thread_num].dataPtr(),
- m_yp[thread_num].dataPtr(),
- m_zp[thread_num].dataPtr(),
- Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(),
- Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(),
- ixyzmin_grid,
- &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2],
- &dx[0], &dx[1], &dx[2],
- &WarpX::nox, &WarpX::noy, &WarpX::noz,
- BL_TO_FORTRAN_ANYD(*exfab),
- BL_TO_FORTRAN_ANYD(*eyfab),
- BL_TO_FORTRAN_ANYD(*ezfab),
- BL_TO_FORTRAN_ANYD(*bxfab),
- BL_TO_FORTRAN_ANYD(*byfab),
- BL_TO_FORTRAN_ANYD(*bzfab),
- &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal,
- &lvect_fieldgathe, &WarpX::field_gathering_algo);
+ FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+ exfab, eyfab, ezfab, bxfab, byfab, bzfab,
+ Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev);
if (np_gather < np)
{
const IntVect& ref_ratio = WarpX::RefRatio(lev-1);
const Box& cbox = amrex::coarsen(box,ref_ratio);
- const std::array<Real,3>& cxyzmin_grid = WarpX::LowerCorner(cbox, lev-1);
- const int* cixyzmin_grid = cbox.loVect();
-
- const FArrayBox* cexfab = &(*cEx)[pti];
- const FArrayBox* ceyfab = &(*cEy)[pti];
- const FArrayBox* cezfab = &(*cEz)[pti];
- const FArrayBox* cbxfab = &(*cBx)[pti];
- const FArrayBox* cbyfab = &(*cBy)[pti];
- const FArrayBox* cbzfab = &(*cBz)[pti];
+ // Data on the grid
+ FArrayBox const* cexfab = &(*cEx)[pti];
+ FArrayBox const* ceyfab = &(*cEy)[pti];
+ FArrayBox const* cezfab = &(*cEz)[pti];
+ FArrayBox const* cbxfab = &(*cBx)[pti];
+ FArrayBox const* cbyfab = &(*cBy)[pti];
+ FArrayBox const* cbzfab = &(*cBz)[pti];
+
if (WarpX::use_fdtd_nci_corr)
{
#if (AMREX_SPACEDIM == 2)
@@ -1494,26 +1242,14 @@ PhysicalParticleContainer::Evolve (int lev,
#endif
}
- long ncrse = np - nfine_gather;
- warpx_geteb_energy_conserving(
- &ncrse,
- m_xp[thread_num].dataPtr()+nfine_gather,
- m_yp[thread_num].dataPtr()+nfine_gather,
- m_zp[thread_num].dataPtr()+nfine_gather,
- Exp.dataPtr()+nfine_gather, Eyp.dataPtr()+nfine_gather, Ezp.dataPtr()+nfine_gather,
- Bxp.dataPtr()+nfine_gather, Byp.dataPtr()+nfine_gather, Bzp.dataPtr()+nfine_gather,
- cixyzmin_grid,
- &cxyzmin_grid[0], &cxyzmin_grid[1], &cxyzmin_grid[2],
- &cdx[0], &cdx[1], &cdx[2],
- &WarpX::nox, &WarpX::noy, &WarpX::noz,
- BL_TO_FORTRAN_ANYD(*cexfab),
- BL_TO_FORTRAN_ANYD(*ceyfab),
- BL_TO_FORTRAN_ANYD(*cezfab),
- BL_TO_FORTRAN_ANYD(*cbxfab),
- BL_TO_FORTRAN_ANYD(*cbyfab),
- BL_TO_FORTRAN_ANYD(*cbzfab),
- &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal,
- &lvect_fieldgathe, &WarpX::field_gathering_algo);
+ // Field gather for particles in gather buffers
+ e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal();
+ FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+ cexfab, ceyfab, cezfab,
+ cbxfab, cbyfab, cbzfab,
+ cEx->nGrow(), e_is_nodal,
+ nfine_gather, np-nfine_gather,
+ thread_num, lev, lev-1);
}
BL_PROFILE_VAR_STOP(blp_pxr_fg);
@@ -1521,10 +1257,10 @@ PhysicalParticleContainer::Evolve (int lev,
//
// Particle Push
//
- BL_PROFILE_VAR_START(blp_pxr_pp);
+ BL_PROFILE_VAR_START(blp_ppc_pp);
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);
+ BL_PROFILE_VAR_STOP(blp_ppc_pp);
//
// Current Deposition
@@ -1561,7 +1297,12 @@ PhysicalParticleContainer::Evolve (int lev,
BL_PROFILE_VAR_STOP(blp_copy);
}
- if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev);
+ if (rho) {
+ DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev);
+ if (has_buffer){
+ DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1);
+ }
+ }
if (cost) {
const Box& tbx = pti.tilebox();
@@ -1742,36 +1483,52 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
Real dt)
{
+ // 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 gi = giv.dataPtr();
+ Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+ Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+ Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+ const Real* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr();
+ const Real* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr();
+ const Real* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr();
+ const Real* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr();
+ const Real* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
+ const Real* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
+
if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
{
- 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");
+ };
}
void
@@ -1800,9 +1557,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
auto& attribs = pti.GetAttribs();
- auto& uxp = attribs[PIdx::ux];
- auto& uyp = attribs[PIdx::uy];
- auto& uzp = attribs[PIdx::uz];
auto& Exp = attribs[PIdx::Ex];
auto& Eyp = attribs[PIdx::Ey];
auto& Ezp = attribs[PIdx::Ez];
@@ -1834,42 +1588,44 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
//
pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]);
- const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev);
- const int* ixyzmin_grid = box.loVect();
-
- const int ll4symtry = false;
- long lvect_fieldgathe = 64;
-
- warpx_geteb_energy_conserving(
- &np,
- m_xp[thread_num].dataPtr(),
- m_yp[thread_num].dataPtr(),
- m_zp[thread_num].dataPtr(),
- Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(),
- Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(),
- ixyzmin_grid,
- &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2],
- &dx[0], &dx[1], &dx[2],
- &WarpX::nox, &WarpX::noy, &WarpX::noz,
- BL_TO_FORTRAN_ANYD(exfab),
- BL_TO_FORTRAN_ANYD(eyfab),
- BL_TO_FORTRAN_ANYD(ezfab),
- BL_TO_FORTRAN_ANYD(bxfab),
- BL_TO_FORTRAN_ANYD(byfab),
- BL_TO_FORTRAN_ANYD(bzfab),
- &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal,
- &lvect_fieldgathe, &WarpX::field_gathering_algo);
-
- warpx_particle_pusher_momenta(&np,
- m_xp[thread_num].dataPtr(),
- m_yp[thread_num].dataPtr(),
- m_zp[thread_num].dataPtr(),
- uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(),
- m_giv[thread_num].dataPtr(),
- Exp.dataPtr(), Eyp.dataPtr(), Ezp.dataPtr(),
- Bxp.dataPtr(), Byp.dataPtr(), Bzp.dataPtr(),
- &this->charge, &this->mass, &dt,
- &WarpX::particle_pusher_algo);
+ int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
+ FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+ &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
+ Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev);
+
+ // This wraps the momentum advance so that inheritors can modify the call.
+ // Extract pointers to the different particle quantities
+ Real* const AMREX_RESTRICT gi = m_giv[thread_num].dataPtr();
+ Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+ Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+ Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+ const Real* const AMREX_RESTRICT Expp = Exp.dataPtr();
+ const Real* const AMREX_RESTRICT Eypp = Eyp.dataPtr();
+ const Real* const AMREX_RESTRICT Ezpp = Ezp.dataPtr();
+ const Real* const AMREX_RESTRICT Bxpp = Bxp.dataPtr();
+ const Real* const AMREX_RESTRICT Bypp = Byp.dataPtr();
+ const Real* const AMREX_RESTRICT Bzpp = Bzp.dataPtr();
+
+ // Loop over the particles and update their momentum
+ const Real q = this->charge;
+ 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],
+ Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
+ amrex::ParallelFor( pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i],
+ Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ }
+ );
+ } else {
+ amrex::Abort("Unknown particle pusher");
+ };
}
}
}
@@ -2034,74 +1790,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
}
}
-int PhysicalParticleContainer::GetRefineFac(const Real x, const Real y, const Real z)
-{
- if (finestLevel() == 0) return 1;
- if (not WarpX::refine_plasma) return 1;
-
- IntVect iv;
- const Geometry& geom = Geom(0);
-
- std::array<Real, 3> offset;
-
-#if ( AMREX_SPACEDIM == 3)
- offset[0] = geom.ProbLo(0);
- offset[1] = geom.ProbLo(1);
- offset[2] = geom.ProbLo(2);
-#elif ( AMREX_SPACEDIM == 2 )
- offset[0] = geom.ProbLo(0);
- offset[1] = 0.0;
- offset[2] = geom.ProbLo(1);
-#endif
-
- AMREX_D_TERM(iv[0]=static_cast<int>(floor((x-offset[0])*geom.InvCellSize(0)));,
- iv[1]=static_cast<int>(floor((y-offset[1])*geom.InvCellSize(1)));,
- iv[2]=static_cast<int>(floor((z-offset[2])*geom.InvCellSize(2))););
-
- iv += geom.Domain().smallEnd();
-
- const int dir = WarpX::moving_window_dir;
-
- IntVect iv2 = iv;
- iv2[dir] = 0;
-
- if ( (*m_refined_injection_mask)(iv2) != -1) return (*m_refined_injection_mask)(iv2);
-
- int ref_fac = 1;
- for (int lev = 0; lev < finestLevel(); ++lev)
- {
- const IntVect rr = m_gdb->refRatio(lev);
- const BoxArray& fine_ba = this->ParticleBoxArray(lev+1);
- const int num_boxes = fine_ba.size();
- Vector<Box> stretched_boxes;
- const int safety_factor = 4;
- for (int i = 0; i < num_boxes; ++i)
- {
- Box bx = fine_ba[i];
- bx.coarsen(ref_fac*rr[dir]);
- bx.setSmall(dir, std::numeric_limits<int>::min()/safety_factor);
- bx.setBig(dir, std::numeric_limits<int>::max()/safety_factor);
- stretched_boxes.push_back(bx);
- }
-
- BoxArray stretched_ba(stretched_boxes.dataPtr(), stretched_boxes.size());
-
- const int num_ghost = 0;
- if ( stretched_ba.intersects(Box(iv, iv), num_ghost) )
- {
- ref_fac *= rr[dir];
- }
- else
- {
- break;
- }
- }
-
- (*m_refined_injection_mask)(iv2) = ref_fac;
-
- return ref_fac;
-}
-
/* \brief Inject particles during the simulation
* \param injection_box: domain where particles should be injected.
*/
@@ -2112,3 +1800,134 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box)
const int lev=0;
AddPlasma(lev, injection_box);
}
+
+/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab,
+ * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp.
+ * \param Exp-Bzp: fields on particles.
+ * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti
+ * \param ngE: number of guard cells for E
+ * \param e_is_nodal: 0 if E is staggered, 1 if E is nodal
+ * \param offset: index of first particle for which fields are gathered
+ * \param np_to_gather: number of particles onto which fields are gathered
+ * \param thread_num: if using OpenMP, thread number
+ * \param lev: level on which particles are located
+ * \param gather_lev: level from which particles gather fields (lev-1) for
+ particles in buffers.
+ */
+void
+PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
+ RealVector& Exp,
+ RealVector& Eyp,
+ RealVector& Ezp,
+ RealVector& Bxp,
+ RealVector& Byp,
+ RealVector& Bzp,
+ FArrayBox const * exfab,
+ FArrayBox const * eyfab,
+ FArrayBox const * ezfab,
+ FArrayBox const * bxfab,
+ FArrayBox const * byfab,
+ FArrayBox const * bzfab,
+ const int ngE, const int e_is_nodal,
+ const long offset,
+ const long np_to_gather,
+ int thread_num,
+ int lev,
+ int gather_lev)
+{
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) ||
+ (gather_lev==(lev )),
+ "Gather buffers only work for lev-1");
+
+ // If no particles, do not do anything
+ if (np_to_gather == 0) return;
+ // Get cell size on gather_lev
+ const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0));
+ // Set staggering shift depending on e_is_nodal
+ const Real stagger_shift = e_is_nodal ? 0.0 : 0.5;
+
+ // Get box from which field is gathered.
+ // If not gathering from the finest level, the box is coarsened.
+ Box box;
+ if (lev == gather_lev) {
+ box = pti.tilebox();
+ } else {
+ const IntVect& ref_ratio = WarpX::RefRatio(gather_lev);
+ box = amrex::coarsen(pti.tilebox(),ref_ratio);
+ }
+
+ // Add guard cells to the box.
+ box.grow(ngE);
+
+ const Array4<const Real>& ex_arr = exfab->array();
+ const Array4<const Real>& ey_arr = eyfab->array();
+ const Array4<const Real>& ez_arr = ezfab->array();
+ const Array4<const Real>& bx_arr = bxfab->array();
+ 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;
+
+ // Lower corner of tile box physical domain
+ const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev);
+
+ const Dim3 lo = lbound(box);
+
+ // Depending on l_lower_in_v and WarpX::nox, call
+ // different versions of template function doGatherShapeN
+ if (WarpX::l_lower_order_in_v){
+ if (WarpX::nox == 1){
+ doGatherShapeN<1,1>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
+ ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+ np_to_gather, dx,
+ xyzmin, lo, stagger_shift);
+ } else if (WarpX::nox == 2){
+ doGatherShapeN<2,1>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
+ ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+ np_to_gather, dx,
+ xyzmin, lo, stagger_shift);
+ } else if (WarpX::nox == 3){
+ doGatherShapeN<3,1>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
+ ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+ np_to_gather, dx,
+ xyzmin, lo, stagger_shift);
+ }
+ } else {
+ if (WarpX::nox == 1){
+ doGatherShapeN<1,0>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
+ ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+ np_to_gather, dx,
+ xyzmin, lo, stagger_shift);
+ } else if (WarpX::nox == 2){
+ doGatherShapeN<2,0>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
+ ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+ np_to_gather, dx,
+ xyzmin, lo, stagger_shift);
+ } else if (WarpX::nox == 3){
+ doGatherShapeN<3,0>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
+ ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+ np_to_gather, dx,
+ xyzmin, lo, stagger_shift);
+ }
+ }
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 08:16:41 +0000