Spatially vary velocity (#2491)

* initial add of GetVelocity and VelocityProperties

* bug fixes and first test implementation

* add documentation to parameters.rst

* Add error if |beta|>=1

* Check for |beta|>=1 in Juttner as well

* update comment on initial_distribution python script

* Doxygenization and floats in python

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Address comments in review

- Add units to analysis_distribution.py constants
- Move a comment on VelocityProperties constructor from implementation to header file

Co-authored-by: Hannah Klion <hannah.klion@gmail.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

13 files changed, 504 insertions, 103 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 03b5dd805..cc28381cf 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -672,9 +672,11 @@ Particle initialization
       Theta is specified by a combination of ``<species_name>.theta_distribution_type``, ``<species_name>.theta``, and ``<species_name>.theta_function(x,y,z)`` (see below).
       For values of :math:`\theta > 0.01`, errors due to ignored relativistic terms exceed 1%.
       Temperatures less than zero are not allowed.
-      It also includes the optional parameter ``<species_name>.beta`` where beta is equal to v/c.
-      The plasma will be initialized to move at bulk velocity beta*c in the
-      ``<species_name>.bulk_vel_dir = (+/-) 'x', 'y', 'z'`` direction. Please leave no whitespace
+      The plasma can be initialized to move at a bulk velocity :math:`\beta = v/c`.
+      The speed is specified by the parameters ``<species_name>.beta_distribution_type``, ``<species_name>.beta``, and ``<species_name>.beta_function(x,y,z)`` (see below).
+      :math:`\beta` can be positive or negative and is limited to the range :math:`-1 < \beta < 1`.
+      The direction of the velocity field is given by ``<species_name>.bulk_vel_dir = (+/-) 'x', 'y', 'z'``, and must be the same across the domain.
+      Please leave no whitespace
       between the sign and the character on input. A direction without a sign will be treated as
       positive. The MB distribution is initialized in the drifting frame by sampling three Gaussian
       distributions in each dimension using, the Box Mueller method, and then the distribution is
@@ -692,10 +694,11 @@ Particle initialization
       The Sobol method used to generate the distribution will not terminate for :math:`\theta \lesssim 0.1`, and the code will abort if it encounters a temperature below that threshold.
       The Maxwell-Boltzmann distribution is recommended for temperatures in the range :math:`0.01 < \theta < 0.1`.
       Errors due to relativistic effects can be expected to approximately between 1% and 10%.
-      It also
-      includes the optional parameter ``<species_name>.beta`` where beta is equal to :math:`\frac{v}{c}`. The plasma
-      will be initialized to move at velocity :math:`\beta \cdot c` in the
-      ``<species_name>.bulk_vel_dir = (+/-) 'x', 'y', 'z'`` direction. Please leave no whitespace
+      The plasma can be initialized to move at a bulk velocity :math:`\beta = v/c`.
+      The speed is specified by the parameters ``<species_name>.beta_distribution_type``, ``<species_name>.beta``, and ``<species_name>.beta_function(x,y,z)`` (see below).
+      :math:`\beta` can be positive or negative and is limited to the range :math:`-1 < \beta < 1`.
+      The direction of the velocity field is given by ``<species_name>.bulk_vel_dir = (+/-) 'x', 'y', 'z'``, and must be the same across the domain.
+      Please leave no whitespace
       between the sign and the character on input. A direction without a sign will be treated as
       positive. The MJ distribution will be initialized in the moving frame using the Sobol method,
       and then the distribution will be transformed to the simulation frame using the flipping method.
@@ -724,6 +727,14 @@ Particle initialization
 
     * If ``parser``, use a spatially-dependent analytic parser function, given by the required parameter ``<species_name>.theta_function(x,y,z)``.
 
+* ``<species_name>.beta_distribution_type`` (`string`) optional (default ``constant``)
+    Only read if ``<species_name>.momentum_distribution_type`` is ``maxwell_boltzmann`` or ``maxwell_juttner``.
+    See documentation for these distributions (above) for constraints on values of beta.
+
+    * If ``constant``, use a constant speed, given by the required float parameter ``<species_name>.beta``.
+
+    * If ``parser``, use a spatially-dependent analytic parser function, given by the required parameter ``<species_name>.beta_function(x,y,z)``.
+
 * ``<species_name>.zinject_plane`` (`float`)
     Only read if  ``<species_name>`` is in ``particles.rigid_injected_species``.
     Injection plane when using the rigid injection method.
diff --git a/Examples/Tests/initial_distribution/analysis_distribution.py b/Examples/Tests/initial_distribution/analysis_distribution.py
index 3da4cf423..d43ced2bb 100755
--- a/Examples/Tests/initial_distribution/analysis_distribution.py
+++ b/Examples/Tests/initial_distribution/analysis_distribution.py
@@ -12,6 +12,8 @@
 # 3 denotes maxwell-juttner distribution.
 # 4 denotes gaussian position distribution.
 # 5 denotes maxwell-juttner distribution w/ spatially varying temperature
+# 6 denotes maxwell-boltzmann distribution w/ constant velocity
+# 7 denotes maxwell-boltzmann distribution w/ spatially-varying velocity
 # The distribution is obtained through reduced diagnostic ParticleHistogram.
 
 import numpy as np
@@ -180,5 +182,80 @@ print('Maxwell-Juttner parser temperature difference:', f5_error)
 
 assert(f5_error < tolerance)
 
+#==============================================
+# maxwell-boltzmann with constant bulk velocity
+#==============================================
+
+# load data
+bin_value_g, bin_data_g = read_reduced_diags_histogram("h6.txt")[2:]
+bin_value_uy, bin_data_uy = read_reduced_diags_histogram("h6uy.txt")[2:]
+
+# Expected values for beta and u = beta*gamma
+beta_const = 0.2
+g_const = 1. / np.sqrt(1. - beta_const * beta_const)
+uy_const = beta_const * g_const
+g_bin_size = 0.004
+g_bin_min = 1.
+uy_bin_size = 0.04
+uy_bin_min = -1.
+V = 8.0 # volume in m^3
+n = 1.0e21 # number density in 1/m^3
+
+f_g = np.zeros_like(bin_value_g)
+i_g = int(np.floor((g_const - g_bin_min) / g_bin_size))
+f_g[i_g] = n * V
+f_peak = np.amax(f_g)
+
+f_uy = np.zeros_like(bin_value_uy)
+i_uy = int(np.floor((-uy_const - uy_bin_min) / uy_bin_size))
+f_uy[i_uy] = n * V
+
+f6_error = np.sum(np.abs(f_g - bin_data_g) + np.abs(f_uy - bin_data_uy)) \
+        / bin_value_g.size / f_peak
+
+print('Maxwell-Boltzmann constant velocity difference:', f6_error)
+
+assert(f6_error < tolerance)
+
+#============================================
+# maxwell-boltzmann with parser bulk velocity
+#============================================
+
+# load data
+bin_value_g, bin_data_g = read_reduced_diags_histogram("h7.txt")[2:]
+bin_value_uy, bin_data_uy_neg = read_reduced_diags_histogram("h7uy_neg.txt")[2:]
+bin_data_uy_pos = read_reduced_diags_histogram("h7uy_pos.txt")[3]
+
+# Expected values for beta and u = beta*gamma
+beta_const = 0.2
+g_const = 1. / np.sqrt(1. - beta_const * beta_const)
+uy_const = beta_const * g_const
+g_bin_size = 0.004
+g_bin_min = 1.
+uy_bin_size = 0.04
+uy_bin_min = -1.
+V = 8.0 # volume in m^3
+n = 1.0e21 # number density in 1/m^3
+
+f_g = np.zeros_like(bin_value_g)
+i_g = int(np.floor((g_const - g_bin_min) / g_bin_size))
+f_g[i_g] = n * V
+f_peak = np.amax(f_g)
+
+f_uy_neg = np.zeros_like(bin_value_uy)
+i_uy_neg = int(np.floor((uy_const - uy_bin_min) / uy_bin_size))
+f_uy_neg[i_uy_neg] = n * V / 2.
+
+f_uy_pos = np.zeros_like(bin_value_uy)
+i_uy_pos = int(np.floor((-uy_const - uy_bin_min) / uy_bin_size))
+f_uy_pos[i_uy_pos] = n * V / 2.
+
+f7_error = np.sum(np.abs(f_g - bin_data_g) + np.abs(f_uy_pos - bin_data_uy_pos) \
+        + np.abs(f_uy_neg - bin_data_uy_neg)) / bin_value_g.size / f_peak
+
+print('Maxwell-Boltzmann parser velocity difference:', f7_error)
+
+assert(f7_error < tolerance)
+
 test_name = filename[:-9] # Could also be os.path.split(os.getcwd())[1]
 checksumAPI.evaluate_checksum(test_name, filename)
diff --git a/Examples/Tests/initial_distribution/inputs b/Examples/Tests/initial_distribution/inputs
index 0b47249df..d301c5751 100644
--- a/Examples/Tests/initial_distribution/inputs
+++ b/Examples/Tests/initial_distribution/inputs
@@ -29,7 +29,7 @@ algo.particle_shape = 1
 #################################
 ############ PLASMA #############
 #################################
-particles.species_names = gaussian maxwell_boltzmann maxwell_juttner beam maxwell_juttner_parser
+particles.species_names = gaussian maxwell_boltzmann maxwell_juttner beam maxwell_juttner_parser velocity_constant velocity_parser
 particles.rigid_injected_species = beam
 
 gaussian.charge                     = -q_e
@@ -93,6 +93,32 @@ maxwell_juttner_parser.momentum_distribution_type = "maxwell_juttner"
 maxwell_juttner_parser.theta_distribution_type    = "parser"
 maxwell_juttner_parser.theta_function(x,y,z)      = "1.0 + heaviside(x,0)"
 
+velocity_constant.charge                     = -q_e
+velocity_constant.mass                       = m_e
+velocity_constant.injection_style            = "NRandomPerCell"
+velocity_constant.num_particles_per_cell     = 1000
+velocity_constant.profile                    = constant
+velocity_constant.density                    = 1.0e21
+velocity_constant.momentum_distribution_type = "maxwell_boltzmann"
+velocity_constant.theta_distribution_type    = "constant"
+velocity_constant.theta                      = 1e-9
+velocity_constant.beta_distribution_type     = "constant"
+velocity_constant.beta                       = 0.2
+velocity_constant.bulk_vel_dir               = -y
+
+velocity_parser.charge                     = -q_e
+velocity_parser.mass                       = m_e
+velocity_parser.injection_style            = "NRandomPerCell"
+velocity_parser.num_particles_per_cell     = 1000
+velocity_parser.profile                    = constant
+velocity_parser.density                    = 1.0e21
+velocity_parser.momentum_distribution_type = "maxwell_boltzmann"
+velocity_parser.theta_distribution_type    = "constant"
+velocity_parser.theta                      = 1e-9
+velocity_parser.beta_distribution_type     = "parser"
+velocity_parser.beta_function(x,y,z)       = "-0.2 + 0.4 * heaviside(z,0)"
+velocity_parser.bulk_vel_dir               = -y
+
 #################################
 ########## DIAGNOSTIC ###########
 #################################
@@ -101,7 +127,9 @@ maxwell_juttner_parser.theta_function(x,y,z)      = "1.0 + heaviside(x,0)"
 # 3 for maxwell-juttner with constant temperature
 # 4 for beam
 # 5 for maxwell-juttner with parser function temperature
-warpx.reduced_diags_names              = h1x h1y h1z h2x h2y h2z h3 h3_filtered h4x h4y h4z bmmntr h5_neg h5_pos
+# 6 for maxwell-boltzmann with constant velocity
+# 7 for maxwell-boltzmann with parser velocity
+warpx.reduced_diags_names              = h1x h1y h1z h2x h2y h2z h3 h3_filtered h4x h4y h4z bmmntr h5_neg h5_pos h6 h6uy h7 h7uy_pos h7uy_neg
 
 h1x.type                                 = ParticleHistogram
 h1x.intervals                            = 1
@@ -223,6 +251,62 @@ h5_pos.bin_max                              = 12.0
 h5_pos.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)"
 h5_pos.filter_function(t,x,y,z,ux,uy,uz)    = "x > 0.0"
 
