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. --- Source/PlasmaInjector.cpp | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 Source/PlasmaInjector.cpp (limited to 'Source/PlasmaInjector.cpp') 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; +} + -- 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/PlasmaInjector.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 ef8962d7d7352dad6fe4bed7885c88edc129349d Mon Sep 17 00:00:00 2001 From: atmyers Date: Tue, 14 Mar 2017 17:17:16 -0700 Subject: some work towards generalizing the initial conditions / plasma injection --- Exec/plasma_acceleration/ParticleProb.cpp | 87 +++++++++---------------------- Exec/plasma_acceleration/inputs | 43 ++++++++------- Source/ParticleContainer.H | 5 ++ Source/ParticleContainer.cpp | 5 +- Source/PhysicalParticleContainer.H | 4 +- Source/PhysicalParticleContainer.cpp | 17 +++--- Source/PlasmaInjector.H | 39 +++++++++----- Source/PlasmaInjector.cpp | 49 +++++++++++------ 8 files changed, 128 insertions(+), 121 deletions(-) (limited to 'Source/PlasmaInjector.cpp') diff --git a/Exec/plasma_acceleration/ParticleProb.cpp b/Exec/plasma_acceleration/ParticleProb.cpp index 7a538abdb..4e0986c5d 100644 --- a/Exec/plasma_acceleration/ParticleProb.cpp +++ b/Exec/plasma_acceleration/ParticleProb.cpp @@ -16,7 +16,7 @@ using namespace amrex; void PhysicalParticleContainer::InitData() { - BL_PROFILE("SPC::InitData()"); + BL_PROFILE("PhysicalParticleContainer::InitData()"); // species_id 0 : the beam // species_id 1 : the electrons of the plasma @@ -33,64 +33,32 @@ PhysicalParticleContainer::InitData() const Geometry& geom = Geom(lev); const Real* dx = geom.CellSize(); - Real weight, gamma, uz; - Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax; - int n_part_per_cell; - { - ParmParse pp("pwfa"); - - // Calculate the particle weight - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e22; - if (species_id == 0) { - pp.query("beam_density", weight); - } else if (species_id == 1) { - pp.query("plasma_density", weight); - } - #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 - - // Calculate the limits between which the particles will be initialized - particle_xmin = particle_ymin = particle_zmin = -2.e-5; - particle_xmax = particle_ymax = particle_zmax = 2.e-5; - if (species_id == 0) { - pp.query("beam_xmin", particle_xmin); - pp.query("beam_xmax", particle_xmax); - pp.query("beam_ymin", particle_ymin); - pp.query("beam_ymax", particle_ymax); - pp.query("beam_zmin", particle_zmin); - pp.query("beam_zmax", particle_zmax); - } else if (species_id == 1) { - pp.query("plasma_xmin", particle_xmin); - pp.query("plasma_xmax", particle_xmax); - pp.query("plasma_ymin", particle_ymin); - pp.query("plasma_ymax", particle_ymax); - pp.query("plasma_zmin", particle_zmin); - pp.query("plasma_zmax", particle_zmax); - } - uz = 0.; - if (species_id == 0) { // relativistic beam - gamma = 1.; - pp.query("beam_gamma", gamma); - uz = std::sqrt( gamma*gamma - 1 ); - } - uz *= PhysConst::c; + Real weight, gamma, uz, scale_fac; + Real particle_xmin, particle_xmax, + particle_ymin, particle_ymax, + particle_zmin, particle_zmax; + 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 + + // Calculate the limits between which the particles will be initialized + uz = 0.; + if (species_id == 0) { // relativistic beam + gamma = plasma_injector->getGamma(); + uz = std::sqrt( gamma*gamma - 1 ); } - + uz *= PhysConst::c; + std::array attribs; attribs.fill(0.0); - attribs[PIdx::w ] = weight; - attribs[PIdx::uz] = uz; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { + 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(); @@ -110,14 +78,11 @@ PhysicalParticleContainer::InitData() Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; #endif - if (x >= particle_xmax || x < particle_xmin || - y >= particle_ymax || y < particle_ymin || - z >= particle_zmax || z < particle_zmin ) - { - continue; - } - else + if (plasma_injector->insideBounds(x, y, z)) { + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + attribs[PIdx::w ] = weight; + attribs[PIdx::uz] = uz; AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); } } diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 91bcd121d..0fad702b9 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -30,30 +30,35 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 + + +# Information about the particle species particles.nspecies = 2 +particles.species_names = beam plasma # The beam -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 +beam.profile = "constant" +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.density = 1e23 # number of particles per m^3 +beam.gamma = 1e9 +beam.num_particles_per_cell = 4 # The plasma -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 +plasma.profile = "constant" +plasma.xmin = -200.e-6 +plasma.xmax = 200.e-6 +plasma.ymin = -200.e-6 +plasma.ymax = 200.e-6 +plasma.zmin = 0.e-6 +plasma.zmax = 200.e-6 +plasma.density = 1e22 # number of particles per m^3 +plasma.gamma = 1 +plasma.num_particles_per_cell = 4 # Moving window warpx.do_moving_window = 1 diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index 934646885..c7eb8a16c 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -71,6 +71,10 @@ public: int nSpecies() {return nspecies;} +protected: + + std::vector species_names; + private: // physical particles (+ laser) @@ -80,5 +84,6 @@ private: // runtime parameters int nspecies = 1; // physical particles only. If WarpX::use_laser, nspecies+1 == allcontainers.size(). + }; #endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index cf219ea42..ff3137246 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -13,7 +13,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) int n = WarpX::use_laser ? nspecies+1 : nspecies; allcontainers.resize(n); for (int i = 0; i < nspecies; ++i) { - allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i)); + allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i, species_names[i])); } if (WarpX::use_laser) { allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1)); @@ -31,6 +31,9 @@ MultiParticleContainer::ReadParameters () pp.query("nspecies", nspecies); BL_ASSERT(nspecies >= 0); + pp.getarr("species_names", species_names); + BL_ASSERT(species_names.size() == nspecies); + initialized = true; } } diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 541de7c98..923d40429 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -10,7 +10,7 @@ class PhysicalParticleContainer : public WarpXParticleContainer { public: - PhysicalParticleContainer (amrex::AmrCore* amr_core, int ispecies); + PhysicalParticleContainer (amrex::AmrCore* amr_core, int ispecies, const std::string& name); virtual ~PhysicalParticleContainer () {} virtual void AllocData () override; @@ -30,6 +30,8 @@ public: protected: + std::string species_name; + std::unique_ptr plasma_injector; }; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 767a41919..c7e98def9 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -8,26 +8,27 @@ using namespace amrex; -PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies) - : WarpXParticleContainer(amr_core, ispecies) +PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies, + const std::string& name) + : WarpXParticleContainer(amr_core, ispecies), + species_name(name) { - std::stringstream ss; - ss << "species" << ispecies; - ParmParse pp(ss.str()); + + ParmParse pp(species_name); 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); + plasma_injector.reset(new ConstantPlasmaInjector(species_id, species_name)); } else if (plasma_profile_s == "double_ramp") { - plasma_injector.reset(new DoubleRampPlasmaInjector); + plasma_injector.reset(new DoubleRampPlasmaInjector(species_id, species_name)); } else if (plasma_profile_s == "custom") { - plasma_injector.reset(new CustomPlasmaInjector); + plasma_injector.reset(new CustomPlasmaInjector(species_id, species_name)); } else { diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H index 5847fb4a6..82a8f0bb8 100644 --- a/Source/PlasmaInjector.H +++ b/Source/PlasmaInjector.H @@ -7,40 +7,51 @@ class PlasmaInjector { public: - PlasmaInjector(); + PlasmaInjector(int ispecies, const std::string& name); - virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) = 0; + virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) = 0; - bool insideBounds(amrex::Real x, amrex::Real y, amrex::Real z); + bool insideBounds(amrex::Real x, amrex::Real y, amrex::Real z); + + int numParticlesPerCell() {return num_particles_per_cell;}; + + amrex::Real getGamma() {return gamma;} protected: - amrex::Real xmin, xmax; - amrex::Real ymin, ymax; - amrex::Real zmin, zmax; + amrex::Real xmin, xmax; + amrex::Real ymin, ymax; + amrex::Real zmin, zmax; + + amrex::Real density; + int num_particles_per_cell; - amrex::Real density; + amrex::Real gamma; + int species_id; + std::string species_name; }; class ConstantPlasmaInjector : public PlasmaInjector { public: - amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); + ConstantPlasmaInjector(int ispecies, const std::string& name); + amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); }; class CustomPlasmaInjector : public PlasmaInjector { public: - amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); + CustomPlasmaInjector(int ispecies, const std::string& name); + amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); }; class DoubleRampPlasmaInjector : public PlasmaInjector { public: - DoubleRampPlasmaInjector(); - amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); - + DoubleRampPlasmaInjector(int ispecies, const std::string& name); + amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); + protected: - amrex::Real plateau_length; - amrex::Real ramp_length; + amrex::Real plateau_length; + amrex::Real ramp_length; }; #endif diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index 5106a9333..6735853e2 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -1,37 +1,51 @@ #include "PlasmaInjector.H" -#include +#include using namespace amrex; -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); +PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) + : species_id(ispecies), species_name(name) +{ + ParmParse pp(species_name); + + 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("num_particles_per_cell", num_particles_per_cell); + pp.get("density", density); + pp.get("gamma", gamma); } 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; + z >= zmax || z < zmin ) return false; + return true; +} + +ConstantPlasmaInjector::ConstantPlasmaInjector(int ispecies, const std::string& name) + : PlasmaInjector(ispecies, name) +{ } Real ConstantPlasmaInjector::getDensity(Real x, Real y, Real z) { return density; } -DoubleRampPlasmaInjector::DoubleRampPlasmaInjector() - : PlasmaInjector() +CustomPlasmaInjector::CustomPlasmaInjector(int ispecies, const std::string& name) + : PlasmaInjector(ispecies, name) { - ParmParse pp("plasma"); +} + +DoubleRampPlasmaInjector::DoubleRampPlasmaInjector(int ispecies, const std::string& name) + : PlasmaInjector(ispecies, name) +{ + ParmParse pp(species_name); pp.get("ramp_length", ramp_length); pp.get("plateau_length", plateau_length); } @@ -40,3 +54,4 @@ Real DoubleRampPlasmaInjector::getDensity(Real x, Real y, Real z) { return density; } + -- cgit v1.2.3 From d96a3110736101738a59b2d5c965ee62cb42b322 Mon Sep 17 00:00:00 2001 From: atmyers Date: Wed, 15 Mar 2017 14:46:12 -0700 Subject: add information about the particle velocities and mass/charge to the particle injector. --- Exec/plasma_acceleration/ParticleProb.cpp | 33 ++++++----------- Exec/plasma_acceleration/inputs | 27 ++++++++++---- Source/PlasmaInjector.H | 11 ++++-- Source/PlasmaInjector.cpp | 59 ++++++++++++++++++++++++++++++- 4 files changed, 99 insertions(+), 31 deletions(-) (limited to 'Source/PlasmaInjector.cpp') diff --git a/Exec/plasma_acceleration/ParticleProb.cpp b/Exec/plasma_acceleration/ParticleProb.cpp index acd1e1a58..010484315 100644 --- a/Exec/plasma_acceleration/ParticleProb.cpp +++ b/Exec/plasma_acceleration/ParticleProb.cpp @@ -18,21 +18,14 @@ PhysicalParticleContainer::InitData() { BL_PROFILE("PhysicalParticleContainer::InitData()"); - // species_id 0 : the beam - // species_id 1 : the electrons of the plasma - // Note: the ions of the plasma are implicitly motionless, and so are not part of the simulation - if (species_id == 0 or species_id == 1) { - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - } else { - amrex::Abort("PhysicalParticleContainer::InitData(): species_id must be 0 or 1"); - } + charge = plasma_injector->getCharge(); + mass = plasma_injector->getMass(); const int lev = 0; const Geometry& geom = Geom(lev); const Real* dx = geom.CellSize(); - Real weight, gamma, uz, scale_fac; + Real scale_fac; int n_part_per_cell = plasma_injector->numParticlesPerCell(); #if BL_SPACEDIM==3 @@ -41,14 +34,6 @@ PhysicalParticleContainer::InitData() scale_face = dx[0]*dx[1]/n_part_per_cell; #endif - // Calculate the limits between which the particles will be initialized - uz = 0.; - if (species_id == 0) { // relativistic beam - gamma = plasma_injector->getGamma(); - uz = std::sqrt( gamma*gamma - 1 ); - } - uz *= PhysConst::c; - std::array attribs; attribs.fill(0.0); for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { @@ -74,14 +59,18 @@ PhysicalParticleContainer::InitData() Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; #endif - if (plasma_injector->insideBounds(x, y, z)) - { + if (plasma_injector->insideBounds(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::uz] = uz; + 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); } - } + } } } } diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 0fad702b9..45832d536 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -31,24 +31,34 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 - # Information about the particle species particles.nspecies = 2 particles.species_names = beam plasma -# The beam -beam.profile = "constant" +# 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.gamma = 1e9 beam.num_particles_per_cell = 4 -# The plasma +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 @@ -56,10 +66,15 @@ 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.gamma = 1 plasma.num_particles_per_cell = 4 +plasma.ux = 0.0 +plasma.uy = 0.0 +plasma.uz = 0.0 + # Moving window warpx.do_moving_window = 1 warpx.moving_window_dir = z diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H index 82a8f0bb8..bcdb6d4dd 100644 --- a/Source/PlasmaInjector.H +++ b/Source/PlasmaInjector.H @@ -15,10 +15,15 @@ public: int numParticlesPerCell() {return num_particles_per_cell;}; - amrex::Real getGamma() {return gamma;} + void getMomentum(amrex::Real* u); + + amrex::Real getCharge() {return charge;} + amrex::Real getMass() {return mass;} protected: + amrex::Real mass, charge; + amrex::Real xmin, xmax; amrex::Real ymin, ymax; amrex::Real zmin, zmax; @@ -26,7 +31,9 @@ protected: amrex::Real density; int num_particles_per_cell; - amrex::Real gamma; + amrex::Real ux; + amrex::Real uy; + amrex::Real uz; int species_id; std::string species_name; diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index 6735853e2..a30b46f4d 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -2,13 +2,58 @@ #include +#include +#include + using namespace amrex; +Real parseChargeName(const std::string& name) { + if (name == "q_e") { + return PhysConst::q_e; + } else { + std::stringstream stringstream; + std::string string; + stringstream << "Charge string '" << name << "' not recognized."; + string = stringstream.str(); + amrex::Abort(string.c_str()); + return 0.0; + } +} + +Real parseChargeString(const std::string& name) { + if(name.substr(0, 1) == "-") + return -1.0 * parseChargeName(name.substr(1, name.size() - 1)); + return parseChargeName(name); +} + +Real parseMassString(const std::string& name) { + if (name == "m_e") { + return PhysConst::m_e; + } else if (name == "m_p"){ + return PhysConst::m_p; + } else { + std::stringstream stringstream; + std::string string; + stringstream << "Mass string '" << name << "' not recognized."; + string = stringstream.str(); + amrex::Abort(string.c_str()); + return 0.0; + } +} + PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) : species_id(ispecies), species_name(name) { ParmParse pp(species_name); + std::string charge_s; + pp.get("charge", charge_s); + charge = parseChargeString(charge_s); + + std::string mass_s; + pp.get("mass", mass_s); + mass = parseMassString(mass_s); + pp.get("xmin", xmin); pp.get("ymin", ymin); pp.get("zmin", zmin); @@ -18,7 +63,19 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.get("num_particles_per_cell", num_particles_per_cell); pp.get("density", density); - pp.get("gamma", gamma); + + ux = 0.; + uy = 0.; + uz = 0.; + pp.query("ux", ux); + pp.query("uy", uy); + pp.query("uz", uz); +} + +void PlasmaInjector::getMomentum(Real* u) { + u[0] = ux * PhysConst::c; + u[1] = uy * PhysConst::c; + u[2] = uz * PhysConst::c; } bool PlasmaInjector::insideBounds(Real x, Real y, Real z) { -- cgit v1.2.3 From 98321e065b39e6c42b432ee85245cb3ec354e416 Mon Sep 17 00:00:00 2001 From: atmyers Date: Wed, 15 Mar 2017 15:32:04 -0700 Subject: convert to lower case when parsing the charge/mass strings. --- Source/PlasmaInjector.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) (limited to 'Source/PlasmaInjector.cpp') diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index a30b46f4d..562388eb4 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -48,10 +48,18 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) std::string charge_s; pp.get("charge", charge_s); + std::transform(charge_s.begin(), + charge_s.end(), + charge_s.begin(), + ::tolower); charge = parseChargeString(charge_s); std::string mass_s; pp.get("mass", mass_s); + std::transform(mass_s.begin(), + mass_s.end(), + mass_s.begin(), + ::tolower); mass = parseMassString(mass_s); pp.get("xmin", xmin); -- cgit v1.2.3 From 1098d924b6ff68ba401a1783c8d82dcde55b6019 Mon Sep 17 00:00:00 2001 From: atmyers Date: Wed, 15 Mar 2017 17:14:11 -0700 Subject: add an additional momentum distribution type and have uniform plasma test use it. --- Exec/Langmuir/inputs | 1 + Exec/Langmuir/inputs.lb | 1 + Exec/Langmuir/inputs.multi.rt | 2 + Exec/Langmuir/inputs.nolb | 1 + Exec/Langmuir/inputs.rt | 1 + Exec/plasma_acceleration/inputs | 14 ++++- Exec/uniform_plasma/ParticleProb.cpp | 95 +--------------------------------- Exec/uniform_plasma/inputs | 25 +++++++-- Exec/uniform_plasma/inputs.rt | 14 +++-- Source/PlasmaInjector.H | 47 +++++++++++++++-- Source/PlasmaInjector.cpp | 99 +++++++++++++++++++++++++++++------- 11 files changed, 176 insertions(+), 124 deletions(-) (limited to 'Source/PlasmaInjector.cpp') diff --git a/Exec/Langmuir/inputs b/Exec/Langmuir/inputs index 6b43860a0..dbbdc4c0c 100644 --- a/Exec/Langmuir/inputs +++ b/Exec/Langmuir/inputs @@ -48,4 +48,5 @@ electrons.zmax = 20.e-6 electrons.num_particles_per_cell = 10 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/Langmuir/inputs.lb b/Exec/Langmuir/inputs.lb index ada715058..349fd28eb 100644 --- a/Exec/Langmuir/inputs.lb +++ b/Exec/Langmuir/inputs.lb @@ -54,4 +54,5 @@ electrons.zmax = 20.e-6 electrons.num_particles_per_cell = 10 electrons.profile = constant electrons.density = 1.e25 # number of electrons per m^3 +electrons.momentum_distribution_type = "constant" electrons.ux = 1.0 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.multi.rt b/Exec/Langmuir/inputs.multi.rt index 259812838..5a2a0c6a2 100644 --- a/Exec/Langmuir/inputs.multi.rt +++ b/Exec/Langmuir/inputs.multi.rt @@ -51,6 +51,7 @@ electrons.zmax = 20.e-6 electrons.num_particles_per_cell = 10 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 electrons.uy = 0. # uy = gamma*beta_y electrons.uz = 0. # uz = gamma*beta_z @@ -66,6 +67,7 @@ positrons.zmax = 20.e-6 positrons.num_particles_per_cell = 10 positrons.profile = constant positrons.density = 1.e25 # number of positrons per m^3 +positrons.momentum_distribution_type = "constant" positrons.ux = 0.01 # ux = gamma*beta_x positrons.uy = 0. # uy = gamma*beta_y positrons.uz = 0. # uz = gamma*beta_z diff --git a/Exec/Langmuir/inputs.nolb b/Exec/Langmuir/inputs.nolb index 028e30985..cc11302ab 100644 --- a/Exec/Langmuir/inputs.nolb +++ b/Exec/Langmuir/inputs.nolb @@ -54,4 +54,5 @@ electrons.zmax = 20.e-6 electrons.num_particles_per_cell = 10 electrons.profile = constant electrons.density = 1.e25 # number of electrons per m^3 +electrons.momentum_distribution_type = "constant" electrons.ux = 1.0 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.rt b/Exec/Langmuir/inputs.rt index 70eed81f8..533da114a 100644 --- a/Exec/Langmuir/inputs.rt +++ b/Exec/Langmuir/inputs.rt @@ -50,4 +50,5 @@ electrons.zmax = 20.e-6 electrons.num_particles_per_cell = 10 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/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 241fc705a..096b4d430 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -35,36 +35,48 @@ warpx.cfl = 1.0 particles.nspecies = 2 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.momentum_distribution_type = "constant" 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 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 +plasma.momentum_distribution_type = "constant" plasma.ux = 0.0 plasma.uy = 0.0 plasma.uz = 0.0 diff --git a/Exec/uniform_plasma/ParticleProb.cpp b/Exec/uniform_plasma/ParticleProb.cpp index cef618f61..0114631c8 100644 --- a/Exec/uniform_plasma/ParticleProb.cpp +++ b/Exec/uniform_plasma/ParticleProb.cpp @@ -17,98 +17,5 @@ using namespace amrex; void PhysicalParticleContainer::InitData() { - BL_PROFILE("PPC::InitData()"); - - // species_id 0 : electrons - // species_id 1 : Hydrogen ions - // Note: the ions of the plasma are implicitly motionless, and so are not part of the simulation - if (species_id == 0) { - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - } else { - charge = PhysConst::q_e; - mass = PhysConst::m_p; - } - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, u_th, ux_m, uy_m, uz_m; - Real particle_shift_x, particle_shift_y, particle_shift_z; - int n_part_per_cell; - { - ParmParse pp("uniform_plasma"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #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 - - u_th = 0.; - ux_m = 0.; - uy_m = 0.; - uz_m = 0.; - pp.query("u_th", u_th); - pp.query("ux_m", ux_m); - pp.query("uy_m", uy_m); - pp.query("uz_m", uz_m); - u_th *= PhysConst::c; - ux_m *= PhysConst::c; - uy_m *= PhysConst::c; - uz_m *= PhysConst::c; - } - - std::array attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - - // Initialize random generator for normal distribution - std::default_random_engine generator; - std::normal_distribution velocity_distribution(0.0, u_th); - - 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_part + #include "AMReX_Real.H" #include "AMReX_ParmParse.H" +class PlasmaMomentumDistribution +{ +public: + virtual void getMomentum(amrex::Real* u) = 0; +}; + +class ConstantMomentumDistribution : public PlasmaMomentumDistribution +{ +public: + ConstantMomentumDistribution(amrex::Real ux, + amrex::Real uy, + amrex::Real uz); + void getMomentum(amrex::Real* u); + +private: + amrex::Real _ux; + amrex::Real _uy; + amrex::Real _uz; +}; + +class GaussianRandomMomentumDistribution : public PlasmaMomentumDistribution +{ +public: + GaussianRandomMomentumDistribution(amrex::Real ux_m, + amrex::Real uy_m, + amrex::Real uz_m, + amrex::Real u_th); + void getMomentum(amrex::Real* u); +private: + amrex::Real _ux_m; + amrex::Real _uy_m; + amrex::Real _uz_m; + amrex::Real _u_th; + + std::normal_distribution momentum_distribution; + std::default_random_engine generator; +}; + class PlasmaInjector { + public: PlasmaInjector(int ispecies, const std::string& name); @@ -31,12 +72,10 @@ protected: amrex::Real density; int num_particles_per_cell; - amrex::Real ux; - amrex::Real uy; - amrex::Real uz; - int species_id; std::string species_name; + + std::unique_ptr mom_dist; }; class ConstantPlasmaInjector : public PlasmaInjector { diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index 562388eb4..cc396d16f 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -41,6 +41,37 @@ Real parseMassString(const std::string& name) { } } +ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, + Real uy, + Real uz) + : _ux(ux), _uy(uy), _uz(uz) +{ +} + +void ConstantMomentumDistribution::getMomentum(Real* u) { + u[0] = _ux; + u[1] = _uy; + u[2] = _uz; +} + +GaussianRandomMomentumDistribution::GaussianRandomMomentumDistribution(Real ux_m, + Real uy_m, + Real uz_m, + Real u_th) + : _ux_m(ux_m), _uy_m(uy_m), _uz_m(uz_m), _u_th(u_th), + momentum_distribution(0.0, u_th) +{ +} + +void GaussianRandomMomentumDistribution::getMomentum(Real* u) { + Real ux_th = momentum_distribution(generator); + Real uy_th = momentum_distribution(generator); + Real uz_th = momentum_distribution(generator); + u[0] = _ux_m + ux_th; + u[1] = _uy_m + uy_th; + u[2] = _uz_m + uz_th; +} + PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) : species_id(ispecies), species_name(name) { @@ -61,29 +92,63 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) mass_s.begin(), ::tolower); mass = parseMassString(mass_s); - - pp.get("xmin", xmin); - pp.get("ymin", ymin); - pp.get("zmin", zmin); - pp.get("xmax", xmax); - pp.get("ymax", ymax); - pp.get("zmax", zmax); + + xmin = std::numeric_limits::lowest(); + ymin = std::numeric_limits::lowest(); + zmin = std::numeric_limits::lowest(); + + xmax = std::numeric_limits::max(); + ymax = std::numeric_limits::max(); + zmax = std::numeric_limits::max(); + + pp.query("xmin", xmin); + pp.query("ymin", ymin); + pp.query("zmin", zmin); + pp.query("xmax", xmax); + pp.query("ymax", ymax); + pp.query("zmax", zmax); pp.get("num_particles_per_cell", num_particles_per_cell); pp.get("density", density); - ux = 0.; - uy = 0.; - uz = 0.; - pp.query("ux", ux); - pp.query("uy", uy); - pp.query("uz", uz); + std::string mom_dist_s; + pp.get("momentum_distribution_type", mom_dist_s); + std::transform(mom_dist_s.begin(), + mom_dist_s.end(), + mom_dist_s.begin(), + ::tolower); + if (mom_dist_s == "constant") { + Real ux = 0.; + Real uy = 0.; + Real uz = 0.; + pp.query("ux", ux); + pp.query("uy", uy); + pp.query("uz", uz); + mom_dist.reset(new ConstantMomentumDistribution(ux, uy, uz)); + } else if (mom_dist_s == "gaussian") { + Real ux_m = 0.; + Real uy_m = 0.; + Real uz_m = 0.; + Real u_th = 0.; + pp.query("ux_m", ux_m); + pp.query("uy_m", uy_m); + pp.query("uz_m", uz_m); + pp.query("u_th", u_th); + mom_dist.reset(new GaussianRandomMomentumDistribution(ux_m, uy_m, uz_m, u_th)); + } else { + std::stringstream stringstream; + std::string string; + stringstream << "Momentum distribution type '" << mom_dist_s << "' not recognized."; + string = stringstream.str(); + amrex::Abort(string.c_str()); + } } void PlasmaInjector::getMomentum(Real* u) { - u[0] = ux * PhysConst::c; - u[1] = uy * PhysConst::c; - u[2] = uz * PhysConst::c; + mom_dist->getMomentum(u); + u[0] *= PhysConst::c; + u[1] *= PhysConst::c; + u[2] *= PhysConst::c; } bool PlasmaInjector::insideBounds(Real x, Real y, Real z) { @@ -118,5 +183,3 @@ DoubleRampPlasmaInjector::DoubleRampPlasmaInjector(int ispecies, const std::stri Real DoubleRampPlasmaInjector::getDensity(Real x, Real y, Real z) { return density; } - - -- cgit v1.2.3 From 737c39407f9f743de8177996acf6ef788dffd214 Mon Sep 17 00:00:00 2001 From: atmyers Date: Thu, 16 Mar 2017 11:14:44 -0700 Subject: Implement a couple of different density types. --- Source/CustomDensityProb.cpp | 5 +- Source/ParticleProb.cpp | 21 -------- Source/PhysicalParticleContainer.H | 4 +- Source/PhysicalParticleContainer.cpp | 23 +-------- Source/PlasmaInjector.H | 91 +++++++++++++++++++++++++---------- Source/PlasmaInjector.cpp | 93 +++++++++++++++++++----------------- 6 files changed, 123 insertions(+), 114 deletions(-) delete mode 100644 Source/ParticleProb.cpp (limited to 'Source/PlasmaInjector.cpp') diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp index 9b7298ee7..11d80e0b5 100644 --- a/Source/CustomDensityProb.cpp +++ b/Source/CustomDensityProb.cpp @@ -4,7 +4,10 @@ #include -amrex::Real CustomPlasmaInjector::getDensity(amrex::Real x, amrex::Real y, amrex::Real z) { +amrex::Real CustomDensityDistribution::getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const +{ amrex::Abort("If running with a custom density profile, you must supply a CustomDensityProb.cpp file"); return 0.0; } diff --git a/Source/ParticleProb.cpp b/Source/ParticleProb.cpp deleted file mode 100644 index 4abec6ea8..000000000 --- a/Source/ParticleProb.cpp +++ /dev/null @@ -1,21 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include - -#include - -#include -#include - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - static_assert(false, - "Each problem must have its own version of PhysicalParticleContainer::InitData()"); -} diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index a399df3fe..e4a5bfdef 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -47,7 +47,9 @@ public: protected: - void InitNPerCell (); + void InitNDiagPerCell (); + void InitNRandomUniformPerCell (); + void InitNRandom (); std::string species_name; std::unique_ptr plasma_injector; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 0e3875ca1..65c59f845 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -2,7 +2,6 @@ #include #include -#include #include #include @@ -13,27 +12,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp : WarpXParticleContainer(amr_core, ispecies), species_name(name) { - - ParmParse pp(species_name); - - 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(species_id, species_name)); - } - - else if (plasma_profile_s == "double_ramp") { - plasma_injector.reset(new DoubleRampPlasmaInjector(species_id, species_name)); - } - - else if (plasma_profile_s == "custom") { - plasma_injector.reset(new CustomPlasmaInjector(species_id, species_name)); - } - - else { - amrex::Abort("Unknown plasma injector type"); - } + plasma_injector.reset(new PlasmaInjector(species_id, species_name)); } void diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H index dfaa264c8..76f11be3e 100644 --- a/Source/PlasmaInjector.H +++ b/Source/PlasmaInjector.H @@ -6,19 +6,71 @@ #include "AMReX_Real.H" #include "AMReX_ParmParse.H" +/// +/// PlasmaDensityDistribution describes how the charge density +/// is set in particle initialization. Subclasses must define a +/// getDensity function that describes the charge density as a +/// function of x, y, and z. +/// +class PlasmaDensityDistribution +{ +public: + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const = 0; +}; + +/// +/// This describes a constant density distribution. +/// +class ConstantDensityDistribution : public PlasmaDensityDistribution +{ +public: + ConstantDensityDistribution(amrex::Real density); + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const override; + +private: + amrex::Real _density; +}; + +/// +/// This describes a custom density distribution. Users can supply +/// a getDensity function. +/// +/// +class CustomDensityDistribution : public PlasmaDensityDistribution +{ +public: + CustomDensityDistribution(amrex::Real density); + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const override; +}; + +/// +/// PlasmaMomentumDistribution describes how the particle momenta +/// are set. Subclasses must define a getMomentum method that fills +/// a u with the 3 components of the particle momentum +/// class PlasmaMomentumDistribution { public: virtual void getMomentum(amrex::Real* u) = 0; }; +/// +/// This is a constant momentum distribution - all particles will +/// have the same ux, uy, and uz +/// class ConstantMomentumDistribution : public PlasmaMomentumDistribution { public: ConstantMomentumDistribution(amrex::Real ux, amrex::Real uy, amrex::Real uz); - void getMomentum(amrex::Real* u); + virtual void getMomentum(amrex::Real* u) override; private: amrex::Real _ux; @@ -26,6 +78,13 @@ private: amrex::Real _uz; }; +/// +/// This is a Gaussian Random momentum distribution. +/// Particles will get random momenta, drawn from a normal. +/// ux_m, ux_y, and ux_z describe the mean components in the x, y, and z +/// directions, while u_th is the standard deviation of the random +/// component. +/// class GaussianRandomMomentumDistribution : public PlasmaMomentumDistribution { public: @@ -33,7 +92,7 @@ public: amrex::Real uy_m, amrex::Real uz_m, amrex::Real u_th); - void getMomentum(amrex::Real* u); + virtual void getMomentum(amrex::Real* u) override; private: amrex::Real _ux_m; amrex::Real _uy_m; @@ -50,7 +109,7 @@ class PlasmaInjector public: PlasmaInjector(int ispecies, const std::string& name); - virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) = 0; + amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z); bool insideBounds(amrex::Real x, amrex::Real y, amrex::Real z); @@ -74,30 +133,10 @@ protected: int species_id; std::string species_name; - - std::unique_ptr mom_dist; -}; - -class ConstantPlasmaInjector : public PlasmaInjector { -public: - ConstantPlasmaInjector(int ispecies, const std::string& name); - amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); -}; - -class CustomPlasmaInjector : public PlasmaInjector { -public: - CustomPlasmaInjector(int ispecies, const std::string& name); - amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); -}; - -class DoubleRampPlasmaInjector : public PlasmaInjector { -public: - DoubleRampPlasmaInjector(int ispecies, const std::string& name); - amrex::Real getDensity( amrex::Real x, amrex::Real y, amrex::Real z); + std::string injection_style; -protected: - amrex::Real plateau_length; - amrex::Real ramp_length; + std::unique_ptr rho_dist; + std::unique_ptr mom_dist; }; #endif diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index cc396d16f..acde1fa69 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -7,15 +7,20 @@ using namespace amrex; +void StringParseAbortMessage(const std::string& var, + const std::string& name) { + std::stringstream stringstream; + std::string string; + stringstream << var << " string '" << name << "' not recognized."; + string = stringstream.str(); + amrex::Abort(string.c_str()); +} + Real parseChargeName(const std::string& name) { if (name == "q_e") { return PhysConst::q_e; } else { - std::stringstream stringstream; - std::string string; - stringstream << "Charge string '" << name << "' not recognized."; - string = stringstream.str(); - amrex::Abort(string.c_str()); + StringParseAbortMessage("Charge", name); return 0.0; } } @@ -32,21 +37,25 @@ Real parseMassString(const std::string& name) { } else if (name == "m_p"){ return PhysConst::m_p; } else { - std::stringstream stringstream; - std::string string; - stringstream << "Mass string '" << name << "' not recognized."; - string = stringstream.str(); - amrex::Abort(string.c_str()); + StringParseAbortMessage("Mass", name); return 0.0; } } +ConstantDensityDistribution::ConstantDensityDistribution(Real density) + : _density(density) +{} + +Real ConstantDensityDistribution::getDensity(Real x, Real y, Real z) const +{ + return _density; +} + ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, Real uy, Real uz) : _ux(ux), _uy(uy), _uz(uz) -{ -} +{} void ConstantMomentumDistribution::getMomentum(Real* u) { u[0] = _ux; @@ -77,6 +86,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) { ParmParse pp(species_name); + // parse charge and mass std::string charge_s; pp.get("charge", charge_s); std::transform(charge_s.begin(), @@ -93,6 +103,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) ::tolower); mass = parseMassString(mass_s); + // parse plasma boundaries xmin = std::numeric_limits::lowest(); ymin = std::numeric_limits::lowest(); zmin = std::numeric_limits::lowest(); @@ -107,10 +118,24 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.query("xmax", xmax); pp.query("ymax", ymax); pp.query("zmax", zmax); - + + // parse density information + std::string rho_dist_s; + pp.get("profile", rho_dist_s); + std::transform(rho_dist_s.begin(), + rho_dist_s.end(), + rho_dist_s.begin(), + ::tolower); + if (rho_dist_s == "constant") { + Real density; + pp.get("density", density); + rho_dist.reset(new ConstantDensityDistribution(density)); + } else { + StringParseAbortMessage("Density distribution type", rho_dist_s); + } pp.get("num_particles_per_cell", num_particles_per_cell); - pp.get("density", density); + // parse momentum information std::string mom_dist_s; pp.get("momentum_distribution_type", mom_dist_s); std::transform(mom_dist_s.begin(), @@ -136,11 +161,15 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.query("u_th", u_th); mom_dist.reset(new GaussianRandomMomentumDistribution(ux_m, uy_m, uz_m, u_th)); } else { - std::stringstream stringstream; - std::string string; - stringstream << "Momentum distribution type '" << mom_dist_s << "' not recognized."; - string = stringstream.str(); - amrex::Abort(string.c_str()); + StringParseAbortMessage("Momentum distribution type", mom_dist_s); + } + + // get injection style + pp.get("injection_style", injection_style); + if (injection_style != "nrandom" or + injection_style != "nrandompercell" or + injection_style != "ndiagpercell") { + StringParseAbortMessage("Injection style", injection_style); } } @@ -158,28 +187,6 @@ bool PlasmaInjector::insideBounds(Real x, Real y, Real z) { return true; } -ConstantPlasmaInjector::ConstantPlasmaInjector(int ispecies, const std::string& name) - : PlasmaInjector(ispecies, name) -{ -} - -Real ConstantPlasmaInjector::getDensity(Real x, Real y, Real z) { - return density; -} - -CustomPlasmaInjector::CustomPlasmaInjector(int ispecies, const std::string& name) - : PlasmaInjector(ispecies, name) -{ -} - -DoubleRampPlasmaInjector::DoubleRampPlasmaInjector(int ispecies, const std::string& name) - : PlasmaInjector(ispecies, name) -{ - ParmParse pp(species_name); - pp.get("ramp_length", ramp_length); - pp.get("plateau_length", plateau_length); -} - -Real DoubleRampPlasmaInjector::getDensity(Real x, Real y, Real z) { - return density; +Real PlasmaInjector::getDensity(Real x, Real y, Real z) { + return rho_dist->getDensity(x, y, z); } -- cgit v1.2.3 From b1e7fa6f54d0811c06165b1728064ebf5069c851 Mon Sep 17 00:00:00 2001 From: atmyers Date: Thu, 16 Mar 2017 13:09:44 -0700 Subject: Plasma injector stuff almost complete. --- Exec/Langmuir/ParticleProb.cpp | 20 ------ Exec/Langmuir/inputs | 4 ++ Exec/Langmuir/inputs.lb | 4 ++ Exec/Langmuir/inputs.multi.rt | 8 +++ Exec/Langmuir/inputs.nolb | 4 ++ Exec/Langmuir/inputs.rt | 4 ++ Exec/laser_injection/ParticleProb.cpp | 99 -------------------------- Exec/plasma_acceleration/CustomDensityProb.cpp | 7 +- Exec/plasma_acceleration/ParticleProb.cpp | 20 ------ Exec/plasma_acceleration/inputs | 6 +- Source/CustomDensityProb.cpp | 6 +- Source/Make.package | 2 +- Source/ParticleContainer.cpp | 6 +- Source/PhysicalParticleContainer.H | 2 +- Source/PhysicalParticleContainer.cpp | 89 ++++++++++++++++++++++- Source/PlasmaInjector.H | 28 +++++--- Source/PlasmaInjector.cpp | 44 ++++++++---- 17 files changed, 177 insertions(+), 176 deletions(-) delete mode 100644 Exec/Langmuir/ParticleProb.cpp delete mode 100644 Exec/laser_injection/ParticleProb.cpp delete mode 100644 Exec/plasma_acceleration/ParticleProb.cpp (limited to 'Source/PlasmaInjector.cpp') diff --git a/Exec/Langmuir/ParticleProb.cpp b/Exec/Langmuir/ParticleProb.cpp deleted file mode 100644 index 806861a19..000000000 --- a/Exec/Langmuir/ParticleProb.cpp +++ /dev/null @@ -1,20 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include - -#include - -#include -#include - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - InitNPerCell(); -} diff --git a/Exec/Langmuir/inputs b/Exec/Langmuir/inputs index dbbdc4c0c..292cd19de 100644 --- a/Exec/Langmuir/inputs +++ b/Exec/Langmuir/inputs @@ -39,14 +39,18 @@ particles.species_names = electrons electrons.charge = -q_e electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + electrons.xmin = -20.e-6 electrons.xmax = 0.e-6 electrons.ymin = -20.e-6 electrons.ymax = 20.e-6 electrons.zmin = -20.e-6 electrons.zmax = 20.e-6 + electrons.num_particles_per_cell = 10 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/Langmuir/inputs.lb b/Exec/Langmuir/inputs.lb index 349fd28eb..5652735b6 100644 --- a/Exec/Langmuir/inputs.lb +++ b/Exec/Langmuir/inputs.lb @@ -45,14 +45,18 @@ particles.species_names = electrons electrons.charge = -q_e electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + electrons.xmin = -20.e-6 electrons.xmax = 0.e-6 electrons.ymin = -20.e-6 electrons.ymax = 20.e-6 electrons.zmin = -20.e-6 electrons.zmax = 20.e-6 + electrons.num_particles_per_cell = 10 electrons.profile = constant electrons.density = 1.e25 # number of electrons per m^3 + electrons.momentum_distribution_type = "constant" electrons.ux = 1.0 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.multi.rt b/Exec/Langmuir/inputs.multi.rt index 5a2a0c6a2..1256652e3 100644 --- a/Exec/Langmuir/inputs.multi.rt +++ b/Exec/Langmuir/inputs.multi.rt @@ -42,15 +42,19 @@ particles.species_names = electrons positrons electrons.charge = -q_e electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + electrons.xmin = -20.e-6 electrons.xmax = 0.e-6 electrons.ymin = -20.e-6 electrons.ymax = 20.e-6 electrons.zmin = -20.e-6 electrons.zmax = 20.e-6 + electrons.num_particles_per_cell = 10 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 electrons.uy = 0. # uy = gamma*beta_y @@ -58,15 +62,19 @@ electrons.uz = 0. # uz = gamma*beta_z positrons.charge = q_e positrons.mass = m_e +positrons.injection_style = "NDiagPerCell" + positrons.xmin = -20.e-6 positrons.xmax = 0.e-6 positrons.ymin = -20.e-6 positrons.ymax = 20.e-6 positrons.zmin = -20.e-6 positrons.zmax = 20.e-6 + positrons.num_particles_per_cell = 10 positrons.profile = constant positrons.density = 1.e25 # number of positrons per m^3 + positrons.momentum_distribution_type = "constant" positrons.ux = 0.01 # ux = gamma*beta_x positrons.uy = 0. # uy = gamma*beta_y diff --git a/Exec/Langmuir/inputs.nolb b/Exec/Langmuir/inputs.nolb index cc11302ab..6ca55b3c7 100644 --- a/Exec/Langmuir/inputs.nolb +++ b/Exec/Langmuir/inputs.nolb @@ -45,14 +45,18 @@ particles.species_names = electrons electrons.charge = -q_e electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + electrons.xmin = -20.e-6 electrons.xmax = 0.e-6 electrons.ymin = -20.e-6 electrons.ymax = 20.e-6 electrons.zmin = -20.e-6 electrons.zmax = 20.e-6 + electrons.num_particles_per_cell = 10 electrons.profile = constant electrons.density = 1.e25 # number of electrons per m^3 + electrons.momentum_distribution_type = "constant" electrons.ux = 1.0 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.rt b/Exec/Langmuir/inputs.rt index 533da114a..feeb48c93 100644 --- a/Exec/Langmuir/inputs.rt +++ b/Exec/Langmuir/inputs.rt @@ -41,14 +41,18 @@ particles.species_names = electrons electrons.charge = -q_e electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + electrons.xmin = -20.e-6 electrons.xmax = 0.e-6 electrons.ymin = -20.e-6 electrons.ymax = 20.e-6 electrons.zmin = -20.e-6 electrons.zmax = 20.e-6 + electrons.num_particles_per_cell = 10 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/laser_injection/ParticleProb.cpp b/Exec/laser_injection/ParticleProb.cpp deleted file mode 100644 index 35376c4a2..000000000 --- a/Exec/laser_injection/ParticleProb.cpp +++ /dev/null @@ -1,99 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include - -#include - -#include -#include -#include - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("PPC::InitData()"); - - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, u_th; - int n_part_per_cell; - { - ParmParse pp("uniform_plasma"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #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 - - u_th = 0.; - pp.query("u_th", u_th); - u_th *= PhysConst::c; - } - - std::array attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - - // Initialize random generator for normal distribution - std::default_random_engine generator; - std::normal_distribution velocity_distribution(0.0, u_th); - std::uniform_real_distribution position_distribution(0.0,1.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_part - -#include - -#include -#include - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - InitNPerCell(); -} diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 096b4d430..978a76048 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -41,6 +41,7 @@ particles.species_names = beam plasma beam.charge = -q_e beam.mass = m_e +beam.injection_style = "NDiagPerCell" beam.xmin = -20.e-6 beam.xmax = 20.e-6 @@ -64,6 +65,7 @@ beam.uz = 1e9 plasma.charge = -q_e plasma.mass = m_e +plasma.injection_style = "NDiagPerCell" plasma.xmin = -200.e-6 plasma.xmax = 200.e-6 @@ -72,8 +74,8 @@ 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.profile = "custom" +plasma.custom_profile_params = 1e22 # number of particles per m^3 plasma.num_particles_per_cell = 4 plasma.momentum_distribution_type = "constant" diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp index 11d80e0b5..1df3d75ad 100644 --- a/Source/CustomDensityProb.cpp +++ b/Source/CustomDensityProb.cpp @@ -4,9 +4,9 @@ #include -amrex::Real CustomDensityDistribution::getDensity(amrex::Real x, - amrex::Real y, - amrex::Real z) const +amrex::Real CustomDensityProfile::getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const { amrex::Abort("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 181f046b2..23968c898 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -6,7 +6,7 @@ CEXE_sources += WarpXWrappers.cpp CEXE_sources += WarpX.cpp WarpXInitData.cpp WarpXEvolve.cpp WarpXIO.cpp WarpXProb.cpp WarpXRegrid.cpp -CEXE_sources += ParticleProb.cpp ParticleIO.cpp +CEXE_sources += ParticleIO.cpp CEXE_sources += ParticleContainer.cpp WarpXParticleContainer.cpp PhysicalParticleContainer.cpp LaserParticleContainer.cpp CEXE_headers += WarpXWrappers.h diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index ff3137246..3ec09ac20 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -31,8 +31,10 @@ MultiParticleContainer::ReadParameters () pp.query("nspecies", nspecies); BL_ASSERT(nspecies >= 0); - pp.getarr("species_names", species_names); - BL_ASSERT(species_names.size() == nspecies); + if (nspecies > 0) { + pp.getarr("species_names", species_names); + BL_ASSERT(species_names.size() == nspecies); + } initialized = true; } diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index e4a5bfdef..98d15e22c 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -49,7 +49,7 @@ protected: void InitNDiagPerCell (); void InitNRandomUniformPerCell (); - void InitNRandom (); + void InitNRandomNormal (); std::string species_name; std::unique_ptr plasma_injector; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 65c59f845..9964f0e11 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -15,6 +15,15 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp plasma_injector.reset(new PlasmaInjector(species_id, species_name)); } +void PhysicalParticleContainer::InitData() { + if (plasma_injector->injection_style == "ndiagpercell") + InitNDiagPerCell(); + else if (plasma_injector->injection_style == "nrandomuniformpercell") + InitNRandomUniformPerCell(); + else if (plasma_injector->injection_style == "nrandomnormal") + InitNRandomNormal(); +} + void PhysicalParticleContainer::AllocData () { @@ -25,8 +34,84 @@ PhysicalParticleContainer::AllocData () } void -PhysicalParticleContainer::InitNPerCell () { - BL_PROFILE("PhysicalParticleContainer::InitNPerCell()"); +PhysicalParticleContainer::InitNRandomUniformPerCell () { + BL_PROFILE("PhysicalParticleContainer::InitNRandomPerCell()"); + + 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); + + // Initialize random generator for normal distribution + std::default_random_engine generator; + std::uniform_real_distribution position_distribution(0.0,1.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::InitNRandomNormal () { + BL_PROFILE("PhysicalParticleContainer::InitNRandomNormal()"); + amrex::Abort("PhysicalParticleContainer::InitNRandomNormal() not implemented yet."); +} + +void +PhysicalParticleContainer::InitNDiagPerCell () { + BL_PROFILE("PhysicalParticleContainer::InitNDiagPerCell()"); charge = plasma_injector->getCharge(); mass = plasma_injector->getMass(); diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H index 76f11be3e..800c3d27d 100644 --- a/Source/PlasmaInjector.H +++ b/Source/PlasmaInjector.H @@ -7,26 +7,28 @@ #include "AMReX_ParmParse.H" /// -/// PlasmaDensityDistribution describes how the charge density +/// PlasmaDensityProfile describes how the charge density /// is set in particle initialization. Subclasses must define a /// getDensity function that describes the charge density as a /// function of x, y, and z. /// -class PlasmaDensityDistribution +class PlasmaDensityProfile { public: virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) const = 0; +protected: + std::string _species_name; }; /// /// This describes a constant density distribution. /// -class ConstantDensityDistribution : public PlasmaDensityDistribution +class ConstantDensityProfile : public PlasmaDensityProfile { public: - ConstantDensityDistribution(amrex::Real density); + ConstantDensityProfile(amrex::Real _density); virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) const override; @@ -37,16 +39,18 @@ private: /// /// This describes a custom density distribution. Users can supply -/// a getDensity function. +/// in their problem directory. /// /// -class CustomDensityDistribution : public PlasmaDensityDistribution +class CustomDensityProfile : public PlasmaDensityProfile { public: - CustomDensityDistribution(amrex::Real density); + CustomDensityProfile(const std::string& species_name); virtual amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z) const override; +private: + std::vector params; }; /// @@ -103,6 +107,11 @@ private: std::default_random_engine generator; }; +/// +/// The PlasmaInjector class parses and stores information about the plasma +/// type used in the particle container. This information is used to create the +/// particles on initialization and whenever the window moves. +/// class PlasmaInjector { @@ -120,6 +129,8 @@ public: amrex::Real getCharge() {return charge;} amrex::Real getMass() {return mass;} + std::string injection_style; + protected: amrex::Real mass, charge; @@ -133,9 +144,8 @@ protected: int species_id; std::string species_name; - std::string injection_style; - std::unique_ptr rho_dist; + std::unique_ptr rho_prof; std::unique_ptr mom_dist; }; diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index acde1fa69..0c5c30c60 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -42,15 +42,21 @@ Real parseMassString(const std::string& name) { } } -ConstantDensityDistribution::ConstantDensityDistribution(Real density) +ConstantDensityProfile::ConstantDensityProfile(Real density) : _density(density) {} -Real ConstantDensityDistribution::getDensity(Real x, Real y, Real z) const +Real ConstantDensityProfile::getDensity(Real x, Real y, Real z) const { return _density; } +CustomDensityProfile::CustomDensityProfile(const std::string& species_name) +{ + ParmParse pp(species_name); + pp.getarr("custom_profile_params", params); +} + ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, Real uy, Real uz) @@ -120,18 +126,20 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.query("zmax", zmax); // parse density information - std::string rho_dist_s; - pp.get("profile", rho_dist_s); - std::transform(rho_dist_s.begin(), - rho_dist_s.end(), - rho_dist_s.begin(), + std::string rho_prof_s; + pp.get("profile", rho_prof_s); + std::transform(rho_prof_s.begin(), + rho_prof_s.end(), + rho_prof_s.begin(), ::tolower); - if (rho_dist_s == "constant") { + if (rho_prof_s == "constant") { Real density; pp.get("density", density); - rho_dist.reset(new ConstantDensityDistribution(density)); + rho_prof.reset(new ConstantDensityProfile(density)); + } else if (rho_prof_s == "custom") { + rho_prof.reset(new CustomDensityProfile(species_name)); } else { - StringParseAbortMessage("Density distribution type", rho_dist_s); + StringParseAbortMessage("Density profile type", rho_prof_s); } pp.get("num_particles_per_cell", num_particles_per_cell); @@ -166,10 +174,16 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) // get injection style pp.get("injection_style", injection_style); - if (injection_style != "nrandom" or - injection_style != "nrandompercell" or - injection_style != "ndiagpercell") { - StringParseAbortMessage("Injection style", injection_style); + std::transform(injection_style.begin(), + injection_style.end(), + injection_style.begin(), + ::tolower); + if (injection_style == "nrandomnormal" or + injection_style == "nrandomuniformpercell" or + injection_style == "ndiagpercell") { + return; + } else { + StringParseAbortMessage("Injection style", injection_style); } } @@ -188,5 +202,5 @@ bool PlasmaInjector::insideBounds(Real x, Real y, Real z) { } Real PlasmaInjector::getDensity(Real x, Real y, Real z) { - return rho_dist->getDensity(x, y, z); + return rho_prof->getDensity(x, y, z); } -- cgit v1.2.3 From 66a4755131d695803cb2b1f7f2b2fe446bc1271e Mon Sep 17 00:00:00 2001 From: atmyers Date: Thu, 16 Mar 2017 16:16:07 -0700 Subject: put these functions into a namespace. --- Source/PlasmaInjector.cpp | 64 ++++++++++++++++++++++++----------------------- 1 file changed, 33 insertions(+), 31 deletions(-) (limited to 'Source/PlasmaInjector.cpp') diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index 0c5c30c60..5592e25f5 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -7,38 +7,40 @@ using namespace amrex; -void StringParseAbortMessage(const std::string& var, - const std::string& name) { - std::stringstream stringstream; - std::string string; - stringstream << var << " string '" << name << "' not recognized."; - string = stringstream.str(); - amrex::Abort(string.c_str()); -} - -Real parseChargeName(const std::string& name) { - if (name == "q_e") { - return PhysConst::q_e; - } else { - StringParseAbortMessage("Charge", name); - return 0.0; +namespace { + void StringParseAbortMessage(const std::string& var, + const std::string& name) { + std::stringstream stringstream; + std::string string; + stringstream << var << " string '" << name << "' not recognized."; + string = stringstream.str(); + amrex::Abort(string.c_str()); } -} - -Real parseChargeString(const std::string& name) { - if(name.substr(0, 1) == "-") - return -1.0 * parseChargeName(name.substr(1, name.size() - 1)); - return parseChargeName(name); -} - -Real parseMassString(const std::string& name) { - if (name == "m_e") { - return PhysConst::m_e; - } else if (name == "m_p"){ - return PhysConst::m_p; - } else { - StringParseAbortMessage("Mass", name); - return 0.0; + + Real parseChargeName(const std::string& name) { + if (name == "q_e") { + return PhysConst::q_e; + } else { + StringParseAbortMessage("Charge", name); + return 0.0; + } + } + + Real parseChargeString(const std::string& name) { + if(name.substr(0, 1) == "-") + return -1.0 * parseChargeName(name.substr(1, name.size() - 1)); + return parseChargeName(name); + } + + Real parseMassString(const std::string& name) { + if (name == "m_e") { + return PhysConst::m_e; + } else if (name == "m_p"){ + return PhysConst::m_p; + } else { + StringParseAbortMessage("Mass", name); + return 0.0; + } } } -- cgit v1.2.3 From 7cbccbeb645642704d5f47a19bba804751dbf6f4 Mon Sep 17 00:00:00 2001 From: atmyers Date: Thu, 16 Mar 2017 16:47:08 -0700 Subject: use std::array to represet momentum. --- Source/PhysicalParticleContainer.cpp | 6 +++--- Source/PlasmaInjector.H | 13 +++++++++---- Source/PlasmaInjector.cpp | 6 +++--- 3 files changed, 15 insertions(+), 10 deletions(-) (limited to 'Source/PlasmaInjector.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 9964f0e11..21fd95ba7 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -89,7 +89,7 @@ PhysicalParticleContainer::InitNRandomUniformPerCell () { if (plasma_injector->insideBounds(x, y, z)) { Real weight; - Real u[3]; + std::array u; weight = plasma_injector->getDensity(x, y, z) * scale_fac; plasma_injector->getMomentum(u); attribs[PIdx::w ] = weight; @@ -156,7 +156,7 @@ PhysicalParticleContainer::InitNDiagPerCell () { if (plasma_injector->insideBounds(x, y, z)) { Real weight; - Real u[3]; + std::array u; weight = plasma_injector->getDensity(x, y, z) * scale_fac; plasma_injector->getMomentum(u); attribs[PIdx::w ] = weight; @@ -322,7 +322,7 @@ PhysicalParticleContainer::AddParticles (int lev, const Box& part_box) Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; #endif Real weight; - Real u[3]; + std::array u; weight = plasma_injector->getDensity(x, y, z) * scale_fac; plasma_injector->getMomentum(u); attribs[PIdx::w ] = weight; diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H index e2f1fe459..ffe59dd2b 100644 --- a/Source/PlasmaInjector.H +++ b/Source/PlasmaInjector.H @@ -2,6 +2,7 @@ #define PLASMA_INJECTOR_H_ #include +#include #include "AMReX_Real.H" #include "AMReX_Array.H" @@ -63,8 +64,9 @@ private: class PlasmaMomentumDistribution { public: + using vec3 = std::array; virtual ~PlasmaMomentumDistribution() {}; - virtual void getMomentum(amrex::Real* u) = 0; + virtual void getMomentum(vec3& u) = 0; }; /// @@ -77,7 +79,7 @@ public: ConstantMomentumDistribution(amrex::Real ux, amrex::Real uy, amrex::Real uz); - virtual void getMomentum(amrex::Real* u) override; + virtual void getMomentum(vec3& u) override; private: amrex::Real _ux; @@ -99,7 +101,7 @@ public: amrex::Real uy_m, amrex::Real uz_m, amrex::Real u_th); - virtual void getMomentum(amrex::Real* u) override; + virtual void getMomentum(vec3& u) override; private: amrex::Real _ux_m; amrex::Real _uy_m; @@ -119,6 +121,9 @@ class PlasmaInjector { public: + + using vec3 = std::array; + PlasmaInjector(int ispecies, const std::string& name); amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z); @@ -127,7 +132,7 @@ public: int numParticlesPerCell() {return num_particles_per_cell;}; - void getMomentum(amrex::Real* u); + void getMomentum(vec3& u); amrex::Real getCharge() {return charge;} amrex::Real getMass() {return mass;} diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp index 5592e25f5..027bb5b57 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/PlasmaInjector.cpp @@ -65,7 +65,7 @@ ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, : _ux(ux), _uy(uy), _uz(uz) {} -void ConstantMomentumDistribution::getMomentum(Real* u) { +void ConstantMomentumDistribution::getMomentum(vec3& u) { u[0] = _ux; u[1] = _uy; u[2] = _uz; @@ -80,7 +80,7 @@ GaussianRandomMomentumDistribution::GaussianRandomMomentumDistribution(Real ux_m { } -void GaussianRandomMomentumDistribution::getMomentum(Real* u) { +void GaussianRandomMomentumDistribution::getMomentum(vec3& u) { Real ux_th = momentum_distribution(generator); Real uy_th = momentum_distribution(generator); Real uz_th = momentum_distribution(generator); @@ -189,7 +189,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) } } -void PlasmaInjector::getMomentum(Real* u) { +void PlasmaInjector::getMomentum(vec3& u) { mom_dist->getMomentum(u); u[0] *= PhysConst::c; u[1] *= PhysConst::c; -- cgit v1.2.3