Implement ability to inject fresh plasma from C++

1 files changed, 91 insertions, 2 deletions

diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index 76c928621..76e5f52e6 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -307,6 +307,95 @@ WarpX::MoveWindow ()
   shiftMF(*Efield[0][1], geom[0], num_shift, dir, Ey_nodal_flag);
   shiftMF(*Efield[0][2], geom[0], num_shift, dir, Ez_nodal_flag);
 
-  // remove particles that are outside of the box
-  mypc->Redistribute(false);  // Redistribute particles
+
+  // optionally fill the new cells with a uniform plasma 
+  // (also level 0 only right now)
+  ParmParse pp("warpx");
+  pp.get("do_plasma_injection", do_plasma_injection);
+  if(do_plasma_injection) {
+
+    // particleBox encloses the cells where we generate particles
+    Box particleBox = geom[0].Domain();
+    int domainLength = particleBox.length(dir);
+    int sign = (num_shift < 0) ? -1 : 1;
+    particleBox.shift(dir, sign*(domainLength - std::abs(num_shift)));
+    particleBox &= geom[0].Domain();
+
+    // get information about the plasma to inject
+    int n_part_per_cell, ispecies;
+    pp.get("injected_plasma_ppc", n_part_per_cell);
+    pp.get("injected_plasma_species", ispecies);
+    Real weight;
+    pp.get("injected_plasma_density", weight);
+
+    const Real* dx  = geom[0].CellSize();
+#if BL_SPACEDIM==3
+    weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell;
+#elif BL_SPACEDIM==2
+    weight *= dx[0]*dx[1]/n_part_per_cell;
+#endif
+
+    WarpXParticleContainer& myspc = mypc->GetParticleContainer(ispecies);
+    const BoxArray& ba = myspc.ParticleBoxArray(0);
+    const DistributionMapping& dm = myspc.ParticleDistributionMap(0); 
+    MultiFab dummy_mf(ba, 1, 0, dm, Fab_noallocate);
+
+    // For each grid, loop only over the cells in the new region
+    for (MFIter mfi(dummy_mf,false); mfi.isValid(); ++mfi) {
+        int gid = mfi.index();
+	Box grid = ba[gid];
+	Box intersectBox = grid & particleBox;
+	if (intersectBox.numPtsOK() == 0) continue;
+	RealBox intersectRealBox { intersectBox, dx, geom[0].ProbLo() };
+
+#if (BL_SPACEDIM == 3)
+	int nx = intersectBox.length(0);
+	int ny = intersectBox.length(1);
+	int nz = intersectBox.length(2);
+#elif (BL_SPACEDIM == 2)
+	int nx = intersectBox.length(0);
+	int ny = 1;
+	int nz = intersectBox.length(1);
+#endif
+
+	for (int k = 0; k < nz; k++) {
+          for (int j = 0; j < ny; j++) {
+            for (int i = 0; i < nx; i++) {
+       	      for (int i_part=0; i_part<n_part_per_cell;i_part++) {
+	        Real particle_shift = (0.5+i_part)/n_part_per_cell;
+#if (BL_SPACEDIM == 3)
+		Real x = intersectRealBox.lo(0) + (i + particle_shift)*dx[0];
+		Real y = intersectRealBox.lo(1) + (j + particle_shift)*dx[1];
+		Real z = intersectRealBox.lo(2) + (k + particle_shift)*dx[2];
+#elif (BL_SPACEDIM == 2)
+		Real x = intersectRealBox.lo(0) + (i + particle_shift)*dx[0];
+		Real y = 0.0;
+		Real z = intersectRealBox.lo(1) + (k + particle_shift)*dx[1];
+#endif   
+
+		int id  = ParticleBase::NextID();
+		int cpu = ParallelDescriptor::MyProc();
+
+		std::vector<Real> pos(3, 0.0);
+#if (BL_SPACEDIM == 3)
+		pos[0] = x;
+		pos[1] = y;
+		pos[2] = z;
+#elif (BL_SPACEDIM == 2)
+		pos[0] = x;
+		pos[1] = z;
+#endif
+		
+		std::vector<Real> attributes(PIdx::nattribs, 0.0);
+		attributes[PIdx::w] = weight;
+		myspc.addOneParticle(id, cpu, pos, attributes);
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // Redistribure (note - this removes particles that are outside of the box)
+  mypc->Redistribute(false);
 }