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 --- Source/ParticleContainer.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'Source/ParticleContainer.cpp') 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; } } -- 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/ParticleContainer.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