aboutsummaryrefslogtreecommitdiff
path: root/Source/Initialization
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Initialization')
-rw-r--r--Source/Initialization/CMakeLists.txt2
-rw-r--r--Source/Initialization/GetVelocity.H88
-rw-r--r--Source/Initialization/GetVelocity.cpp20
-rw-r--r--Source/Initialization/InjectorMomentum.H63
-rw-r--r--Source/Initialization/Make.package2
-rw-r--r--Source/Initialization/PlasmaInjector.H2
-rw-r--r--Source/Initialization/PlasmaInjector.cpp61
-rw-r--r--Source/Initialization/VelocityProperties.H54
-rw-r--r--Source/Initialization/VelocityProperties.cpp71
9 files changed, 285 insertions, 78 deletions
diff --git a/Source/Initialization/CMakeLists.txt b/Source/Initialization/CMakeLists.txt
index 74e8c9488..99cca7764 100644
--- a/Source/Initialization/CMakeLists.txt
+++ b/Source/Initialization/CMakeLists.txt
@@ -6,5 +6,7 @@ target_sources(WarpX
PlasmaInjector.cpp
WarpXInitData.cpp
TemperatureProperties.cpp
+ VelocityProperties.cpp
GetTemperature.cpp
+ GetVelocity.cpp
)
diff --git a/Source/Initialization/GetVelocity.H b/Source/Initialization/GetVelocity.H
new file mode 100644
index 000000000..7b1d5f213
--- /dev/null
+++ b/Source/Initialization/GetVelocity.H
@@ -0,0 +1,88 @@
+/* Copyright 2021 Hannah Klion
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef GET_VELOCITY_H_
+#define GET_VELOCITY_H_
+
+#include "VelocityProperties.H"
+
+/** Get velocity at a point on the grid
+ *
+ * Class to get velocity at a point on the grid, either constant (m_velocity)
+ * or a spatially varying value computed using the parser function (m_velocity_parser).
+ * It also stores the direction of the velocity field. It provides the velocity information
+ * held by the VelocityProperties instance passed to the constructor.
+ */
+struct GetVelocity
+{
+ /* Type of velocity initialization */
+ VelocityInitType m_type;
+
+ /* Velocity direction */
+ int m_dir; //! Index x=0, y=1, z=2
+ int m_sign_dir; //! Sign of the velocity direction positive=1, negative=-1
+
+ /** Constant velocity value, if m_type == VelConstantValue */
+ amrex::Real m_velocity;
+ /** Velocity parser function, if m_type == VelParserFunction */
+ amrex::ParserExecutor<3> m_velocity_parser;
+
+ /**
+ * \brief Construct the functor with information provided by vel
+ *
+ * \param[in] vel: const reference to the VelocityProperties object that will be used to
+ * populate the functor
+ */
+ GetVelocity (VelocityProperties const& vel) noexcept;
+
+ /**
+ * \brief Functor call. Returns the value of velocity at the location (x,y,z)
+ *
+ * \param[in] x x-coordinate of given location
+ * \param[in] y y-coordinate of given location
+ * \param[in] z z-cooridnate of given location
+ *
+ *\return value of velocity at (x,y,z).
+ * m_velocity if m_type is VelConstantValue
+ * m_velocity_parser(x,y,z) if m_type is VelParserFunction
+ */
+ AMREX_GPU_HOST_DEVICE
+ amrex::Real operator() (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
+ {
+ switch (m_type)
+ {
+ case (VelConstantValue):
+ {
+ return m_sign_dir * m_velocity;
+ }
+ case (VelParserFunction):
+ {
+ return m_sign_dir * m_velocity_parser(x,y,z);
+ }
+ default:
+ {
+ amrex::Abort("Get initial velocity: unknown type");
+ return 0.0;
+ }
+ }
+ }
+
+ /**
+ * \brief Returns the index of the direction of the bulk velocity
+ *
+ *\return index of direction of velocity.
+ * 0: x
+ * 1: y
+ * 2: z
+ */
+ AMREX_GPU_HOST_DEVICE
+ int direction () const noexcept
+ {
+ return m_dir;
+ }
+};
+#endif
diff --git a/Source/Initialization/GetVelocity.cpp b/Source/Initialization/GetVelocity.cpp
new file mode 100644
index 000000000..2d2a342a7
--- /dev/null
+++ b/Source/Initialization/GetVelocity.cpp
@@ -0,0 +1,20 @@
+/* Copyright 2021 Hannah Klion
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "GetVelocity.H"
+
+GetVelocity::GetVelocity (VelocityProperties const& vel) noexcept {
+ m_type = vel.m_type;
+ m_dir = vel.m_dir;
+ m_sign_dir = vel.m_sign_dir;
+ if (m_type == VelConstantValue) {
+ m_velocity = vel.m_velocity;
+ }
+ else if (m_type == VelParserFunction) {
+ m_velocity_parser = vel.m_ptr_velocity_parser->compile<3>();
+ }
+}
diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H
index 794f62eeb..ea6c961b9 100644
--- a/Source/Initialization/InjectorMomentum.H
+++ b/Source/Initialization/InjectorMomentum.H
@@ -10,7 +10,9 @@
#include "CustomMomentumProb.H"
#include "GetTemperature.H"
+#include "GetVelocity.H"
#include "TemperatureProperties.H"
+#include "VelocityProperties.H"
#include "Utils/WarpXConst.H"
#include <AMReX.H>
@@ -138,11 +140,10 @@ private:
struct InjectorMomentumBoltzmann
{
// Constructor whose inputs are:
- // a pointer to the initial temperature container temperature,
- // boost velocity/c beta,
- // and boost direction dir respectively.
- InjectorMomentumBoltzmann(GetTemperature const& t, amrex::Real b, int d) noexcept
- : dir(d), beta(b), temperature(t)
+ // a reference to the initial temperature container t,
+ // a reference to the initial velocity container b
+ InjectorMomentumBoltzmann(GetTemperature const& t, GetVelocity const& b) noexcept
+ : velocity(b), temperature(t)
{}
AMREX_GPU_HOST_DEVICE
@@ -156,7 +157,14 @@ struct InjectorMomentumBoltzmann
if (theta < 0) {
amrex::Abort("Negative temperature parameter theta encountered, which is not allowed");
}
+ // Calculate local velocity and abort if |beta|>=1
+ amrex::Real const beta = velocity(x,y,z);
+ if (beta <= -1 || beta >= 1) {
+ amrex::Abort("beta = v/c magnitude greater than or equal to 1");
+ }
+ // Calculate the value of vave from the local temperature
amrex::Real const vave = std::sqrt(amrex::Real(2) * theta);
+ int const dir = velocity.direction();
amrex::Real x1, x2, gamma;
amrex::Real u[3];
x1 = amrex::Random(engine);
@@ -197,19 +205,20 @@ struct InjectorMomentumBoltzmann
AMREX_GPU_HOST_DEVICE
amrex::XDim3
- getBulkMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/) const noexcept
+ getBulkMomentum (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
{
using namespace amrex;
Real u[3];
for (int idim = 0; idim < 3; ++idim) u[idim] = 0.0_rt;
+ const Real beta = velocity(x,y,z);
+ int const dir = velocity.direction();
const Real gamma = static_cast<amrex::Real>(1./sqrt(1+beta*beta));
u[dir] = gamma*beta;
return XDim3 {u[0],u[1],u[2]};
}
private:
- int dir;
- amrex::Real beta;
+ GetVelocity velocity;
GetTemperature temperature;
};
@@ -219,11 +228,10 @@ private:
struct InjectorMomentumJuttner
{
// Constructor whose inputs are:
- // the temperature parameter theta,
- // boost velocity/c beta,
- // and boost direction dir respectively.
- InjectorMomentumJuttner(GetTemperature const& t, amrex::Real b, int d) noexcept
- : dir(d), beta(b), temperature(t)
+ // a reference to the initial temperature container t,
+ // a reference to the initial velocity container b
+ InjectorMomentumJuttner(GetTemperature const& t, GetVelocity const& b) noexcept
+ : velocity(b), temperature(t)
{}
AMREX_GPU_HOST_DEVICE
@@ -241,6 +249,12 @@ struct InjectorMomentumJuttner
if (theta < 0.1) {
amrex::Abort("Temeprature parameter theta is less than minimum 0.1 allowed for Maxwell-Juttner");
}
+ // Calculate local velocity and abort if |beta|>=1
+ amrex::Real const beta = velocity(x,y,z);
+ if (beta <= -1 || beta >= 1) {
+ amrex::Abort("beta = v/c magnitude greater than or equal to 1");
+ }
+ int const dir = velocity.direction();
x1 = static_cast<amrex::Real>(0.);
gamma = static_cast<amrex::Real>(0.);
u[dir] = static_cast<amrex::Real>(0.);
@@ -291,19 +305,20 @@ struct InjectorMomentumJuttner
AMREX_GPU_HOST_DEVICE
amrex::XDim3
- getBulkMomentum (amrex::Real /*x*/, amrex::Real /*y*/, amrex::Real /*z*/) const noexcept
+ getBulkMomentum (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept
{
using namespace amrex;
Real u[3];
for (int idim = 0; idim < 3; ++idim) u[idim] = 0.0_rt;
+ Real const beta = velocity(x,y,z);
+ int const dir = velocity.direction();
const Real gamma = static_cast<Real>(1./sqrt(1.+beta*beta));
u[dir] = gamma*beta;
return XDim3 {u[0],u[1],u[2]};
}
private:
- int dir;
- amrex::Real beta;
+ GetVelocity velocity;
GetTemperature temperature;
};
@@ -412,16 +427,16 @@ struct InjectorMomentum
{ }
InjectorMomentum (InjectorMomentumBoltzmann* t,
- GetTemperature const& temperature, amrex::Real beta, int dir)
+ GetTemperature const& temperature, GetVelocity const& velocity)
: type(Type::boltzmann),
- object(t, temperature, beta, dir)
+ object(t, temperature, velocity)
{ }
// This constructor stores a InjectorMomentumJuttner in union object.
InjectorMomentum (InjectorMomentumJuttner* t,
- GetTemperature const& temperature, amrex::Real beta, int dir)
+ GetTemperature const& temperature, GetVelocity const& velocity)
: type(Type::juttner),
- object(t, temperature, beta, dir)
+ object(t, temperature, velocity)
{ }
// This constructor stores a InjectorMomentumCustom in union object.
@@ -570,11 +585,11 @@ private:
int a_flux_normal_axis, int a_flux_direction) noexcept
: gaussianflux(a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th,a_flux_normal_axis,a_flux_direction) {}
Object (InjectorMomentumBoltzmann*,
- GetTemperature const& t, amrex::Real b, int dir) noexcept
- : boltzmann(t,b,dir) {}
+ GetTemperature const& t, GetVelocity const& b) noexcept
+ : boltzmann(t,b) {}
Object (InjectorMomentumJuttner*,
- GetTemperature const& t, amrex::Real b, int dir) noexcept
- : juttner(t,b,dir) {}
+ GetTemperature const& t, GetVelocity const& b) noexcept
+ : juttner(t,b) {}
Object (InjectorMomentumRadialExpansion*,
amrex::Real u_over_r) noexcept
: radial_expansion(u_over_r) {}
diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package
index 73075d858..13eed4963 100644
--- a/Source/Initialization/Make.package
+++ b/Source/Initialization/Make.package
@@ -5,5 +5,7 @@ CEXE_sources += InjectorDensity.cpp
CEXE_sources += InjectorMomentum.cpp
CEXE_sources += TemperatureProperties.cpp
CEXE_sources += GetTemperature.cpp
+CEXE_sources += VelocityProperties.cpp
+CEXE_sources += GetVelocity.cpp
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Initialization
diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H
index 3139444ac..b3b7fa4fd 100644
--- a/Source/Initialization/PlasmaInjector.H
+++ b/Source/Initialization/PlasmaInjector.H
@@ -12,6 +12,7 @@
#include "InjectorDensity.H"
#include "InjectorMomentum.H"
#include "TemperatureProperties.H"
+#include "VelocityProperties.H"
#include "Particles/SpeciesPhysicalProperties.H"
#include "InjectorPosition_fwd.H"
@@ -153,6 +154,7 @@ protected:
// Keep a pointer to TemperatureProperties to ensure the lifetime of the
// contained Parser
std::unique_ptr<TemperatureProperties> h_mom_temp;
+ std::unique_ptr<VelocityProperties> h_mom_vel;
void parseDensity (amrex::ParmParse& pp);
void parseMomentum (amrex::ParmParse& pp);
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index fab5d3d58..e6f33ccd5 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -10,6 +10,7 @@
#include "PlasmaInjector.H"
#include "Initialization/GetTemperature.H"
+#include "Initialization/GetVelocity.H"
#include "Initialization/InjectorDensity.H"
#include "Initialization/InjectorMomentum.H"
#include "Initialization/InjectorPosition.H"
@@ -555,67 +556,19 @@ void PlasmaInjector::parseMomentum (amrex::ParmParse& pp)
ux_m, uy_m, uz_m, ux_th, uy_th, uz_th,
flux_normal_axis, flux_direction));
} else if (mom_dist_s == "maxwell_boltzmann"){
- amrex::Real beta = 0._rt;
- int dir = 0;
- std::string direction = "x";
h_mom_temp = std::make_unique<TemperatureProperties>(pp);
GetTemperature getTemp(*h_mom_temp.get());
- queryWithParser(pp, "beta", beta);
- if(beta < 0){
- amrex::Abort("Please enter a positive beta value. Drift direction is set with <s_name>.bulk_vel_dir = 'x' or '+x', '-x', 'y' or '+y', etc.");
- }
- pp.query("bulk_vel_dir", direction);
- if(direction[0] == '-'){
- beta = -beta;
- }
- if((direction == "x" || direction[1] == 'x') ||
- (direction == "X" || direction[1] == 'X')){
- dir = 0;
- } else if ((direction == "y" || direction[1] == 'y') ||
- (direction == "Y" || direction[1] == 'Y')){
- dir = 1;
- } else if ((direction == "z" || direction[1] == 'z') ||
- (direction == "Z" || direction[1] == 'Z')){
- dir = 2;
- } else{
- std::stringstream stringstream;
- stringstream << "Cannot interpret <s_name>.bulk_vel_dir input '" << direction << "'. Please enter +/- x, y, or z with no whitespace between the sign and other character.";
- direction = stringstream.str();
- amrex::Abort(direction.c_str());
- }
+ h_mom_vel = std::make_unique<VelocityProperties>(pp);
+ GetVelocity getVel(*h_mom_vel.get());
// Construct InjectorMomentum with InjectorMomentumBoltzmann.
- h_inj_mom.reset(new InjectorMomentum((InjectorMomentumBoltzmann*)nullptr, getTemp, beta, dir));
+ h_inj_mom.reset(new InjectorMomentum((InjectorMomentumBoltzmann*)nullptr, getTemp, getVel));
} else if (mom_dist_s == "maxwell_juttner"){
- amrex::Real beta = 0._rt;
- int dir = 0;
h_mom_temp = std::make_unique<TemperatureProperties>(pp);
GetTemperature getTemp(*h_mom_temp.get());
- std::string direction = "x";
- queryWithParser(pp, "beta", beta);
- if(beta < 0){
- amrex::Abort("Please enter a positive beta value. Drift direction is set with <s_name>.bulk_vel_dir = 'x' or '+x', '-x', 'y' or '+y', etc.");
- }
- pp.query("bulk_vel_dir", direction);
- if(direction[0] == '-'){
- beta = -beta;
- }
- if((direction == "x" || direction[1] == 'x') ||
- (direction == "X" || direction[1] == 'X')){
- dir = 0;
- } else if ((direction == "y" || direction[1] == 'y') ||
- (direction == "Y" || direction[1] == 'Y')){
- dir = 1;
- } else if ((direction == "z" || direction[1] == 'z') ||
- (direction == "Z" || direction[1] == 'Z')){
- dir = 2;
- } else{
- std::stringstream stringstream;
- stringstream << "Cannot interpret <s_name>.bulk_vel_dir input '" << direction << "'. Please enter +/- x, y, or z with no whitespace between the sign and other character.";
- direction = stringstream.str();
- amrex::Abort(direction.c_str());
- }
+ h_mom_vel = std::make_unique<VelocityProperties>(pp);
+ GetVelocity getVel(*h_mom_vel.get());
// Construct InjectorMomentum with InjectorMomentumJuttner.
- h_inj_mom.reset(new InjectorMomentum((InjectorMomentumJuttner*)nullptr, getTemp, beta, dir));
+ h_inj_mom.reset(new InjectorMomentum((InjectorMomentumJuttner*)nullptr, getTemp, getVel));
} else if (mom_dist_s == "radial_expansion") {
amrex::Real u_over_r = 0._rt;
queryWithParser(pp, "u_over_r", u_over_r);
diff --git a/Source/Initialization/VelocityProperties.H b/Source/Initialization/VelocityProperties.H
new file mode 100644
index 000000000..e9e6362bf
--- /dev/null
+++ b/Source/Initialization/VelocityProperties.H
@@ -0,0 +1,54 @@
+/* Copyright 2021 Hannah Klion
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef VELOCITY_PROPERTIES_H_
+#define VELOCITY_PROPERTIES_H_
+
+#include <Utils/WarpXUtil.H>
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_Parser.H>
+#include <AMReX_REAL.H>
+
+/* Type of velocity initialization. Used by VelocityProperties and GetVelocity. */
+enum VelocityInitType {VelConstantValue, VelParserFunction};
+
+/**
+ * \brief Struct to store velocity properties, for use in momentum initialization.
+ *
+ * Reads in and stores velocity used to initialize the Maxwell-Boltzmann and Maxwell-Juttner
+ * momentum distributions in InjectorMomentum. The information is read from the parameters of
+ * the species being initialized, and will be accessed by GetVelocity.
+ */
+struct VelocityProperties
+{
+ /**
+ * \brief Read runtime parameters to populate constant or spatially-varying velocity
+ * information
+ *
+ * Construct VelocityProperties based on the passed parameters.
+ * If velocity is a constant, store value. If a parser, make and
+ * store the parser function
+ *
+ * \param[in] pp: Reference to the parameter parser object for the species being initialized
+ */
+ VelocityProperties (amrex::ParmParse& pp);
+
+ /* Type of velocity initialization */
+ VelocityInitType m_type;
+
+ /* Velocity direction */
+ int m_dir; // Index x=0, y=1, z=2
+ int m_sign_dir; // Sign of the velocity direction positive=1, negative=-1
+
+ /* Constant velocity value, if m_type == VelConstantValue */
+ amrex::Real m_velocity;
+ /* Storage of the parser function, if m_type == VelParserFunction */
+ std::unique_ptr<amrex::Parser> m_ptr_velocity_parser;
+};
+
+#endif
diff --git a/Source/Initialization/VelocityProperties.cpp b/Source/Initialization/VelocityProperties.cpp
new file mode 100644
index 000000000..6a0e74ecf
--- /dev/null
+++ b/Source/Initialization/VelocityProperties.cpp
@@ -0,0 +1,71 @@
+/* Copyright 2021 Hannah Klion
+ *
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "VelocityProperties.H"
+
+VelocityProperties::VelocityProperties (amrex::ParmParse& pp) {
+ // Set defaults
+ std::string vel_dist_s = "constant";
+ std::string vel_dir_s = "x";
+ m_velocity = 0;
+
+ pp.query("bulk_vel_dir", vel_dir_s);
+ if(vel_dir_s[0] == '-'){
+ m_sign_dir = -1;
+ }
+ else {
+ m_sign_dir = 1;
+ }
+
+ if ((vel_dir_s == "x" || vel_dir_s[1] == 'x') ||
+ (vel_dir_s == "X" || vel_dir_s[1] == 'X')){
+ m_dir = 0;
+ }
+ else if ((vel_dir_s == "y" || vel_dir_s[1] == 'y') ||
+ (vel_dir_s == "Y" || vel_dir_s[1] == 'Y')){
+ m_dir = 1;
+ }
+ else if ((vel_dir_s == "z" || vel_dir_s[1] == 'z') ||
+ (vel_dir_s == "Z" || vel_dir_s[1] == 'Z')) {
+ m_dir = 2;
+ }
+ else {
+ std::stringstream stringstream;
+ stringstream << "Cannot interpret <s_name>.bulk_vel_dir input '" << vel_dir_s <<
+ "'. Please enter +/- x, y, or z with no whitespace between the sign and" <<
+ " other character.";
+ vel_dir_s = stringstream.str();
+ amrex::Abort(vel_dir_s.c_str());
+ }
+
+ pp.query("beta_distribution_type", vel_dist_s);
+ if (vel_dist_s == "constant") {
+ queryWithParser(pp, "beta", m_velocity);
+ m_type = VelConstantValue;
+ if (m_velocity >= 1 || m_velocity <= -1) {
+ std::stringstream stringstream;
+ stringstream << "Magnitude of velocity beta = " << m_velocity <<
+ " is greater than or equal to 1";
+ amrex::Abort(stringstream.str().c_str());
+ }
+ }
+ else if (vel_dist_s == "parser") {
+ std::string str_beta_function;
+ Store_parserString(pp, "beta_function(x,y,z)", str_beta_function);
+ m_ptr_velocity_parser =
+ std::make_unique<amrex::Parser>(makeParser(str_beta_function,{"x","y","z"}));
+ m_type = VelParserFunction;
+ }
+ else {
+ std::stringstream stringstream;
+ std::string string;
+ stringstream << "Velocity distribution type '" << vel_dist_s << "' not recognized." << std::endl;
+ string = stringstream.str();
+ amrex::Abort(string.c_str());
+ }
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 04:31:36 +0000