WarpX
SharedDepositionUtils.H
Go to the documentation of this file.
1 /* Copyright 2022 Noah Kaplan, Andrew Myers, Phil Miller
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef WARPX_SHAREDDEPOSITIONUTILS_H_
8 #define WARPX_SHAREDDEPOSITIONUTILS_H_
9 
11 #include "Particles/ShapeFactors.H"
13 #include "Utils/WarpXConst.H"
14 #ifdef WARPX_DIM_RZ
15 # include "Utils/WarpX_Complex.H"
16 #endif
17 
18 #include <AMReX.H>
19 
20 /*
21  * \brief gets the maximum width, height, or length of a tilebox. In number of cells.
22  * \param nCells : Number of cells in the direction to be considered
23  * \param tilesize : The 1D tilesize in the direction to be considered
24  */
26 int getMaxTboxAlongDim (int nCells, int tilesize){
27  int maxTilesize = 0;
28  const int nTiles = nCells / tilesize;
29  const int remainder = nCells % tilesize;
30  maxTilesize = tilesize + int(std::ceil((amrex::Real) remainder / nTiles));
31  return maxTilesize;
32 }
33 
34 /*
35  * \brief atomically add the values from the local deposition buffer back to the global array.
36  * \param bx : Box defining the index space of the local buffer
37  * \param global : The global array
38  * \param local : The local array
39  */
40 #if defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)
42 void addLocalToGlobal (const amrex::Box& bx,
43  const amrex::Array4<amrex::Real>& global,
44  const amrex::Array4<amrex::Real>& local) noexcept
45 {
46  using namespace amrex::literals;
47 
48  const auto lo = amrex::lbound(bx);
49  const auto len = amrex::length(bx);
50  for (int icell = threadIdx.x; icell < bx.numPts(); icell += blockDim.x)
51  {
52  int k = icell / (len.x*len.y);
53  int j = (icell - k*(len.x*len.y)) / len.x;
54  int i = (icell - k*(len.x*len.y)) - j*len.x;
55  i += lo.x;
56  j += lo.y;
57  k += lo.z;
58  if (amrex::Math::abs(local(i, j, k)) > 0.0_rt) {
59  amrex::Gpu::Atomic::AddNoRet( &global(i, j, k), local(i, j, k));
60  }
61  }
62 }
63 #endif
64 
65 #if defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)
66 template <int depos_order>
68 void depositComponent (const GetParticlePosition<PIdx>& GetPosition,
69  const amrex::ParticleReal * const wp,
70  const amrex::ParticleReal * const uxp,
71  const amrex::ParticleReal * const uyp,
72  const amrex::ParticleReal * const uzp,
73  const int* ion_lev,
74  amrex::Array4<amrex::Real> const& j_buff,
75  amrex::IntVect const j_type,
76  const amrex::Real relative_time,
77  const std::array<amrex::Real,3>& dx,
78  const std::array<amrex::Real,3>& xyzmin,
79  const amrex::Dim3 lo,
80  const amrex::Real q,
81  const int n_rz_azimuthal_modes,
82  const unsigned int ip,
83  const int zdir, const int NODE, const int CELL, const int dir)
84 {
85  using namespace amrex::literals;
86 
87 #if !defined(WARPX_DIM_RZ)
88  amrex::ignore_unused(n_rz_azimuthal_modes);
89 #endif
90 
91  // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
92  // (do_ionization=1)
93  const bool do_ionization = ion_lev;
94  const amrex::Real dzi = 1.0_rt/dx[2];
95 #if defined(WARPX_DIM_1D_Z)
96  const amrex::Real invvol = dzi;
97 #endif
98 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
99  const amrex::Real dxi = 1.0_rt/dx[0];
100  const amrex::Real invvol = dxi*dzi;
101 #elif defined(WARPX_DIM_3D)
102  const amrex::Real dxi = 1.0_rt/dx[0];
103  const amrex::Real dyi = 1.0_rt/dx[1];
104  const amrex::Real invvol = dxi*dyi*dzi;
105 #endif
106 
107 #if (AMREX_SPACEDIM >= 2)
108  const amrex::Real xmin = xyzmin[0];
109 #endif
110 #if defined(WARPX_DIM_3D)
111  const amrex::Real ymin = xyzmin[1];
112 #endif
113  const amrex::Real zmin = xyzmin[2];
114 
115  const amrex::Real clightsq = 1.0_rt/PhysConst::c/PhysConst::c;
116 
117  // --- Get particle quantities
118  const amrex::Real gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
119  + uyp[ip]*uyp[ip]*clightsq
120  + uzp[ip]*uzp[ip]*clightsq);
121  amrex::Real wq = q*wp[ip];
122  if (do_ionization){
123  wq *= ion_lev[ip];
124  }
125 
126  amrex::ParticleReal xp, yp, zp;
127  GetPosition(ip, xp, yp, zp);
128 
129  const amrex::Real vx = uxp[ip]*gaminv;
130  const amrex::Real vy = uyp[ip]*gaminv;
131  const amrex::Real vz = uzp[ip]*gaminv;
132  // pcurrent is the particle current in the deposited direction
133 #if defined(WARPX_DIM_RZ)
134  // In RZ, wqx is actually wqr, and wqy is wqtheta
135  // Convert to cylindrical at the mid point
136  const amrex::Real xpmid = xp + relative_time*vx;
137  const amrex::Real ypmid = yp + relative_time*vy;
138  const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
139  amrex::Real costheta;
140  amrex::Real sintheta;
141  if (rpmid > 0._rt) {
142  costheta = xpmid/rpmid;
143  sintheta = ypmid/rpmid;
144  } else {
145  costheta = 1._rt;
146  sintheta = 0._rt;
147  }
148  const Complex xy0 = Complex{costheta, sintheta};
149  const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
150  const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
151 #else
152  const amrex::Real wqx = wq*invvol*vx;
153  const amrex::Real wqy = wq*invvol*vy;
154 #endif
155  const amrex::Real wqz = wq*invvol*vz;
156 
157  amrex::Real pcurrent = 0.0;
158  if (dir == 0) {
159  pcurrent = wqx;
160  } else if (dir == 1) {
161  pcurrent = wqy;
162  } else if (dir == 2) {
163  pcurrent = wqz;
164  }
165 
166  // --- Compute shape factors
167  Compute_shape_factor< depos_order > const compute_shape_factor;
168 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D)
169 
170  // x direction
171  // Get particle position after 1/2 push back in position
172 #if defined(WARPX_DIM_RZ)
173  // Keep these double to avoid bug in single precision
174  const double xmid = (rpmid - xmin)*dxi;
175 #else
176  const double xmid = ((xp - xmin) + relative_time*vx)*dxi;
177 #endif
178  // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
179  // sx_j[xyz] shape factor along x for the centering of each current
180  // There are only two possible centerings, node or cell centered, so at most only two shape factor
181  // arrays will be needed.
182  // Keep these double to avoid bug in single precision
183  double sx_node[depos_order + 1] = {0.};
184  double sx_cell[depos_order + 1] = {0.};
185  int j_node = 0;
186  int j_cell = 0;
187  if (j_type[0] == NODE) {
188  j_node = compute_shape_factor(sx_node, xmid);
189  }
190  if (j_type[0] == CELL) {
191  j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
192  }
193 
194  amrex::Real sx_j[depos_order + 1] = {0._rt};
195  for (int ix=0; ix<=depos_order; ix++)
196  {
197  sx_j[ix] = ((j_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
198  }
199 
200  int const j_j = ((j_type[0] == NODE) ? j_node : j_cell);
201 #endif //AMREX_SPACEDIM >= 2
202 
203 #if defined(WARPX_DIM_3D)
204  // y direction
205  // Keep these double to avoid bug in single precision
206  const double ymid = ((yp - ymin) + relative_time*vy)*dyi;
207  double sy_node[depos_order + 1] = {0.};
208  double sy_cell[depos_order + 1] = {0.};
209  int k_node = 0;
210  int k_cell = 0;
211  if (j_type[1] == NODE) {
212  k_node = compute_shape_factor(sy_node, ymid);
213  }
214  if (j_type[1] == CELL) {
215  k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
216  }
217  amrex::Real sy_j[depos_order + 1] = {0._rt};
218  for (int iy=0; iy<=depos_order; iy++)
219  {
220  sy_j[iy] = ((j_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
221  }
222  int const k_j = ((j_type[1] == NODE) ? k_node : k_cell);
223 #endif
224 
225  // z direction
226  // Keep these double to avoid bug in single precision
227  const double zmid = ((zp - zmin) + relative_time*vz)*dzi;
228  double sz_node[depos_order + 1] = {0.};
229  double sz_cell[depos_order + 1] = {0.};
230  int l_node = 0;
231  int l_cell = 0;
232  if (j_type[zdir] == NODE) {
233  l_node = compute_shape_factor(sz_node, zmid);
234  }
235  if (j_type[zdir] == CELL) {
236  l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
237  }
238  amrex::Real sz_j[depos_order + 1] = {0._rt};
239  for (int iz=0; iz<=depos_order; iz++)
240  {
241  sz_j[iz] = ((j_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
242  }
243  int const l_j = ((j_type[zdir] == NODE) ? l_node : l_cell);
244 
245  // Deposit current into j_buff
246 #if defined(WARPX_DIM_1D_Z)
247  for (int iz=0; iz<=depos_order; iz++){
249  &j_buff(lo.x+l_j+iz, 0, 0, 0),
250  sz_j[iz]*pcurrent);
251  }
252 #endif
253 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
254  for (int iz=0; iz<=depos_order; iz++){
255  for (int ix=0; ix<=depos_order; ix++){
257  &j_buff(lo.x+j_j+ix, lo.y+l_j+iz, 0, 0),
258  sx_j[ix]*sz_j[iz]*pcurrent);
259 #if defined(WARPX_DIM_RZ)
260  Complex xy = xy0; // Note that xy is equal to e^{i m theta}
261  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
262  // The factor 2 on the weighting comes from the normalization of the modes
263  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());
264  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());
265  xy = xy*xy0;
266  }
267 #endif
268  }
269  }
270 #elif defined(WARPX_DIM_3D)
271  for (int iz=0; iz<=depos_order; iz++){
272  for (int iy=0; iy<=depos_order; iy++){
273  for (int ix=0; ix<=depos_order; ix++){
275  &j_buff(lo.x+j_j+ix, lo.y+k_j+iy, lo.z+l_j+iz),
276  sx_j[ix]*sy_j[iy]*sz_j[iz]*pcurrent);
277  }
278  }
279  }
280 #endif
281 }
282 #endif
283 
284 #endif // WARPX_SHAREDDEPOSITIONUTILS_H_
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE int getMaxTboxAlongDim(int nCells, int tilesize)
Definition: SharedDepositionUtils.H:26
NODE
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
constexpr int iz
constexpr int iy
constexpr int ix
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 length(Array4< T > const &a) noexcept
const int[]
i
Definition: check_interp_points_and_weights.py:174
tuple dx
lab frame
Definition: stencil.py:429
xmin
Definition: stencil.py:424
Definition: ShapeFactors.H:29
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T real() const noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T imag() const noexcept