diff options
| -rw-r--r-- | GNUmakefile | 4 | ||||
| -rw-r--r-- | Source/Initialization/InjectorMomentum.H | 72 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 22 |
3 files changed, 95 insertions, 3 deletions
diff --git a/GNUmakefile b/GNUmakefile index ef86420f6..5cd8f99e4 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -5,8 +5,8 @@ OPENBC_HOME ?= ../openbc_poisson DEBUG = FALSE #DEBUG = TRUE -#DIM = 2 -DIM = 3 +DIM = 2 +#DIM = 3 #QED = TRUE #QED_TABLE_GEN = TRUE 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; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 32711cd31..300c41f7b 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_boltzmann"){ + 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 InjectorMomentumBoltzmann. + inj_mom.reset(new InjectorMomentum((InjectorMomentumBoltzmann*)nullptr, theta, beta, dir)); } else if (mom_dist_s == "maxwell_juttner"){ Real beta = 0.; Real theta = 10.; |