From c2351d612918c865a54447f5fe24c69cf09ddccf Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 14 Feb 2017 13:20:09 -0800 Subject: A function for sampling the EM field along a slice. --- Source/WarpX.cpp | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) (limited to 'Source/WarpX.cpp') diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 598ba9ce7..6ed613fa0 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -288,3 +288,100 @@ WarpX::Copy(MultiFab& dstmf, int dcomp, int ncomp, const MultiFab& srcmf, int sc dfab.SetBoxType(dst_typ); } } + +Box +WarpX::getIndexBox(const RealBox& real_box) { + BL_ASSERT(max_level == 0); + + IntVect slice_lo, slice_hi; + + D_TERM(slice_lo[0]=floor((real_box.lo(0) - geom[0].ProbLo(0))/geom[0].CellSize(0));, + slice_lo[1]=floor((real_box.lo(1) - geom[0].ProbLo(1))/geom[0].CellSize(1));, + slice_lo[2]=floor((real_box.lo(2) - geom[0].ProbLo(2))/geom[0].CellSize(2));); + + D_TERM(slice_hi[0]=floor((real_box.hi(0) - geom[0].ProbLo(0))/geom[0].CellSize(0));, + slice_hi[1]=floor((real_box.hi(1) - geom[0].ProbLo(1))/geom[0].CellSize(1));, + slice_hi[2]=floor((real_box.hi(2) - geom[0].ProbLo(2))/geom[0].CellSize(2));); + + return Box(slice_lo, slice_hi) & geom[0].Domain(); +} + +void +WarpX::fillSlice(Real z_coord) { + BL_ASSERT(max_level == 0); + + // Get our slice and convert to index space + RealBox real_slice = geom[0].ProbDomain(); + real_slice.setLo(2, z_coord); + real_slice.setHi(2, z_coord); + Box slice_box = getIndexBox(real_slice); + + // define the multifab that stores slice + BoxArray ba = Efield[0][0]->boxArray(); + const IndexType correct_typ(IntVect::TheZeroVector()); + ba.convert(correct_typ); + const DistributionMapping& dm = dmap[0]; + std::vector< std::pair > isects; + ba.intersections(slice_box, isects, false, 0); + Array boxes; + Array procs; + for (int i = 0; i < isects.size(); ++i) { + procs.push_back(dm[isects[i].first]); + boxes.push_back(isects[i].second); + } + procs.push_back(ParallelDescriptor::MyProc()); + BoxArray slice_ba(&boxes[0], boxes.size()); + DistributionMapping slice_dmap(procs); + MultiFab slice(slice_ba, 6, 0, slice_dmap); + + const MultiFab* mfs[6]; + mfs[0] = Efield[0][0].get(); + mfs[1] = Efield[0][1].get(); + mfs[2] = Efield[0][2].get(); + mfs[3] = Bfield[0][0].get(); + mfs[4] = Bfield[0][1].get(); + mfs[5] = Bfield[0][2].get(); + + IntVect flags[6]; + flags[0] = WarpX::Ex_nodal_flag; + flags[1] = WarpX::Ey_nodal_flag; + flags[2] = WarpX::Ez_nodal_flag; + flags[3] = WarpX::Bx_nodal_flag; + flags[4] = WarpX::By_nodal_flag; + flags[5] = WarpX::Bz_nodal_flag; + + // Fill the slice with sampled data + const Real* dx = geom[0].CellSize(); + const MultiFab& mf = *(Efield[0][0]); + const IntVect& flag = Ex_nodal_flag; + const Geometry& g = geom[0]; + for (MFIter mfi(slice); mfi.isValid(); ++mfi) { + int slice_gid = mfi.index(); + Box grid = slice_ba[slice_gid]; + int xlo = grid.smallEnd(0), ylo = grid.smallEnd(1), zlo = grid.smallEnd(2); + int xhi = grid.bigEnd(0), yhi = grid.bigEnd(1), zhi = grid.bigEnd(2); + + IntVect iv = grid.smallEnd(); + ba.intersections(Box(iv, iv), isects, true, 0); + BL_ASSERT(!isects.empty()); + int full_gid = isects[0].first; + + for (int k = zlo; k <= zhi; k++) { + for (int j = ylo; j <= yhi; j++) { + for (int i = xlo; i <= xhi; i++) { + for (int comp = 0; comp < 6; comp++) { + Real x = g.ProbLo(0) + i*dx[0]; + Real y = g.ProbLo(1) + j*dx[1]; + Real z = z_coord; + + D_TERM(iv[0]=floor((x - g.ProbLo(0) + 0.5*g.CellSize(0)*flags[comp][0])/g.CellSize(0));, + iv[1]=floor((y - g.ProbLo(1) + 0.5*g.CellSize(1)*flags[comp][1])/g.CellSize(1));, + iv[2]=floor((z - g.ProbLo(2) + 0.5*g.CellSize(2)*flags[comp][2])/g.CellSize(2));); + + slice[slice_gid](IntVect(i, j, k), comp) = (*mfs[comp])[full_gid](iv); + } + } + } + } + } +} -- cgit v1.2.3 From 8a08f34d13d22edb431c10d9a7ea3fab17f1c773 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 14 Feb 2017 15:45:23 -0800 Subject: Another function for sampling at an arbitrary set of points. --- Source/WarpX.H | 10 +++++++-- Source/WarpX.cpp | 63 ++++++++++++++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 67 insertions(+), 6 deletions(-) (limited to 'Source/WarpX.cpp') diff --git a/Source/WarpX.H b/Source/WarpX.H index 5da070ae3..06a3be8df 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -32,7 +32,12 @@ public: MultiParticleContainer& GetPartContainer () { return *mypc; } - void fillSlice(Real z_coord); + void fillSlice(Real z_coord) const; + + void sampleAtPoints(const Array& x, + const Array& y, + const Array& z, + Array >& result) const; static void FillBoundary (MultiFab& mf, const Geometry& geom, const IntVect& nodalflag); @@ -57,7 +62,7 @@ public: protected: virtual void ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/) override {} - Box getIndexBox(const RealBox& real_box); + Box getIndexBox(const RealBox& real_box) const; private: @@ -100,6 +105,7 @@ private: // Particle container std::unique_ptr mypc; + // Fields: First array for level, second for direction Array > > current; Array > > Efield; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6ed613fa0..65a52f9bc 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -290,7 +290,7 @@ WarpX::Copy(MultiFab& dstmf, int dcomp, int ncomp, const MultiFab& srcmf, int sc } Box -WarpX::getIndexBox(const RealBox& real_box) { +WarpX::getIndexBox(const RealBox& real_box) const { BL_ASSERT(max_level == 0); IntVect slice_lo, slice_hi; @@ -307,7 +307,7 @@ WarpX::getIndexBox(const RealBox& real_box) { } void -WarpX::fillSlice(Real z_coord) { +WarpX::fillSlice(Real z_coord) const { BL_ASSERT(max_level == 0); // Get our slice and convert to index space @@ -352,8 +352,6 @@ WarpX::fillSlice(Real z_coord) { // Fill the slice with sampled data const Real* dx = geom[0].CellSize(); - const MultiFab& mf = *(Efield[0][0]); - const IntVect& flag = Ex_nodal_flag; const Geometry& g = geom[0]; for (MFIter mfi(slice); mfi.isValid(); ++mfi) { int slice_gid = mfi.index(); @@ -385,3 +383,60 @@ WarpX::fillSlice(Real z_coord) { } } } + +void WarpX::sampleAtPoints(const Array& x, + const Array& y, + const Array& z, + Array >& result) const { + BL_ASSERT((x.size() == y.size()) and (y.size() == z.size())); + BL_ASSERT(max_level == 0); + + const MultiFab* mfs[6]; + mfs[0] = Efield[0][0].get(); + mfs[1] = Efield[0][1].get(); + mfs[2] = Efield[0][2].get(); + mfs[3] = Bfield[0][0].get(); + mfs[4] = Bfield[0][1].get(); + mfs[5] = Bfield[0][2].get(); + + IntVect flags[6]; + flags[0] = WarpX::Ex_nodal_flag; + flags[1] = WarpX::Ey_nodal_flag; + flags[2] = WarpX::Ez_nodal_flag; + flags[3] = WarpX::Bx_nodal_flag; + flags[4] = WarpX::By_nodal_flag; + flags[5] = WarpX::Bz_nodal_flag; + + const unsigned npoints = x.size(); + const int ncomps = 6; + result.resize(ncomps); + for (int i = 0; i < ncomps; i++) { + result[i].resize(npoints, 0); + } + + BoxArray ba = Efield[0][0]->boxArray(); + const IndexType correct_typ(IntVect::TheZeroVector()); + ba.convert(correct_typ); + std::vector< std::pair > isects; + + IntVect iv; + const Geometry& g = geom[0]; + for (int i = 0; i < npoints; ++i) { + for (int comp = 0; comp < ncomps; comp++) { + D_TERM(iv[0]=floor((x[i] - g.ProbLo(0) + 0.5*g.CellSize(0)*flags[comp][0])/g.CellSize(0));, + iv[1]=floor((y[i] - g.ProbLo(1) + 0.5*g.CellSize(1)*flags[comp][1])/g.CellSize(1));, + iv[2]=floor((z[i] - g.ProbLo(2) + 0.5*g.CellSize(2)*flags[comp][2])/g.CellSize(2));); + + ba.intersections(Box(iv, iv), isects, true, 0); + const int grid = isects[0].first; + const int who = dmap[0][grid]; + if (who == ParallelDescriptor::MyProc()) { + result[comp][i] = (*mfs[comp])[grid](iv); + } + } + } + + for (int i = 0; i < ncomps; i++) { + ParallelDescriptor::ReduceRealSum(result[i].dataPtr(), result[i].size()); + } +} -- cgit v1.2.3 From 845796c42a090b97662f5804d3fb13bf4a35a297 Mon Sep 17 00:00:00 2001 From: atmyers Date: Thu, 23 Feb 2017 15:27:52 -0800 Subject: Initial commit of plasma injector class. --- Exec/plasma_acceleration/CustomDensityProb.cpp | 8 +++++ Source/CustomDensityProb.cpp | 9 +++++ Source/Make.package | 3 ++ Source/PlasmaInjector.H | 46 ++++++++++++++++++++++++++ Source/PlasmaInjector.cpp | 40 ++++++++++++++++++++++ Source/WarpX.H | 3 +- Source/WarpX.cpp | 22 ++++++++++++ 7 files changed, 130 insertions(+), 1 deletion(-) create mode 100644 Exec/plasma_acceleration/CustomDensityProb.cpp create mode 100644 Source/CustomDensityProb.cpp create mode 100644 Source/PlasmaInjector.H create mode 100644 Source/PlasmaInjector.cpp (limited to 'Source/WarpX.cpp') diff --git a/Exec/plasma_acceleration/CustomDensityProb.cpp b/Exec/plasma_acceleration/CustomDensityProb.cpp new file mode 100644 index 000000000..62e0ae2aa --- /dev/null +++ b/Exec/plasma_acceleration/CustomDensityProb.cpp @@ -0,0 +1,8 @@ +#include + +#include + +Real CustomPlasmaInjector::getDensity(Real x, Real y, Real z) { + std::cout << "Using custom getDensity" << std::endl; + return 0.0; +} diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp new file mode 100644 index 000000000..54e5e0fad --- /dev/null +++ b/Source/CustomDensityProb.cpp @@ -0,0 +1,9 @@ +#include + +#include + +Real CustomPlasmaInjector::getDensity(Real x, Real y, Real z) { + static_assert(false, + "If running with a custom density profile, you must supply a CustomDensityProb.cpp file"); + return 0.0; +} diff --git a/Source/Make.package b/Source/Make.package index bf0d4d0c3..c846c3fa3 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -15,6 +15,9 @@ CEXE_headers += WarpX.H ParticleIterator.H WarpX_f.H WarpXConst.H CEXE_headers += ParticleContainer.H WarpXParticleContainer.H PhysicalParticleContainer.H LaserParticleContainer.H +CEXE_headers += PlasmaInjector.H +CEXE_sources += PlasmaInjector.cpp CustomDensityProb.cpp + F90EXE_sources += WarpX_f.F90 WarpX_picsar.F90 WarpX_laser.F90 ifeq ($(USE_OPENBC_POISSON),TRUE) diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H new file mode 100644 index 000000000..42b8dd79f --- /dev/null +++ b/Source/PlasmaInjector.H @@ -0,0 +1,46 @@ +#ifndef PLASMA_INJECTOR_H_ +#define PLASMA_INJECTOR_H_ + +#include "Real.H" +#include "ParmParse.H" + +class PlasmaInjector +{ +public: + PlasmaInjector(); + + virtual Real getDensity(Real x, Real y, Real z) = 0; + + bool insideBounds(Real x, Real y, Real z); + +protected: + + Real xmin, xmax; + Real ymin, ymax; + Real zmin, zmax; + + Real density; + +}; + +class ConstantPlasmaInjector : public PlasmaInjector { +public: + Real getDensity( Real x, Real y, Real z); +}; + +class CustomPlasmaInjector : public PlasmaInjector { +public: + Real getDensity( Real x, Real y, Real z); +}; + +class DoubleRampPlasmaInjector : public PlasmaInjector { +public: + DoubleRampPlasmaInjector(); + Real getDensity( Real x, Real y, Real z); + +protected: + Real plateau_length; + Real ramp_length; +}; + +#endif diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp new file mode 100644 index 000000000..3f4d4d336 --- /dev/null +++ b/Source/PlasmaInjector.cpp @@ -0,0 +1,40 @@ +#include "PlasmaInjector.H" + +#include + +PlasmaInjector::PlasmaInjector() { + ParmParse pp("plasma"); + + pp.get("xmin", xmin); + pp.get("ymin", ymin); + pp.get("zmin", zmin); + pp.get("xmax", xmax); + pp.get("ymax", ymax); + pp.get("zmax", zmax); + + pp.get("density", density); +} + +bool PlasmaInjector::insideBounds(Real x, Real y, Real z) { + if (x >= xmax || x < xmin || + y >= ymax || y < ymin || + z >= zmax || z < zmin ) return true; + return false; +} + +Real ConstantPlasmaInjector::getDensity(Real x, Real y, Real z) { + return density; +} + +DoubleRampPlasmaInjector::DoubleRampPlasmaInjector() + : PlasmaInjector() +{ + ParmParse pp("plasma"); + pp.get("ramp_length", ramp_length); + pp.get("plateau_length", plateau_length); +} + +Real DoubleRampPlasmaInjector::getDensity(Real x, Real y, Real z) { + return density; +} + diff --git a/Source/WarpX.H b/Source/WarpX.H index 9f15073c8..f1fd3ee7d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -10,7 +10,7 @@ #include #include - +#include #include class WarpX @@ -133,6 +133,7 @@ private: Array injected_plasma_ppc; Array injected_plasma_species; Array injected_plasma_density; + std::unique_ptr plasma_injector; // Other runtime parameters int verbose = 1; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 8ff55e65b..8a7f4e214 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -176,6 +176,28 @@ WarpX::ReadParameters () pp.query("use_laser", use_laser); } + { + ParmParse pp("plasma"); + std::string plasma_profile_s; + pp.get("profile", plasma_profile_s); + std::transform(plasma_profile_s.begin(), plasma_profile_s.end(), plasma_profile_s.begin(), ::tolower); + if (plasma_profile_s == "constant") { + plasma_injector.reset(new ConstantPlasmaInjector); + } + + else if (plasma_profile_s == "double_ramp") { + plasma_injector.reset(new DoubleRampPlasmaInjector); + } + + else if (plasma_profile_s == "custom") { + plasma_injector.reset(new CustomPlasmaInjector); + } + + else { + BoxLib::Abort("Unknown plasma injector type"); + } + } + { ParmParse pp("interpolation"); pp.query("nox", nox); -- cgit v1.2.3 From a2ae29034bc0626a833cf1ac9fe672f2dd35561b Mon Sep 17 00:00:00 2001 From: atmyers Date: Tue, 14 Mar 2017 12:13:14 -0700 Subject: Update plasma injection stuff to work with amrex. --- Exec/plasma_acceleration/CustomDensityProb.cpp | 2 ++ Source/PlasmaInjector.H | 26 +++++++++++++------------- Source/PlasmaInjector.cpp | 2 ++ Source/WarpX.H | 12 ++++++------ Source/WarpX.cpp | 4 ++-- 5 files changed, 25 insertions(+), 21 deletions(-) (limited to 'Source/WarpX.cpp') diff --git a/Exec/plasma_acceleration/CustomDensityProb.cpp b/Exec/plasma_acceleration/CustomDensityProb.cpp index 62e0ae2aa..247803999 100644 --- a/Exec/plasma_acceleration/CustomDensityProb.cpp +++ b/Exec/plasma_acceleration/CustomDensityProb.cpp @@ -2,6 +2,8 @@ #include +using namespace amrex; + Real CustomPlasmaInjector::getDensity(Real x, Real y, Real z) { std::cout << "Using custom getDensity" << std::endl; return 0.0; diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H index 42b8dd79f..5847fb4a6 100644 --- a/Source/PlasmaInjector.H +++ b/Source/PlasmaInjector.H @@ -1,46 +1,46 @@ #ifndef PLASMA_INJECTOR_H_ #define PLASMA_INJECTOR_H_ -#include "Real.H" -#include "ParmParse.H" +#include "AMReX_Real.H" +#include "AMReX_ParmParse.H" class PlasmaInjector { public: PlasmaInjector(); - virtual Real getDensity(Real x, Real y, Real z) = 0; + virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) = 0; - bool insideBounds(Real x, Real y, Real z); + bool insideBounds(amrex::Real x, amrex::Real y, amrex::Real z); protected: - Real xmin, xmax; - Real ymin, ymax; - Real zmin, zmax; + amrex::Real xmin, xmax; + amrex::Real ymin, ymax; + amrex::Real zmin, zmax; - Real density; + amrex::Real density; }; class ConstantPlasmaInjector : public PlasmaInjector { public: - Real getDensity( Real x, Real y, Real z); + amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); }; class CustomPlasmaInjector : public PlasmaInjector { public: - Real getDensity( Real x, Real y, Real z); + amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); }; class DoubleRampPlasmaInjector : public PlasmaInjector { public: DoubleRampPlasmaInjector(); - Real getDensity( Real x, Real y, Real z); + amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); protected: - Real plateau_length; - Real ramp_length; + amrex::Real plateau_length; + amrex::Real ramp_length; }; #endif diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index 3f4d4d336..5106a9333 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -2,6 +2,8 @@ #include +using namespace amrex; + PlasmaInjector::PlasmaInjector() { ParmParse pp("plasma"); diff --git a/Source/WarpX.H b/Source/WarpX.H index 4c35f0b2b..152ab577c 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -35,12 +35,12 @@ public: MultiParticleContainer& GetPartContainer () { return *mypc; } - void fillSlice(Real z_coord) const; + void fillSlice(amrex::Real z_coord) const; - void sampleAtPoints(const Array& x, - const Array& y, - const Array& z, - Array >& result) const; + void sampleAtPoints(const amrex::Array& x, + const amrex::Array& y, + const amrex::Array& z, + amrex::Array >& result) const; static void FillBoundary (amrex::MultiFab& mf, const amrex::Geometry& geom, const amrex::IntVect& nodalflag); @@ -119,7 +119,7 @@ protected: //! Delete level data. Called by AmrCore::regrid. virtual void ClearLevel (int lev) override; - Box getIndexBox(const RealBox& real_box) const; + amrex::Box getIndexBox(const amrex::RealBox& real_box) const; private: diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 419db0b11..152daa545 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -197,7 +197,7 @@ WarpX::ReadParameters () } else { - BoxLib::Abort("Unknown plasma injector type"); + amrex::Abort("Unknown plasma injector type"); } } @@ -375,7 +375,7 @@ WarpX::fillSlice(Real z_coord) const { procs.push_back(ParallelDescriptor::MyProc()); BoxArray slice_ba(&boxes[0], boxes.size()); DistributionMapping slice_dmap(procs); - MultiFab slice(slice_ba, 6, 0, slice_dmap); + MultiFab slice(slice_ba, slice_dmap, 6, 0); const MultiFab* mfs[6]; mfs[0] = Efield[0][0].get(); -- cgit v1.2.3 From a44b32cad0d8df3decdb4f0ac1dfc66364236136 Mon Sep 17 00:00:00 2001 From: atmyers Date: Tue, 14 Mar 2017 15:28:53 -0700 Subject: Move the plasma injector stuff to particle container. --- Exec/plasma_acceleration/inputs | 37 +++++++++++++++++++----------------- Source/PhysicalParticleContainer.H | 5 +++++ Source/PhysicalParticleContainer.cpp | 24 +++++++++++++++++++++++ Source/WarpX.H | 2 -- Source/WarpX.cpp | 22 --------------------- 5 files changed, 49 insertions(+), 41 deletions(-) (limited to 'Source/WarpX.cpp') diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 26ade8abf..91bcd121d 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -32,25 +32,28 @@ algo.particle_pusher = 0 warpx.cfl = 1.0 particles.nspecies = 2 -pwfa.num_particles_per_cell = 4 - # The beam -pwfa.beam_xmin = -20.e-6 -pwfa.beam_xmax = 20.e-6 -pwfa.beam_ymin = -20.e-6 -pwfa.beam_ymax = 20.e-6 -pwfa.beam_zmin = -150.e-6 -pwfa.beam_zmax = -100.e-6 -pwfa.beam_density = 1e23 # number of electrons per m^3 -pwfa.beam_gamma = 1e9 +species0.profile = "constant" +species0.xmin = -20.e-6 +species0.xmax = 20.e-6 +species0.ymin = -20.e-6 +species0.ymax = 20.e-6 +species0.zmin = -150.e-6 +species0.zmax = -100.e-6 +species0.density = 1e23 # number of electrons per m^3 +species0.gamma = 1e9 +species0.num_particles_per_cell = 4 + # The plasma -pwfa.plasma_xmin = -200.e-6 -pwfa.plasma_xmax = 200.e-6 -pwfa.plasma_ymin = -200.e-6 -pwfa.plasma_ymax = 200.e-6 -pwfa.plasma_zmin = 0.e-6 -pwfa.plasma_zmax = 200.e-6 -pwfa.plasma_density = 1e22 # number of electrons per m^3 +species1.profile = "constant" +species1.xmin = -200.e-6 +species1.xmax = 200.e-6 +species1.ymin = -200.e-6 +species1.ymax = 200.e-6 +species1.zmin = 0.e-6 +species1.zmax = 200.e-6 +species1.density = 1e22 # number of electrons per m^3 +species1.num_particles_per_cell = 4 # Moving window warpx.do_moving_window = 1 diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index c6d13185e..541de7c98 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -3,6 +3,7 @@ #include +#include #include class PhysicalParticleContainer @@ -26,6 +27,10 @@ public: amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::Real t, amrex::Real dt) override; virtual void PostRestart () override {} + +protected: + + std::unique_ptr plasma_injector; }; #endif diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 8c229d41f..767a41919 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -1,6 +1,8 @@ #include +#include #include +#include #include #include @@ -9,6 +11,28 @@ using namespace amrex; PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies) : WarpXParticleContainer(amr_core, ispecies) { + std::stringstream ss; + ss << "species" << ispecies; + ParmParse pp(ss.str()); + + std::string plasma_profile_s; + pp.get("profile", plasma_profile_s); + std::transform(plasma_profile_s.begin(), plasma_profile_s.end(), plasma_profile_s.begin(), ::tolower); + if (plasma_profile_s == "constant") { + plasma_injector.reset(new ConstantPlasmaInjector); + } + + else if (plasma_profile_s == "double_ramp") { + plasma_injector.reset(new DoubleRampPlasmaInjector); + } + + else if (plasma_profile_s == "custom") { + plasma_injector.reset(new CustomPlasmaInjector); + } + + else { + amrex::Abort("Unknown plasma injector type"); + } } void diff --git a/Source/WarpX.H b/Source/WarpX.H index 152ab577c..c198990dc 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -13,7 +13,6 @@ #include #include -#include class WarpX : public amrex::AmrCore @@ -183,7 +182,6 @@ private: amrex::Array injected_plasma_ppc; amrex::Array injected_plasma_species; amrex::Array injected_plasma_density; - std::unique_ptr plasma_injector; // Other runtime parameters int verbose = 1; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 152daa545..7ac23307b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -179,28 +179,6 @@ WarpX::ReadParameters () pp.query("use_laser", use_laser); } - { - ParmParse pp("plasma"); - std::string plasma_profile_s; - pp.get("profile", plasma_profile_s); - std::transform(plasma_profile_s.begin(), plasma_profile_s.end(), plasma_profile_s.begin(), ::tolower); - if (plasma_profile_s == "constant") { - plasma_injector.reset(new ConstantPlasmaInjector); - } - - else if (plasma_profile_s == "double_ramp") { - plasma_injector.reset(new DoubleRampPlasmaInjector); - } - - else if (plasma_profile_s == "custom") { - plasma_injector.reset(new CustomPlasmaInjector); - } - - else { - amrex::Abort("Unknown plasma injector type"); - } - } - { ParmParse pp("interpolation"); pp.query("nox", nox); -- cgit v1.2.3 From 8b7bbe4b1bff5d6176970cbb0405640f702706b2 Mon Sep 17 00:00:00 2001 From: atmyers Date: Wed, 15 Mar 2017 15:35:14 -0700 Subject: refactor particle initialization code, both for initial conditions and when the window moves. --- Exec/plasma_acceleration/ParticleProb.cpp | 58 +-------------- Exec/plasma_acceleration/inputs | 8 -- Source/PhysicalParticleContainer.H | 32 ++++++-- Source/PhysicalParticleContainer.cpp | 120 ++++++++++++++++++++++++++++++ Source/WarpX.H | 2 - Source/WarpX.cpp | 6 -- Source/WarpXEvolve.cpp | 15 +--- Source/WarpXParticleContainer.H | 3 - Source/WarpXParticleContainer.cpp | 43 ----------- 9 files changed, 149 insertions(+), 138 deletions(-) (limited to 'Source/WarpX.cpp') diff --git a/Exec/plasma_acceleration/ParticleProb.cpp b/Exec/plasma_acceleration/ParticleProb.cpp index 010484315..806861a19 100644 --- a/Exec/plasma_acceleration/ParticleProb.cpp +++ b/Exec/plasma_acceleration/ParticleProb.cpp @@ -16,61 +16,5 @@ using namespace amrex; void PhysicalParticleContainer::InitData() { - BL_PROFILE("PhysicalParticleContainer::InitData()"); - - charge = plasma_injector->getCharge(); - mass = plasma_injector->getMass(); - - const int lev = 0; - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real scale_fac; - int n_part_per_cell = plasma_injector->numParticlesPerCell(); - -#if BL_SPACEDIM==3 - scale_fac = dx[0]*dx[1]*dx[2]/n_part_per_cell; -#elif BL_SPACEDIM==2 - scale_face = dx[0]*dx[1]/n_part_per_cell; -#endif - - std::array attribs; - attribs.fill(0.0); - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_partinsideBounds(x, y, z)) { - Real weight; - Real u[3]; - weight = plasma_injector->getDensity(x, y, z) * scale_fac; - plasma_injector->getMomentum(u); - attribs[PIdx::w ] = weight; - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } + InitNPerCell(); } diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 45832d536..241fc705a 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -38,27 +38,22 @@ particles.species_names = beam plasma # The beam species information beam.charge = -q_e beam.mass = m_e - beam.xmin = -20.e-6 beam.xmax = 20.e-6 beam.ymin = -20.e-6 beam.ymax = 20.e-6 beam.zmin = -150.e-6 beam.zmax = -100.e-6 - beam.profile = "constant" beam.density = 1e23 # number of particles per m^3 beam.num_particles_per_cell = 4 - beam.ux = 0.0 beam.uy = 0.0 beam.uz = 1e9 - # The plasma species information plasma.charge = -q_e plasma.mass = m_e - plasma.profile = "constant" plasma.xmin = -200.e-6 plasma.xmax = 200.e-6 @@ -66,7 +61,6 @@ plasma.ymin = -200.e-6 plasma.ymax = 200.e-6 plasma.zmin = 0.e-6 plasma.zmax = 200.e-6 - plasma.profile = "constant" plasma.density = 1e22 # number of particles per m^3 plasma.num_particles_per_cell = 4 @@ -84,5 +78,3 @@ warpx.moving_window_v = 1.0 # in units of the speed of light 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/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 923d40429..a399df3fe 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -10,7 +10,9 @@ class PhysicalParticleContainer : public WarpXParticleContainer { public: - PhysicalParticleContainer (amrex::AmrCore* amr_core, int ispecies, const std::string& name); + PhysicalParticleContainer (amrex::AmrCore* amr_core, + int ispecies, + const std::string& name); virtual ~PhysicalParticleContainer () {} virtual void AllocData () override; @@ -18,20 +20,36 @@ public: virtual void InitData () override; virtual void FieldGather(int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) override; virtual void Evolve (int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, - amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::Real t, amrex::Real dt) override; + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz, + amrex::MultiFab& jx, + amrex::MultiFab& jy, + amrex::MultiFab& jz, + amrex::Real t, + amrex::Real dt) override; virtual void PostRestart () override {} + // Inject 'ppc' particles per cell in Box 'part_box' + void AddParticles (int lev, const amrex::Box& part_box); + protected: - std::string species_name; + void InitNPerCell (); + std::string species_name; std::unique_ptr plasma_injector; }; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index c7e98def9..0e3875ca1 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -45,6 +45,68 @@ PhysicalParticleContainer::AllocData () resizeData(); } +void +PhysicalParticleContainer::InitNPerCell () { + BL_PROFILE("PhysicalParticleContainer::InitNPerCell()"); + + charge = plasma_injector->getCharge(); + mass = plasma_injector->getMass(); + + const int lev = 0; + const Geometry& geom = Geom(lev); + const Real* dx = geom.CellSize(); + + Real scale_fac; + int n_part_per_cell = plasma_injector->numParticlesPerCell(); + +#if BL_SPACEDIM==3 + scale_fac = dx[0]*dx[1]*dx[2]/n_part_per_cell; +#elif BL_SPACEDIM==2 + scale_fac = dx[0]*dx[1]/n_part_per_cell; +#endif + + std::array attribs; + attribs.fill(0.0); + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { + const Box& tile_box = mfi.tilebox(); + RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + const auto& boxlo = tile_box.smallEnd(); + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) + { + for (int i_part=0; i_partinsideBounds(x, y, z)) { + Real weight; + Real u[3]; + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + plasma_injector->getMomentum(u); + attribs[PIdx::w ] = weight; + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + } + } + } + } +} + + void PhysicalParticleContainer::FieldGather (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, @@ -152,6 +214,64 @@ PhysicalParticleContainer::FieldGather (int lev, } } +void +PhysicalParticleContainer::AddParticles (int lev, const Box& part_box) +{ + const Geometry& geom = Geom(lev); + const Real* dx = geom.CellSize(); + + Real scale_fac; + int ppc = plasma_injector->numParticlesPerCell(); + +#if BL_SPACEDIM==3 + scale_fac = dx[0]*dx[1]*dx[2]/ppc; +#elif BL_SPACEDIM==2 + scale_fac = dx[0]*dx[1]/ppc; +#endif + + std::array attribs; + attribs.fill(0.0); + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + const Box& intersectBox = tile_box & part_box; + if (intersectBox.ok()) + { + RealBox real_box { intersectBox, dx, geom.ProbLo() }; + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + const IntVect& boxlo = tile_box.smallEnd(); + for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) + { + for (int i_part=0; i_part < ppc; i_part++) + { + Real particle_shift = (0.5+i_part)/ppc; +#if (BL_SPACEDIM == 3) + Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; + Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; + Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; +#elif (BL_SPACEDIM == 2) + Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; + Real y = 0.0; + Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; +#endif + Real weight; + Real u[3]; + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + plasma_injector->getMomentum(u); + attribs[PIdx::w ] = weight; + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + } + } + } + } +} + void PhysicalParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, diff --git a/Source/WarpX.H b/Source/WarpX.H index c198990dc..135a0f690 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -179,9 +179,7 @@ private: // Plasma injection parameters int do_plasma_injection = 0; int num_injected_species = -1; - amrex::Array injected_plasma_ppc; amrex::Array injected_plasma_species; - amrex::Array injected_plasma_density; // Other runtime parameters int verbose = 1; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 7ac23307b..bcdbc709b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -165,15 +165,9 @@ WarpX::ReadParameters () 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); } pp.query("use_laser", use_laser); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 5e5154554..875a0392d 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -262,20 +262,11 @@ WarpX::InjectPlasma (int num_shift, int dir) const Real* dx = geom[lev].CellSize(); - for (int i = 0; i < num_injected_species; ++i) - { - int ppc = injected_plasma_ppc[i]; - Real density = injected_plasma_density[i]; -#if BL_SPACEDIM==3 - Real weight = density * dx[0]*dx[1]*dx[2]/ppc; -#elif BL_SPACEDIM==2 - Real weight = density * dx[0]*dx[1]/ppc; -#endif - + for (int i = 0; i < num_injected_species; ++i) { int ispecies = injected_plasma_species[i]; WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - - pc.AddParticles(lev, particleBox, weight, ppc); + PhysicalParticleContainer* ppc = (PhysicalParticleContainer*) &pc; + ppc->AddParticles(lev, particleBox); } } } diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index b2f308ac5..6fc8b1415 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -67,9 +67,6 @@ public: amrex::Real x, amrex::Real y, amrex::Real z, const std::array& attribs); - // Inject 'ppc' particles per cell in Box 'part_box' - void AddParticles (int lev, const amrex::Box& part_box, amrex::Real weight, int ppc); - void ReadHeader (std::istream& is); void WriteHeader (std::ostream& os) const; diff --git a/Source/WarpXParticleContainer.cpp b/Source/WarpXParticleContainer.cpp index 76860ee1e..25b92a712 100644 --- a/Source/WarpXParticleContainer.cpp +++ b/Source/WarpXParticleContainer.cpp @@ -65,49 +65,6 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, particle_tile.push_back(attribs); } -void -WarpXParticleContainer::AddParticles (int lev, const Box& part_box, Real weight, int ppc) -{ - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - std::array attribs; - attribs.fill(0.0); - attribs[PIdx::w] = weight; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - const Box& intersectBox = tile_box & part_box; - if (intersectBox.ok()) - { - RealBox real_box { intersectBox, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part < ppc; i_part++) - { - Real particle_shift = (0.5+i_part)/ppc; -#if (BL_SPACEDIM == 3) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} - void WarpXParticleContainer::AddNParticles (int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, -- cgit v1.2.3 From 94793352a37e0130d8f73b1773ecc773bc678286 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 21 Mar 2017 13:32:03 -0700 Subject: Fixing some tests. --- Exec/Langmuir/inputs.rt | 2 +- Exec/uniform_plasma/inputs.rt | 1 + Regression/WarpX-tests.ini | 8 ++++---- Source/WarpX.H | 2 ++ Source/WarpX.cpp | 2 ++ 5 files changed, 10 insertions(+), 5 deletions(-) (limited to 'Source/WarpX.cpp') diff --git a/Exec/Langmuir/inputs.rt b/Exec/Langmuir/inputs.rt index feeb48c93..62e62928b 100644 --- a/Exec/Langmuir/inputs.rt +++ b/Exec/Langmuir/inputs.rt @@ -55,4 +55,4 @@ electrons.profile = constant electrons.density = 1.e25 # number of electrons per m^3 electrons.momentum_distribution_type = "constant" -electrons.ux = 0.01 # ux = gamma*beta_x + diff --git a/Exec/uniform_plasma/inputs.rt b/Exec/uniform_plasma/inputs.rt index c3ae94312..bf6bb127a 100644 --- a/Exec/uniform_plasma/inputs.rt +++ b/Exec/uniform_plasma/inputs.rt @@ -36,6 +36,7 @@ particles.species_names = electrons electrons.charge = -q_e electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" electrons.num_particles_per_cell = 10 electrons.profile = constant electrons.density = 1.e25 # number of electrons per m^3 diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 02de34fda..ed7786184 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -65,7 +65,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.ux=0.01 langmuirwave.particle_xmax=0.e-6 +runtime_params = electrons.ux=0.01 electrons.particle_xmax=0.e-6 [Langmuir_x] buildDir = Exec/Langmuir @@ -78,7 +78,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.ux=0.01 langmuirwave.particle_xmax=0.e-6 +runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 [Langmuir_y] buildDir = Exec/Langmuir @@ -91,7 +91,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.uy=0.01 langmuirwave.particle_ymax=0.e-6 +runtime_params = electrons.uy=0.01 electrons.ymax=0.e-6 [Langmuir_z] buildDir = Exec/Langmuir @@ -104,7 +104,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.uz=0.01 langmuirwave.particle_zmax=0.e-6 +runtime_params = electrons.uz=0.01 electrons.zmax=0.e-6 [Langmuir_multi] buildDir = Exec/Langmuir diff --git a/Source/WarpX.H b/Source/WarpX.H index 135a0f690..382a6f80c 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -34,12 +34,14 @@ public: MultiParticleContainer& GetPartContainer () { return *mypc; } +#if (BL_SPACEDIM == 3) void fillSlice(amrex::Real z_coord) const; void sampleAtPoints(const amrex::Array& x, const amrex::Array& y, const amrex::Array& z, amrex::Array >& result) const; +#endif static void FillBoundary (amrex::MultiFab& mf, const amrex::Geometry& geom, const amrex::IntVect& nodalflag); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index bcdbc709b..226777868 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -321,6 +321,7 @@ WarpX::getIndexBox(const RealBox& real_box) const { return Box(slice_lo, slice_hi) & geom[0].Domain(); } +#if (BL_SPACEDIM == 3) void WarpX::fillSlice(Real z_coord) const { BL_ASSERT(max_level == 0); @@ -455,3 +456,4 @@ void WarpX::sampleAtPoints(const Array& x, ParallelDescriptor::ReduceRealSum(result[i].dataPtr(), result[i].size()); } } +#endif -- cgit v1.2.3