Merged in boosted_injection (pull request #44)

Generalize plasma injection, for boosted-frame simulations

1 files changed, 93 insertions, 23 deletions

diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp
index c3bad5fa7..e5f6b0c82 100644
--- a/Source/PhysicalParticleContainer.cpp
+++ b/Source/PhysicalParticleContainer.cpp
@@ -71,7 +71,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
 }
 
 void
-PhysicalParticleContainer::AddParticles (int lev, Box part_box)
+PhysicalParticleContainer::AddParticles (int lev)
 {
     BL_PROFILE("PhysicalParticleContainer::AddParticles()");
 
@@ -100,20 +100,38 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box)
         return;
     }
 
-    if ( not plasma_injector->doInjection() ) return;
+    if ( plasma_injector->doInjection() ) {
+        AddPlasma( lev );
+    }
+}
+
+/**
+ * Create new macroparticles for this species, with a fixed
+ * number of particles per cell (in the cells of `part_realbox`).
+ * The new particles are only created inside the intersection of `part_realbox`
+ * with the local grid for the current proc.
+ * @param lev the index of the refinement level
+ * @param part_realbox the box in which new particles should be created
+ * (this box should correspond to an integer number of cells in each direction,
+ * but its boundaries need not be aligned with the actual cells of the simulation)
+ */
+void
+PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox )
+{
 
+    // If no part_realbox is provided, initialize particles in the whole domain
     const Geometry& geom = Geom(lev);
-    if (!part_box.ok()) part_box = geom.Domain();
+    if (!part_realbox.ok()) part_realbox = geom.ProbDomain();
 
     int num_ppc = plasma_injector->num_particles_per_cell;
 
-    const std::array<Real,3>& dx = WarpX::CellSize(lev);
+    const Real* dx = geom.CellSize();
 
     Real scale_fac;
 #if BL_SPACEDIM==3
     scale_fac = dx[0]*dx[1]*dx[2]/num_ppc;
 #elif BL_SPACEDIM==2
-    scale_fac = dx[0]*dx[2]/num_ppc;
+    scale_fac = dx[0]*dx[1]/num_ppc;
 #endif
 
 #ifdef _OPENMP
@@ -132,38 +150,90 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box)
         std::array<Real,PIdx::nattribs> attribs;
         attribs.fill(0.0);
 
+        // Loop through the tiles
         for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
-            const Box& tile_box  = mfi.tilebox();
-            const Box& intersectBox = tile_box & part_box;
-            if (!intersectBox.ok()) continue;
 
-            const std::array<Real, 3>& tile_corner =
-                WarpX::LowerCorner(intersectBox, lev);
+            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<BL_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);
+            }
+            if (no_overlap == 1) continue; // Go to the next tile
 
             const int grid_id = mfi.index();
             const int tile_id = mfi.LocalTileIndex();
 
-            const auto& boxlo = intersectBox.smallEnd();
-            for (IntVect iv = intersectBox.smallEnd();
-                 iv <= intersectBox.bigEnd(); intersectBox.next(iv)) {
+            // 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)) {
                 for (int i_part=0; i_part<num_ppc;i_part++) {
                     std::array<Real, 3> r;
                     plasma_injector->getPositionUnitBox(r, i_part);
 #if ( BL_SPACEDIM == 3 )
-                    Real x = tile_corner[0] + (iv[0]-boxlo[0] + r[0])*dx[0];
-                    Real y = tile_corner[1] + (iv[1]-boxlo[1] + r[1])*dx[1];
-                    Real z = tile_corner[2] + (iv[2]-boxlo[2] + r[2])*dx[2];
+                    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 ( BL_SPACEDIM == 2 )
-                    Real x = tile_corner[0] + (iv[0]-boxlo[0] + r[0])*dx[0];
-                    Real y = 0.;
-                    Real z = tile_corner[2] + (iv[1]-boxlo[1] + r[2])*dx[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 ( BL_SPACEDIM == 3 )
+                    if(!tile_realbox.contains( RealVect{x, y, z} )) continue;
+#elif ( BL_SPACEDIM == 2 )
+                    if(!tile_realbox.contains( RealVect{x, z} )) continue;
+#endif
+
                     if (plasma_injector->insideBounds(x, y, z)) {
-                        Real weight;
+                        Real dens;
                         std::array<Real, 3> u;
-                        plasma_injector->getMomentum(u, x, y, z);
-                        weight = plasma_injector->getDensity(x, y, z) * scale_fac;
-                        attribs[PIdx::w ] = weight;
+                        if (WarpX::gamma_boost == 1.){
+                            // Lab-frame simulation
+                            plasma_injector->getMomentum(u, x, y, z);
+                            dens = plasma_injector->getDensity(x, y, z);
+                        } else {
+                            // Boosted-frame simulation: call `getMomentum`
+                            // and `getDensity` with lab-frame Parameters
+                            // (Assumes that the plasma has a low velocity,
+                            // and that the boost is along z)
+                            Real t = WarpX::GetInstance().gett_new(lev);
+                            Real v_boost = WarpX::beta_boost*PhysConst::c;
+                            Real z_lab = WarpX::gamma_boost*( z - v_boost*t );
+                            plasma_injector->getMomentum(u, x, y, z_lab);
+                            dens = plasma_injector->getDensity(x, y, z);
+                            // Perform Lorentz transform
+                            // (Assumes that the plasma has a low velocity)
+                            u[2] = WarpX::gamma_boost * ( u[2] - v_boost );
+                            dens = WarpX::gamma_boost * dens;
+                        }
+                        attribs[PIdx::w ] = dens * scale_fac;
                         attribs[PIdx::ux] = u[0];
                         attribs[PIdx::uy] = u[1];
                         attribs[PIdx::uz] = u[2];