From f717fec86e8dae31660890499ff4b1a518a5024b Mon Sep 17 00:00:00 2001 From: atmyers Date: Sun, 22 Jan 2017 11:19:44 -0800 Subject: Implement ability to inject fresh plasma from C++ --- Source/WarpXEvolve.cpp | 93 ++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 91 insertions(+), 2 deletions(-) (limited to 'Source/WarpXEvolve.cpp') 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 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 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); } -- cgit v1.2.3 From 75691dcfa0eec0488567cdeb14d85113c2401586 Mon Sep 17 00:00:00 2001 From: atmyers Date: Sun, 22 Jan 2017 11:23:29 -0800 Subject: Fix typo in comment. --- Source/WarpXEvolve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 76e5f52e6..98ff18983 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -396,6 +396,6 @@ WarpX::MoveWindow () } } - // Redistribure (note - this removes particles that are outside of the box) + // Redistribute (note - this removes particles that are outside of the box) mypc->Redistribute(false); } -- cgit v1.2.3 From cddcef197c594362bb3ddf3b1abaf750c64dbc71 Mon Sep 17 00:00:00 2001 From: atmyers Date: Mon, 23 Jan 2017 10:56:18 -0800 Subject: Fixing this check to be assertion safe in respose to Weiqun's comment. --- Source/WarpXEvolve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 98ff18983..a3905b2ce 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -345,7 +345,7 @@ WarpX::MoveWindow () int gid = mfi.index(); Box grid = ba[gid]; Box intersectBox = grid & particleBox; - if (intersectBox.numPtsOK() == 0) continue; + if (intersectBox.isEmpty()) continue; RealBox intersectRealBox { intersectBox, dx, geom[0].ProbLo() }; #if (BL_SPACEDIM == 3) -- cgit v1.2.3 From 05e1438b61322cfd831e9dea8040852fd6759aaf Mon Sep 17 00:00:00 2001 From: atmyers Date: Mon, 23 Jan 2017 11:25:04 -0800 Subject: Only ParmParse the plasma injection parameters at initialization. --- Source/WarpX.H | 8 ++++++-- Source/WarpX.cpp | 7 +++++++ Source/WarpXEvolve.cpp | 12 +++--------- 3 files changed, 16 insertions(+), 11 deletions(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/WarpX.H b/Source/WarpX.H index 0113413d1..691930d1e 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -105,10 +105,14 @@ private: int moving_window_dir = 0; Real moving_window_x; Real moving_window_v; - int do_plasma_injection = 0; - // runtime parameters + // Plasma injection parameters + int do_plasma_injection = 0; + int injected_plasma_ppc; + int injected_plasma_species; + Real injected_plasma_density; + // Other runtime parameters int verbose = 1; int max_step = std::numeric_limits::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index bb6509786..fbfc063ac 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -133,6 +133,13 @@ WarpX::ReadParameters () pp.get("moving_window_v", moving_window_v); } + pp.query("do_plasma_injection", do_plasma_injection); + if (do_plasma_injection) { + pp.get("injected_plasma_ppc", injected_plasma_ppc); + pp.get("injected_plasma_species", injected_plasma_species); + pp.get("injected_plasma_density", injected_plasma_density); + } + moving_window_x = geom[0].ProbLo(moving_window_dir); moving_window_v = 0.0; pp.query("moving_window_v", moving_window_v); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index a3905b2ce..3f1001a42 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -307,11 +307,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); - // 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 @@ -321,12 +318,9 @@ WarpX::MoveWindow () 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); + int n_part_per_cell = injected_plasma_ppc; + int ispecies = injected_plasma_species; + Real weight = injected_plasma_density; const Real* dx = geom[0].CellSize(); #if BL_SPACEDIM==3 -- cgit v1.2.3 From ee98feb72d6df609facf18bf448a92ee5e7023a2 Mon Sep 17 00:00:00 2001 From: atmyers Date: Mon, 23 Jan 2017 17:09:03 -0800 Subject: Add ability to inject several species at once. --- Exec/plasma_acceleration/inputs | 7 +- Source/WarpX.H | 7 +- Source/WarpX.cpp | 13 +++- Source/WarpXEvolve.cpp | 150 +++++++++++++++++++++------------------- 4 files changed, 96 insertions(+), 81 deletions(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 28d169e92..ef8729aaa 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -59,6 +59,7 @@ warpx.moving_window_v = 1.0 # in units of the speed of light # Particle Injection warpx.do_plasma_injection = 1 -warpx.injected_plasma_species = 1 -warpx.injected_plasma_density = 1e22 -warpx.injected_plasma_ppc = 4 +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 7e21285bf..be290c6f2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -109,9 +109,10 @@ private: // Plasma injection parameters int do_plasma_injection = 0; - int injected_plasma_ppc; - int injected_plasma_species; - Real injected_plasma_density; + int num_injected_species; + Array injected_plasma_ppc; + Array injected_plasma_species; + Array injected_plasma_density; // Other runtime parameters int verbose = 1; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index fbfc063ac..08d432866 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -135,9 +135,16 @@ WarpX::ReadParameters () pp.query("do_plasma_injection", do_plasma_injection); if (do_plasma_injection) { - pp.get("injected_plasma_ppc", injected_plasma_ppc); - pp.get("injected_plasma_species", injected_plasma_species); - pp.get("injected_plasma_density", injected_plasma_density); + 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); 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; iGetParticleContainer(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_partGetParticleContainer(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 pos(3, 0.0); + std::vector 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 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 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; iRedistribute(false); -- cgit v1.2.3