+h5_pos.type                                 = ParticleHistogram
+h5_pos.intervals                            = 1
+h5_pos.path                                 = "./"
+h5_pos.species                              = maxwell_juttner_parser
+h5_pos.bin_number                           = 50
+h5_pos.bin_min                              = 1.0
+h5_pos.bin_max                              = 12.0
+h5_pos.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)"
+h5_pos.filter_function(t,x,y,z,ux,uy,uz)    = "x > 0.0"
+
+h6.type                                   = ParticleHistogram
+h6.intervals                              = 1
+h6.path                                   = "./"
+h6.species                                = velocity_constant
+h6.bin_number                             = 50
+h6.bin_min                                = 1.0
+h6.bin_max                                = 1.2
+h6.histogram_function(t,x,y,z,ux,uy,uz)   = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)"
+
+h6uy.type                                   = ParticleHistogram
+h6uy.intervals                              = 1
+h6uy.path                                   = "./"
+h6uy.species                                = velocity_constant
+h6uy.bin_number                             = 50
+h6uy.bin_min                                = -1.0
+h6uy.bin_max                                = 1.0
+h6uy.histogram_function(t,x,y,z,ux,uy,uz)   = "uy"
+
+h7.type                                   = ParticleHistogram
+h7.intervals                              = 1
+h7.path                                   = "./"
+h7.species                                = velocity_parser
+h7.bin_number                             = 50
+h7.bin_min                                = 1.0
+h7.bin_max                                = 1.2
+h7.histogram_function(t,x,y,z,ux,uy,uz)   = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)"
+
+h7uy_pos.type                                 = ParticleHistogram
+h7uy_pos.intervals                            = 1
+h7uy_pos.path                                 = "./"
+h7uy_pos.species                              = velocity_parser
+h7uy_pos.bin_number                           = 50
+h7uy_pos.bin_min                              = -1.0
+h7uy_pos.bin_max                              = 1.0
+h7uy_pos.histogram_function(t,x,y,z,ux,uy,uz) = "uy"
+h7uy_pos.filter_function(t,x,y,z,ux,uy,uz)    = "z > 0.0"
+
+h7uy_neg.type                                 = ParticleHistogram
+h7uy_neg.intervals                            = 1
+h7uy_neg.path                                 = "./"
+h7uy_neg.species                              = velocity_parser
+h7uy_neg.bin_number                           = 50
+h7uy_neg.bin_min                              = -1.0
+h7uy_neg.bin_max                              = 1.0
+h7uy_neg.histogram_function(t,x,y,z,ux,uy,uz) = "uy"
+h7uy_neg.filter_function(t,x,y,z,ux,uy,uz)    = "z < 0.0"
 
 # our little beam monitor
 bmmntr.type                              = BeamRelevant
