diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Initialization/InjectorFlux.H | 146 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.H | 9 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 45 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 20 |
4 files changed, 207 insertions, 13 deletions
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; |