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/CustomDensityProb.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 Source/CustomDensityProb.cpp (limited to 'Source/CustomDensityProb.cpp') diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp new file mode 100644 index 000000000..54e5e0fad --- /dev/null +++ b/Source/CustomDensityProb.cpp @@ -0,0 +1,9 @@ +#include + +#include + +Real CustomPlasmaInjector::getDensity(Real x, Real y, Real z) { + static_assert(false, + "If running with a custom density profile, you must supply a CustomDensityProb.cpp file"); + return 0.0; +} -- cgit v1.2.3 From ee779da1c75fde34e4fa88e0a4a61623b411c1ee Mon Sep 17 00:00:00 2001 From: atmyers Date: Wed, 15 Mar 2017 12:33:09 -0700 Subject: fix default CustomPlasmaInjector --- Source/CustomDensityProb.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) (limited to 'Source/CustomDensityProb.cpp') diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp index 54e5e0fad..9b7298ee7 100644 --- a/Source/CustomDensityProb.cpp +++ b/Source/CustomDensityProb.cpp @@ -2,8 +2,9 @@ #include -Real CustomPlasmaInjector::getDensity(Real x, Real y, Real z) { - static_assert(false, - "If running with a custom density profile, you must supply a CustomDensityProb.cpp file"); - return 0.0; +#include + +amrex::Real CustomPlasmaInjector::getDensity(amrex::Real x, amrex::Real y, amrex::Real z) { + amrex::Abort("If running with a custom density profile, you must supply a CustomDensityProb.cpp file"); + return 0.0; } -- 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/CustomDensityProb.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/CustomDensityProb.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