diff --git a/Regression/Checksum/benchmarks_json/initial_distribution.json b/Regression/Checksum/benchmarks_json/initial_distribution.json
index 1d3004fe6..90deea3ff 100644
--- a/Regression/Checksum/benchmarks_json/initial_distribution.json
+++ b/Regression/Checksum/benchmarks_json/initial_distribution.json
@@ -2,8 +2,8 @@
   "beam": {
     "particle_cpu": 0.0,
     "particle_id": 869553169926.0,
-    "particle_momentum_x": 1.0642030918892651e-16,
-    "particle_momentum_y": 1.0642283537845086e-16,
+    "particle_momentum_x": 1.0642030918892653e-16,
+    "particle_momentum_y": 1.0642283537844856e-16,
     "particle_momentum_z": 1.0642015613690699e-16,
     "particle_position_x": 97548.74915484959,
     "particle_position_y": 97419.42550664608,
@@ -14,7 +14,7 @@
     "particle_cpu": 0.0,
     "particle_id": 131072256000.0,
     "particle_momentum_x": 1.1148694854418736e-18,
-    "particle_momentum_y": 1.1139504847071663e-18,
+    "particle_momentum_y": 1.1139504847073219e-18,
     "particle_momentum_z": 1.11450954619821e-18,
     "particle_position_x": 256080.8791703085,
     "particle_position_y": 255990.74838430836,
@@ -22,21 +22,21 @@
     "particle_weight": 8e+21
   },
   "lev=0": {
-    "Bx": 2.295513267966755e-12,
-    "By": 2.2621088131755748e-12,
-    "Bz": 2.3259320950698525e-12,
-    "Ex": 103489.99666142697,
-    "Ey": 102764.83503312292,
-    "Ez": 109758.37177938089,
-    "jx": 190321750907.71967,
-    "jy": 188988153118.15323,
-    "jz": 201849513650.8026
+    "Bx": 1.4041930126107817e-11,
+    "By": 2.2620996239033253e-12,
+    "Bz": 2.345593514787453e-12,
+    "Ex": 103489.94137386356,
+    "Ey": 2681711.7522189138,
+    "Ez": 109758.3778080657,
+    "jx": 190321649231.9422,
+    "jy": 4931762417404.209,
+    "jz": 201849524737.76642
   },
   "maxwell_boltzmann": {
     "particle_cpu": 0.0,
     "particle_id": 393216256000.0,
-    "particle_momentum_x": 1.117492093974982e-18,
-    "particle_momentum_y": 1.1136299905181956e-18,
+    "particle_momentum_x": 1.1174920939749821e-18,
+    "particle_momentum_y": 1.1136299905175387e-18,
     "particle_momentum_z": 1.114502805331281e-18,
     "particle_position_x": 255953.32756914443,
     "particle_position_y": 256048.5390456779,
@@ -47,7 +47,7 @@
     "particle_cpu": 0.0,
     "particle_id": 655360256000.0,
     "particle_momentum_x": 2.21325802726127e-16,
-    "particle_momentum_y": 2.2167623692374665e-16,
+    "particle_momentum_y": 2.2167623692374692e-16,
     "particle_momentum_z": 2.2153818320139063e-16,
     "particle_position_x": 255924.63477257418,
     "particle_position_y": 255995.2040137927,
@@ -58,11 +58,33 @@
     "particle_cpu": 0.0,
     "particle_id": 1167591168000.0,
     "particle_momentum_x": 3.239944229468624e-16,
-    "particle_momentum_y": 3.2438215468627077e-16,
+    "particle_momentum_y": 3.24382154686272e-16,
     "particle_momentum_z": 3.240434151669339e-16,
     "particle_position_x": 256015.2596849932,
     "particle_position_y": 255991.72140724873,
     "particle_position_z": 255927.88992890544,
     "particle_weight": 8e+21
+  },
+  "velocity_constant": {
+    "particle_cpu": 0.0,
+    "particle_id": 1429735168000.0,
+    "particle_momentum_x": 3.5298815518958025e-21,
+    "particle_momentum_y": 2.8541316836401376e-17,
+    "particle_momentum_z": 3.529131393024785e-21,
+    "particle_position_x": 256013.66100185944,
+    "particle_position_y": 255987.4495912913,
+    "particle_position_z": 255995.3840195281,
+    "particle_weight": 8e+21
+  },
+  "velocity_parser": {
+    "particle_cpu": 0.0,
+    "particle_id": 1691879168000.0,
+    "particle_momentum_x": 3.52626487830545e-21,
+    "particle_momentum_y": 2.85413258278335e-17,
+    "particle_momentum_z": 3.531621710284665e-21,
+    "particle_position_x": 256049.20358254094,
+    "particle_position_y": 255976.85796854852,
+    "particle_position_z": 256083.54377426513,
+    "particle_weight": 8e+21
   }
 }
\ No newline at end of file
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());
+    }
+}