aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/CustomDensityProb.cpp13
-rw-r--r--Source/Make.package5
-rw-r--r--Source/ParticleContainer.H5
-rw-r--r--Source/ParticleContainer.cpp7
-rw-r--r--Source/ParticleProb.cpp21
-rw-r--r--Source/PhysicalParticleContainer.H39
-rw-r--r--Source/PhysicalParticleContainer.cpp169
-rw-r--r--Source/PlasmaInjector.H161
-rw-r--r--Source/PlasmaInjector.cpp208
-rw-r--r--Source/WarpX.H15
-rw-r--r--Source/WarpX.cpp160
-rw-r--r--Source/WarpXEvolve.cpp15
-rw-r--r--Source/WarpXParticleContainer.H3
-rw-r--r--Source/WarpXParticleContainer.cpp43
-rw-r--r--Source/WarpXProb.cpp11
-rw-r--r--Source/WarpXRegrid.cpp2
-rw-r--r--Source/WarpX_picsar.F9072
17 files changed, 772 insertions, 177 deletions
diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp
new file mode 100644
index 000000000..1df3d75ad
--- /dev/null
+++ b/Source/CustomDensityProb.cpp
@@ -0,0 +1,13 @@
+#include <PlasmaInjector.H>
+
+#include <iostream>
+
+#include <AMReX.H>
+
+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 e1ea6aa29..64265957b 100644
--- a/Source/Make.package
+++ b/Source/Make.package
@@ -6,7 +6,7 @@ CEXE_sources += WarpXWrappers.cpp WarpX_py.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 WarpX_py.H
@@ -15,6 +15,9 @@ CEXE_headers += WarpX.H WarpX_f.H WarpXConst.H
CEXE_headers += ParticleContainer.H WarpXParticleContainer.H PhysicalParticleContainer.H LaserParticleContainer.H
+CEXE_headers += PlasmaInjector.H
+CEXE_sources += PlasmaInjector.cpp CustomDensityProb.cpp
+
F90EXE_sources += WarpX_f.F90 WarpX_picsar.F90 WarpX_laser.F90
ifeq ($(USE_OPENBC_POISSON),TRUE)
diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H
index 8d5c4b7e0..20946e6e5 100644
--- a/Source/ParticleContainer.H
+++ b/Source/ParticleContainer.H
@@ -73,6 +73,10 @@ public:
int nSpecies() {return nspecies;}
+protected:
+
+ std::vector<std::string> species_names;
+
private:
// physical particles (+ laser)
@@ -82,5 +86,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 b27534feb..07cc8de09 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,11 @@ MultiParticleContainer::ReadParameters ()
pp.query("nspecies", nspecies);
BL_ASSERT(nspecies >= 0);
+ if (nspecies > 0) {
+ pp.getarr("species_names", species_names);
+ BL_ASSERT(species_names.size() == nspecies);
+ }
+
initialized = true;
}
}
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 <cmath>
-
-#include <AMReX_BLProfiler.H>
-
-#include <ParticleContainer.H>
-#include <WarpXConst.H>
-
-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 c6d13185e..40a8a6068 100644
--- a/Source/PhysicalParticleContainer.H
+++ b/Source/PhysicalParticleContainer.H
@@ -3,13 +3,16 @@
#include <map>
+#include <PlasmaInjector.H>
#include <WarpXParticleContainer.H>
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;
@@ -17,15 +20,39 @@ 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 particles in Box 'part_box'
+ void AddParticles (int lev, amrex::Box part_box = amrex::Box());
+
+protected:
+
+ 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<PlasmaInjector> plasma_injector;
};
#endif
diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp
index 9e2532403..d499f7dfb 100644
--- a/Source/PhysicalParticleContainer.cpp
+++ b/Source/PhysicalParticleContainer.cpp
@@ -1,4 +1,5 @@
#include <limits>
+#include <sstream>
#include <ParticleContainer.H>
#include <WarpX_f.H>
@@ -6,9 +7,16 @@
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)
{
+ plasma_injector.reset(new PlasmaInjector(species_id, species_name));
+}
+
+void PhysicalParticleContainer::InitData() {
+ AddParticles(0); // Note - only one level right now
}
void
@@ -21,6 +29,152 @@ PhysicalParticleContainer::AllocData ()
}
void
+PhysicalParticleContainer::AddNRandomUniformPerCell (int lev, Box part_box) {
+ BL_PROFILE("PhysicalParticleContainer::AddNRandomPerCell()");
+
+ charge = plasma_injector->getCharge();
+ mass = plasma_injector->getMass();
+
+ 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();
+
+#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<Real,PIdx::nattribs> attribs;
+ attribs.fill(0.0);
+
+ // Initialize random generator for normal distribution
+ std::default_random_engine generator;
+ std::uniform_real_distribution<Real> position_distribution(0.0,1.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()) continue;
+
+ RealBox tile_real_box { intersectBox, dx, geom.ProbLo() };
+
+ const int grid_id = mfi.index();
+ const int tile_id = mfi.LocalTileIndex();
+
+ const auto& boxlo = intersectBox.smallEnd();
+ for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv))
+ {
+ for (int i_part=0; i_part<n_part_per_cell;i_part++)
+ {
+ // Randomly generate the positions (uniformly inside each cell)
+ Real particle_shift_x = position_distribution(generator);
+ Real particle_shift_y = position_distribution(generator);
+ Real particle_shift_z = position_distribution(generator);
+
+#if (BL_SPACEDIM == 3)
+ Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0];
+ Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_y)*dx[1];
+ Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift_z)*dx[2];
+#elif (BL_SPACEDIM == 2)
+ Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0];
+ Real y = 0.0;
+ Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_z)*dx[1];
+#endif
+
+ if (plasma_injector->insideBounds(x, y, z)) {
+ Real weight;
+ std::array<Real, 3> 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);
+ }
+ }
+ }
+ }
+}
+
+void
+PhysicalParticleContainer::AddNRandomNormal (int lev, Box part_box) {
+ BL_PROFILE("PhysicalParticleContainer::AddNRandomNormal()");
+ amrex::Abort("PhysicalParticleContainer::AddNRandomNormal() not implemented yet.");
+}
+
+void
+PhysicalParticleContainer::AddNDiagPerCell (int lev, Box part_box) {
+ BL_PROFILE("PhysicalParticleContainer::AddNDiagPerCell()");
+
+ charge = plasma_injector->getCharge();
+ mass = plasma_injector->getMass();
+
+ 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();
+
+#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<Real,PIdx::nattribs> 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()) continue;
+
+ RealBox tile_real_box { intersectBox, dx, geom.ProbLo() };
+
+ const int grid_id = mfi.index();
+ const int tile_id = mfi.LocalTileIndex();
+
+ const auto& boxlo = intersectBox.smallEnd();
+ for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv))
+ {
+ for (int i_part=0; i_part<n_part_per_cell;i_part++)
+ {
+ Real particle_shift = (0.5+i_part)/n_part_per_cell;
+#if (BL_SPACEDIM == 3)
+ Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0];
+ Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1];
+ Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2];
+#elif (BL_SPACEDIM == 2)
+ Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0];
+ Real y = 0.0;
+ Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1];
+#endif
+
+ if (plasma_injector->insideBounds(x, y, z)) {
+ Real weight;
+ std::array<Real, 3> 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);
+ }
+ }
+ }
+ }
+}
+
+
+void
PhysicalParticleContainer::FieldGather (int lev,
const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
@@ -128,6 +282,17 @@ PhysicalParticleContainer::FieldGather (int lev,
}
void
+PhysicalParticleContainer::AddParticles (int lev, Box part_box)
+{
+ 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
PhysicalParticleContainer::Evolve (int lev,
const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz,
diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H
new file mode 100644
index 000000000..c346f79da
--- /dev/null
+++ b/Source/PlasmaInjector.H
@@ -0,0 +1,161 @@
+#ifndef PLASMA_INJECTOR_H_
+#define PLASMA_INJECTOR_H_
+
+#include <random>
+#include <array>
+
+#include "AMReX_REAL.H"
+#include "AMReX_Array.H"
+#include "AMReX_ParmParse.H"
+
+///
+/// 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 PlasmaDensityProfile
+{
+public:
+ virtual ~PlasmaDensityProfile() {};
+ 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 ConstantDensityProfile : public PlasmaDensityProfile
+{
+public:
+ ConstantDensityProfile(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
+/// in their problem directory.
+///
+///
+class CustomDensityProfile : public PlasmaDensityProfile
+{
+public:
+ CustomDensityProfile(const std::string& species_name);
+ virtual amrex::Real getDensity(amrex::Real x,
+ amrex::Real y,
+ amrex::Real z) const override;
+private:
+ amrex::Array<amrex::Real> params;
+};
+
+///
+/// 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:
+ using vec3 = std::array<amrex::Real, 3>;
+ virtual ~PlasmaMomentumDistribution() {};
+ virtual void getMomentum(vec3& 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);
+ virtual void getMomentum(vec3& u) override;
+
+private:
+ amrex::Real _ux;
+ amrex::Real _uy;
+ 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:
+ GaussianRandomMomentumDistribution(amrex::Real ux_m,
+ amrex::Real uy_m,
+ amrex::Real uz_m,
+ amrex::Real u_th);
+ virtual void getMomentum(vec3& u) override;
+private:
+ amrex::Real _ux_m;
+ amrex::Real _uy_m;
+ amrex::Real _uz_m;
+ amrex::Real _u_th;
+
+ // xxxxx These random number generators are not thread safe.
+ std::normal_distribution<amrex::Real> momentum_distribution;
+ 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
+{
+
+public:
+
+ using vec3 = std::array<amrex::Real, 3>;
+
+ PlasmaInjector(int ispecies, const std::string& name);
+
+ amrex::Real getDensity(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;};
+
+ void getMomentum(vec3& u);
+
+ amrex::Real getCharge() {return charge;}
+ amrex::Real getMass() {return mass;}
+
+ std::string injection_style;
+
+protected:
+
+ amrex::Real mass, charge;
+
+ amrex::Real xmin, xmax;
+ amrex::Real ymin, ymax;
+ amrex::Real zmin, zmax;
+
+ amrex::Real density;
+ int num_particles_per_cell;
+
+ int species_id;
+ std::string species_name;
+
+ std::unique_ptr<PlasmaDensityProfile> rho_prof;
+ std::unique_ptr<PlasmaMomentumDistribution> mom_dist;
+};
+
+#endif
diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp
new file mode 100644
index 000000000..027bb5b57
--- /dev/null
+++ b/Source/PlasmaInjector.cpp
@@ -0,0 +1,208 @@
+#include "PlasmaInjector.H"
+
+#include <sstream>
+
+#include <WarpXConst.H>
+#include <AMReX.H>
+
+using namespace amrex;
+
+namespace {
+ void StringParseAbortMessage(const std::string& var,
+ const std::string& name) {
+ std::stringstream stringstream;
+ std::string string;
+ stringstream << var << " string '" << name << "' not recognized.";
+ string = stringstream.str();
+ amrex::Abort(string.c_str());
+ }
+
+ Real parseChargeName(const std::string& name) {
+ if (name == "q_e") {
+ return PhysConst::q_e;
+ } else {
+ StringParseAbortMessage("Charge", name);
+ return 0.0;
+ }
+ }
+
+ Real parseChargeString(const std::string& name) {
+ if(name.substr(0, 1) == "-")
+ return -1.0 * parseChargeName(name.substr(1, name.size() - 1));
+ return parseChargeName(name);
+ }
+
+ Real parseMassString(const std::string& name) {
+ if (name == "m_e") {
+ return PhysConst::m_e;
+ } else if (name == "m_p"){
+ return PhysConst::m_p;
+ } else {
+ StringParseAbortMessage("Mass", name);
+ return 0.0;
+ }
+ }
+}
+
+ConstantDensityProfile::ConstantDensityProfile(Real density)
+ : _density(density)
+{}
+
+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)
+ : _ux(ux), _uy(uy), _uz(uz)
+{}
+
+void ConstantMomentumDistribution::getMomentum(vec3& u) {
+ u[0] = _ux;
+ u[1] = _uy;
+ u[2] = _uz;
+}
+
+GaussianRandomMomentumDistribution::GaussianRandomMomentumDistribution(Real ux_m,
+ Real uy_m,
+ Real uz_m,
+ Real u_th)
+ : _ux_m(ux_m), _uy_m(uy_m), _uz_m(uz_m), _u_th(u_th),
+ momentum_distribution(0.0, u_th)
+{
+}
+
+void GaussianRandomMomentumDistribution::getMomentum(vec3& u) {
+ Real ux_th = momentum_distribution(generator);
+ Real uy_th = momentum_distribution(generator);
+ Real uz_th = momentum_distribution(generator);
+ u[0] = _ux_m + ux_th;
+ u[1] = _uy_m + uy_th;
+ u[2] = _uz_m + uz_th;
+}
+
+PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name)
+ : species_id(ispecies), species_name(name)
+{
+ ParmParse pp(species_name);
+
+ // parse charge and mass
+ std::string charge_s;
+ pp.get("charge", charge_s);
+ std::transform(charge_s.begin(),
+ charge_s.end(),
+ charge_s.begin(),
+ ::tolower);
+ charge = parseChargeString(charge_s);
+
+ std::string mass_s;
+ pp.get("mass", mass_s);
+ std::transform(mass_s.begin(),
+ mass_s.end(),
+ mass_s.begin(),
+ ::tolower);
+ mass = parseMassString(mass_s);
+
+ // parse plasma boundaries
+ xmin = std::numeric_limits<amrex::Real>::lowest();
+ ymin = std::numeric_limits<amrex::Real>::lowest();
+ zmin = std::numeric_limits<amrex::Real>::lowest();
+
+ xmax = std::numeric_limits<amrex::Real>::max();
+ ymax = std::numeric_limits<amrex::Real>::max();
+ zmax = std::numeric_limits<amrex::Real>::max();
+
+ pp.query("xmin", xmin);
+ pp.query("ymin", ymin);
+ pp.query("zmin", zmin);
+ pp.query("xmax", xmax);
+ pp.query("ymax", ymax);
+ pp.query("zmax", zmax);
+
+ // parse density information
+ 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_prof_s == "constant") {
+ Real density;
+ pp.get("density", density);
+ rho_prof.reset(new ConstantDensityProfile(density));
+ } else if (rho_prof_s == "custom") {
+ rho_prof.reset(new CustomDensityProfile(species_name));
+ } else {
+ StringParseAbortMessage("Density profile type", rho_prof_s);
+ }
+ pp.get("num_particles_per_cell", num_particles_per_cell);
+
+ // parse momentum information
+ std::string mom_dist_s;
+ pp.get("momentum_distribution_type", mom_dist_s);
+ std::transform(mom_dist_s.begin(),
+ mom_dist_s.end(),
+ mom_dist_s.begin(),
+ ::tolower);
+ if (mom_dist_s == "constant") {
+ Real ux = 0.;
+ Real uy = 0.;
+ Real uz = 0.;
+ pp.query("ux", ux);
+ pp.query("uy", uy);
+ pp.query("uz", uz);
+ mom_dist.reset(new ConstantMomentumDistribution(ux, uy, uz));
+ } else if (mom_dist_s == "gaussian") {
+ Real ux_m = 0.;
+ Real uy_m = 0.;
+ Real uz_m = 0.;
+ Real u_th = 0.;
+ pp.query("ux_m", ux_m);
+ pp.query("uy_m", uy_m);
+ pp.query("uz_m", uz_m);
+ pp.query("u_th", u_th);
+ mom_dist.reset(new GaussianRandomMomentumDistribution(ux_m, uy_m, uz_m, u_th));
+ } else {
+ StringParseAbortMessage("Momentum distribution type", mom_dist_s);
+ }
+
+ // get injection style
+ pp.get("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);
+ }
+}
+
+void PlasmaInjector::getMomentum(vec3& u) {
+ mom_dist->getMomentum(u);
+ u[0] *= PhysConst::c;
+ u[1] *= PhysConst::c;
+ u[2] *= PhysConst::c;
+}
+
+bool PlasmaInjector::insideBounds(Real x, Real y, Real z) {
+ if (x >= xmax || x < xmin ||
+ y >= ymax || y < ymin ||
+ z >= zmax || z < zmin ) return false;
+ return true;
+}
+
+Real PlasmaInjector::getDensity(Real x, Real y, Real z) {
+ return rho_prof->getDensity(x, y, z);
+}
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 0067361a1..32128c90e 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -36,6 +36,15 @@ public:
MultiParticleContainer& GetPartContainer () { return *mypc; }
+#if (BL_SPACEDIM == 3)
+ void fillSlice(amrex::Real z_coord) const;
+
+ void sampleAtPoints(const amrex::Array<amrex::Real>& x,
+ const amrex::Array<amrex::Real>& y,
+ const amrex::Array<amrex::Real>& z,
+ amrex::Array<amrex::Array<amrex::Real> >& result) const;
+#endif
+
static void FillBoundary (amrex::MultiFab& mf, const amrex::Geometry& geom, const amrex::IntVect& nodalflag);
static void shiftMF(amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift,
@@ -113,6 +122,8 @@ protected:
//! Delete level data. Called by AmrCore::regrid.
virtual void ClearLevel (int lev) override;
+ amrex::Box getIndexBox(const amrex::RealBox& real_box) const;
+
private:
// Singleton is used when the code is run from python
@@ -142,6 +153,7 @@ private:
void WriteWarpXHeader(const std::string& name) const;
void WriteJobInfo (const std::string& dir) const;
+// static void WriteBuildInfo (std::ostream& os);
static void PackPlotDataPtrs (amrex::Array<const amrex::MultiFab*>& pmf,
const amrex::Array<std::unique_ptr<amrex::MultiFab> >& data);
static void Copy (amrex::MultiFab& dstmf, int dcomp, int ncomp, const amrex::MultiFab& srcmf, int scomp);
@@ -157,6 +169,7 @@ private:
// Particle container
std::unique_ptr<MultiParticleContainer> mypc;
+
// Fields: First array for level, second for direction
amrex::Array<amrex::Array< std::unique_ptr<amrex::MultiFab> > > current;
amrex::Array<amrex::Array< std::unique_ptr<amrex::MultiFab> > > Efield;
@@ -171,9 +184,7 @@ private:
// Plasma injection parameters
int do_plasma_injection = 0;
int num_injected_species = -1;
- amrex::Array<int> injected_plasma_ppc;
amrex::Array<int> injected_plasma_species;
- amrex::Array<amrex::Real> injected_plasma_density;
// Other runtime parameters
int verbose = 1;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 052b5ed09..13c46ec39 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);
@@ -312,3 +306,157 @@ WarpX::Copy(MultiFab& dstmf, int dcomp, int ncomp, const MultiFab& srcmf, int sc
dfab.SetBoxType(dst_typ);
}
}
+
+Box
+WarpX::getIndexBox(const RealBox& real_box) const {
+ BL_ASSERT(max_level == 0);
+
+ IntVect slice_lo, slice_hi;
+
+ D_TERM(slice_lo[0]=floor((real_box.lo(0) - geom[0].ProbLo(0))/geom[0].CellSize(0));,
+ slice_lo[1]=floor((real_box.lo(1) - geom[0].ProbLo(1))/geom[0].CellSize(1));,
+ slice_lo[2]=floor((real_box.lo(2) - geom[0].ProbLo(2))/geom[0].CellSize(2)););
+
+ D_TERM(slice_hi[0]=floor((real_box.hi(0) - geom[0].ProbLo(0))/geom[0].CellSize(0));,
+ slice_hi[1]=floor((real_box.hi(1) - geom[0].ProbLo(1))/geom[0].CellSize(1));,
+ slice_hi[2]=floor((real_box.hi(2) - geom[0].ProbLo(2))/geom[0].CellSize(2)););
+
+ return Box(slice_lo, slice_hi) & geom[0].Domain();
+}
+
+#if (BL_SPACEDIM == 3)
+void
+WarpX::fillSlice(Real z_coord) const {
+ BL_ASSERT(max_level == 0);
+
+ // Get our slice and convert to index space
+ RealBox real_slice = geom[0].ProbDomain();
+ real_slice.setLo(2, z_coord);
+ real_slice.setHi(2, z_coord);
+ Box slice_box = getIndexBox(real_slice);
+
+ // define the multifab that stores slice
+ BoxArray ba = Efield[0][0]->boxArray();
+ const IndexType correct_typ(IntVect::TheZeroVector());
+ ba.convert(correct_typ);
+ const DistributionMapping& dm = dmap[0];
+ std::vector< std::pair<int, Box> > isects;
+ ba.intersections(slice_box, isects, false, 0);
+ Array<Box> boxes;
+ Array<int> procs;
+ for (int i = 0; i < isects.size(); ++i) {
+ procs.push_back(dm[isects[i].first]);
+ boxes.push_back(isects[i].second);
+ }
+ procs.push_back(ParallelDescriptor::MyProc());
+ BoxArray slice_ba(&boxes[0], boxes.size());
+ DistributionMapping slice_dmap(procs);
+ MultiFab slice(slice_ba, slice_dmap, 6, 0);
+
+ const MultiFab* mfs[6];
+ mfs[0] = Efield[0][0].get();
+ mfs[1] = Efield[0][1].get();
+ mfs[2] = Efield[0][2].get();
+ mfs[3] = Bfield[0][0].get();
+ mfs[4] = Bfield[0][1].get();
+ mfs[5] = Bfield[0][2].get();
+
+ IntVect flags[6];
+ flags[0] = WarpX::Ex_nodal_flag;
+ flags[1] = WarpX::Ey_nodal_flag;
+ flags[2] = WarpX::Ez_nodal_flag;
+ flags[3] = WarpX::Bx_nodal_flag;
+ flags[4] = WarpX::By_nodal_flag;
+ flags[5] = WarpX::Bz_nodal_flag;
+
+ // Fill the slice with sampled data
+ const Real* dx = geom[0].CellSize();
+ const Geometry& g = geom[0];
+ for (MFIter mfi(slice); mfi.isValid(); ++mfi) {
+ int slice_gid = mfi.index();
+ Box grid = slice_ba[slice_gid];
+ int xlo = grid.smallEnd(0), ylo = grid.smallEnd(1), zlo = grid.smallEnd(2);
+ int xhi = grid.bigEnd(0), yhi = grid.bigEnd(1), zhi = grid.bigEnd(2);
+
+ IntVect iv = grid.smallEnd();
+ ba.intersections(Box(iv, iv), isects, true, 0);
+ BL_ASSERT(!isects.empty());
+ int full_gid = isects[0].first;
+
+ for (int k = zlo; k <= zhi; k++) {
+ for (int j = ylo; j <= yhi; j++) {
+ for (int i = xlo; i <= xhi; i++) {
+ for (int comp = 0; comp < 6; comp++) {
+ Real x = g.ProbLo(0) + i*dx[0];
+ Real y = g.ProbLo(1) + j*dx[1];
+ Real z = z_coord;
+
+ D_TERM(iv[0]=floor((x - g.ProbLo(0) + 0.5*g.CellSize(0)*flags[comp][0])/g.CellSize(0));,
+ iv[1]=floor((y - g.ProbLo(1) + 0.5*g.CellSize(1)*flags[comp][1])/g.CellSize(1));,
+ iv[2]=floor((z - g.ProbLo(2) + 0.5*g.CellSize(2)*flags[comp][2])/g.CellSize(2)););
+
+ slice[slice_gid](IntVect(i, j, k), comp) = (*mfs[comp])[full_gid](iv);
+ }
+ }
+ }
+ }
+ }
+}
+
+void WarpX::sampleAtPoints(const Array<Real>& x,
+ const Array<Real>& y,
+ const Array<Real>& z,
+ Array<Array<Real> >& result) const {
+ BL_ASSERT((x.size() == y.size()) and (y.size() == z.size()));
+ BL_ASSERT(max_level == 0);
+
+ const MultiFab* mfs[6];
+ mfs[0] = Efield[0][0].get();
+ mfs[1] = Efield[0][1].get();
+ mfs[2] = Efield[0][2].get();
+ mfs[3] = Bfield[0][0].get();
+ mfs[4] = Bfield[0][1].get();
+ mfs[5] = Bfield[0][2].get();
+
+ IntVect flags[6];
+ flags[0] = WarpX::Ex_nodal_flag;
+ flags[1] = WarpX::Ey_nodal_flag;
+ flags[2] = WarpX::Ez_nodal_flag;
+ flags[3] = WarpX::Bx_nodal_flag;
+ flags[4] = WarpX::By_nodal_flag;
+ flags[5] = WarpX::Bz_nodal_flag;
+
+ const unsigned npoints = x.size();
+ const int ncomps = 6;
+ result.resize(ncomps);
+ for (int i = 0; i < ncomps; i++) {
+ result[i].resize(npoints, 0);
+ }
+
+ BoxArray ba = Efield[0][0]->boxArray();
+ const IndexType correct_typ(IntVect::TheZeroVector());
+ ba.convert(correct_typ);
+ std::vector< std::pair<int, Box> > isects;
+
+ IntVect iv;
+ const Geometry& g = geom[0];
+ for (int i = 0; i < npoints; ++i) {
+ for (int comp = 0; comp < ncomps; comp++) {
+ D_TERM(iv[0]=floor((x[i] - g.ProbLo(0) + 0.5*g.CellSize(0)*flags[comp][0])/g.CellSize(0));,
+ iv[1]=floor((y[i] - g.ProbLo(1) + 0.5*g.CellSize(1)*flags[comp][1])/g.CellSize(1));,
+ iv[2]=floor((z[i] - g.ProbLo(2) + 0.5*g.CellSize(2)*flags[comp][2])/g.CellSize(2)););
+
+ ba.intersections(Box(iv, iv), isects, true, 0);
+ const int grid = isects[0].first;
+ const int who = dmap[0][grid];
+ if (who == ParallelDescriptor::MyProc()) {
+ result[comp][i] = (*mfs[comp])[grid](iv);
+ }
+ }
+ }
+
+ for (int i = 0; i < ncomps; i++) {
+ ParallelDescriptor::ReduceRealSum(result[i].dataPtr(), result[i].size());
+ }
+}
+#endif
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index a0852ae90..3e96213a8 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -289,20 +289,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);
+ auto& ppc = dynamic_cast<PhysicalParticleContainer&>(pc);
+ ppc.AddParticles(lev, particleBox);
}
}
}
diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H
index ac0ebfd8c..12152c286 100644
--- a/Source/WarpXParticleContainer.H
+++ b/Source/WarpXParticleContainer.H
@@ -69,9 +69,6 @@ public:
amrex::Real x, amrex::Real y, amrex::Real z,
const std::array<amrex::Real,PIdx::nattribs>& 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 2965824ff..31e4702ed 100644
--- a/Source/WarpXParticleContainer.cpp
+++ b/Source/WarpXParticleContainer.cpp
@@ -66,49 +66,6 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile,
}
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<Real,PIdx::nattribs> 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,
int nattr, const Real* attr, int uniqueparticles)
diff --git a/Source/WarpXProb.cpp b/Source/WarpXProb.cpp
index c90b0f6a8..4b11f75f1 100644
--- a/Source/WarpXProb.cpp
+++ b/Source/WarpXProb.cpp
@@ -1,9 +1,3 @@
-
-//
-// Each problem must have its own version of WarpX::InitLevelData(int lev)
-// to initialize the field data on this level.
-//
-
#include <WarpX.H>
using namespace amrex;
@@ -11,14 +5,9 @@ using namespace amrex;
void
WarpX::InitLevelData (int lev)
{
- static_assert(false,
- "Each problem must have its own version of WarpX::InitLevelData(int lev)");
-
-#if 0
for (int i = 0; i < 3; ++i) {
current[lev][i]->setVal(0.0);
Efield[lev][i]->setVal(0.0);
Bfield[lev][i]->setVal(0.0);
}
-#endif
}
diff --git a/Source/WarpXRegrid.cpp b/Source/WarpXRegrid.cpp
index 32508b9ba..d3eced1e6 100644
--- a/Source/WarpXRegrid.cpp
+++ b/Source/WarpXRegrid.cpp
@@ -217,7 +217,7 @@ WarpX::getCostCountDM (const Array<long>& cost, const BoxArray& ba)
{
DistributionMapping res;
int nprocs = ParallelDescriptor::NProcs();
- const int factor = 1.5; // A process can get up to 'factor' times of the average number of boxes.
+ const Real factor = 1.5; // A process can get up to 'factor' times of the average number of boxes.
int nmax = (cost.size()+nprocs-1) / nprocs * factor;
Real eff;
res.KnapSackProcessorMap(cost, nprocs, &eff, true, nmax);
diff --git a/Source/WarpX_picsar.F90 b/Source/WarpX_picsar.F90
index d3d7be39e..6dfb16250 100644
--- a/Source/WarpX_picsar.F90
+++ b/Source/WarpX_picsar.F90
@@ -1,6 +1,7 @@
#if (BL_SPACEDIM == 3)
#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb3d_energy_conserving
+#define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic
#elif (BL_SPACEDIM == 2)
@@ -232,74 +233,9 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n
! Dimension 3
#if (BL_SPACEDIM==3)
-
- SELECT CASE(current_depo_algo)
-
- ! Scalar classical current deposition subroutines
- CASE(3)
-
- IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN
- CALL depose_jxjyjz_scalar_1_1_1(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard)
- ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN
- CALL depose_jxjyjz_scalar_2_2_2(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard)
- ELSE IF ((nox.eq.3).and.(noy.eq.3).and.(noz.eq.3)) THEN
- CALL depose_jxjyjz_scalar_3_3_3(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard)
- ELSE
- CALL pxr_depose_jxjyjz_esirkepov_n(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
- ENDIF
-
- ! Optimized classical current deposition
- CASE(2)
-
- IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN
- CALL depose_jxjyjz_vecHVv2_1_1_1(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard)
- ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN
- CALL depose_jxjyjz_vecHVv2_2_2_2(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard)
- ELSE IF ((nox.eq.3).and.(noy.eq.3).and.(noz.eq.3)) THEN
- CALL depose_jxjyjz_vecHVv3_3_3_3(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard)
- ELSE
- CALL pxr_depose_jxjyjz_esirkepov_n(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
- ENDIF
-
- ! Esirkepov non optimized
- CASE(1)
-
- CALL pxr_depose_jxjyjz_esirkepov_n(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
-
- ! Optimized Esirkepov
- CASE DEFAULT
-
- IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN
- CALL depose_jxjyjz_esirkepov_1_1_1(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
- ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN
- CALL depose_jxjyjz_esirkepov_2_2_2(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
- ELSE IF ((nox.eq.3).and.(noy.eq.3).and.(noz.eq.3)) THEN
- CALL depose_jxjyjz_esirkepov_3_3_3(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
- ELSE
- CALL pxr_depose_jxjyjz_esirkepov_n(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
- dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
- nox,noy,noz,.TRUE._c_long,.FALSE._c_long)
- ENDIF
-
- END SELECT
+ CALL WRPX_PXR_CURRENT_DEPOSITION(jx,jy,jz,np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin, &
+ dt,dx,dy,dz,nx,ny,nz,nxguard,nyguard,nzguard, &
+ nox,noy,noz,current_depo_algo)
! Dimension 2
#elif (BL_SPACEDIM==2)
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-13 13:33:05 +0000