From a44b32cad0d8df3decdb4f0ac1dfc66364236136 Mon Sep 17 00:00:00 2001 From: atmyers Date: Tue, 14 Mar 2017 15:28:53 -0700 Subject: Move the plasma injector stuff to particle container. --- Source/PhysicalParticleContainer.cpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 8c229d41f..767a41919 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -1,6 +1,8 @@ #include +#include #include +#include #include #include @@ -9,6 +11,28 @@ using namespace amrex; PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies) : WarpXParticleContainer(amr_core, ispecies) { + std::stringstream ss; + ss << "species" << ispecies; + ParmParse pp(ss.str()); + + std::string plasma_profile_s; + pp.get("profile", plasma_profile_s); + std::transform(plasma_profile_s.begin(), plasma_profile_s.end(), plasma_profile_s.begin(), ::tolower); + if (plasma_profile_s == "constant") { + plasma_injector.reset(new ConstantPlasmaInjector); + } + + else if (plasma_profile_s == "double_ramp") { + plasma_injector.reset(new DoubleRampPlasmaInjector); + } + + else if (plasma_profile_s == "custom") { + plasma_injector.reset(new CustomPlasmaInjector); + } + + else { + amrex::Abort("Unknown plasma injector type"); + } } void -- 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/PhysicalParticleContainer.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 8b7bbe4b1bff5d6176970cbb0405640f702706b2 Mon Sep 17 00:00:00 2001 From: atmyers Date: Wed, 15 Mar 2017 15:35:14 -0700 Subject: refactor particle initialization code, both for initial conditions and when the window moves. --- Exec/plasma_acceleration/ParticleProb.cpp | 58 +-------------- Exec/plasma_acceleration/inputs | 8 -- Source/PhysicalParticleContainer.H | 32 ++++++-- Source/PhysicalParticleContainer.cpp | 120 ++++++++++++++++++++++++++++++ Source/WarpX.H | 2 - Source/WarpX.cpp | 6 -- Source/WarpXEvolve.cpp | 15 +--- Source/WarpXParticleContainer.H | 3 - Source/WarpXParticleContainer.cpp | 43 ----------- 9 files changed, 149 insertions(+), 138 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Exec/plasma_acceleration/ParticleProb.cpp b/Exec/plasma_acceleration/ParticleProb.cpp index 010484315..806861a19 100644 --- a/Exec/plasma_acceleration/ParticleProb.cpp +++ b/Exec/plasma_acceleration/ParticleProb.cpp @@ -16,61 +16,5 @@ using namespace amrex; void PhysicalParticleContainer::InitData() { - BL_PROFILE("PhysicalParticleContainer::InitData()"); - - charge = plasma_injector->getCharge(); - mass = plasma_injector->getMass(); - - const int lev = 0; - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real scale_fac; - int n_part_per_cell = plasma_injector->numParticlesPerCell(); - -#if BL_SPACEDIM==3 - scale_fac = dx[0]*dx[1]*dx[2]/n_part_per_cell; -#elif BL_SPACEDIM==2 - scale_face = dx[0]*dx[1]/n_part_per_cell; -#endif - - std::array attribs; - attribs.fill(0.0); - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_partinsideBounds(x, y, z)) { - Real weight; - Real u[3]; - weight = plasma_injector->getDensity(x, y, z) * scale_fac; - plasma_injector->getMomentum(u); - attribs[PIdx::w ] = weight; - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } + InitNPerCell(); } diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 45832d536..241fc705a 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -38,27 +38,22 @@ particles.species_names = beam plasma # The beam species information beam.charge = -q_e beam.mass = m_e - beam.xmin = -20.e-6 beam.xmax = 20.e-6 beam.ymin = -20.e-6 beam.ymax = 20.e-6 beam.zmin = -150.e-6 beam.zmax = -100.e-6 - beam.profile = "constant" beam.density = 1e23 # number of particles per m^3 beam.num_particles_per_cell = 4 - beam.ux = 0.0 beam.uy = 0.0 beam.uz = 1e9 - # The plasma species information plasma.charge = -q_e plasma.mass = m_e - plasma.profile = "constant" plasma.xmin = -200.e-6 plasma.xmax = 200.e-6 @@ -66,7 +61,6 @@ plasma.ymin = -200.e-6 plasma.ymax = 200.e-6 plasma.zmin = 0.e-6 plasma.zmax = 200.e-6 - plasma.profile = "constant" plasma.density = 1e22 # number of particles per m^3 plasma.num_particles_per_cell = 4 @@ -84,5 +78,3 @@ warpx.moving_window_v = 1.0 # in units of the speed of light warpx.do_plasma_injection = 1 warpx.num_injected_species = 1 warpx.injected_plasma_species = 1 0 -warpx.injected_plasma_density = 1e22 1e22 -warpx.injected_plasma_ppc = 4 4 diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 923d40429..a399df3fe 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -10,7 +10,9 @@ class PhysicalParticleContainer : public WarpXParticleContainer { public: - PhysicalParticleContainer (amrex::AmrCore* amr_core, int ispecies, const std::string& name); + PhysicalParticleContainer (amrex::AmrCore* amr_core, + int ispecies, + const std::string& name); virtual ~PhysicalParticleContainer () {} virtual void AllocData () override; @@ -18,20 +20,36 @@ public: virtual void InitData () override; virtual void FieldGather(int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) override; virtual void Evolve (int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, - amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::Real t, amrex::Real dt) override; + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz, + amrex::MultiFab& jx, + amrex::MultiFab& jy, + amrex::MultiFab& jz, + amrex::Real t, + amrex::Real dt) override; virtual void PostRestart () override {} + // Inject 'ppc' particles per cell in Box 'part_box' + void AddParticles (int lev, const amrex::Box& part_box); + protected: - std::string species_name; + void InitNPerCell (); + std::string species_name; std::unique_ptr plasma_injector; }; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index c7e98def9..0e3875ca1 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -45,6 +45,68 @@ PhysicalParticleContainer::AllocData () resizeData(); } +void +PhysicalParticleContainer::InitNPerCell () { + BL_PROFILE("PhysicalParticleContainer::InitNPerCell()"); + + charge = plasma_injector->getCharge(); + mass = plasma_injector->getMass(); + + const int lev = 0; + const Geometry& geom = Geom(lev); + const Real* dx = geom.CellSize(); + + Real scale_fac; + int n_part_per_cell = plasma_injector->numParticlesPerCell(); + +#if BL_SPACEDIM==3 + scale_fac = dx[0]*dx[1]*dx[2]/n_part_per_cell; +#elif BL_SPACEDIM==2 + scale_fac = dx[0]*dx[1]/n_part_per_cell; +#endif + + std::array attribs; + attribs.fill(0.0); + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { + const Box& tile_box = mfi.tilebox(); + RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + const auto& boxlo = tile_box.smallEnd(); + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) + { + for (int i_part=0; i_partinsideBounds(x, y, z)) { + Real weight; + Real u[3]; + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + plasma_injector->getMomentum(u); + attribs[PIdx::w ] = weight; + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + } + } + } + } +} + + void PhysicalParticleContainer::FieldGather (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, @@ -152,6 +214,64 @@ PhysicalParticleContainer::FieldGather (int lev, } } +void +PhysicalParticleContainer::AddParticles (int lev, const Box& part_box) +{ + const Geometry& geom = Geom(lev); + const Real* dx = geom.CellSize(); + + Real scale_fac; + int ppc = plasma_injector->numParticlesPerCell(); + +#if BL_SPACEDIM==3 + scale_fac = dx[0]*dx[1]*dx[2]/ppc; +#elif BL_SPACEDIM==2 + scale_fac = dx[0]*dx[1]/ppc; +#endif + + std::array attribs; + attribs.fill(0.0); + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + const Box& intersectBox = tile_box & part_box; + if (intersectBox.ok()) + { + RealBox real_box { intersectBox, dx, geom.ProbLo() }; + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + const IntVect& boxlo = tile_box.smallEnd(); + for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) + { + for (int i_part=0; i_part < ppc; i_part++) + { + Real particle_shift = (0.5+i_part)/ppc; +#if (BL_SPACEDIM == 3) + Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; + Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; + Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; +#elif (BL_SPACEDIM == 2) + Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; + Real y = 0.0; + Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; +#endif + Real weight; + Real u[3]; + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + plasma_injector->getMomentum(u); + attribs[PIdx::w ] = weight; + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + } + } + } + } +} + void PhysicalParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, diff --git a/Source/WarpX.H b/Source/WarpX.H index c198990dc..135a0f690 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -179,9 +179,7 @@ private: // Plasma injection parameters int do_plasma_injection = 0; int num_injected_species = -1; - amrex::Array injected_plasma_ppc; amrex::Array injected_plasma_species; - amrex::Array injected_plasma_density; // Other runtime parameters int verbose = 1; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 7ac23307b..bcdbc709b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -165,15 +165,9 @@ WarpX::ReadParameters () pp.query("do_plasma_injection", do_plasma_injection); if (do_plasma_injection) { pp.get("num_injected_species", num_injected_species); - injected_plasma_ppc.resize(num_injected_species); - pp.getarr("injected_plasma_ppc", injected_plasma_ppc, - 0, num_injected_species); injected_plasma_species.resize(num_injected_species); pp.getarr("injected_plasma_species", injected_plasma_species, 0, num_injected_species); - injected_plasma_density.resize(num_injected_species); - pp.getarr("injected_plasma_density", injected_plasma_density, - 0, num_injected_species); } pp.query("use_laser", use_laser); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 5e5154554..875a0392d 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -262,20 +262,11 @@ WarpX::InjectPlasma (int num_shift, int dir) const Real* dx = geom[lev].CellSize(); - for (int i = 0; i < num_injected_species; ++i) - { - int ppc = injected_plasma_ppc[i]; - Real density = injected_plasma_density[i]; -#if BL_SPACEDIM==3 - Real weight = density * dx[0]*dx[1]*dx[2]/ppc; -#elif BL_SPACEDIM==2 - Real weight = density * dx[0]*dx[1]/ppc; -#endif - + for (int i = 0; i < num_injected_species; ++i) { int ispecies = injected_plasma_species[i]; WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - - pc.AddParticles(lev, particleBox, weight, ppc); + PhysicalParticleContainer* ppc = (PhysicalParticleContainer*) &pc; + ppc->AddParticles(lev, particleBox); } } } diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index b2f308ac5..6fc8b1415 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -67,9 +67,6 @@ public: amrex::Real x, amrex::Real y, amrex::Real z, const std::array& attribs); - // Inject 'ppc' particles per cell in Box 'part_box' - void AddParticles (int lev, const amrex::Box& part_box, amrex::Real weight, int ppc); - void ReadHeader (std::istream& is); void WriteHeader (std::ostream& os) const; diff --git a/Source/WarpXParticleContainer.cpp b/Source/WarpXParticleContainer.cpp index 76860ee1e..25b92a712 100644 --- a/Source/WarpXParticleContainer.cpp +++ b/Source/WarpXParticleContainer.cpp @@ -65,49 +65,6 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, particle_tile.push_back(attribs); } -void -WarpXParticleContainer::AddParticles (int lev, const Box& part_box, Real weight, int ppc) -{ - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - std::array attribs; - attribs.fill(0.0); - attribs[PIdx::w] = weight; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - const Box& intersectBox = tile_box & part_box; - if (intersectBox.ok()) - { - RealBox real_box { intersectBox, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part < ppc; i_part++) - { - Real particle_shift = (0.5+i_part)/ppc; -#if (BL_SPACEDIM == 3) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} - void WarpXParticleContainer::AddNParticles (int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, -- cgit v1.2.3 From 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/PhysicalParticleContainer.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/PhysicalParticleContainer.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 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/PhysicalParticleContainer.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 From 1439e2a8138ffdd81a5901aa2713f9c29dfab68d Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 21 Mar 2017 13:06:30 -0700 Subject: Use the plasma injector when the window moves as well as on initialization. --- Source/PhysicalParticleContainer.H | 10 ++-- Source/PhysicalParticleContainer.cpp | 108 +++++++++++------------------------ 2 files changed, 37 insertions(+), 81 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 98d15e22c..40a8a6068 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -42,14 +42,14 @@ public: virtual void PostRestart () override {} - // Inject 'ppc' particles per cell in Box 'part_box' - void AddParticles (int lev, const amrex::Box& part_box); + // Inject particles in Box 'part_box' + void AddParticles (int lev, amrex::Box part_box = amrex::Box()); protected: - void InitNDiagPerCell (); - void InitNRandomUniformPerCell (); - void InitNRandomNormal (); + void AddNDiagPerCell (int lev, amrex::Box part_box = amrex::Box()); + void AddNRandomUniformPerCell (int lev, amrex::Box part_box = amrex::Box()); + void AddNRandomNormal (int lev, amrex::Box part_box = amrex::Box()); std::string species_name; std::unique_ptr plasma_injector; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 21fd95ba7..bc42b36a7 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -16,12 +16,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } 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(); + AddParticles(0); // Note - only one level right now } void @@ -34,15 +29,16 @@ PhysicalParticleContainer::AllocData () } void -PhysicalParticleContainer::InitNRandomUniformPerCell () { - BL_PROFILE("PhysicalParticleContainer::InitNRandomPerCell()"); +PhysicalParticleContainer::AddNRandomUniformPerCell (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNRandomPerCell()"); charge = plasma_injector->getCharge(); mass = plasma_injector->getMass(); - const int lev = 0; const Geometry& geom = Geom(lev); const Real* dx = geom.CellSize(); + + if (!part_box.ok()) part_box = geom.Domain(); Real scale_fac; int n_part_per_cell = plasma_injector->numParticlesPerCell(); @@ -62,13 +58,16 @@ PhysicalParticleContainer::InitNRandomUniformPerCell () { for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; + const Box& intersectBox = tile_box & part_box; + if (!intersectBox.ok()) continue; + + RealBox tile_real_box { intersectBox, 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)) + const auto& boxlo = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) { for (int i_part=0; i_partgetCharge(); mass = plasma_injector->getMass(); - const int lev = 0; const Geometry& geom = Geom(lev); const Real* dx = geom.CellSize(); + if (!part_box.ok()) part_box = geom.Domain(); + Real scale_fac; int n_part_per_cell = plasma_injector->numParticlesPerCell(); @@ -133,13 +133,16 @@ PhysicalParticleContainer::InitNDiagPerCell () { attribs.fill(0.0); for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; + const Box& intersectBox = tile_box & part_box; + if (!intersectBox.ok()) continue; + + RealBox tile_real_box { intersectBox, 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)) + const auto& boxlo = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) { for (int i_part=0; i_partinsideBounds(x, y, z)) { Real weight; std::array u; @@ -279,61 +282,14 @@ PhysicalParticleContainer::FieldGather (int lev, } void -PhysicalParticleContainer::AddParticles (int lev, const Box& part_box) +PhysicalParticleContainer::AddParticles (int lev, Box part_box) { - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real scale_fac; - int ppc = plasma_injector->numParticlesPerCell(); - -#if BL_SPACEDIM==3 - scale_fac = dx[0]*dx[1]*dx[2]/ppc; -#elif BL_SPACEDIM==2 - scale_fac = dx[0]*dx[1]/ppc; -#endif - - std::array attribs; - attribs.fill(0.0); - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - const Box& intersectBox = tile_box & part_box; - if (intersectBox.ok()) - { - RealBox real_box { intersectBox, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part < ppc; i_part++) - { - Real particle_shift = (0.5+i_part)/ppc; -#if (BL_SPACEDIM == 3) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - Real weight; - std::array u; - 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); - } - } - } - } + if (plasma_injector->injection_style == "ndiagpercell") + AddNDiagPerCell(lev, part_box); + else if (plasma_injector->injection_style == "nrandomuniformpercell") + AddNRandomUniformPerCell(lev, part_box); + else if (plasma_injector->injection_style == "nrandomnormal") + AddNRandomNormal(lev, part_box); } void -- cgit v1.2.3