diff options
| -rw-r--r-- | Exec/plasma_acceleration/inputs | 7 | ||||
| -rw-r--r-- | Source/WarpX.H | 9 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 16 | ||||
| -rw-r--r-- | Source/WarpXEvolve.cpp | 93 |
4 files changed, 120 insertions, 5 deletions
diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 951ea4e85..ef8729aaa 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -56,3 +56,10 @@ pwfa.plasma_density = 1e22 # number of electrons per m^3 warpx.do_moving_window = 1 warpx.moving_window_dir = 2 # 0, 1, 2, for x, y, z warpx.moving_window_v = 1.0 # in units of the speed of light + +# Particle Injection +warpx.do_plasma_injection = 1 +warpx.num_injected_species = 1 +warpx.injected_plasma_species = 1 0 +warpx.injected_plasma_density = 1e22 1e22 +warpx.injected_plasma_ppc = 4 4 diff --git a/Source/WarpX.H b/Source/WarpX.H index b03b4e9fc..be290c6f2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -77,6 +77,7 @@ private: void ComputeDt (); void MoveWindow (); + void InjectPlasma(int num_shift, int dir); void WriteWarpXHeader(const std::string& name) const; void WriteCheckPointFile() const; @@ -106,8 +107,14 @@ private: Real moving_window_x; Real moving_window_v; - // runtime parameters + // Plasma injection parameters + int do_plasma_injection = 0; + int num_injected_species; + Array<int> injected_plasma_ppc; + Array<int> injected_plasma_species; + Array<Real> injected_plasma_density; + // Other runtime parameters int verbose = 1; int max_step = std::numeric_limits<int>::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index bb6509786..00453fa5c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -133,6 +133,20 @@ WarpX::ReadParameters () pp.get("moving_window_v", moving_window_v); } + pp.query("do_plasma_injection", do_plasma_injection); + if (do_plasma_injection) { + pp.get("num_injected_species", num_injected_species); + injected_plasma_ppc.resize(num_injected_species); + pp.getarr("injected_plasma_ppc", injected_plasma_ppc, + 0, num_injected_species); + injected_plasma_species.resize(num_injected_species); + pp.getarr("injected_plasma_species", injected_plasma_species, + 0, num_injected_species); + injected_plasma_density.resize(num_injected_species); + pp.getarr("injected_plasma_density", injected_plasma_density, + 0, num_injected_species); + } + moving_window_x = geom[0].ProbLo(moving_window_dir); moving_window_v = 0.0; pp.query("moving_window_v", moving_window_v); @@ -238,8 +252,6 @@ WarpX::shiftMF(MultiFab& mf, const Geometry& geom, int num_shift, mf[mfi].copy(tmpmf[mfi], srcBox, 0, dstBox, 0, mf.nComp()); mf[mfi].SetBoxType(dst_typ); } - - mf.FillBoundary(geom.periodicity()); } void diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 76c928621..975d54d6a 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -272,6 +272,93 @@ WarpX::ComputeDt () } void +WarpX::InjectPlasma(int num_shift, int dir) { + + 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 a dummy mf to loop over + const Real* dx = geom[0].CellSize(); + 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 + for (MFIter mfi(dummy_mf,false); mfi.isValid(); ++mfi) { + int gid = mfi.index(); + Box grid = ba[gid]; + Box intersectBox = grid & particleBox; + if (intersectBox.isEmpty()) 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 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]; +#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); + + 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 + attributes[PIdx::w] = weight; + myspc->addOneParticle(id, cpu, pos, attributes); + } + } + } + } + } + } + } +} + +void WarpX::MoveWindow () { @@ -307,6 +394,8 @@ 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 + InjectPlasma(num_shift, dir); + + // Redistribute (note - this removes particles that are outside of the box) + mypc->Redistribute(false); } |