diff options
| -rw-r--r-- | Docs/source/running_cpp/parameters.rst | 9 | ||||
| -rw-r--r-- | Source/Initialization/InjectorMomentum.H | 92 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 22 |
3 files changed, 122 insertions, 1 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index e5727e376..70cf7c199 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -266,6 +266,15 @@ Particle initialization well as standard deviations along each direction ``<species_name>.ux_th``, ``<species_name>.uy_th`` and ``<species_name>.uz_th``. + * ``maxwell_juttner``: Maxwell-Juttner distribution for high temperature plasma. This mode + requires a dimensionless temperature parameter ``<species_name>.theta``, where theta is equal + to kb*T/(m*c^2), where kb is the Boltzmann constant, and m is the mass of the species. It also + includes the optional parameter ``<species_name>.beta`` where beta is equal to v/c. The plasma + will be initialized to move at velocity beta*c in the ``<species_name>.direction = 'x', 'y', 'z'``, + direction. 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. + Both the Sobol and the flipping method can be found in Zenitani 2015 (Phys. Plasmas 22, 042116). + * ``radial_expansion``: momentum depends on the radial coordinate linearly. This requires additional parameter ``u_over_r`` which is the slope. diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H index ba0d26fc4..f5819a908 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -47,6 +47,81 @@ private: amrex::Real m_ux_th, m_uy_th, m_uz_th; }; +// 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). +struct InjectorMomentumJuttner +{ + // Constructor whose inputs are: + // the temperature parameter theta, + // boost velocity/c beta, + // and boost direction dir respectively. + InjectorMomentumJuttner(amrex::Real t, amrex::Real b, int d) noexcept + : theta(t), beta(b), dir(d) + {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + // Sobol method for sampling MJ Speeds, + // from Zenitani 2015 (Phys. Plasmas 22, 042116). + amrex::Real x1, x2, gamma; + amrex::Real u [3]; + x1 = 0.; + gamma = 0.; + u[dir] = 0.; + // This condition is equation 10 in Zenitani, + // though x1 is defined differently. + while(u[dir]-gamma <= x1) + { + u[dir] = -theta* + std::log(amrex::Random()*amrex::Random()*amrex::Random()); + gamma = std::sqrt(1+std::pow(u[dir],2)); + x1 = theta*std::log(amrex::Random()); + } + // The following code samples a random unit vector + // and multiplies the result by speed u[dir]. + x1 = amrex::Random(); + x2 = amrex::Random(); + // Direction dir is an input parameter that sets the boost direction: + // 'x' -> d = 0, 'y' -> d = 1, 'z' -> d = 2. + u[(dir+1)%3] = 2*u[dir]*std::sqrt(x1*(1-x1))*std::sin(2*M_PI*x2); + u[(dir+2)%3] = 2*u[dir]*std::sqrt(x1*(1-x1))*std::cos(2*M_PI*x2); + // The value of dir is the boost direction to be transformed. + u[dir] = u[dir]*(2*x1-1); + x1 = amrex::Random(); + // 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. + 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, theta; +}; + + // struct whose getMomentum returns momentum for 1 particle, for // radial expansion struct InjectorMomentumRadialExpansion @@ -124,6 +199,13 @@ struct InjectorMomentum object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th) { } + // This constructor stores a InjectorMomentumJuttner in union object. + InjectorMomentum (InjectorMomentumJuttner* t, + amrex::Real theta, amrex::Real beta, int dir) + : type(Type::juttner), + object(t, theta, beta, dir) + { } + // This constructor stores a InjectorMomentumCustom in union object. InjectorMomentum (InjectorMomentumCustom* t, std::string const& a_species_name) @@ -165,6 +247,10 @@ struct InjectorMomentum { return object.gaussian.getMomentum(x,y,z); } + case Type::juttner: + { + return object.juttner.getMomentum(x,y,z); + } case Type::constant: { return object.constant.getMomentum(x,y,z); @@ -186,7 +272,7 @@ struct InjectorMomentum } private: - enum struct Type { constant, custom, gaussian, radial_expansion, parser }; + enum struct Type { constant, custom, gaussian, juttner, radial_expansion, parser}; Type type; // An instance of union Object constructs and stores any one of @@ -204,6 +290,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 (InjectorMomentumJuttner*, + amrex::Real t, amrex::Real b, int dir) noexcept + : juttner(t,b,dir) {} Object (InjectorMomentumRadialExpansion*, amrex::Real u_over_r) noexcept : radial_expansion(u_over_r) {} @@ -215,6 +304,7 @@ private: InjectorMomentumConstant constant; InjectorMomentumCustom custom; InjectorMomentumGaussian gaussian; + InjectorMomentumJuttner juttner; InjectorMomentumRadialExpansion radial_expansion; InjectorMomentumParser parser; }; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index af04c6b31..32711cd31 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -269,6 +269,28 @@ void PlasmaInjector::parseMomentum (ParmParse& pp) // 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 == "maxwell_juttner"){ + Real beta = 0.; + Real theta = 10.; + int dir = 0; + std::string direction = "x"; + pp.query("beta", beta); + pp.query("theta", theta); + pp.query("direction", direction); + if(direction == "x" || direction == "X"){ + dir = 0; + } else if (direction == "y" || direction == "Y"){ + dir = 1; + } else if (direction == "z" || direction == "Z"){ + dir = 2; + } else{ + std::stringstream stringstream; + stringstream << "Direction " << direction << " is not recognzied. Please enter x, y, or z."; + direction = stringstream.str(); + amrex::Abort(direction.c_str()); + } + // Construct InjectorMomentum with InjectorMomentumJuttner. + inj_mom.reset(new InjectorMomentum((InjectorMomentumJuttner*)nullptr, theta, beta, dir)); } else if (mom_dist_s == "radial_expansion") { Real u_over_r = 0.; pp.query("u_over_r", u_over_r); |