aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Docs/source/usage/parameters.rst12
-rw-r--r--Examples/Tests/flux_injection/inputs_3d8
-rw-r--r--Examples/Tests/flux_injection/inputs_rz4
-rw-r--r--Python/pywarpx/picmi.py6
-rw-r--r--Source/Initialization/InjectorFlux.H146
-rw-r--r--Source/Initialization/PlasmaInjector.H9
-rw-r--r--Source/Initialization/PlasmaInjector.cpp45
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp20
8 files changed, 227 insertions, 23 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index f7eba5583..0a12dc82d 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -725,8 +725,8 @@ Particle initialization
For more information on the `openPMD format <https://github.com/openPMD>`__ and how to build WarpX with it, please visit :ref:`the install section <install-developers>`.
* ``NFluxPerCell``: Continuously inject a flux of macroparticles from a planar surface.
- The density specified by the density profile is interpreted to have the units of #/m^2/s.
This requires the additional parameters:
+ ``<species_name>.flux_profile`` (see the description of this parameter further below)
``<species_name>.surface_flux_pos`` (`double`, location of the injection plane [meter])
``<species_name>.flux_normal_axis`` (`x`, `y`, or `z` for 3D, `x` or `z` for 2D, or `r`, `t`, or `z` for RZ. When `flux_normal_axis` is `r` or `t`, the `x` and `y` components of the user-specified momentum distribution are interpreted as the `r` and `t` components respectively)
``<species_name>.flux_direction`` (`-1` or `+1`, direction of flux relative to the plane)
@@ -802,6 +802,16 @@ Particle initialization
``electrons.density_function(x,y,z) = "n0+n0*x**2*1.e12"`` where ``n0`` is a
user-defined constant, see above. WARNING: where ``density_function(x,y,z)`` is close to zero, particles will still be injected between ``xmin`` and ``xmax`` etc., with a null weight. This is undesirable because it results in useless computing. To avoid this, see option ``density_min`` below.
+* ``<species_name>.flux_profile`` (`string`)
+ Defines the expression of the flux, when using ``<species_name>.injection_style=NFluxPerCell``
+
+ * ``constant``: Constant flux. This requires the additional parameter ``<species_name>.flux``.
+ i.e., the injection flux in :math:`m^{-2}.s^{-1}`.
+
+ * ``parse_flux_function``: the flux is given by a function in the input file.
+ It requires the additional argument ``<species_name>.flux_function(x,y,z,t)``, which is a
+ mathematical expression for the flux of the species.
+
* ``<species_name>.density_min`` (`float`) optional (default `0.`)
Minimum plasma density. No particle is injected where the density is below this value.
diff --git a/Examples/Tests/flux_injection/inputs_3d b/Examples/Tests/flux_injection/inputs_3d
index e3de9dc36..7d277e7c2 100644
--- a/Examples/Tests/flux_injection/inputs_3d
+++ b/Examples/Tests/flux_injection/inputs_3d
@@ -38,8 +38,8 @@ electron.num_particles_per_cell = 100
electron.surface_flux_pos = -4.
electron.flux_normal_axis = y
electron.flux_direction = +1
-electron.profile = constant
-electron.density = 1.
+electron.flux_profile = parse_flux_function
+electron.flux_function(x,y,z,t) = "1."
electron.momentum_distribution_type = gaussianflux
electron.ux_th = 0.1
electron.uy_th = 0.1
@@ -53,8 +53,8 @@ proton.num_particles_per_cell = 100
proton.surface_flux_pos = 4.
proton.flux_normal_axis = x
proton.flux_direction = -1
-proton.profile = constant
-proton.density = 1.
+proton.flux_profile = constant
+proton.flux = 1.
proton.momentum_distribution_type = gaussianflux
proton.ux_th = 0.1
proton.ux_m = 0.05
diff --git a/Examples/Tests/flux_injection/inputs_rz b/Examples/Tests/flux_injection/inputs_rz
index 4d63e6701..7770ab1da 100644
--- a/Examples/Tests/flux_injection/inputs_rz
+++ b/Examples/Tests/flux_injection/inputs_rz
@@ -49,8 +49,8 @@ electron.xmin = 1.5
electron.xmax = 1.9
electron.zmin = -1.0
electron.zmax = 1.0
-electron.profile = constant
-electron.density = 1.
+electron.flux_profile = constant
+electron.flux = 1.
electron.momentum_distribution_type = gaussianflux
electron.uy_th = 0.
electron.uy_m = 10.
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index efa14badd..a607010c7 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -401,10 +401,10 @@ class UniformFluxDistribution(picmistandard.PICMI_UniformFluxDistribution, Densi
self.set_mangle_dict()
self.set_species_attributes(species, layout)
- species.profile = "constant"
- species.density = self.flux
+ species.flux_profile = "constant"
+ species.flux = self.flux
if density_scale is not None:
- species.density *= density_scale
+ species.flux *= density_scale
species.flux_normal_axis = self.flux_normal_axis
species.surface_flux_pos = self.surface_flux_position
species.flux_direction = self.flux_direction
diff --git a/Source/Initialization/InjectorFlux.H b/Source/Initialization/InjectorFlux.H
new file mode 100644
index 000000000..adfe75b3e
--- /dev/null
+++ b/Source/Initialization/InjectorFlux.H
@@ -0,0 +1,146 @@
+/* Copyright 2023 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef INJECTOR_FLUX_H_
+#define INJECTOR_FLUX_H_
+
+#include "Utils/WarpXConst.H"
+
+#include <AMReX.H>
+#include <AMReX_Array.H>
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_Math.H>
+#include <AMReX_Parser.H>
+#include <AMReX_REAL.H>
+
+#include <cmath>
+#include <string>
+
+// struct whose getFlux returns constant flux.
+struct InjectorFluxConstant
+{
+ InjectorFluxConstant (amrex::Real a_flux) noexcept : m_flux(a_flux) {}
+
+ AMREX_GPU_HOST_DEVICE
+ amrex::Real
+ getFlux (amrex::Real, amrex::Real, amrex::Real, amrex::Real) const noexcept
+ {
+ return m_flux;
+ }
+
+private:
+ amrex::Real m_flux;
+};
+
+// struct whose getFlux returns local flux computed from parser.
+struct InjectorFluxParser
+{
+ InjectorFluxParser (amrex::ParserExecutor<4> const& a_parser) noexcept
+ : m_parser(a_parser) {}
+
+ AMREX_GPU_HOST_DEVICE
+ amrex::Real
+ getFlux (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t) const noexcept
+ {
+ return m_parser(x,y,z,t);
+ }
+
+ amrex::ParserExecutor<4> m_parser;
+};
+
+// Base struct for flux injector.
+// InjectorFlux contains a union (called Object) that holds any one
+// instance of:
+// - InjectorFluxConstant : to generate constant flux;
+// - InjectorFluxParser : to generate flux from parser;
+// The choice is made at runtime, depending in the constructor called.
+// This mimics virtual functions.
+struct InjectorFlux
+{
+ // This constructor stores a InjectorFluxConstant in union object.
+ InjectorFlux (InjectorFluxConstant* t, amrex::Real a_flux)
+ : type(Type::constant),
+ object(t,a_flux)
+ { }
+
+ // This constructor stores a InjectorFluxParser in union object.
+ InjectorFlux (InjectorFluxParser* t, amrex::ParserExecutor<4> const& a_parser)
+ : type(Type::parser),
+ object(t,a_parser)
+ { }
+
+ // Explicitly prevent the compiler from generating copy constructors
+ // and copy assignment operators.
+ InjectorFlux (InjectorFlux const&) = delete;
+ InjectorFlux (InjectorFlux&&) = delete;
+ void operator= (InjectorFlux const&) = delete;
+ void operator= (InjectorFlux &&) = delete;
+
+ void clear ()
+ {
+ switch (type)
+ {
+ case Type::constant:
+ case Type::parser:
+ {
+ break;
+ }
+ }
+ }
+
+ // call getFlux from the object stored in the union
+ // (the union is called Object, and the instance is called object).
+ AMREX_GPU_HOST_DEVICE
+ amrex::Real
+ getFlux (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t) const noexcept
+ {
+ switch (type)
+ {
+ case Type::parser:
+ {
+ return object.parser.getFlux(x,y,z,t);
+ }
+ case Type::constant:
+ {
+ return object.constant.getFlux(x,y,z,t);
+ }
+ default:
+ {
+ amrex::Abort("InjectorFlux: unknown type");
+ return 0.0;
+ }
+ }
+ }
+
+private:
+ enum struct Type { constant, parser };
+ Type type;
+
+ // An instance of union Object constructs and stores any one of
+ // the objects declared (constant or parser).
+ union Object {
+ Object (InjectorFluxConstant*, amrex::Real a_flux) noexcept
+ : constant(a_flux) {}
+ Object (InjectorFluxParser*, amrex::ParserExecutor<4> const& a_parser) noexcept
+ : parser(a_parser) {}
+ InjectorFluxConstant constant;
+ InjectorFluxParser parser;
+ };
+ Object object;
+};
+
+// In order for InjectorFlux to be trivially copyable, its destructor
+// must be trivial. So we have to rely on a custom deleter for unique_ptr.
+struct InjectorFluxDeleter {
+ void operator () (InjectorFlux* p) const {
+ if (p) {
+ p->clear();
+ delete p;
+ }
+ }
+};
+
+#endif
diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H
index c841b6f86..0f33aa062 100644
--- a/Source/Initialization/PlasmaInjector.H
+++ b/Source/Initialization/PlasmaInjector.H
@@ -10,6 +10,7 @@
#define PLASMA_INJECTOR_H_
#include "InjectorDensity.H"
+#include "InjectorFlux.H"
#include "InjectorMomentum.H"
#include "TemperatureProperties.H"
#include "VelocityProperties.H"
@@ -119,6 +120,7 @@ public:
bool radially_weighted = true;
std::string str_density_function;
+ std::string str_flux_function;
std::string str_momentum_function_ux;
std::string str_momentum_function_uy;
std::string str_momentum_function_uz;
@@ -131,6 +133,7 @@ public:
InjectorPosition* getInjectorPosition ();
InjectorDensity* getInjectorDensity ();
+ InjectorFlux* getInjectorFlux ();
InjectorMomentum* getInjectorMomentum ();
protected:
@@ -140,6 +143,7 @@ protected:
PhysicalSpecies physical_species = PhysicalSpecies::unspecified;
amrex::Real density;
+ amrex::Real flux;
int species_id;
std::string species_name;
@@ -151,6 +155,10 @@ protected:
InjectorDensity* d_inj_rho = nullptr;
std::unique_ptr<amrex::Parser> density_parser;
+ std::unique_ptr<InjectorFlux,InjectorFluxDeleter> h_inj_flux;
+ InjectorFlux* d_inj_flux = nullptr;
+ std::unique_ptr<amrex::Parser> flux_parser;
+
std::unique_ptr<InjectorMomentum,InjectorMomentumDeleter> h_inj_mom;
InjectorMomentum* d_inj_mom = nullptr;
std::unique_ptr<amrex::Parser> ux_parser;
@@ -163,6 +171,7 @@ protected:
std::unique_ptr<VelocityProperties> h_mom_vel;
void parseDensity (const amrex::ParmParse& pp);
+ void parseFlux (const amrex::ParmParse& pp);
void parseMomentum (const amrex::ParmParse& pp);
};
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index 22c72f580..865e531e7 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -342,7 +342,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
(InjectorPositionRandomPlane*)nullptr,
xmin, xmax, ymin, ymax, zmin, zmax,
flux_normal_axis);
- parseDensity(pp_species_name);
+ parseFlux(pp_species_name);
parseMomentum(pp_species_name);
} else if (injection_style == "nuniformpercell") {
// Note that for RZ, three numbers are expected, r, theta, and z.
@@ -504,6 +504,16 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
#endif
}
+ if (h_inj_flux) {
+#ifdef AMREX_USE_GPU
+ d_inj_flux = static_cast<InjectorFlux*>
+ (amrex::The_Arena()->alloc(sizeof(InjectorFlux)));
+ amrex::Gpu::htod_memcpy_async(d_inj_flux, h_inj_flux.get(), sizeof(InjectorFlux));
+#else
+ d_inj_flux = h_inj_flux.get();
+#endif
+ }
+
if (h_inj_mom) {
#ifdef AMREX_USE_GPU
d_inj_mom = static_cast<InjectorMomentum*>
@@ -569,6 +579,33 @@ void PlasmaInjector::parseDensity (const amrex::ParmParse& pp)
}
}
+// Depending on injection type at runtime, initialize inj_flux
+// so that inj_flux->getFlux calls
+// InjectorFlux[Constant or Parser or etc.].getFlux.
+void PlasmaInjector::parseFlux (const amrex::ParmParse& pp)
+{
+ // parse flux information
+ std::string flux_prof_s;
+ pp.get("flux_profile", flux_prof_s);
+ std::transform(flux_prof_s.begin(), flux_prof_s.end(),
+ flux_prof_s.begin(), ::tolower);
+ if (flux_prof_s == "constant") {
+ utils::parser::getWithParser(pp, "flux", flux);
+ // Construct InjectorFlux with InjectorFluxConstant.
+ h_inj_flux.reset(new InjectorFlux((InjectorFluxConstant*)nullptr, flux));
+ } else if (flux_prof_s == "parse_flux_function") {
+ utils::parser::Store_parserString(
+ pp, "flux_function(x,y,z,t)", str_flux_function);
+ // Construct InjectorFlux with InjectorFluxParser.
+ flux_parser = std::make_unique<amrex::Parser>(
+ utils::parser::makeParser(str_flux_function,{"x","y","z","t"}));
+ h_inj_flux.reset(new InjectorFlux((InjectorFluxParser*)nullptr,
+ flux_parser->compile<4>()));
+ } else {
+ StringParseAbortMessage("Flux profile type", flux_prof_s);
+ }
+}
+
// Depending on injection type at runtime, initialize inj_mom
// so that inj_mom->getMomentum calls
// InjectorMomentum[Constant or Gaussian or etc.].getMomentum.
@@ -731,6 +768,12 @@ PlasmaInjector::getInjectorDensity ()
return d_inj_rho;
}
+InjectorFlux*
+PlasmaInjector::getInjectorFlux ()
+{
+ return d_inj_flux;
+}
+
InjectorMomentum*
PlasmaInjector::getInjectorMomentum ()
{
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 685f000d7..0793094ff 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -1475,10 +1475,8 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
}
InjectorPosition* inj_pos = plasma_injector->getInjectorPosition();
- InjectorDensity* inj_rho = plasma_injector->getInjectorDensity();
+ InjectorFlux* inj_flux = plasma_injector->getInjectorFlux();
InjectorMomentum* inj_mom = plasma_injector->getInjectorMomentum();
- const amrex::Real density_min = plasma_injector->density_min;
- const amrex::Real density_max = plasma_injector->density_max;
constexpr int level_zero = 0;
const amrex::Real t = WarpX::GetInstance().gett_new(level_zero);
@@ -1794,7 +1792,7 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
#ifdef WARPX_DIM_RZ
// Conversion from cylindrical to Cartesian coordinates
// Replace the x and y, setting an angle theta.
- // These x and y are used to get the momentum and density
+ // These x and y are used to get the momentum and flux
Real theta;
if (nmodes == 1 && rz_random_theta) {
// With only 1 mode, the angle doesn't matter so
@@ -1821,14 +1819,12 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
pu.y = sin_theta*ur + cos_theta*ut;
}
#endif
- Real dens = inj_rho->getDensity(ppos.x, ppos.y, ppos.z);
- // Remove particle if density below threshold
- if ( dens < density_min ){
+ Real flux = inj_flux->getFlux(ppos.x, ppos.y, ppos.z, t);
+ // Remove particle if flux is negative or 0
+ if ( flux <=0 ){
p.id() = -1;
continue;
}
- // Cut density if above threshold
- dens = amrex::min(dens, density_max);
if (loc_do_field_ionization) {
p_ion_level[ip] = loc_ionization_initial_level;
@@ -1854,12 +1850,12 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
#ifdef WARPX_DIM_RZ
// The particle weight is proportional to the user-specified
- // flux (denoted as `dens` here) and the emission surface within
+ // flux and the emission surface within
// one cell (captured partially by `scale_fac`).
// For cylindrical emission (flux_normal_axis==0
// or flux_normal_axis==2), the emission surface depends on
// the radius ; thus, the calculation is finalized here
- Real t_weight = dens * scale_fac * dt;
+ Real t_weight = flux * scale_fac * dt;
if (loc_flux_normal_axis != 1) {
if (radially_weighted) {
t_weight *= 2._rt*MathConst::pi*radial_position;
@@ -1872,7 +1868,7 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
}
const Real weight = t_weight;
#else
- const Real weight = dens * scale_fac * dt;
+ const Real weight = flux * scale_fac * dt;
#endif
pa[PIdx::w ][ip] = weight;
pa[PIdx::ux][ip] = pu.x;
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 14:54:57 +0000