WarpX
CurrentDeposition.H
Go to the documentation of this file.
1 /* Copyright 2019 Axel Huebl, David Grote, Maxence Thevenet
2  * Remi Lehe, Weiqun Zhang, Michael Rowan
3  *
4  * This file is part of WarpX.
5  *
6  * License: BSD-3-Clause-LBNL
7  */
8 #ifndef CURRENTDEPOSITION_H_
9 #define CURRENTDEPOSITION_H_
10 
12 #include "Particles/ShapeFactors.H"
13 #include "Utils/WarpX_Complex.H"
14 
15 #include <AMReX.H>
16 #include <AMReX_Array4.H>
17 #include <AMReX_REAL.H>
18 
19 using namespace amrex::literals;
20 
40 template <int depos_order>
41 void doDepositionShapeN(const GetParticlePosition& GetPosition,
42  const amrex::ParticleReal * const wp,
43  const amrex::ParticleReal * const uxp,
44  const amrex::ParticleReal * const uyp,
45  const amrex::ParticleReal * const uzp,
46  const int * const ion_lev,
47  amrex::FArrayBox& jx_fab,
48  amrex::FArrayBox& jy_fab,
49  amrex::FArrayBox& jz_fab,
50  const long np_to_depose, const amrex::Real dt,
51  const std::array<amrex::Real,3>& dx,
52  const std::array<amrex::Real,3>& xyzmin,
53  const amrex::Dim3 lo,
54  const amrex::Real q,
55  const long n_rz_azimuthal_modes)
56 {
57 #if !defined(WARPX_DIM_RZ)
58  amrex::ignore_unused(n_rz_azimuthal_modes);
59 #endif
60 
61  // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
62  // (do_ionization=1)
63  const bool do_ionization = ion_lev;
64  const amrex::Real dxi = 1.0/dx[0];
65  const amrex::Real dzi = 1.0/dx[2];
66 #if !(defined WARPX_DIM_RZ)
67  const amrex::Real dts2dx = 0.5*dt*dxi;
68 #endif
69  const amrex::Real dts2dz = 0.5*dt*dzi;
70 #if (AMREX_SPACEDIM == 2)
71  const amrex::Real invvol = dxi*dzi;
72 #elif (defined WARPX_DIM_3D)
73  const amrex::Real dyi = 1.0/dx[1];
74  const amrex::Real dts2dy = 0.5*dt*dyi;
75  const amrex::Real invvol = dxi*dyi*dzi;
76 #endif
77 
78  const amrex::Real xmin = xyzmin[0];
79 #if (defined WARPX_DIM_3D)
80  const amrex::Real ymin = xyzmin[1];
81 #endif
82  const amrex::Real zmin = xyzmin[2];
83 
84  const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c;
85 
86  amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
87  amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
88  amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
89  amrex::IntVect const jx_type = jx_fab.box().type();
90  amrex::IntVect const jy_type = jy_fab.box().type();
91  amrex::IntVect const jz_type = jz_fab.box().type();
92 
93  constexpr int zdir = (AMREX_SPACEDIM - 1);
94  constexpr int NODE = amrex::IndexType::NODE;
95  constexpr int CELL = amrex::IndexType::CELL;
96 
97  // Loop over particles and deposit into jx_fab, jy_fab and jz_fab
98  amrex::ParallelFor(
99  np_to_depose,
100  [=] AMREX_GPU_DEVICE (long ip) {
101  // --- Get particle quantities
102  const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
103  + uyp[ip]*uyp[ip]*clightsq
104  + uzp[ip]*uzp[ip]*clightsq);
105  amrex::Real wq = q*wp[ip];
106  if (do_ionization){
107  wq *= ion_lev[ip];
108  }
109 
110  amrex::ParticleReal xp, yp, zp;
111  GetPosition(ip, xp, yp, zp);
112 
113  const amrex::Real vx = uxp[ip]*gaminv;
114  const amrex::Real vy = uyp[ip]*gaminv;
115  const amrex::Real vz = uzp[ip]*gaminv;
116  // wqx, wqy wqz are particle current in each direction
117 #if (defined WARPX_DIM_RZ)
118  // In RZ, wqx is actually wqr, and wqy is wqtheta
119  // Convert to cylinderical at the mid point
120  const amrex::Real xpmid = xp - 0.5*dt*vx;
121  const amrex::Real ypmid = yp - 0.5*dt*vy;
122  const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
123  amrex::Real costheta;
124  amrex::Real sintheta;
125  if (rpmid > 0.) {
126  costheta = xpmid/rpmid;
127  sintheta = ypmid/rpmid;
128  } else {
129  costheta = 1.;
130  sintheta = 0.;
131  }
132  const Complex xy0 = Complex{costheta, sintheta};
133  const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
134  const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
135 #else
136  const amrex::Real wqx = wq*invvol*vx;
137  const amrex::Real wqy = wq*invvol*vy;
138 #endif
139  const amrex::Real wqz = wq*invvol*vz;
140 
141  // --- Compute shape factors
142  // x direction
143  // Get particle position after 1/2 push back in position
144 #if (defined WARPX_DIM_RZ)
145  // Keep these double to avoid bug in single precision
146  const double xmid = (rpmid - xmin)*dxi;
147 #else
148  const double xmid = (xp - xmin)*dxi - dts2dx*vx;
149 #endif
150  // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
151  // sx_j[xyz] shape factor along x for the centering of each current
152  // There are only two possible centerings, node or cell centered, so at most only two shape factor
153  // arrays will be needed.
154  // Keep these double to avoid bug in single precision
155  double sx_node[depos_order + 1];
156  double sx_cell[depos_order + 1];
157  int j_node = 0;
158  int j_cell = 0;
159  Compute_shape_factor< depos_order > const compute_shape_factor;
160  if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) {
161  j_node = compute_shape_factor(sx_node, xmid);
162  }
163  if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
164  j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
165  }
166 
167  amrex::Real sx_jx[depos_order + 1] = {0.};
168  amrex::Real sx_jy[depos_order + 1] = {0.};
169  amrex::Real sx_jz[depos_order + 1] = {0.};
170  for (int ix=0; ix<=depos_order; ix++)
171  {
172  sx_jx[ix] = ((jx_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
173  sx_jy[ix] = ((jy_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
174  sx_jz[ix] = ((jz_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
175  }
176 
177  int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell);
178  int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell);
179  int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell);
180 
181 #if (defined WARPX_DIM_3D)
182  // y direction
183  // Keep these double to avoid bug in single precision
184  const double ymid = (yp - ymin)*dyi - dts2dy*vy;
185  double sy_node[depos_order + 1];
186  double sy_cell[depos_order + 1];
187  int k_node = 0;
188  int k_cell = 0;
189  if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) {
190  k_node = compute_shape_factor(sy_node, ymid);
191  }
192  if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
193  k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
194  }
195  amrex::Real sy_jx[depos_order + 1] = {0.};
196  amrex::Real sy_jy[depos_order + 1] = {0.};
197  amrex::Real sy_jz[depos_order + 1] = {0.};
198  for (int iy=0; iy<=depos_order; iy++)
199  {
200  sy_jx[iy] = ((jx_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
201  sy_jy[iy] = ((jy_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
202  sy_jz[iy] = ((jz_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
203  }
204  int const k_jx = ((jx_type[1] == NODE) ? k_node : k_cell);
205  int const k_jy = ((jy_type[1] == NODE) ? k_node : k_cell);
206  int const k_jz = ((jz_type[1] == NODE) ? k_node : k_cell);
207 #endif
208 
209  // z direction
210  // Keep these double to avoid bug in single precision
211  const double zmid = (zp - zmin)*dzi - dts2dz*vz;
212  double sz_node[depos_order + 1];
213  double sz_cell[depos_order + 1];
214  int l_node = 0;
215  int l_cell = 0;
216  if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) {
217  l_node = compute_shape_factor(sz_node, zmid);
218  }
219  if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
220  l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
221  }
222  amrex::Real sz_jx[depos_order + 1] = {0.};
223  amrex::Real sz_jy[depos_order + 1] = {0.};
224  amrex::Real sz_jz[depos_order + 1] = {0.};
225  for (int iz=0; iz<=depos_order; iz++)
226  {
227  sz_jx[iz] = ((jx_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
228  sz_jy[iz] = ((jy_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
229  sz_jz[iz] = ((jz_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
230  }
231  int const l_jx = ((jx_type[zdir] == NODE) ? l_node : l_cell);
232  int const l_jy = ((jy_type[zdir] == NODE) ? l_node : l_cell);
233  int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell);
234 
235  // Deposit current into jx_arr, jy_arr and jz_arr
236 #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
237  for (int iz=0; iz<=depos_order; iz++){
238  for (int ix=0; ix<=depos_order; ix++){
239  amrex::Gpu::Atomic::Add(
240  &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
241  sx_jx[ix]*sz_jx[iz]*wqx);
242  amrex::Gpu::Atomic::Add(
243  &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
244  sx_jy[ix]*sz_jy[iz]*wqy);
245  amrex::Gpu::Atomic::Add(
246  &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
247  sx_jz[ix]*sz_jz[iz]*wqz);
248 #if (defined WARPX_DIM_RZ)
249  Complex xy = xy0; // Note that xy is equal to e^{i m theta}
250  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
251  // The factor 2 on the weighting comes from the normalization of the modes
252  amrex::Gpu::Atomic::Add( &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 2*imode-1), 2.*sx_jx[ix]*sz_jx[iz]*wqx*xy.real());
253  amrex::Gpu::Atomic::Add( &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 2*imode ), 2.*sx_jx[ix]*sz_jx[iz]*wqx*xy.imag());
254  amrex::Gpu::Atomic::Add( &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 2*imode-1), 2.*sx_jy[ix]*sz_jy[iz]*wqy*xy.real());
255  amrex::Gpu::Atomic::Add( &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 2*imode ), 2.*sx_jy[ix]*sz_jy[iz]*wqy*xy.imag());
256  amrex::Gpu::Atomic::Add( &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 2*imode-1), 2.*sx_jz[ix]*sz_jz[iz]*wqz*xy.real());
257  amrex::Gpu::Atomic::Add( &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 2*imode ), 2.*sx_jz[ix]*sz_jz[iz]*wqz*xy.imag());
258  xy = xy*xy0;
259  }
260 #endif
261  }
262  }
263 #elif (defined WARPX_DIM_3D)
264  for (int iz=0; iz<=depos_order; iz++){
265  for (int iy=0; iy<=depos_order; iy++){
266  for (int ix=0; ix<=depos_order; ix++){
267  amrex::Gpu::Atomic::Add(
268  &jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz),
269  sx_jx[ix]*sy_jx[iy]*sz_jx[iz]*wqx);
270  amrex::Gpu::Atomic::Add(
271  &jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz),
272  sx_jy[ix]*sy_jy[iy]*sz_jy[iz]*wqy);
273  amrex::Gpu::Atomic::Add(
274  &jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz),
275  sx_jz[ix]*sy_jz[iy]*sz_jz[iz]*wqz);
276  }
277  }
278  }
279 #endif
280  }
281  );
282 }
283 
305 template <int depos_order>
307  const amrex::ParticleReal * const wp,
308  const amrex::ParticleReal * const uxp,
309  const amrex::ParticleReal * const uyp,
310  const amrex::ParticleReal * const uzp,
311  const int * ion_lev,
312  const amrex::Array4<amrex::Real>& Jx_arr,
313  const amrex::Array4<amrex::Real>& Jy_arr,
314  const amrex::Array4<amrex::Real>& Jz_arr,
315  const long np_to_depose,
316  const amrex::Real dt,
317  const std::array<amrex::Real,3>& dx,
318  const std::array<amrex::Real, 3> xyzmin,
319  const amrex::Dim3 lo,
320  const amrex::Real q,
321  const long n_rz_azimuthal_modes)
322 {
323  using namespace amrex;
324 #if !defined(WARPX_DIM_RZ)
325  ignore_unused(n_rz_azimuthal_modes);
326 #endif
327 
328  // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
329  // (do_ionization=1)
330  bool const do_ionization = ion_lev;
331  Real const dxi = 1.0_rt / dx[0];
332 #if !(defined WARPX_DIM_RZ)
333  Real const dtsdx0 = dt*dxi;
334 #endif
335  Real const xmin = xyzmin[0];
336 #if (defined WARPX_DIM_3D)
337  Real const dyi = 1.0_rt / dx[1];
338  Real const dtsdy0 = dt*dyi;
339  Real const ymin = xyzmin[1];
340 #endif
341  Real const dzi = 1.0_rt / dx[2];
342  Real const dtsdz0 = dt*dzi;
343  Real const zmin = xyzmin[2];
344 
345 #if (defined WARPX_DIM_3D)
346  Real const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]);
347  Real const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]);
348  Real const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]);
349 #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
350  Real const invdtdx = 1.0_rt / (dt*dx[2]);
351  Real const invdtdz = 1.0_rt / (dt*dx[0]);
352  Real const invvol = 1.0_rt / (dx[0]*dx[2]);
353 #endif
354 
355 #if (defined WARPX_DIM_RZ)
356  Complex const I = Complex{0._rt, 1._rt};
357 #endif
358 
359  Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c );
360 
361  // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
362  amrex::ParallelFor(
363  np_to_depose,
364  [=] AMREX_GPU_DEVICE (long const ip) {
365 
366  // --- Get particle quantities
367  Real const gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
368  + uyp[ip]*uyp[ip]*clightsq
369  + uzp[ip]*uzp[ip]*clightsq);
370 
371  // wqx, wqy wqz are particle current in each direction
372  Real wq = q*wp[ip];
373  if (do_ionization){
374  wq *= ion_lev[ip];
375  }
376 
377  ParticleReal xp, yp, zp;
378  GetPosition(ip, xp, yp, zp);
379 
380  Real const wqx = wq*invdtdx;
381 #if (defined WARPX_DIM_3D)
382  Real const wqy = wq*invdtdy;
383 #endif
384  Real const wqz = wq*invdtdz;
385 
386  // computes current and old position in grid units
387 #if (defined WARPX_DIM_RZ)
388  Real const xp_mid = xp - 0.5_rt * dt*uxp[ip]*gaminv;
389  Real const yp_mid = yp - 0.5_rt * dt*uyp[ip]*gaminv;
390  Real const xp_old = xp - dt*uxp[ip]*gaminv;
391  Real const yp_old = yp - dt*uyp[ip]*gaminv;
392  Real const rp_new = std::sqrt(xp*xp
393  + yp*yp);
394  Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
395  Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
396  Real costheta_new, sintheta_new;
397  if (rp_new > 0._rt) {
398  costheta_new = xp/rp_new;
399  sintheta_new = yp/rp_new;
400  } else {
401  costheta_new = 1._rt;
402  sintheta_new = 0._rt;
403  }
404  amrex::Real costheta_mid, sintheta_mid;
405  if (rp_mid > 0._rt) {
406  costheta_mid = xp_mid/rp_mid;
407  sintheta_mid = yp_mid/rp_mid;
408  } else {
409  costheta_mid = 1._rt;
410  sintheta_mid = 0._rt;
411  }
412  amrex::Real costheta_old, sintheta_old;
413  if (rp_old > 0._rt) {
414  costheta_old = xp_old/rp_old;
415  sintheta_old = yp_old/rp_old;
416  } else {
417  costheta_old = 1._rt;
418  sintheta_old = 0._rt;
419  }
420  const Complex xy_new0 = Complex{costheta_new, sintheta_new};
421  const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
422  const Complex xy_old0 = Complex{costheta_old, sintheta_old};
423  // Keep these double to avoid bug in single precision
424  double const x_new = (rp_new - xmin)*dxi;
425  double const x_old = (rp_old - xmin)*dxi;
426 #else
427  // Keep these double to avoid bug in single precision
428  double const x_new = (xp - xmin)*dxi;
429  double const x_old = x_new - dtsdx0*uxp[ip]*gaminv;
430 #endif
431 #if (defined WARPX_DIM_3D)
432  // Keep these double to avoid bug in single precision
433  double const y_new = (yp - ymin)*dyi;
434  double const y_old = y_new - dtsdy0*uyp[ip]*gaminv;
435 #endif
436  // Keep these double to avoid bug in single precision
437  double const z_new = (zp - zmin)*dzi;
438  double const z_old = z_new - dtsdz0*uzp[ip]*gaminv;
439 
440 #if (defined WARPX_DIM_RZ)
441  Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
442 #elif (defined WARPX_DIM_XZ)
443  Real const vy = uyp[ip]*gaminv;
444 #endif
445 
446  // Shape factor arrays
447  // Note that there are extra values above and below
448  // to possibly hold the factor for the old particle
449  // which can be at a different grid location.
450  // Keep these double to avoid bug in single precision
451  double sx_new[depos_order + 3] = {0.};
452  double sx_old[depos_order + 3] = {0.};
453 #if (defined WARPX_DIM_3D)
454  // Keep these double to avoid bug in single precision
455  double sy_new[depos_order + 3] = {0.};
456  double sy_old[depos_order + 3] = {0.};
457 #endif
458  // Keep these double to avoid bug in single precision
459  double sz_new[depos_order + 3] = {0.};
460  double sz_old[depos_order + 3] = {0.};
461 
462  // --- Compute shape factors
463  // Compute shape factors for position as they are now and at old positions
464  // [ijk]_new: leftmost grid point that the particle touches
465  Compute_shape_factor< depos_order > compute_shape_factor;
466  Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor;
467 
468  const int i_new = compute_shape_factor(sx_new+1, x_new);
469  const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
470 #if (defined WARPX_DIM_3D)
471  const int j_new = compute_shape_factor(sy_new+1, y_new);
472  const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
473 #endif
474  const int k_new = compute_shape_factor(sz_new+1, z_new);
475  const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new);
476 
477  // computes min/max positions of current contributions
478  int dil = 1, diu = 1;
479  if (i_old < i_new) dil = 0;
480  if (i_old > i_new) diu = 0;
481 #if (defined WARPX_DIM_3D)
482  int djl = 1, dju = 1;
483  if (j_old < j_new) djl = 0;
484  if (j_old > j_new) dju = 0;
485 #endif
486  int dkl = 1, dku = 1;
487  if (k_old < k_new) dkl = 0;
488  if (k_old > k_new) dku = 0;
489 
490 #if (defined WARPX_DIM_3D)
491 
492  for (int k=dkl; k<=depos_order+2-dku; k++) {
493  for (int j=djl; j<=depos_order+2-dju; j++) {
494  amrex::Real sdxi = 0._rt;
495  for (int i=dil; i<=depos_order+1-diu; i++) {
496  sdxi += wqx*(sx_old[i] - sx_new[i])*((sy_new[j] + 0.5_rt*(sy_old[j] - sy_new[j]))*sz_new[k] +
497  (0.5_rt*sy_new[j] + 1._rt/3._rt*(sy_old[j] - sy_new[j]))*(sz_old[k] - sz_new[k]));
498  amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdxi);
499  }
500  }
501  }
502  for (int k=dkl; k<=depos_order+2-dku; k++) {
503  for (int i=dil; i<=depos_order+2-diu; i++) {
504  amrex::Real sdyj = 0._rt;
505  for (int j=djl; j<=depos_order+1-dju; j++) {
506  sdyj += wqy*(sy_old[j] - sy_new[j])*((sz_new[k] + 0.5_rt*(sz_old[k] - sz_new[k]))*sx_new[i] +
507  (0.5_rt*sz_new[k] + 1._rt/3._rt*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
508  amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdyj);
509  }
510  }
511  }
512  for (int j=djl; j<=depos_order+2-dju; j++) {
513  for (int i=dil; i<=depos_order+2-diu; i++) {
514  amrex::Real sdzk = 0._rt;
515  for (int k=dkl; k<=depos_order+1-dku; k++) {
516  sdzk += wqz*(sz_old[k] - sz_new[k])*((sx_new[i] + 0.5_rt*(sx_old[i] - sx_new[i]))*sy_new[j] +
517  (0.5_rt*sx_new[i] + 1._rt/3._rt*(sx_old[i] - sx_new[i]))*(sy_old[j] - sy_new[j]));
518  amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdzk);
519  }
520  }
521  }
522 
523 #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
524 
525  for (int k=dkl; k<=depos_order+2-dku; k++) {
526  amrex::Real sdxi = 0._rt;
527  for (int i=dil; i<=depos_order+1-diu; i++) {
528  sdxi += wqx*(sx_old[i] - sx_new[i])*(sz_new[k] + 0.5_rt*(sz_old[k] - sz_new[k]));
529  amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi);
530 #if (defined WARPX_DIM_RZ)
531  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
532  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
533  // The factor 2 comes from the normalization of the modes
534  const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
535  amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real());
536  amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
537  xy_mid = xy_mid*xy_mid0;
538  }
539 #endif
540  }
541  }
542  for (int k=dkl; k<=depos_order+2-dku; k++) {
543  for (int i=dil; i<=depos_order+2-diu; i++) {
544  Real const sdyj = wq*vy*invvol*((sz_new[k] + 0.5_rt * (sz_old[k] - sz_new[k]))*sx_new[i] +
545  (0.5_rt * sz_new[k] + 1._rt / 3._rt *(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
546  amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
547 #if (defined WARPX_DIM_RZ)
548  Complex xy_new = xy_new0;
549  Complex xy_mid = xy_mid0;
550  Complex xy_old = xy_old0;
551  // Throughout the following loop, xy_ takes the value e^{i m theta_}
552  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
553  // The factor 2 comes from the normalization of the modes
554  // The minus sign comes from the different convention with respect to Davidson et al.
555  const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode*
556  (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old));
557  amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real());
558  amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
559  xy_new = xy_new*xy_new0;
560  xy_mid = xy_mid*xy_mid0;
561  xy_old = xy_old*xy_old0;
562  }
563 #endif
564  }
565  }
566  for (int i=dil; i<=depos_order+2-diu; i++) {
567  Real sdzk = 0._rt;
568  for (int k=dkl; k<=depos_order+1-dku; k++) {
569  sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5_rt * (sx_old[i] - sx_new[i]));
570  amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
571 #if (defined WARPX_DIM_RZ)
572  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
573  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
574  // The factor 2 comes from the normalization of the modes
575  const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
576  amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real());
577  amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
578  xy_mid = xy_mid*xy_mid0;
579  }
580 #endif
581  }
582  }
583 
584 #endif
585  }
586  );
587 }
588 
614 template <int depos_order>
616  const amrex::ParticleReal* const wp,
617  const amrex::ParticleReal* const uxp,
618  const amrex::ParticleReal* const uyp,
619  const amrex::ParticleReal* const uzp,
620  const int* const ion_lev,
621  amrex::FArrayBox& jx_fab,
622  amrex::FArrayBox& jy_fab,
623  amrex::FArrayBox& jz_fab,
624  const long np_to_depose,
625  const amrex::Real dt,
626  const std::array<amrex::Real,3>& dx,
627  const std::array<amrex::Real,3>& xyzmin,
628  const amrex::Dim3 lo,
629  const amrex::Real q,
630  const long n_rz_azimuthal_modes)
631 {
632 #if (defined WARPX_DIM_RZ)
633  amrex::ignore_unused(GetPosition,
634  wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab,
635  np_to_depose, dt, dx, xyzmin, lo, q, n_rz_azimuthal_modes);
636  amrex::Abort("Vay deposition not implemented in RZ geometry");
637 #endif
638 
639 #if !(defined WARPX_DIM_RZ)
640  amrex::ignore_unused(n_rz_azimuthal_modes);
641 
642  // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1
643  const bool do_ionization = ion_lev;
644 
645  // Inverse cell volume in each direction
646  const amrex::Real dxi = 1._rt / dx[0];
647  const amrex::Real dzi = 1._rt / dx[2];
648 #if (defined WARPX_DIM_3D)
649  const amrex::Real dyi = 1._rt / dx[1];
650 #endif
651 
652  // Inverse of time step
653  const amrex::Real invdt = 1._rt / dt;
654 
655  // Total inverse cell volume
656 #if (defined WARPX_DIM_XZ)
657  const amrex::Real invvol = dxi * dzi;
658 #elif (defined WARPX_DIM_3D)
659  const amrex::Real invvol = dxi * dyi * dzi;
660 #endif
661 
662  // Lower bound of physical domain in each direction
663  const amrex::Real xmin = xyzmin[0];
664  const amrex::Real zmin = xyzmin[2];
665 #if (defined WARPX_DIM_3D)
666  const amrex::Real ymin = xyzmin[1];
667 #endif
668 
669  // Auxiliary constants
670 #if (defined WARPX_DIM_3D)
671  const amrex::Real onethird = 1._rt / 3._rt;
672  const amrex::Real onesixth = 1._rt / 6._rt;
673 #endif
674 
675  // Inverse of light speed squared
676  const amrex::Real invcsq = 1._rt / (PhysConst::c * PhysConst::c);
677 
678  // Arrays where D will be stored
679  amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
680  amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
681  amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
682 
683  // Loop over particles and deposit (Dx,Dy,Dz) into jx_fab, jy_fab and jz_fab
684  amrex::ParallelFor(np_to_depose, [=] AMREX_GPU_DEVICE (long ip)
685  {
686  // Inverse of Lorentz factor gamma
687  const amrex::Real invgam = 1._rt / std::sqrt(1._rt + uxp[ip] * uxp[ip] * invcsq
688  + uyp[ip] * uyp[ip] * invcsq
689  + uzp[ip] * uzp[ip] * invcsq);
690  // Product of particle charges and weights
691  amrex::Real wq = q * wp[ip];
692  if (do_ionization) wq *= ion_lev[ip];
693 
694  // Current particle positions (in physical units)
695  amrex::ParticleReal xp, yp, zp;
696  GetPosition(ip, xp, yp, zp);
697 
698  // Particle velocities
699  const amrex::Real vx = uxp[ip] * invgam;
700  const amrex::Real vy = uyp[ip] * invgam;
701  const amrex::Real vz = uzp[ip] * invgam;
702 
703  // Particle current densities
704 #if (defined WARPX_DIM_XZ)
705  const amrex::Real wqy = wq * vy * invvol;
706 #endif
707 
708  // Current and old particle positions in grid units
709  // Keep these double to avoid bug in single precision.
710  double const x_new = (xp - xmin) * dxi;
711  double const x_old = x_new - vx * dt * dxi;
712 #if (defined WARPX_DIM_3D)
713  // Keep these double to avoid bug in single precision.
714  double const y_new = (yp - ymin) * dyi;
715  double const y_old = y_new - vy * dt * dyi;
716 #endif
717  // Keep these double to avoid bug in single precision.
718  double const z_new = (zp - zmin) * dzi;
719  double const z_old = z_new - vz * dt * dzi;
720 
721  // Shape factor arrays for current and old positions (nodal)
722  // Keep these double to avoid bug in single precision.
723  double sx_new[depos_order+1] = {0.};
724  double sx_old[depos_order+1] = {0.};
725 #if (defined WARPX_DIM_3D)
726  // Keep these double to avoid bug in single precision.
727  double sy_new[depos_order+1] = {0.};
728  double sy_old[depos_order+1] = {0.};
729 #endif
730  // Keep these double to avoid bug in single precision.
731  double sz_new[depos_order+1] = {0.};
732  double sz_old[depos_order+1] = {0.};
733 
734  // Compute shape factors for current positions
735 
736  // i_new leftmost grid point in x that the particle touches
737  // sx_new shape factor along x for the centering of each current
738  Compute_shape_factor< depos_order > const compute_shape_factor;
739  const int i_new = compute_shape_factor(sx_new, x_new);
740 #if (defined WARPX_DIM_3D)
741  // j_new leftmost grid point in y that the particle touches
742  // sy_new shape factor along y for the centering of each current
743  const int j_new = compute_shape_factor(sy_new, y_new);
744 #endif
745  // k_new leftmost grid point in z that the particle touches
746  // sz_new shape factor along z for the centering of each current
747  const int k_new = compute_shape_factor(sz_new, z_new);
748 
749  // Compute shape factors for old positions
750 
751  // i_old leftmost grid point in x that the particle touches
752  // sx_old shape factor along x for the centering of each current
753  const int i_old = compute_shape_factor(sx_old, x_old);
754 #if (defined WARPX_DIM_3D)
755  // j_old leftmost grid point in y that the particle touches
756  // sy_old shape factor along y for the centering of each current
757  const int j_old = compute_shape_factor(sy_old, y_old);
758 #endif
759  // k_old leftmost grid point in z that the particle touches
760  // sz_old shape factor along z for the centering of each current
761  const int k_old = compute_shape_factor(sz_old, z_old);
762 
763  // Deposit current into jx_arr, jy_arr and jz_arr
764 #if (defined WARPX_DIM_XZ)
765 
766  for (int k=0; k<=depos_order; k++) {
767  for (int i=0; i<=depos_order; i++) {
768 
769  // Re-casting sx_new and sz_new from double to amrex::Real so that
770  // Atomic::Add has consistent types in its argument
771  amrex::Real const sxn_szn = sx_new[i] * sz_new[k];
772  amrex::Real const sxo_szn = sx_old[i] * sz_new[k];
773  amrex::Real const sxn_szo = sx_new[i] * sz_old[k];
774  amrex::Real const sxo_szo = sx_old[i] * sz_old[k];
775 
776  // Jx
777  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
778  wq * invvol * invdt * 0.5_rt * sxn_szn);
779 
780  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0),
781  - wq * invvol * invdt * 0.5_rt * sxo_szn);
782 
783  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0),
784  wq * invvol * invdt * 0.5_rt * sxn_szo);
785 
786  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
787  - wq * invvol * invdt * 0.5_rt * sxo_szo);
788 
789  // Jy
790  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
791  wqy * 0.25_rt * sxn_szn);
792 
793  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0),
794  wqy * 0.25_rt * sxn_szo);
795 
796  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0),
797  wqy * 0.25_rt * sxo_szn);
798 
799  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
800  wqy * 0.25_rt * sxo_szo);
801 
802  // Jz
803  amrex::Gpu::Atomic::Add(&jz_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
804  wq * invvol * invdt * 0.5_rt * sxn_szn);
805 
806  amrex::Gpu::Atomic::Add(&jz_arr(lo.x+i_new+i,lo.y+k_old+k,0,0),
807  - wq * invvol * invdt * 0.5_rt * sxn_szo);
808 
809  amrex::Gpu::Atomic::Add(&jz_arr(lo.x+i_old+i,lo.y+k_new+k,0,0),
810  wq * invvol * invdt * 0.5_rt * sxo_szn);
811 
812  amrex::Gpu::Atomic::Add(&jz_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
813  - wq * invvol * invdt * 0.5_rt * sxo_szo);
814  }
815  }
816 
817 #elif (defined WARPX_DIM_3D)
818 
819  for (int k=0; k<=depos_order; k++) {
820  for (int j=0; j<=depos_order; j++) {
821 
822  amrex::Real const syn_szn = sy_new[j] * sz_new[k];
823  amrex::Real const syo_szn = sy_old[j] * sz_new[k];
824  amrex::Real const syn_szo = sy_new[j] * sz_old[k];
825  amrex::Real const syo_szo = sy_old[j] * sz_old[k];
826 
827  for (int i=0; i<=depos_order; i++) {
828 
829  amrex::Real const sxn_syn_szn = sx_new[i] * syn_szn;
830  amrex::Real const sxo_syn_szn = sx_old[i] * syn_szn;
831  amrex::Real const sxn_syo_szn = sx_new[i] * syo_szn;
832  amrex::Real const sxo_syo_szn = sx_old[i] * syo_szn;
833  amrex::Real const sxn_syn_szo = sx_new[i] * syn_szo;
834  amrex::Real const sxo_syn_szo = sx_old[i] * syn_szo;
835  amrex::Real const sxn_syo_szo = sx_new[i] * syo_szo;
836  amrex::Real const sxo_syo_szo = sx_old[i] * syo_szo;
837 
838  // Jx
839  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k),
840  wq * invvol * invdt * onethird * sxn_syn_szn);
841 
842  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k),
843  - wq * invvol * invdt * onethird * sxo_syn_szn);
844 
845  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k),
846  wq * invvol * invdt * onesixth * sxn_syo_szn);
847 
848  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_old + j,lo.z + k_new + k),
849  - wq * invvol * invdt * onesixth * sxo_syo_szn);
850 
851  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k),
852  wq * invvol * invdt * onesixth * sxn_syn_szo);
853 
854  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k),
855  - wq * invvol * invdt * onesixth * sxo_syn_szo);
856 
857  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k),
858  wq * invvol * invdt * onethird * sxn_syo_szo);
859 
860  amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k),
861  - wq * invvol * invdt * onethird * sxo_syo_szo);
862 
863  // Jy
864  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k),
865  wq * invvol * invdt * onethird * sxn_syn_szn);
866 
867  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k),
868  - wq * invvol * invdt * onethird * sxn_syo_szn);
869 
870  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k),
871  wq * invvol * invdt * onesixth * sxo_syn_szn);
872 
873  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k),
874  - wq * invvol * invdt * onesixth * sxo_syo_szn);
875 
876  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k),
877  wq * invvol * invdt * onesixth * sxn_syn_szo);
878 
879  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k),
880  - wq * invvol * invdt * onesixth * sxn_syo_szo);
881 
882  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k),
883  wq * invvol * invdt * onethird * sxo_syn_szo);
884 
885  amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k),
886  - wq * invvol * invdt * onethird * sxo_syo_szo);
887 
888  // Jz
889  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k),
890  wq * invvol * invdt * onethird * sxn_syn_szn);
891 
892  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k),
893  - wq * invvol * invdt * onethird * sxn_syn_szo);
894 
895  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k),
896  wq * invvol * invdt * onesixth * sxo_syn_szn);
897 
898  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k),
899  - wq * invvol * invdt * onesixth * sxo_syn_szo);
900 
901  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k),
902  wq * invvol * invdt * onesixth * sxn_syo_szn);
903 
904  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k),
905  - wq * invvol * invdt * onesixth * sxn_syo_szo);
906 
907  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k),
908  wq * invvol * invdt * onethird * sxo_syo_szn);
909 
910  amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k),
911  - wq * invvol * invdt * onethird * sxo_syo_szo);
912  }
913  }
914  }
915 #endif
916  } );
917 #endif // #if !(defined WARPX_DIM_RZ)
918 }
919 #endif // CURRENTDEPOSITION_H_
amrex::GpuComplex< amrex::Real > Complex
Definition: WarpX_Complex.H:22
Definition: wp_parser.tab.cpp:115
void doVayDepositionShapeN(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 *const ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, const long np_to_depose, const amrex::Real dt, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const long n_rz_azimuthal_modes)
Vay current deposition (Vay et al, 2013) for thread thread_num: deposit D in real space and store the...
Definition: CurrentDeposition.H:615
int dx
Definition: compute_domain.py:35
Definition: ShapeFactors.H:118
void doDepositionShapeN(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 *const ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, const long np_to_depose, const amrex::Real dt, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const long n_rz_azimuthal_modes)
Current Deposition for thread thread_num /param GetPosition : A functor for returning the particle po...
Definition: CurrentDeposition.H:41
Definition: ShapeFactors.H:22
i
Definition: check_interp_points_and_weights.py:171
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel...
Definition: GetAndSetPosition.H:25
Definition: PML.H:52
void doEsirkepovDepositionShapeN(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, const amrex::Array4< amrex::Real > &Jx_arr, const amrex::Array4< amrex::Real > &Jy_arr, const amrex::Array4< amrex::Real > &Jz_arr, const long np_to_depose, const amrex::Real dt, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > xyzmin, const amrex::Dim3 lo, const amrex::Real q, const long n_rz_azimuthal_modes)
Esirkepov Current Deposition for thread thread_num.
Definition: CurrentDeposition.H:306