Add ability to inject several species at once.

1 files changed, 78 insertions, 72 deletions

diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index 3f1001a42..975d54d6a 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -272,45 +272,10 @@ WarpX::ComputeDt ()
 }
 
 void
-WarpX::MoveWindow ()
-{
-  
-  if (do_moving_window == 0) return;
-
-  // compute the number of cells to shift
-  int dir = moving_window_dir;
-  Real new_lo[BL_SPACEDIM];
-  Real new_hi[BL_SPACEDIM];    
-  const Real* current_lo = geom[0].ProbLo();
-  const Real* current_hi = geom[0].ProbHi();
-  const Real* dx = geom[0].CellSize();
-  moving_window_x += moving_window_v * dt[0];
-  int num_shift = (moving_window_x - current_lo[dir]) / dx[dir];
-  
-  if (num_shift == 0) return;
-  
-  // update the problem domain
-  for (int i=0; i<BL_SPACEDIM; i++) {
-    new_lo[i] = current_lo[i];
-    new_hi[i] = current_hi[i];
-  }
-  new_lo[dir] = current_lo[dir] + num_shift * dx[dir];
-  new_hi[dir] = current_hi[dir] + num_shift * dx[dir];
-  RealBox new_box(new_lo, new_hi);
-  geom[0].ProbDomain(new_box);
-  
-  // shift the mesh fields (Note - only on level 0 for now)
-  shiftMF(*Bfield[0][0], geom[0], num_shift, dir, Bx_nodal_flag);
-  shiftMF(*Bfield[0][1], geom[0], num_shift, dir, By_nodal_flag);
-  shiftMF(*Bfield[0][2], geom[0], num_shift, dir, Bz_nodal_flag);
-  shiftMF(*Efield[0][0], geom[0], num_shift, dir, Ex_nodal_flag);
-  shiftMF(*Efield[0][1], geom[0], num_shift, dir, Ey_nodal_flag);
-  shiftMF(*Efield[0][2], geom[0], num_shift, dir, Ez_nodal_flag);
+WarpX::InjectPlasma(int num_shift, int dir) {
 
-  // optionally fill the new cells with a uniform plasma 
-  // (also level 0 only right now)
   if(do_plasma_injection) {
-
+    
     // particleBox encloses the cells where we generate particles
     Box particleBox = geom[0].Domain();
     int domainLength = particleBox.length(dir);
@@ -318,20 +283,11 @@ WarpX::MoveWindow ()
     particleBox.shift(dir, sign*(domainLength - std::abs(num_shift)));
     particleBox &= geom[0].Domain();
 
-    int n_part_per_cell = injected_plasma_ppc;
-    int ispecies = injected_plasma_species;
-    Real weight = injected_plasma_density;
-
+    // get a dummy mf to loop over
     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); 
+    WarpXParticleContainer* myspc = &(mypc->GetParticleContainer(0));
+    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
@@ -355,40 +311,90 @@ WarpX::MoveWindow ()
 	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;
+	      for (int ispec=0; ispec < num_injected_species; ispec++) {
+		int ispecies = injected_plasma_species[ispec];
+		myspc = &(mypc->GetParticleContainer(ispecies));
+		for (int i_part=0; i_part < injected_plasma_ppc[ispec]; i_part++) {
+		  Real particle_shift = (0.5+i_part)/injected_plasma_ppc[ispec];
 #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];
+		  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];
+		  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();
+		  int id  = ParticleBase::NextID();
+		  int cpu = ParallelDescriptor::MyProc();
 
-		std::vector<Real> pos(3, 0.0);
+		  std::vector<Real> pos(3, 0.0);
 #if (BL_SPACEDIM == 3)
-		pos[0] = x;
-		pos[1] = y;
-		pos[2] = z;
+		  pos[0] = x;
+		  pos[1] = y;
+		  pos[2] = z;
 #elif (BL_SPACEDIM == 2)
-		pos[0] = x;
-		pos[1] = z;
+		  pos[0] = x;
+		  pos[1] = z;
+#endif
+
+		  std::vector<Real> attributes(PIdx::nattribs, 0.0);
+
+		  Real weight = injected_plasma_density[ispec];
+#if BL_SPACEDIM==3
+		  weight *= dx[0]*dx[1]*dx[2]/injected_plasma_ppc[ispec];
+#elif BL_SPACEDIM==2
+		  weight *= dx[0]*dx[1]/injected_plasma_ppc[ispec];
 #endif
-		
-		std::vector<Real> attributes(PIdx::nattribs, 0.0);
-		attributes[PIdx::w] = weight;
-		myspc.addOneParticle(id, cpu, pos, attributes);
-            }
-          }
+		  attributes[PIdx::w] = weight;
+		  myspc->addOneParticle(id, cpu, pos, attributes);
+	      }
+	    }
+	  }
         }
       }
     }
   }
+}
+
+void
+WarpX::MoveWindow ()
+{
+  
+  if (do_moving_window == 0) return;
+
+  // compute the number of cells to shift
+  int dir = moving_window_dir;
+  Real new_lo[BL_SPACEDIM];
+  Real new_hi[BL_SPACEDIM];    
+  const Real* current_lo = geom[0].ProbLo();
+  const Real* current_hi = geom[0].ProbHi();
+  const Real* dx = geom[0].CellSize();
+  moving_window_x += moving_window_v * dt[0];
+  int num_shift = (moving_window_x - current_lo[dir]) / dx[dir];
+  
+  if (num_shift == 0) return;
+  
+  // update the problem domain
+  for (int i=0; i<BL_SPACEDIM; i++) {
+    new_lo[i] = current_lo[i];
+    new_hi[i] = current_hi[i];
+  }
+  new_lo[dir] = current_lo[dir] + num_shift * dx[dir];
+  new_hi[dir] = current_hi[dir] + num_shift * dx[dir];
+  RealBox new_box(new_lo, new_hi);
+  geom[0].ProbDomain(new_box);
+  
+  // shift the mesh fields (Note - only on level 0 for now)
+  shiftMF(*Bfield[0][0], geom[0], num_shift, dir, Bx_nodal_flag);
+  shiftMF(*Bfield[0][1], geom[0], num_shift, dir, By_nodal_flag);
+  shiftMF(*Bfield[0][2], geom[0], num_shift, dir, Bz_nodal_flag);
+  shiftMF(*Efield[0][0], geom[0], num_shift, dir, Ex_nodal_flag);
+  shiftMF(*Efield[0][1], geom[0], num_shift, dir, Ey_nodal_flag);
+  shiftMF(*Efield[0][2], geom[0], num_shift, dir, Ez_nodal_flag);
+
+  InjectPlasma(num_shift, dir);
 
   // Redistribute (note - this removes particles that are outside of the box)
   mypc->Redistribute(false);