diff options
Diffstat (limited to 'Source/Initialization')
| -rw-r--r-- | Source/Initialization/CustomDensityProb.H | 49 | ||||
| -rw-r--r-- | Source/Initialization/CustomDensityProb.cpp | 12 | ||||
| -rw-r--r-- | Source/Initialization/CustomMomentumProb.H | 30 | ||||
| -rw-r--r-- | Source/Initialization/CustomMomentumProb.cpp | 14 | ||||
| -rw-r--r-- | Source/Initialization/InjectorDensity.H | 202 | ||||
| -rw-r--r-- | Source/Initialization/InjectorDensity.cpp | 77 | ||||
| -rw-r--r-- | Source/Initialization/InjectorMomentum.H | 223 | ||||
| -rw-r--r-- | Source/Initialization/InjectorMomentum.cpp | 40 | ||||
| -rw-r--r-- | Source/Initialization/InjectorPosition.H | 146 | ||||
| -rw-r--r-- | Source/Initialization/Make.package | 15 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.H | 293 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 345 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaProfiles.cpp | 41 |
13 files changed, 929 insertions, 558 deletions
diff --git a/Source/Initialization/CustomDensityProb.H b/Source/Initialization/CustomDensityProb.H new file mode 100644 index 000000000..b00830e6c --- /dev/null +++ b/Source/Initialization/CustomDensityProb.H @@ -0,0 +1,49 @@ +#ifndef CUSTOM_DENSITY_PROB_H_ +#define CUSTOM_DENSITY_PROB_H_ + +#include <AMReX_ParmParse.H> +#include <AMReX_Arena.H> +#include <AMReX_Gpu.H> +#include <AMReX_Dim3.H> + +// An example of Custom Density Profile + +// struct whose getDensity returns density at a given position computed from +// a custom function, with runtime input parameters. +struct InjectorDensityCustom +{ + InjectorDensityCustom (std::string const& species_name) + : p(nullptr) + { + // Read parameters for custom density profile from file, and + // store them in managed memory. + amrex::ParmParse pp(species_name); + std::vector<amrex::Real> v; + pp.getarr("custom_profile_params", v); + p = static_cast<amrex::Real*> + (amrex::The_Managed_Arena()->alloc(sizeof(amrex::Real)*v.size())); + for (int i = 0; i < static_cast<int>(v.size()); ++i) { + p[i] = v[i]; + } + } + + // Return density at given position, using user-defined parameters + // stored in p. + AMREX_GPU_HOST_DEVICE + amrex::Real + getDensity (amrex::Real, amrex::Real, amrex::Real) const noexcept + { + return p[0]; + } + + // Note that we are not allowed to have non-trivial destructor. + // So we rely on clear() to free memory. + void clear () { + amrex::The_Managed_Arena()->free(p); + } + +private: + amrex::Real* p; +}; + +#endif diff --git a/Source/Initialization/CustomDensityProb.cpp b/Source/Initialization/CustomDensityProb.cpp deleted file mode 100644 index 3efcb13c5..000000000 --- a/Source/Initialization/CustomDensityProb.cpp +++ /dev/null @@ -1,12 +0,0 @@ -#include <PlasmaInjector.H> - -#include <iostream> - -using namespace amrex; - -/// -/// This "custom" density profile just does constant -/// -Real CustomDensityProfile::getDensity(Real x, Real y, Real z) const { - return params[0]; -} diff --git a/Source/Initialization/CustomMomentumProb.H b/Source/Initialization/CustomMomentumProb.H new file mode 100644 index 000000000..f8bc29a05 --- /dev/null +++ b/Source/Initialization/CustomMomentumProb.H @@ -0,0 +1,30 @@ +#ifndef CUSTOM_MOMENTUM_PROB_H +#define CUSTOM_MOMENTUM_PROB_H + +#include <AMReX_ParmParse.H> +#include <AMReX_Gpu.H> +#include <AMReX_Arena.H> +#include <AMReX_Dim3.H> + +// An example of Custom Momentum Profile + +// struct whose getDensity returns momentum at a given position computed from +// a custom function. +struct InjectorMomentumCustom +{ + InjectorMomentumCustom (std::string const& /*a_species_name*/) {} + + // Return momentum at given position (illustration: momentum=0). + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real, amrex::Real, amrex::Real) const noexcept + { + return {0., 0., 0.}; + } + + // Note that we are not allowed to have non-trivial destructor. + // So we rely on clear() to free memory if needed. + void clear () { } +}; + +#endif diff --git a/Source/Initialization/CustomMomentumProb.cpp b/Source/Initialization/CustomMomentumProb.cpp deleted file mode 100644 index fa21252d0..000000000 --- a/Source/Initialization/CustomMomentumProb.cpp +++ /dev/null @@ -1,14 +0,0 @@ -#include <PlasmaInjector.H> - -#include <iostream> - -using namespace amrex; - -/// -/// This "custom" momentum distribution just does 0 momentum -/// -void CustomMomentumDistribution::getMomentum(vec3& u, Real x, Real y, Real z) { - u[0] = 0; - u[1] = 0; - u[2] = 0; -} diff --git a/Source/Initialization/InjectorDensity.H b/Source/Initialization/InjectorDensity.H new file mode 100644 index 000000000..b7f5c26eb --- /dev/null +++ b/Source/Initialization/InjectorDensity.H @@ -0,0 +1,202 @@ +#ifndef INJECTOR_DENSITY_H_ +#define INJECTOR_DENSITY_H_ + +#include <AMReX_Gpu.H> +#include <AMReX_Dim3.H> +#include <GpuParser.H> +#include <CustomDensityProb.H> +#include <WarpXConst.H> + +// struct whose getDensity returns constant density. +struct InjectorDensityConstant +{ + InjectorDensityConstant (amrex::Real a_rho) noexcept : m_rho(a_rho) {} + + AMREX_GPU_HOST_DEVICE + amrex::Real + getDensity (amrex::Real, amrex::Real, amrex::Real) const noexcept + { + return m_rho; + } + +private: + amrex::Real m_rho; +}; + +// struct whose getDensity returns local density computed from parser. +struct InjectorDensityParser +{ + InjectorDensityParser (WarpXParser const& a_parser) noexcept + : m_parser(a_parser) {} + + AMREX_GPU_HOST_DEVICE + amrex::Real + getDensity (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + return m_parser(x,y,z); + } + + // InjectorDensityParser constructs this GpuParser from WarpXParser. + GpuParser m_parser; +}; + +// struct whose getDensity returns local density computed from predefined profile. +struct InjectorDensityPredefined +{ + InjectorDensityPredefined (std::string const& a_species_name) noexcept; + + void clear (); + + AMREX_GPU_HOST_DEVICE + amrex::Real + getDensity (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + // Choices for profile are: + // - parabolic_channel + switch (profile) + { + case Profile::parabolic_channel: + { + amrex::Real z_start = p[0]; + amrex::Real ramp_up = p[1]; + amrex::Real plateau = p[2]; + amrex::Real ramp_down = p[3]; + amrex::Real rc = p[4]; + amrex::Real n0 = p[5]; + amrex::Real n; + amrex::Real kp = PhysConst::q_e/PhysConst::c + *std::sqrt( n0/(PhysConst::m_e*PhysConst::ep0) ); + + if ((z-z_start)>=0 and + (z-z_start)<ramp_up ) { + n = (z-z_start)/ramp_up; + } else if ((z-z_start)>=ramp_up and + (z-z_start)< ramp_up+plateau ) { + n = 1.; + } else if ((z-z_start)>=ramp_up+plateau and + (z-z_start)< ramp_up+plateau+ramp_down) { + n = 1.-((z-z_start)-ramp_up-plateau)/ramp_down; + } else { + n = 0.; + } + n *= n0*(1.+4.*(x*x+y*y)/(kp*kp*rc*rc*rc*rc)); + return n; + } + default: + amrex::Abort("InjectorDensityPredefined: how did we get here?"); + return 0.0; + } + } + +private: + enum struct Profile { null, parabolic_channel }; + Profile profile; + amrex::Real* p; +}; + +// Base struct for density injector. +// InjectorDensity contains a union (called Object) that holds any one +// instance of: +// - InjectorDensityConstant : to generate constant density; +// - InjectorDensityParser : to generate density from parser; +// - InjectorDensityCustom : to generate density from custom profile; +// - InjectorDensityPredefined: to generate density from predefined profile; +// The choice is made at runtime, depending in the constructor called. +// This mimics virtual functions, except the struct is stored in managed memory +// and member functions are made __host__ __device__ to run on CPU and GPU. +// This struct inherits from amrex::Gpu::Managed to provide new and delete +// operators in managed memory when running on GPU. Nothing special on CPU. +struct InjectorDensity + : public amrex::Gpu::Managed +{ + // This constructor stores a InjectorDensityConstant in union object. + InjectorDensity (InjectorDensityConstant* t, amrex::Real a_rho) + : type(Type::constant), + object(t,a_rho) + { } + + // This constructor stores a InjectorDensityParser in union object. + InjectorDensity (InjectorDensityParser* t, WarpXParser const& a_parser) + : type(Type::parser), + object(t,a_parser) + { } + + // This constructor stores a InjectorDensityCustom in union object. + InjectorDensity (InjectorDensityCustom* t, std::string const& a_species_name) + : type(Type::custom), + object(t,a_species_name) + { } + + // This constructor stores a InjectorDensityPredefined in union object. + InjectorDensity (InjectorDensityPredefined* t, std::string const& a_species_name) + : type(Type::predefined), + object(t,a_species_name) + { } + + // Explicitly prevent the compiler from generating copy constructors + // and copy assignment operators. + InjectorDensity (InjectorDensity const&) = delete; + InjectorDensity (InjectorDensity&&) = delete; + void operator= (InjectorDensity const&) = delete; + void operator= (InjectorDensity &&) = delete; + + ~InjectorDensity (); + + std::size_t sharedMemoryNeeded () const noexcept; + + // call getDensity from the object stored in the union + // (the union is called Object, and the instance is called object). + AMREX_GPU_HOST_DEVICE + amrex::Real + getDensity (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + switch (type) + { + case Type::parser: + { + return object.parser.getDensity(x,y,z); + } + case Type::constant: + { + return object.constant.getDensity(x,y,z); + } + case Type::custom: + { + return object.custom.getDensity(x,y,z); + } + case Type::predefined: + { + return object.predefined.getDensity(x,y,z); + } + default: + { + amrex::Abort("InjectorDensity: unknown type"); + return 0.0; + } + } + } + +private: + enum struct Type { constant, custom, predefined, parser }; + Type type; + + // An instance of union Object constructs and stores any one of + // the objects declared (constant or parser or custom or predefined). + union Object { + Object (InjectorDensityConstant*, amrex::Real a_rho) noexcept + : constant(a_rho) {} + Object (InjectorDensityParser*, WarpXParser const& a_parser) noexcept + : parser(a_parser) {} + Object (InjectorDensityCustom*, std::string const& a_species_name) noexcept + : custom(a_species_name) {} + Object (InjectorDensityPredefined*, std::string const& a_species_name) noexcept + : predefined(a_species_name) {} + InjectorDensityConstant constant; + InjectorDensityParser parser; + InjectorDensityCustom custom; + InjectorDensityPredefined predefined; + }; + Object object; +}; + +#endif diff --git a/Source/Initialization/InjectorDensity.cpp b/Source/Initialization/InjectorDensity.cpp new file mode 100644 index 000000000..7fed85b75 --- /dev/null +++ b/Source/Initialization/InjectorDensity.cpp @@ -0,0 +1,77 @@ +#include <PlasmaInjector.H> + +using namespace amrex; + +InjectorDensity::~InjectorDensity () +{ + switch (type) + { + case Type::parser: + { + object.parser.m_parser.clear(); + break; + } + case Type::custom: + { + object.custom.clear(); + break; + } + case Type::predefined: + { + object.predefined.clear(); + break; + } + } +} + +// Compute the amount of memory needed in GPU Shared Memory. +std::size_t +InjectorDensity::sharedMemoryNeeded () const noexcept +{ + switch (type) + { + case Type::parser: + { + // For parser injector, the 3D position of each particle + // is stored in shared memory. + return amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3; + } + default: + return 0; + } +} + +InjectorDensityPredefined::InjectorDensityPredefined ( + std::string const& a_species_name) noexcept + : profile(Profile::null) +{ + ParmParse pp(a_species_name); + + std::vector<amrex::Real> v; + // Read parameters for the predefined plasma profile, + // and store them in managed memory + pp.getarr("predefined_profile_params", v); + p = static_cast<amrex::Real*> + (amrex::The_Managed_Arena()->alloc(sizeof(amrex::Real)*v.size())); + for (int i = 0; i < static_cast<int>(v.size()); ++i) { + p[i] = v[i]; + } + + // Parse predefined profile name, and update member variable profile. + std::string which_profile_s; + pp.query("predefined_profile_name", which_profile_s); + std::transform(which_profile_s.begin(), which_profile_s.end(), + which_profile_s.begin(), ::tolower); + if (which_profile_s == "parabolic_channel"){ + profile = Profile::parabolic_channel; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(v.size() > 6, + "InjectorDensityPredefined::parabolic_channel: not enough parameters"); + } +} + +// Note that we are not allowed to have non-trivial destructor. +// So we rely on clear() to free memory. +void InjectorDensityPredefined::clear () +{ + amrex::The_Managed_Arena()->free(p); +} diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H new file mode 100644 index 000000000..399ee7759 --- /dev/null +++ b/Source/Initialization/InjectorMomentum.H @@ -0,0 +1,223 @@ +#ifndef INJECTOR_MOMENTUM_H_ +#define INJECTOR_MOMENTUM_H_ + +#include <AMReX_Gpu.H> +#include <AMReX_Dim3.H> +#include <GpuParser.H> +#include <CustomMomentumProb.H> + +// struct whose getMomentum returns constant momentum. +struct InjectorMomentumConstant +{ + InjectorMomentumConstant (amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) noexcept + : m_ux(a_ux), m_uy(a_uy), m_uz(a_uz) {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real, amrex::Real, amrex::Real) const noexcept + { + return amrex::XDim3{m_ux,m_uy,m_uz}; + } +private: + amrex::Real m_ux, m_uy, m_uz; +}; + +// struct whose getMomentum returns momentum for 1 particle, from random +// gaussian distribution. +struct InjectorMomentumGaussian +{ + InjectorMomentumGaussian (amrex::Real a_ux_m, amrex::Real a_uy_m, + amrex::Real a_uz_m, amrex::Real a_ux_th, + amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept + : m_ux_m(a_ux_m), m_uy_m(a_uy_m), m_uz_m(a_uz_m), + m_ux_th(a_ux_th), m_uy_th(a_uy_th), m_uz_th(a_uz_th) + {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + return amrex::XDim3{amrex::RandomNormal(m_ux_m, m_ux_th), + amrex::RandomNormal(m_uy_m, m_uy_th), + amrex::RandomNormal(m_uz_m, m_uz_th)}; + } +private: + amrex::Real m_ux_m, m_uy_m, m_uz_m; + amrex::Real m_ux_th, m_uy_th, m_uz_th; +}; + +// struct whose getMomentum returns momentum for 1 particle, for +// radial expansion +struct InjectorMomentumRadialExpansion +{ + InjectorMomentumRadialExpansion (amrex::Real a_u_over_r) noexcept + : u_over_r(a_u_over_r) + {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + return {x*u_over_r, y*u_over_r, z*u_over_r}; + } + +private: + amrex::Real u_over_r; +}; + +// struct whose getMomentumm returns local momentum computed from parser. +struct InjectorMomentumParser +{ + InjectorMomentumParser (WarpXParser const& a_ux_parser, + WarpXParser const& a_uy_parser, + WarpXParser const& a_uz_parser) noexcept + : m_ux_parser(a_ux_parser), m_uy_parser(a_uy_parser), + m_uz_parser(a_uz_parser) {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + return amrex::XDim3{m_ux_parser(x,y,z),m_uy_parser(x,y,z),m_uz_parser(x,y,z)}; + } + + GpuParser m_ux_parser, m_uy_parser, m_uz_parser; +}; + +// Base struct for momentum injector. +// InjectorMomentum contains a union (called Object) that holds any one +// instance of: +// - InjectorMomentumConstant : to generate constant density; +// - InjectorMomentumGaussian : to generate gaussian distribution; +// - InjectorMomentumRadialExpansion: to generate radial expansion; +// - InjectorMomentumParser : to generate momentum from parser; +// The choice is made at runtime, depending in the constructor called. +// This mimics virtual functions, except the struct is stored in managed memory +// and member functions are made __host__ __device__ to run on CPU and GPU. +// This struct inherits from amrex::Gpu::Managed to provide new and delete +// operators in managed memory when running on GPU. Nothing special on CPU. +struct InjectorMomentum + : public amrex::Gpu::Managed +{ + // This constructor stores a InjectorMomentumConstant in union object. + InjectorMomentum (InjectorMomentumConstant* t, + amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) + : type(Type::constant), + object(t, a_ux, a_uy, a_uz) + { } + + // This constructor stores a InjectorMomentumParser in union object. + InjectorMomentum (InjectorMomentumParser* t, + WarpXParser const& a_ux_parser, + WarpXParser const& a_uy_parser, + WarpXParser const& a_uz_parser) + : type(Type::parser), + object(t, a_ux_parser, a_uy_parser, a_uz_parser) + { } + + // This constructor stores a InjectorMomentumGaussian in union object. + InjectorMomentum (InjectorMomentumGaussian* t, + amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m, + amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th) + : type(Type::gaussian), + object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th) + { } + + // This constructor stores a InjectorMomentumCustom in union object. + InjectorMomentum (InjectorMomentumCustom* t, + std::string const& a_species_name) + : type(Type::custom), + object(t, a_species_name) + { } + + // This constructor stores a InjectorMomentumRadialExpansion in union object. + InjectorMomentum (InjectorMomentumRadialExpansion* t, + amrex::Real u_over_r) + : type(Type::radial_expansion), + object(t, u_over_r) + { } + + // Explicitly prevent the compiler from generating copy constructors + // and copy assignment operators. + InjectorMomentum (InjectorMomentum const&) = delete; + InjectorMomentum (InjectorMomentum&&) = delete; + void operator= (InjectorMomentum const&) = delete; + void operator= (InjectorMomentum &&) = delete; + + ~InjectorMomentum (); + + std::size_t sharedMemoryNeeded () const noexcept; + + // call getMomentum from the object stored in the union + // (the union is called Object, and the instance is called object). + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + switch (type) + { + case Type::parser: + { + return object.parser.getMomentum(x,y,z); + } + case Type::gaussian: + { + return object.gaussian.getMomentum(x,y,z); + } + case Type::constant: + { + return object.constant.getMomentum(x,y,z); + } + case Type::radial_expansion: + { + return object.radial_expansion.getMomentum(x,y,z); + } + case Type::custom: + { + return object.custom.getMomentum(x,y,z); + } + default: + { + amrex::Abort("InjectorMomentum: unknown type"); + return {0.0,0.0,0.0}; + } + } + } + +private: + enum struct Type { constant, custom, gaussian, radial_expansion, parser }; + Type type; + + // An instance of union Object constructs and stores any one of + // the objects declared (constant or custom or gaussian or + // radial_expansion or parser). + union Object { + Object (InjectorMomentumConstant*, + amrex::Real a_ux, amrex::Real a_uy, amrex::Real a_uz) noexcept + : constant(a_ux,a_uy,a_uz) {} + Object (InjectorMomentumCustom*, + std::string const& a_species_name) noexcept + : custom(a_species_name) {} + Object (InjectorMomentumGaussian*, + amrex::Real a_ux_m, amrex::Real a_uy_m, + amrex::Real a_uz_m, amrex::Real a_ux_th, + amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept + : gaussian(a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th) {} + Object (InjectorMomentumRadialExpansion*, + amrex::Real u_over_r) noexcept + : radial_expansion(u_over_r) {} + Object (InjectorMomentumParser*, + WarpXParser const& a_ux_parser, + WarpXParser const& a_uy_parser, + WarpXParser const& a_uz_parser) noexcept + : parser(a_ux_parser, a_uy_parser, a_uz_parser) {} + InjectorMomentumConstant constant; + InjectorMomentumCustom custom; + InjectorMomentumGaussian gaussian; + InjectorMomentumRadialExpansion radial_expansion; + InjectorMomentumParser parser; + }; + Object object; +}; + +#endif diff --git a/Source/Initialization/InjectorMomentum.cpp b/Source/Initialization/InjectorMomentum.cpp new file mode 100644 index 000000000..a197b5bef --- /dev/null +++ b/Source/Initialization/InjectorMomentum.cpp @@ -0,0 +1,40 @@ +#include <PlasmaInjector.H> + +using namespace amrex; + +InjectorMomentum::~InjectorMomentum () +{ + switch (type) + { + case Type::parser: + { + object.parser.m_ux_parser.clear(); + object.parser.m_uy_parser.clear(); + object.parser.m_uz_parser.clear(); + break; + } + case Type::custom: + { + object.custom.clear(); + break; + } + } +} + +// Compute the amount of memory needed in GPU Shared Memory. +std::size_t +InjectorMomentum::sharedMemoryNeeded () const noexcept +{ + switch (type) + { + case Type::parser: + { + // For parser injector, the 3D position of each particle + // is stored in shared memory. + return amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3; + } + default: + return 0; + } +} + diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H new file mode 100644 index 000000000..19bb092dd --- /dev/null +++ b/Source/Initialization/InjectorPosition.H @@ -0,0 +1,146 @@ +#ifndef INJECTOR_POSITION_H_ +#define INJECTOR_POSITION_H_ + +#include <AMReX_Gpu.H> +#include <AMReX_Dim3.H> +#include <AMReX_Utility.H> + +// struct whose getPositionUnitBox returns x, y and z for a particle with +// random distribution inside a unit cell. +struct InjectorPositionRandom +{ + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + { + return amrex::XDim3{amrex::Random(), amrex::Random(), amrex::Random()}; + } +}; + +// struct whose getPositionUnitBox returns x, y and z for a particle with +// regular distribution inside a unit cell. +struct InjectorPositionRegular +{ + InjectorPositionRegular (amrex::Dim3 const& a_ppc) noexcept : ppc(a_ppc) {} + + // i_part: particle number within the cell, required to evenly space + // particles within the cell. + // ref_fac: the number of particles evenly-spaced within a cell + // is a_ppc*(ref_fac**AMREX_SPACEDIM). + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + { + int nx = ref_fac*ppc.x; + int ny = ref_fac*ppc.y; +#if (AMREX_SPACEDIM == 3) + int nz = ref_fac*ppc.z; +#else + int nz = 1; +#endif + int ix_part = i_part/(ny*nz); // written this way backward compatibility + int iz_part = (i_part-ix_part*(ny*nz)) / ny; + int iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; + return amrex::XDim3{(0.5+ix_part)/nx, (0.5+iy_part)/ny, (0.5+iz_part) / nz}; + } +private: + amrex::Dim3 ppc; +}; + +// Base struct for position injector. +// InjectorPosition contains a union (called Object) that holds any one +// instance of: +// - InjectorPositionRandom : to generate random distribution; +// - InjectorPositionRegular: to generate regular distribution. +// The choice is made at runtime, depending in the constructor called. +// This mimics virtual functions, except the struct is stored in managed memory +// and member functions are made __host__ __device__ to run on CPU and GPU. +// This struct inherits from amrex::Gpu::Managed to provide new and delete +// operators in managed memory when running on GPU. Nothing special on CPU. +struct InjectorPosition + : public amrex::Gpu::Managed +{ + // This constructor stores a InjectorPositionRandom in union object. + InjectorPosition (InjectorPositionRandom* t, + amrex::Real a_xmin, amrex::Real a_xmax, + amrex::Real a_ymin, amrex::Real a_ymax, + amrex::Real a_zmin, amrex::Real a_zmax) + : type(Type::random), + object(t), + xmin(a_xmin), xmax(a_xmax), + ymin(a_ymin), ymax(a_ymax), + zmin(a_zmin), zmax(a_zmax) + { } + + // This constructor stores a InjectorPositionRegular in union object. + InjectorPosition (InjectorPositionRegular* t, + amrex::Real a_xmin, amrex::Real a_xmax, + amrex::Real a_ymin, amrex::Real a_ymax, + amrex::Real a_zmin, amrex::Real a_zmax, + amrex::Dim3 const& a_ppc) + : type(Type::regular), + object(t, a_ppc), + xmin(a_xmin), xmax(a_xmax), + ymin(a_ymin), ymax(a_ymax), + zmin(a_zmin), zmax(a_zmax) + { } + + // Explicitly prevent the compiler from generating copy constructors + // and copy assignment operators. + InjectorPosition (InjectorPosition const&) = delete; + InjectorPosition (InjectorPosition&&) = delete; + void operator= (InjectorPosition const&) = delete; + void operator= (InjectorPosition &&) = delete; + + std::size_t sharedMemoryNeeded () const noexcept { return 0; } + + // call getPositionUnitBox from the object stored in the union + // (the union is called Object, and the instance is called object). + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + { + switch (type) + { + case Type::regular: + { + return object.regular.getPositionUnitBox(i_part, ref_fac); + } + default: + { + return object.random.getPositionUnitBox(i_part, ref_fac); + } + }; + } + + // bool: whether position specified is within bounds. + AMREX_GPU_HOST_DEVICE + bool + insideBounds (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + return (x < xmax and x >= xmin and + y < ymax and y >= ymin and + z < zmax and z >= zmin); + } + +private: + enum struct Type { random, regular }; + Type type; + + // An instance of union Object constructs and stores any one of + // the objects declared (random or regular). + union Object { + Object (InjectorPositionRandom*) noexcept : random() {} + Object (InjectorPositionRegular*, amrex::Dim3 const& a_ppc) noexcept + : regular(a_ppc) {} + InjectorPositionRandom random; + InjectorPositionRegular regular; + }; + Object object; + + amrex::Real xmin, xmax; + amrex::Real ymin, ymax; + amrex::Real zmin, zmax; +}; + +#endif diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index edcf402c9..2c6458b6d 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -1,9 +1,18 @@ -CEXE_sources += CustomDensityProb.cpp -CEXE_sources += PlasmaProfiles.cpp CEXE_sources += WarpXInitData.cpp -CEXE_sources += CustomMomentumProb.cpp + CEXE_sources += PlasmaInjector.cpp CEXE_headers += PlasmaInjector.H +CEXE_headers += InjectorPosition.H + +CEXE_headers += InjectorDensity.H +CEXE_sources += InjectorDensity.cpp + +CEXE_headers += InjectorMomentum.H +CEXE_sources += InjectorMomentum.cpp + +CEXE_headers += CustomDensityProb.H +CEXE_headers += CustomMomentumProb.H + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Initialization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Initialization diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index f998e217e..f7e86bff5 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -1,250 +1,16 @@ #ifndef PLASMA_INJECTOR_H_ #define PLASMA_INJECTOR_H_ -#include <array> +#include <InjectorPosition.H> +#include <InjectorDensity.H> +#include <InjectorMomentum.H> -#include "AMReX_REAL.H" +#include <array> #include <AMReX_Vector.H> #include <WarpXConst.H> #include <WarpXParser.H> -#include "AMReX_ParmParse.H" -#include "AMReX_Utility.H" - -enum class predefined_profile_flag { Null, parabolic_channel }; - -/// -/// 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::Vector<amrex::Real> params; -}; - -/// -/// This describes predefined density distributions. -/// -class PredefinedDensityProfile : public PlasmaDensityProfile -{ -public: - PredefinedDensityProfile(const std::string& species_name); - virtual amrex::Real getDensity(amrex::Real x, - amrex::Real y, - amrex::Real z) const override; - amrex::Real ParabolicChannel(amrex::Real x, - amrex::Real y, - amrex::Real z) const; -private: - predefined_profile_flag which_profile = predefined_profile_flag::Null; - amrex::Vector<amrex::Real> params; -}; - -/// -/// This describes a density function parsed in the input file. -/// -class ParseDensityProfile : public PlasmaDensityProfile -{ -public: - ParseDensityProfile(const std::string _parse_density_function); - virtual amrex::Real getDensity(amrex::Real x, - amrex::Real y, - amrex::Real z) const override; -private: - std::string _parse_density_function; - WarpXParser parser_density; -}; - -/// -/// 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, amrex::Real x, amrex::Real y, amrex::Real z) = 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, amrex::Real x, amrex::Real y, amrex::Real z) override; - -private: - amrex::Real _ux; - amrex::Real _uy; - amrex::Real _uz; -}; - -/// -/// This describes a custom momentum distribution. Users can supply -/// in their problem directory. -/// -/// -class CustomMomentumDistribution : public PlasmaMomentumDistribution -{ -public: - CustomMomentumDistribution(const std::string& species_name); - virtual void getMomentum(vec3& u, amrex::Real x, amrex::Real y, amrex::Real z) override; - -private: - amrex::Vector<amrex::Real> params; -}; - - -/// -/// 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 ux_th, - amrex::Real uy_th, - amrex::Real uz_th); - virtual void getMomentum(vec3& u, amrex::Real x, amrex::Real y, amrex::Real z) override; -private: - amrex::Real _ux_m; - amrex::Real _uy_m; - amrex::Real _uz_m; - amrex::Real _ux_th; - amrex::Real _uy_th; - amrex::Real _uz_th; -}; - -/// -/// This is a radially expanding momentum distribution -/// Particles will have a radial momentum proportional to their -/// radius, with proportionality constant u_over_r -class RadialExpansionMomentumDistribution : public PlasmaMomentumDistribution -{ -public: - RadialExpansionMomentumDistribution( amrex::Real u_over_r ); - virtual void getMomentum(vec3& u, amrex::Real x, amrex::Real y, amrex::Real z) override; -private: - amrex::Real _u_over_r; -}; - -/// -/// This describes a momentum distribution function parsed in the input file. -/// -class ParseMomentumFunction : public PlasmaMomentumDistribution -{ -public: - ParseMomentumFunction(const std::string _parse_momentum_function_ux, - const std::string _parse_momentum_function_uy, - const std::string _parse_momentum_function_uz); - virtual void getMomentum(vec3& u, - amrex::Real x, - amrex::Real y, - amrex::Real z) override; -private: - std::string _parse_momentum_function_ux; - std::string _parse_momentum_function_uy; - std::string _parse_momentum_function_uz; - WarpXParser parser_ux; - WarpXParser parser_uy; - WarpXParser parser_uz; -}; - - -/// -/// PlasmaParticlePosition describes how particles are initialized -/// into each cell box. Subclasses must define a -/// getPositionUnitBox function that returns the position of -/// particle number i_part in a unitary box. -/// -class PlasmaParticlePosition{ -public: - using vec3 = std::array<amrex::Real, 3>; - virtual ~PlasmaParticlePosition() {}; - virtual void getPositionUnitBox(vec3& r, int i_part, int ref_fac=1) = 0; -}; - -/// -/// Particles are initialized with a random uniform -/// distribution inside each cell -/// -class RandomPosition : public PlasmaParticlePosition{ -public: - RandomPosition(int num_particles_per_cell); - virtual void getPositionUnitBox(vec3& r, int i_part, int ref_fac=1) override; -private: - amrex::Real _x; - amrex::Real _y; - amrex::Real _z; - int _num_particles_per_cell; -}; - -/// -/// Particles are regularly distributed inside each cell. The user provides -/// a 3d (resp. 2d) vector num_particles_per_cell_each_dim that contains -/// the number of particles per cell along each dimension. -/// -class RegularPosition : public PlasmaParticlePosition{ -public: - RegularPosition(const amrex::Vector<int>& num_particles_per_cell_each_dim); - virtual void getPositionUnitBox(vec3& r, int i_part, int ref_fac=1) override; -private: - amrex::Real _x; - amrex::Real _y; - amrex::Real _z; - amrex::Vector<int> _num_particles_per_cell_each_dim; -}; +#include <AMReX_ParmParse.H> +#include <AMReX_Utility.H> /// /// The PlasmaInjector class parses and stores information about the plasma @@ -256,28 +22,23 @@ class PlasmaInjector public: - using vec3 = std::array<amrex::Real, 3>; - - PlasmaInjector(); - - PlasmaInjector(int ispecies, const std::string& name); + PlasmaInjector (); - amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z); + PlasmaInjector (int ispecies, const std::string& name); - bool insideBounds(amrex::Real x, amrex::Real y, amrex::Real z); + bool insideBounds (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept; int num_particles_per_cell; amrex::Vector<int> num_particles_per_cell_each_dim; - void getMomentum(vec3& u, amrex::Real x, amrex::Real y, amrex::Real z); + // gamma * beta + amrex::XDim3 getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept; - void getPositionUnitBox(vec3& r, int i_part, int ref_fac=1); + amrex::Real getCharge () {return charge;} + amrex::Real getMass () {return mass;} - amrex::Real getCharge() {return charge;} - amrex::Real getMass() {return mass;} - - bool doInjection() { return part_pos != NULL;} + bool doInjection () const noexcept { return inj_pos != NULL;} bool add_single_particle = false; amrex::Vector<amrex::Real> single_particle_pos; @@ -305,6 +66,21 @@ public: amrex::Real xmin, xmax; amrex::Real ymin, ymax; amrex::Real zmin, zmax; + amrex::Real density_min = 0; + amrex::Real density_max = std::numeric_limits<amrex::Real>::max(); + + InjectorPosition* getInjectorPosition (); + InjectorDensity* getInjectorDensity (); + InjectorMomentum* getInjectorMomentum (); + + // When running on GPU, injector for position, momentum and density store + // particle 3D positions in shared memory IF using the parser. + std::size_t + sharedMemoryNeeded () const noexcept { + return amrex::max(inj_pos->sharedMemoryNeeded(), + inj_rho->sharedMemoryNeeded(), + inj_mom->sharedMemoryNeeded()); + } protected: @@ -315,13 +91,12 @@ protected: int species_id; std::string species_name; - std::unique_ptr<PlasmaDensityProfile> rho_prof; - std::unique_ptr<PlasmaMomentumDistribution> mom_dist; - std::unique_ptr<PlasmaParticlePosition> part_pos; - - void parseDensity(amrex::ParmParse pp); - void parseMomentum(amrex::ParmParse pp); + std::unique_ptr<InjectorPosition> inj_pos; + std::unique_ptr<InjectorDensity > inj_rho; + std::unique_ptr<InjectorMomentum> inj_mom; + void parseDensity (amrex::ParmParse& pp); + void parseMomentum (amrex::ParmParse& pp); }; #endif diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index f9642d1b6..280f9e58a 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -55,192 +55,34 @@ namespace { } } -ConstantDensityProfile::ConstantDensityProfile(Real density) - : _density(density) -{} +PlasmaInjector::PlasmaInjector () {} -Real ConstantDensityProfile::getDensity(Real x, Real y, Real z) const -{ - return _density; -} - -CustomDensityProfile::CustomDensityProfile(const std::string& species_name) +PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) + : species_id(ispecies), species_name(name) { ParmParse pp(species_name); - pp.getarr("custom_profile_params", params); -} -PredefinedDensityProfile::PredefinedDensityProfile(const std::string& species_name) -{ - ParmParse pp(species_name); - std::string which_profile_s; - pp.getarr("predefined_profile_params", params); - pp.query("predefined_profile_name", which_profile_s); - if (which_profile_s == "parabolic_channel"){ - which_profile = predefined_profile_flag::parabolic_channel; - } -} - -ParseDensityProfile::ParseDensityProfile(std::string parse_density_function) - : _parse_density_function(parse_density_function) -{ - parser_density.define(parse_density_function); - parser_density.registerVariables({"x","y","z"}); - - ParmParse pp("my_constants"); - std::set<std::string> symbols = parser_density.symbols(); - symbols.erase("x"); - symbols.erase("y"); - symbols.erase("z"); // after removing variables, we are left with constants - for (auto it = symbols.begin(); it != symbols.end(); ) { - Real v; - if (pp.query(it->c_str(), v)) { - parser_density.setConstant(*it, v); - it = symbols.erase(it); - } else { - ++it; - } - } - for (auto const& s : symbols) { // make sure there no unknown symbols - amrex::Abort("ParseDensityProfile: Unknown symbol "+s); - } -} - -Real ParseDensityProfile::getDensity(Real x, Real y, Real z) const -{ - return parser_density.eval(x,y,z); -} - -ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, - Real uy, - Real uz) - : _ux(ux), _uy(uy), _uz(uz) -{} - -void ConstantMomentumDistribution::getMomentum(vec3& u, Real x, Real y, Real z) { - u[0] = _ux; - u[1] = _uy; - u[2] = _uz; -} - -CustomMomentumDistribution::CustomMomentumDistribution(const std::string& species_name) -{ - ParmParse pp(species_name); - pp.getarr("custom_momentum_params", params); -} - -GaussianRandomMomentumDistribution::GaussianRandomMomentumDistribution(Real ux_m, - Real uy_m, - Real uz_m, - Real ux_th, - Real uy_th, - Real uz_th) - : _ux_m(ux_m), _uy_m(uy_m), _uz_m(uz_m), _ux_th(ux_th), _uy_th(uy_th), _uz_th(uz_th) -{ -} - -void GaussianRandomMomentumDistribution::getMomentum(vec3& u, Real x, Real y, Real z) { - Real ux_th = amrex::RandomNormal(0.0, _ux_th); - Real uy_th = amrex::RandomNormal(0.0, _uy_th); - Real uz_th = amrex::RandomNormal(0.0, _uz_th); - - u[0] = _ux_m + ux_th; - u[1] = _uy_m + uy_th; - u[2] = _uz_m + uz_th; -} -RadialExpansionMomentumDistribution::RadialExpansionMomentumDistribution(Real u_over_r) : _u_over_r( u_over_r ) -{ -} - -void RadialExpansionMomentumDistribution::getMomentum(vec3& u, Real x, Real y, Real z) { - u[0] = _u_over_r * x; - u[1] = _u_over_r * y; - u[2] = _u_over_r * z; -} - -ParseMomentumFunction::ParseMomentumFunction(std::string parse_momentum_function_ux, - std::string parse_momentum_function_uy, - std::string parse_momentum_function_uz) - : _parse_momentum_function_ux(parse_momentum_function_ux), - _parse_momentum_function_uy(parse_momentum_function_uy), - _parse_momentum_function_uz(parse_momentum_function_uz) -{ - parser_ux.define(parse_momentum_function_ux); - parser_uy.define(parse_momentum_function_uy); - parser_uz.define(parse_momentum_function_uz); - - amrex::Array<std::reference_wrapper<WarpXParser>,3> parsers{parser_ux, parser_uy, parser_uz}; - ParmParse pp("my_constants"); - for (auto& p : parsers) { - auto& parser = p.get(); - parser.registerVariables({"x","y","z"}); - std::set<std::string> symbols = parser.symbols(); - symbols.erase("x"); - symbols.erase("y"); - symbols.erase("z"); // after removing variables, we are left with constants - for (auto it = symbols.begin(); it != symbols.end(); ) { - Real v; - if (pp.query(it->c_str(), v)) { - parser.setConstant(*it, v); - it = symbols.erase(it); - } else { - ++it; - } - } - for (auto const& s : symbols) { // make sure there no unknown symbols - amrex::Abort("ParseMomentumFunction: Unknown symbol "+s); - } - } -} - -void ParseMomentumFunction::getMomentum(vec3& u, Real x, Real y, Real z) -{ - u[0] = parser_ux.eval(x,y,z); - u[1] = parser_uy.eval(x,y,z); - u[2] = parser_uz.eval(x,y,z); -} - -RandomPosition::RandomPosition(int num_particles_per_cell): - _num_particles_per_cell(num_particles_per_cell) -{} - -void RandomPosition::getPositionUnitBox(vec3& r, int i_part, int ref_fac){ - r[0] = amrex::Random(); - r[1] = amrex::Random(); - r[2] = amrex::Random(); -} - -RegularPosition::RegularPosition(const amrex::Vector<int>& num_particles_per_cell_each_dim) - : _num_particles_per_cell_each_dim(num_particles_per_cell_each_dim) -{} + pp.query("radially_weighted", radially_weighted); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(radially_weighted, "ERROR: Only radially_weighted=true is supported"); -void RegularPosition::getPositionUnitBox(vec3& r, int i_part, int ref_fac) -{ - int nx = ref_fac*_num_particles_per_cell_each_dim[0]; - int ny = ref_fac*_num_particles_per_cell_each_dim[1]; -#if AMREX_SPACEDIM == 3 - int nz = ref_fac*_num_particles_per_cell_each_dim[2]; -#else - int nz = 1; -#endif - - int ix_part = i_part/(ny * nz); - int iy_part = (i_part % (ny * nz)) % ny; - int iz_part = (i_part % (ny * nz)) / ny; + // parse plasma boundaries + xmin = std::numeric_limits<amrex::Real>::lowest(); + ymin = std::numeric_limits<amrex::Real>::lowest(); + zmin = std::numeric_limits<amrex::Real>::lowest(); - r[0] = (0.5+ix_part)/nx; - r[1] = (0.5+iy_part)/ny; - r[2] = (0.5+iz_part)/nz; -} + xmax = std::numeric_limits<amrex::Real>::max(); + ymax = std::numeric_limits<amrex::Real>::max(); + zmax = std::numeric_limits<amrex::Real>::max(); -PlasmaInjector::PlasmaInjector(){ - part_pos = NULL; -} + pp.query("xmin", xmin); + pp.query("ymin", ymin); + pp.query("zmin", zmin); + pp.query("xmax", xmax); + pp.query("ymax", ymax); + pp.query("zmax", zmax); -PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) - : species_id(ispecies), species_name(name) -{ - ParmParse pp(species_name); + pp.query("density_min", density_min); + pp.query("density_max", density_max); // parse charge and mass std::string charge_s; @@ -290,9 +132,14 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) gaussian_beam = true; parseMomentum(pp); } + // Depending on injection type at runtime, initialize inj_pos + // so that inj_pos->getPositionUnitBox calls + // InjectorPosition[Random or Regular].getPositionUnitBox. else if (part_pos_s == "nrandompercell") { pp.query("num_particles_per_cell", num_particles_per_cell); - part_pos.reset(new RandomPosition(num_particles_per_cell)); + // Construct InjectorPosition with InjectorPositionRandom. + inj_pos.reset(new InjectorPosition((InjectorPositionRandom*)nullptr, + xmin, xmax, ymin, ymax, zmin, zmax)); parseDensity(pp); parseMomentum(pp); } else if (part_pos_s == "nuniformpercell") { @@ -301,7 +148,12 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) #if ( AMREX_SPACEDIM == 2 ) num_particles_per_cell_each_dim[2] = 1; #endif - part_pos.reset(new RegularPosition(num_particles_per_cell_each_dim)); + // Construct InjectorPosition from InjectorPositionRegular. + inj_pos.reset(new InjectorPosition((InjectorPositionRegular*)nullptr, + xmin, xmax, ymin, ymax, zmin, zmax, + Dim3{num_particles_per_cell_each_dim[0], + num_particles_per_cell_each_dim[1], + num_particles_per_cell_each_dim[2]})); num_particles_per_cell = num_particles_per_cell_each_dim[0] * num_particles_per_cell_each_dim[1] * num_particles_per_cell_each_dim[2]; @@ -310,52 +162,71 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) } else { StringParseAbortMessage("Injection style", part_pos_s); } +} - pp.query("radially_weighted", radially_weighted); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(radially_weighted, "ERROR: Only radially_weighted=true is supported"); - - // 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(); +namespace { +WarpXParser makeParser (std::string const& parse_function) +{ + WarpXParser parser(parse_function); + parser.registerVariables({"x","y","z"}); - pp.query("xmin", xmin); - pp.query("ymin", ymin); - pp.query("zmin", zmin); - pp.query("xmax", xmax); - pp.query("ymax", ymax); - pp.query("zmax", zmax); + ParmParse pp("my_constants"); + std::set<std::string> symbols = parser.symbols(); + symbols.erase("x"); + symbols.erase("y"); + symbols.erase("z"); // after removing variables, we are left with constants + for (auto it = symbols.begin(); it != symbols.end(); ) { + Real v; + if (pp.query(it->c_str(), v)) { + parser.setConstant(*it, v); + it = symbols.erase(it); + } else { + ++it; + } + } + for (auto const& s : symbols) { // make sure there no unknown symbols + amrex::Abort("PlasmaInjector::makeParser: Unknown symbol "+s); + } + return parser; +} } -void PlasmaInjector::parseDensity(ParmParse pp){ +// Depending on injection type at runtime, initialize inj_rho +// so that inj_rho->getDensity calls +// InjectorPosition[Constant or Custom or etc.].getDensity. +void PlasmaInjector::parseDensity (ParmParse& pp) +{ // 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); + std::transform(rho_prof_s.begin(), rho_prof_s.end(), + rho_prof_s.begin(), ::tolower); if (rho_prof_s == "constant") { pp.get("density", density); - rho_prof.reset(new ConstantDensityProfile(density)); + // Construct InjectorDensity with InjectorDensityConstant. + inj_rho.reset(new InjectorDensity((InjectorDensityConstant*)nullptr, density)); } else if (rho_prof_s == "custom") { - rho_prof.reset(new CustomDensityProfile(species_name)); + // Construct InjectorDensity with InjectorDensityCustom. + inj_rho.reset(new InjectorDensity((InjectorDensityCustom*)nullptr, species_name)); } else if (rho_prof_s == "predefined") { - rho_prof.reset(new PredefinedDensityProfile(species_name)); + // Construct InjectorDensity with InjectorDensityPredefined. + inj_rho.reset(new InjectorDensity((InjectorDensityPredefined*)nullptr,species_name)); } else if (rho_prof_s == "parse_density_function") { pp.get("density_function(x,y,z)", str_density_function); - rho_prof.reset(new ParseDensityProfile(str_density_function)); + // Construct InjectorDensity with InjectorDensityParser. + inj_rho.reset(new InjectorDensity((InjectorDensityParser*)nullptr, + makeParser(str_density_function))); } else { StringParseAbortMessage("Density profile type", rho_prof_s); } } -void PlasmaInjector::parseMomentum(ParmParse pp){ +// Depending on injection type at runtime, initialize inj_mom +// so that inj_mom->getMomentum calls +// InjectorMomentum[Constant or Custom or etc.].getMomentum. +void PlasmaInjector::parseMomentum (ParmParse& pp) +{ // parse momentum information std::string mom_dist_s; pp.get("momentum_distribution_type", mom_dist_s); @@ -370,9 +241,11 @@ void PlasmaInjector::parseMomentum(ParmParse pp){ pp.query("ux", ux); pp.query("uy", uy); pp.query("uz", uz); - mom_dist.reset(new ConstantMomentumDistribution(ux, uy, uz)); + // Construct InjectorMomentum with InjectorMomentumConstant. + inj_mom.reset(new InjectorMomentum((InjectorMomentumConstant*)nullptr, ux,uy, uz)); } else if (mom_dist_s == "custom") { - mom_dist.reset(new CustomMomentumDistribution(species_name)); + // Construct InjectorMomentum with InjectorMomentumCustom. + inj_mom.reset(new InjectorMomentum((InjectorMomentumCustom*)nullptr, species_name)); } else if (mom_dist_s == "gaussian") { Real ux_m = 0.; Real uy_m = 0.; @@ -386,42 +259,56 @@ void PlasmaInjector::parseMomentum(ParmParse pp){ pp.query("ux_th", ux_th); pp.query("uy_th", uy_th); pp.query("uz_th", uz_th); - mom_dist.reset(new GaussianRandomMomentumDistribution(ux_m, uy_m, uz_m, - ux_th, uy_th, uz_th)); + // Construct InjectorMomentum with InjectorMomentumGaussian. + inj_mom.reset(new InjectorMomentum((InjectorMomentumGaussian*)nullptr, + ux_m, uy_m, uz_m, ux_th, uy_th, uz_th)); } else if (mom_dist_s == "radial_expansion") { Real u_over_r = 0.; pp.query("u_over_r", u_over_r); - mom_dist.reset(new RadialExpansionMomentumDistribution(u_over_r)); + // Construct InjectorMomentum with InjectorMomentumRadialExpansion. + inj_mom.reset(new InjectorMomentum + ((InjectorMomentumRadialExpansion*)nullptr, u_over_r)); } else if (mom_dist_s == "parse_momentum_function") { pp.get("momentum_function_ux(x,y,z)", str_momentum_function_ux); pp.get("momentum_function_uy(x,y,z)", str_momentum_function_uy); pp.get("momentum_function_uz(x,y,z)", str_momentum_function_uz); - mom_dist.reset(new ParseMomentumFunction(str_momentum_function_ux, - str_momentum_function_uy, - str_momentum_function_uz)); + // Construct InjectorMomentum with InjectorMomentumParser. + inj_mom.reset(new InjectorMomentum((InjectorMomentumParser*)nullptr, + makeParser(str_momentum_function_ux), + makeParser(str_momentum_function_uy), + makeParser(str_momentum_function_uz))); } else { StringParseAbortMessage("Momentum distribution type", mom_dist_s); } } -void PlasmaInjector::getPositionUnitBox(vec3& r, int i_part, int ref_fac) { - return part_pos->getPositionUnitBox(r, i_part, ref_fac); +XDim3 PlasmaInjector::getMomentum (Real x, Real y, Real z) const noexcept +{ + return inj_mom->getMomentum(x, y, z); // gamma*beta +} + +bool PlasmaInjector::insideBounds (Real x, Real y, Real z) const noexcept +{ + return (x < xmax and x >= xmin and + y < ymax and y >= ymin and + z < zmax and z >= zmin); } -void PlasmaInjector::getMomentum(vec3& u, Real x, Real y, Real z) { - mom_dist->getMomentum(u, x, y, z); - u[0] *= PhysConst::c; - u[1] *= PhysConst::c; - u[2] *= PhysConst::c; +InjectorPosition* +PlasmaInjector::getInjectorPosition () +{ + return inj_pos.get(); } -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; +InjectorDensity* +PlasmaInjector::getInjectorDensity () +{ + return inj_rho.get(); } -Real PlasmaInjector::getDensity(Real x, Real y, Real z) { - return rho_prof->getDensity(x, y, z); +InjectorMomentum* +PlasmaInjector::getInjectorMomentum () +{ + return inj_mom.get(); } + diff --git a/Source/Initialization/PlasmaProfiles.cpp b/Source/Initialization/PlasmaProfiles.cpp deleted file mode 100644 index d9d207f7e..000000000 --- a/Source/Initialization/PlasmaProfiles.cpp +++ /dev/null @@ -1,41 +0,0 @@ -#include <PlasmaInjector.H> -#include <cmath> -#include <iostream> -#include <WarpXConst.H> - -using namespace amrex; - -Real PredefinedDensityProfile::getDensity(Real x, Real y, Real z) const { - Real n; - if ( which_profile == predefined_profile_flag::parabolic_channel ) { - n = ParabolicChannel(x,y,z); - } - return n; -} - -/// -/// plateau between linear upramp and downramp, and parab transverse profile -/// -Real PredefinedDensityProfile::ParabolicChannel(Real x, Real y, Real z) const { - // params = [z_start ramp_up plateau ramp_down rc n0] - Real z_start = params[0]; - Real ramp_up = params[1]; - Real plateau = params[2]; - Real ramp_down = params[3]; - Real rc = params[4]; - Real n0 = params[5]; - Real n; - Real kp = PhysConst::q_e/PhysConst::c*sqrt( n0/(PhysConst::m_e*PhysConst::ep0) ); - - if ((z-z_start)>=0 and (z-z_start)<ramp_up ) { - n = (z-z_start)/ramp_up; - } else if ((z-z_start)>=ramp_up and (z-z_start)<ramp_up+plateau ) { - n = 1; - } else if ((z-z_start)>=ramp_up+plateau and (z-z_start)<ramp_up+plateau+ramp_down) { - n = 1-((z-z_start)-ramp_up-plateau)/ramp_down; - } else { - n = 0; - } - n *= n0*(1+4*(x*x+y*y)/(kp*kp*std::pow(rc,4))); - return n; -} |