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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
|
/* Copyright 2022 Noah Kaplan, Andrew Myers, Phil Miller
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef SHAREDDEPOSITIONUTILS_H_
#define SHAREDDEPOSITIONUTILS_H_
#include "Particles/Pusher/GetAndSetPosition.H"
#include "Particles/ShapeFactors.H"
#include "Utils/WarpXAlgorithmSelection.H"
#include "Utils/WarpXConst.H"
#ifdef WARPX_DIM_RZ
# include "Utils/WarpX_Complex.H"
#endif
#include <AMReX.H>
/*
* \brief gets the maximum width, height, or length of a tilebox. In number of cells.
* \param nCells : Number of cells in the direction to be considered
* \param tilesize : The 1D tilesize in the direction to be considered
*/
AMREX_FORCE_INLINE
int getMaxTboxAlongDim (int nCells, int tilesize){
int maxTilesize = 0;
const int nTiles = nCells / tilesize;
const int remainder = nCells % tilesize;
maxTilesize = tilesize + int(std::ceil((amrex::Real) remainder / nTiles));
return maxTilesize;
}
/*
* \brief atomically add the values from the local deposition buffer back to the global array.
* \param bx : Box defining the index space of the local buffer
* \param global : The global array
* \param local : The local array
*/
#if defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void addLocalToGlobal (const amrex::Box& bx,
const amrex::Array4<amrex::Real>& global,
const amrex::Array4<amrex::Real>& local) noexcept
{
using namespace amrex::literals;
const auto lo = amrex::lbound(bx);
const auto len = amrex::length(bx);
for (int icell = threadIdx.x; icell < bx.numPts(); icell += blockDim.x)
{
int k = icell / (len.x*len.y);
int j = (icell - k*(len.x*len.y)) / len.x;
int i = (icell - k*(len.x*len.y)) - j*len.x;
i += lo.x;
j += lo.y;
k += lo.z;
if (amrex::Math::abs(local(i, j, k)) > 0.0_rt) {
amrex::Gpu::Atomic::AddNoRet( &global(i, j, k), local(i, j, k));
}
}
}
#endif
#if defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)
template <int depos_order>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void depositComponent (const GetParticlePosition& GetPosition,
const amrex::ParticleReal * const wp,
const amrex::ParticleReal * const uxp,
const amrex::ParticleReal * const uyp,
const amrex::ParticleReal * const uzp,
const int* ion_lev,
amrex::Array4<amrex::Real> const& j_buff,
amrex::IntVect const j_type,
const amrex::Real relative_time,
const std::array<amrex::Real,3>& dx,
const std::array<amrex::Real,3>& xyzmin,
const amrex::Dim3 lo,
const amrex::Real q,
const int n_rz_azimuthal_modes,
const unsigned int ip,
const int zdir, const int NODE, const int CELL, const int dir)
{
using namespace amrex::literals;
#if !defined(WARPX_DIM_RZ)
amrex::ignore_unused(n_rz_azimuthal_modes);
#endif
// Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
// (do_ionization=1)
const bool do_ionization = ion_lev;
const amrex::Real dzi = 1.0_rt/dx[2];
#if defined(WARPX_DIM_1D_Z)
const amrex::Real invvol = dzi;
#endif
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
const amrex::Real dxi = 1.0_rt/dx[0];
const amrex::Real invvol = dxi*dzi;
#elif defined(WARPX_DIM_3D)
const amrex::Real dxi = 1.0_rt/dx[0];
const amrex::Real dyi = 1.0_rt/dx[1];
const amrex::Real invvol = dxi*dyi*dzi;
#endif
#if (AMREX_SPACEDIM >= 2)
const amrex::Real xmin = xyzmin[0];
#endif
#if defined(WARPX_DIM_3D)
const amrex::Real ymin = xyzmin[1];
#endif
const amrex::Real zmin = xyzmin[2];
const amrex::Real clightsq = 1.0_rt/PhysConst::c/PhysConst::c;
// --- Get particle quantities
const amrex::Real gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
+ uyp[ip]*uyp[ip]*clightsq
+ uzp[ip]*uzp[ip]*clightsq);
amrex::Real wq = q*wp[ip];
if (do_ionization){
wq *= ion_lev[ip];
}
amrex::ParticleReal xp, yp, zp;
GetPosition(ip, xp, yp, zp);
const amrex::Real vx = uxp[ip]*gaminv;
const amrex::Real vy = uyp[ip]*gaminv;
const amrex::Real vz = uzp[ip]*gaminv;
// pcurrent is the particle current in the deposited direction
#if defined(WARPX_DIM_RZ)
// In RZ, wqx is actually wqr, and wqy is wqtheta
// Convert to cylinderical at the mid point
const amrex::Real xpmid = xp + relative_time*vx;
const amrex::Real ypmid = yp + relative_time*vy;
const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
amrex::Real costheta;
amrex::Real sintheta;
if (rpmid > 0._rt) {
costheta = xpmid/rpmid;
sintheta = ypmid/rpmid;
} else {
costheta = 1._rt;
sintheta = 0._rt;
}
const Complex xy0 = Complex{costheta, sintheta};
const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
#else
const amrex::Real wqx = wq*invvol*vx;
const amrex::Real wqy = wq*invvol*vy;
#endif
const amrex::Real wqz = wq*invvol*vz;
amrex::Real pcurrent = 0.0;
if (dir == 0) {
pcurrent = wqx;
} else if (dir == 1) {
pcurrent = wqy;
} else if (dir == 2) {
pcurrent = wqz;
}
// --- Compute shape factors
Compute_shape_factor< depos_order > const compute_shape_factor;
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D)
// x direction
// Get particle position after 1/2 push back in position
#if defined(WARPX_DIM_RZ)
// Keep these double to avoid bug in single precision
const double xmid = (rpmid - xmin)*dxi;
#else
const double xmid = ((xp - xmin) + relative_time*vx)*dxi;
#endif
// j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
// sx_j[xyz] shape factor along x for the centering of each current
// There are only two possible centerings, node or cell centered, so at most only two shape factor
// arrays will be needed.
// Keep these double to avoid bug in single precision
double sx_node[depos_order + 1] = {0.};
double sx_cell[depos_order + 1] = {0.};
int j_node = 0;
int j_cell = 0;
if (j_type[0] == NODE) {
j_node = compute_shape_factor(sx_node, xmid);
}
if (j_type[0] == CELL) {
j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
}
amrex::Real sx_j[depos_order + 1] = {0._rt};
for (int ix=0; ix<=depos_order; ix++)
{
sx_j[ix] = ((j_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
}
int const j_j = ((j_type[0] == NODE) ? j_node : j_cell);
#endif //AMREX_SPACEDIM >= 2
#if defined(WARPX_DIM_3D)
// y direction
// Keep these double to avoid bug in single precision
const double ymid = ((yp - ymin) + relative_time*vy)*dyi;
double sy_node[depos_order + 1] = {0.};
double sy_cell[depos_order + 1] = {0.};
int k_node = 0;
int k_cell = 0;
if (j_type[1] == NODE) {
k_node = compute_shape_factor(sy_node, ymid);
}
if (j_type[1] == CELL) {
k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
}
amrex::Real sy_j[depos_order + 1] = {0._rt};
for (int iy=0; iy<=depos_order; iy++)
{
sy_j[iy] = ((j_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
}
int const k_j = ((j_type[1] == NODE) ? k_node : k_cell);
#endif
// z direction
// Keep these double to avoid bug in single precision
const double zmid = ((zp - zmin) + relative_time*vz)*dzi;
double sz_node[depos_order + 1] = {0.};
double sz_cell[depos_order + 1] = {0.};
int l_node = 0;
int l_cell = 0;
if (j_type[zdir] == NODE) {
l_node = compute_shape_factor(sz_node, zmid);
}
if (j_type[zdir] == CELL) {
l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
}
amrex::Real sz_j[depos_order + 1] = {0._rt};
for (int iz=0; iz<=depos_order; iz++)
{
sz_j[iz] = ((j_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
}
int const l_j = ((j_type[zdir] == NODE) ? l_node : l_cell);
// Deposit current into j_buff
#if defined(WARPX_DIM_1D_Z)
for (int iz=0; iz<=depos_order; iz++){
amrex::Gpu::Atomic::AddNoRet(
&j_buff(lo.x+l_j+iz, 0, 0, 0),
sz_j[iz]*pcurrent);
}
#endif
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
for (int iz=0; iz<=depos_order; iz++){
for (int ix=0; ix<=depos_order; ix++){
amrex::Gpu::Atomic::AddNoRet(
&j_buff(lo.x+j_j+ix, lo.y+l_j+iz, 0, 0),
sx_j[ix]*sz_j[iz]*pcurrent);
#if defined(WARPX_DIM_RZ)
Complex xy = xy0; // Note that xy is equal to e^{i m theta}
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 on the weighting comes from the normalization of the modes
amrex::Gpu::Atomic::AddNoRet( &j_buff(lo.x+j_j+ix, lo.y+l_j+iz, 0, 2*imode-1), 2._rt*sx_j[ix]*sz_j[iz]*wqx*xy.real());
amrex::Gpu::Atomic::AddNoRet( &j_buff(lo.x+j_j+ix, lo.y+l_j+iz, 0, 2*imode ), 2._rt*sx_j[ix]*sz_j[iz]*wqx*xy.imag());
xy = xy*xy0;
}
#endif
}
}
#elif defined(WARPX_DIM_3D)
for (int iz=0; iz<=depos_order; iz++){
for (int iy=0; iy<=depos_order; iy++){
for (int ix=0; ix<=depos_order; ix++){
amrex::Gpu::Atomic::AddNoRet(
&j_buff(lo.x+j_j+ix, lo.y+k_j+iy, lo.z+l_j+iz),
sx_j[ix]*sy_j[iy]*sz_j[iz]*pcurrent);
}
}
}
#endif
}
#endif
#endif // SHAREDDEPOSITIONUTILS_H_
|
