diff options
| -rw-r--r-- | Docs/source/running_cpp/parameters.rst | 9 | ||||
| -rw-r--r-- | Source/Initialization/InjectorMomentum.H | 62 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 22 |
3 files changed, 92 insertions, 1 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index e5727e376..34d469988 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 method and the flipping method can be found in Zenitani 2015. + * ``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..9b92dea74 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -47,6 +47,51 @@ 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. +struct InjectorMomentumJuttner +{ + InjectorMomentumJuttner(amrex::Real theta, amrex::Real beta, int dir) noexcept + : t(theta), b(beta), G(1/std::sqrt(1-pow(beta,2))), d(dir) + {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real, amrex::Real, amrex::Real) const noexcept + { + amrex::Real x1, x2; + amrex::Real g, ux; + amrex::Real u [3]; + x1 = 0.; + g = 0.; + u[d] = 0.; + while( u[d]-g <= x1) + { + u[d] = -t*std::log(amrex::Random()*amrex::Random()*amrex::Random()); + g = std::sqrt(1+std::pow(u[d],2)); + x1 = t*std::log(amrex::Random()); + } + x1 = amrex::Random(); + x2 = amrex::Random(); + ux = u[d]*(2*x1-1); + if(-b*ux/g>x2) + { + ux = -ux; + } + ux = G*(ux+g*b); + x2 = amrex::Random(); + u[(d+1)%3] = 2*u[d]*std::sqrt(x1*(1-x1))*std::cos(2*M_PI*x2); + u[(d+2)%3] = 2*u[d]*std::sqrt(x1*(1-x1))*std::cos(2*M_PI*x2); + u[d] = ux; + return amrex::XDim3 {u[0],u[1],u[2]}; + } + +private: + int d; + amrex::Real b, t, G; +}; + + // struct whose getMomentum returns momentum for 1 particle, for // radial expansion struct InjectorMomentumRadialExpansion @@ -124,6 +169,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 +217,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 +242,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 +260,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 +274,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..417959ab3 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"){ + dir = 0; + }else if(direction == "y"){ + dir = 1; + }else if(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); |