1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
/* Copyright 2019-2020 Neil Zaim, Yinjian Zhao
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef WARPX_PARTICLE_UTILS_H_
#define WARPX_PARTICLE_UTILS_H_
#include "Particles/WarpXParticleContainer.H"
#include "Utils/WarpXConst.H"
#include <AMReX_DenseBins.H>
#include <AMReX_Particles.H>
#include <AMReX_BaseFwd.H>
namespace ParticleUtils {
/**
* \brief Find the particles and count the particles that are in each cell. More specifically
* this function returns an amrex::DenseBins object containing an offset array and a permutation
* array which can be used to loop over all the cells in a tile and apply an algorithm to
* particles of a given species present in each cell.
* Note that this does *not* rearrange particle arrays.
*
* @param[in] lev the index of the refinement level.
* @param[in] mfi the MultiFAB iterator.
* @param[in] ptile the particle tile.
*/
amrex::DenseBins<WarpXParticleContainer::ParticleType>
findParticlesInEachCell( int const lev, amrex::MFIter const& mfi,
WarpXParticleContainer::ParticleTileType const& ptile);
/**
* \brief Return (relativistic) particle energy given velocity and mass.
* Note the use of `double` since this calculation is prone to error with
* single precision.
*
* @param[in] u2 square of particle speed (i.e. u dot u where u = gamma*v)
* @param[in] mass particle mass
* @param[out] energy particle energy in eV
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void getEnergy ( amrex::ParticleReal const u2, double const mass,
double& energy )
{
using std::sqrt;
using namespace amrex::literals;
constexpr auto c2 = PhysConst::c * PhysConst::c;
energy = mass * u2 / (sqrt(1.0_rt + u2 / c2) + 1.0_rt) / PhysConst::q_e;
}
/**
* \brief Return (relativistic) collision energy assuming the target (with
* mass M) is stationary and the projectile is approaching with the
* the given speed and mass m. Note the use of `double` since this
* calculation is prone to error with single precision.
*
* @param[in] u2 square of particle speed (i.e. u dot u where u = gamma*v)
* @param[in] m, M mass of projectile and target, respectively
* @param[out] energy particle energy in eV
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void getCollisionEnergy ( amrex::ParticleReal const u2, double const m,
double const M, double& gamma, double& energy )
{
using std::sqrt;
using namespace amrex::literals;
constexpr auto c2 = PhysConst::c * PhysConst::c;
gamma = sqrt(1.0_rt + u2 / c2);
energy = (
2.0_rt * m * M * u2 / (gamma + 1.0_rt)
/ (M + m + sqrt(m*m + M*M + 2.0_rt * m * M * gamma))
) / PhysConst::q_e;
}
/**
* \brief Perform a Lorentz transformation of the given velocity
* to a frame moving with velocity (Vx, Vy, Vz) relative to the present one.
*
* @param[in/out] ux,uy,uz components of velocity vector in the current
frame - importantly these quantities are gamma * velocity
* @param[in] Vx,Vy,Vz velocity of the new frame relative to the current one,
NOT gamma*velocity!
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void doLorentzTransform ( amrex::ParticleReal& ux, amrex::ParticleReal& uy,
amrex::ParticleReal& uz,
amrex::ParticleReal const Vx, amrex::ParticleReal const Vy,
amrex::ParticleReal const Vz )
{
// precompute repeatedly used quantities
constexpr auto c2 = PhysConst::c * PhysConst::c;
const auto V2 = (Vx*Vx + Vy*Vy + Vz*Vz);
const auto gamma_V = 1.0_prt / sqrt(1.0_prt - V2 / c2);
const auto gamma_u = sqrt(1.0_prt + (ux*ux + uy*uy + uz*uz) / c2);
// copy velocity vector values
const auto vx = ux;
const auto vy = uy;
const auto vz = uz;
ux = vx * (1.0_prt + (gamma_V - 1.0_prt) * Vx*Vx/V2)
+ vy * (gamma_V - 1.0_prt) * Vx*Vy/V2
+ vz * (gamma_V - 1.0_prt) * Vx*Vz/V2
- gamma_V * Vx * gamma_u;
uy = vy * (1.0_prt + (gamma_V - 1.0_prt) * Vy*Vy/V2)
+ vx * (gamma_V - 1.0_prt) * Vx*Vy/V2
+ vz * (gamma_V - 1.0_prt) * Vy*Vz/V2
- gamma_V * Vy * gamma_u;
uz = vz * (1.0_prt + (gamma_V - 1.0_prt) * Vz*Vz/V2)
+ vx * (gamma_V - 1.0_prt) * Vx*Vz/V2
+ vy * (gamma_V - 1.0_prt) * Vy*Vz/V2
- gamma_V * Vz * gamma_u;
}
/**
* \brief Generate random unit vector in 3 dimensions
* https://mathworld.wolfram.com/SpherePointPicking.html
*
* @param[out] x x-component of resulting random vector
* @param[out] y y-component of resulting random vector
* @param[out] z z-component of resulting random vector
* @param[in] engine the random-engine
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void getRandomVector ( amrex::ParticleReal& x, amrex::ParticleReal& y,
amrex::ParticleReal& z, amrex::RandomEngine const& engine )
{
using std::sqrt;
using std::cos;
using std::sin;
using namespace amrex::literals;
auto const theta = amrex::Random(engine) * 2.0_prt * MathConst::pi;
z = 2.0_prt * amrex::Random(engine) - 1.0_prt;
auto const xy = sqrt(1_prt - z*z);
x = xy * cos(theta);
y = xy * sin(theta);
}
/** \brief Function to perform scattering of a particle that results in a
* random velocity vector with given magnitude. This is used in isotropic
* collision events.
*
* @param[in/out] ux, uy, uz colliding particle's velocity
* @param[in] vp velocity magnitude of the colliding particle after collision.
* @param[in] engine the random-engine
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void RandomizeVelocity ( amrex::ParticleReal& ux, amrex::ParticleReal& uy,
amrex::ParticleReal& uz,
const amrex::ParticleReal vp,
amrex::RandomEngine const& engine )
{
amrex::ParticleReal x, y, z;
// generate random unit vector for the new velocity direction
getRandomVector(x, y, z, engine);
// scale new vector to have the desired magnitude
ux = x * vp;
uy = y * vp;
uz = z * vp;
}
}
#endif // WARPX_PARTICLE_UTILS_H_
|
