1 files changed, 71 insertions, 1 deletions

diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H
index f5819a908..5bbad36c3 100644
--- a/Source/Initialization/InjectorMomentum.H
+++ b/Source/Initialization/InjectorMomentum.H
@@ -47,6 +47,62 @@ private:
     amrex::Real m_ux_th, m_uy_th, m_uz_th;
 };
 
+// struct whose getMomentum returns momentum for 1 particle with relativistic
+// drift velocity beta, from the Maxwell-Boltzmann distribution.
+struct InjectorMomentumBoltzmann
+{
+    InjectorMomentumBoltzmann(amrex::Real t, amrex::Real b, int d) noexcept
+        : vave(std::sqrt(2*t)), beta(b), dir(d)
+        {}
+
+    AMREX_GPU_HOST_DEVICE
+    amrex::XDim3
+    getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
+    {
+        amrex::Real x1, x2, gamma;
+        amrex::Real u[3];
+        vave = std::sqrt(2*t);
+        x1 = amrex::Random();
+        x2 = amrex::Random();
+        // Each value of sqrt(-log(x1))*sin(2*pi*x2) is a sample from a Gaussian
+        // distribution with sigma = vave.
+        u[(dir+1)%3] = vave* std::sqrt(-std::log(x1)) *std::sin(2*M_PI*x2);
+        u[(dir+2)%3] = vave* std::sqrt(-std::log(x1)) *std::cos(2*M_PI*x2);
+        u[dir] = vave*std::sqrt(-std::log(amrex::Random()))*
+            std::sin(2*M_PI*amrex::Random());
+        gamma = std::sqrt(std::pow(u[0],2)+std::pow(u[1],2)+std::pow(u[2],2));
+        gamma = np.sqrt(1+gamma**2);
+        // The following condition is equtaion 32 in Zenitani, called
+        // The flipping method. It transforms the intergral: d3x' -> d3x
+        // where d3x' is the volume element for positions in the boosted frame.
+        // The particle positions and densities can be initialized in the
+        // simulation frame with this method.
+        // The flipping method can similarly transform any
+        // symmetric distribution from one reference frame to another moving at
+        // a relative velocity of beta.
+        // An equivalent alternative to this method native to WarpX
+        // would be to initialize the particle positions and densities in the
+        // frame moving at speed beta, and then perform a Lorentz transform
+        // on their positions and MJ sampled velocities to the simulation frame.
+        x1 = amrex::Random();
+        if(-beta*u[dir]/gamma > x1)
+        {
+          u[dir] = -u[dir];
+        }
+        // This Lorentz transform is equation 17 in Zenitani.
+        // It transforms the integral d3u' -> d3u
+        // where d3u' is the volume element for momentum in the boosted frame.
+        u[dir] = 1/std::sqrt(1-pow(beta,2))*(u[dir]+gamma*beta);
+        // Note that if beta = 0 then the flipping method and Lorentz transform
+        // have no effect on the u[dir] direction.
+        return amrex::XDim3 {u[0],u[1],u[2]};
+    }
+
+private:
+    int dir;
+    amrex::Real beta, vave;
+};
+
 // struct whose getMomentum returns momentum for 1 particle with relativistc
 // drift velocity beta, from the Maxwell-Juttner distribution. Method is from
 // Zenitani 2015 (Phys. Plasmas 22, 042116).
@@ -198,6 +254,12 @@ struct InjectorMomentum
         : type(Type::gaussian),
           object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th)
     { }
+    
+    InjectorMomentum (InjectorMomentumBoltzmann* t,
+                       amrex::Real theta, amrex::Real beta, int dir)
+         : type(Type::boltzmann),
+           object(t, theta, beta, dir)
+    { }
 
      // This constructor stores a InjectorMomentumJuttner in union object.
      InjectorMomentum (InjectorMomentumJuttner* t,
@@ -247,6 +309,10 @@ struct InjectorMomentum
         {
             return object.gaussian.getMomentum(x,y,z);
         }
+        case Type::boltzmann:
+        {
+            return object.boltzmann.getMomentum(x,y,z);
+        }
         case Type::juttner:
         {
             return object.juttner.getMomentum(x,y,z);
@@ -272,7 +338,7 @@ struct InjectorMomentum
     }
 
 private:
-    enum struct Type { constant, custom, gaussian, juttner, radial_expansion, parser};
+    enum struct Type { constant, custom, gaussian, boltzmann, juttner, radial_expansion, parser};
     Type type;
 
     // An instance of union Object constructs and stores any one of
@@ -290,6 +356,9 @@ private:
                 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 (InjectorMomentumBoltzmann*,
+                amrex::Real t, amrex::Real b, int dir) noexcept
+            : boltzmann(t,b,dir) {}
         Object (InjectorMomentumJuttner*,
                 amrex::Real t, amrex::Real b, int dir) noexcept
             : juttner(t,b,dir) {}
@@ -304,6 +373,7 @@ private:
         InjectorMomentumConstant constant;
         InjectorMomentumCustom   custom;
         InjectorMomentumGaussian gaussian;
+        InjectorMomentumBoltzmann boltzmann;
         InjectorMomentumJuttner juttner;
         InjectorMomentumRadialExpansion radial_expansion;
         InjectorMomentumParser   parser;