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 
11 #include "Parallelization/KernelTimer.H"
12 #include "Particles/Pusher/GetAndSetPosition.H"
13 #include "Particles/ShapeFactors.H"
14 #include "Utils/WarpXAlgorithmSelection.H"
15 #include "Utils/WarpXConst.H"
16 #ifdef WARPX_DIM_RZ
17 # include "Utils/WarpX_Complex.H"
18 #endif
19 
20 #include <AMReX.H>
21 #include <AMReX_Arena.H>
22 #include <AMReX_Array4.H>
23 #include <AMReX_REAL.H>
24 
25 using namespace amrex::literals;
26 
52 template <int depos_order>
53 void doDepositionShapeN(const GetParticlePosition& GetPosition,
54  const amrex::ParticleReal * const wp,
55  const amrex::ParticleReal * const uxp,
56  const amrex::ParticleReal * const uyp,
57  const amrex::ParticleReal * const uzp,
58  const int * const ion_lev,
59  amrex::FArrayBox& jx_fab,
60  amrex::FArrayBox& jy_fab,
61  amrex::FArrayBox& jz_fab,
62  const long np_to_depose,
63  const amrex::Real relative_time,
64  const std::array<amrex::Real,3>& dx,
65  const std::array<amrex::Real,3>& xyzmin,
66  const amrex::Dim3 lo,
67  const amrex::Real q,
68  const int n_rz_azimuthal_modes,
69  amrex::Real* cost,
70  const long load_balance_costs_update_algo)
71 {
72 #if !defined(WARPX_DIM_RZ)
73  amrex::ignore_unused(n_rz_azimuthal_modes);
74 #endif
75 
76 #if !defined(AMREX_USE_GPU)
77  amrex::ignore_unused(cost, load_balance_costs_update_algo);
78 #endif
79 
80  // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
81  // (do_ionization=1)
82  const bool do_ionization = ion_lev;
83  const amrex::Real dzi = 1.0_rt/dx[2];
84 #if defined(WARPX_DIM_1D_Z)
85  const amrex::Real invvol = dzi;
86 #endif
87 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
88  const amrex::Real dxi = 1.0_rt/dx[0];
89  const amrex::Real invvol = dxi*dzi;
90 #elif defined(WARPX_DIM_3D)
91  const amrex::Real dxi = 1.0_rt/dx[0];
92  const amrex::Real dyi = 1.0_rt/dx[1];
93  const amrex::Real invvol = dxi*dyi*dzi;
94 #endif
95 
96 #if (AMREX_SPACEDIM >= 2)
97  const amrex::Real xmin = xyzmin[0];
98 #endif
99 #if defined(WARPX_DIM_3D)
100  const amrex::Real ymin = xyzmin[1];
101 #endif
102  const amrex::Real zmin = xyzmin[2];
103 
104  const amrex::Real clightsq = 1.0_rt/PhysConst::c/PhysConst::c;
105 
106  amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
107  amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
108  amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
109  amrex::IntVect const jx_type = jx_fab.box().type();
110  amrex::IntVect const jy_type = jy_fab.box().type();
111  amrex::IntVect const jz_type = jz_fab.box().type();
112 
113  constexpr int zdir = WARPX_ZINDEX;
114  constexpr int NODE = amrex::IndexType::NODE;
115  constexpr int CELL = amrex::IndexType::CELL;
116 
117  // Loop over particles and deposit into jx_fab, jy_fab and jz_fab
118 #if defined(WARPX_USE_GPUCLOCK)
119  amrex::Real* cost_real = nullptr;
120  if( load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::GpuClock) {
121  cost_real = (amrex::Real *) amrex::The_Managed_Arena()->alloc(sizeof(amrex::Real));
122  *cost_real = 0._rt;
123  }
124 #endif
125  amrex::ParallelFor(
126  np_to_depose,
127  [=] AMREX_GPU_DEVICE (long ip) {
128 #if defined(WARPX_USE_GPUCLOCK)
129  KernelTimer kernelTimer(cost && load_balance_costs_update_algo
130  == LoadBalanceCostsUpdateAlgo::GpuClock, cost_real);
131 #endif
132 
133  // --- Get particle quantities
134  const amrex::Real gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
135  + uyp[ip]*uyp[ip]*clightsq
136  + uzp[ip]*uzp[ip]*clightsq);
137  amrex::Real wq = q*wp[ip];
138  if (do_ionization){
139  wq *= ion_lev[ip];
140  }
141 
142  amrex::ParticleReal xp, yp, zp;
143  GetPosition(ip, xp, yp, zp);
144 
145  const amrex::Real vx = uxp[ip]*gaminv;
146  const amrex::Real vy = uyp[ip]*gaminv;
147  const amrex::Real vz = uzp[ip]*gaminv;
148  // wqx, wqy wqz are particle current in each direction
149 #if defined(WARPX_DIM_RZ)
150  // In RZ, wqx is actually wqr, and wqy is wqtheta
151  // Convert to cylinderical at the mid point
152  const amrex::Real xpmid = xp + relative_time*vx;
153  const amrex::Real ypmid = yp + relative_time*vy;
154  const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
155  amrex::Real costheta;
156  amrex::Real sintheta;
157  if (rpmid > 0._rt) {
158  costheta = xpmid/rpmid;
159  sintheta = ypmid/rpmid;
160  } else {
161  costheta = 1._rt;
162  sintheta = 0._rt;
163  }
164  const Complex xy0 = Complex{costheta, sintheta};
165  const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
166  const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
167 #else
168  const amrex::Real wqx = wq*invvol*vx;
169  const amrex::Real wqy = wq*invvol*vy;
170 #endif
171  const amrex::Real wqz = wq*invvol*vz;
172 
173  // --- Compute shape factors
174  Compute_shape_factor< depos_order > const compute_shape_factor;
175 #if (AMREX_SPACEDIM >= 2)
176  // x direction
177  // Get particle position after 1/2 push back in position
178 #if defined(WARPX_DIM_RZ)
179  // Keep these double to avoid bug in single precision
180  const double xmid = (rpmid - xmin)*dxi;
181 #else
182  const double xmid = ((xp - xmin) + relative_time*vx)*dxi;
183 #endif
184  // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
185  // sx_j[xyz] shape factor along x for the centering of each current
186  // There are only two possible centerings, node or cell centered, so at most only two shape factor
187  // arrays will be needed.
188  // Keep these double to avoid bug in single precision
189  double sx_node[depos_order + 1] = {0.};
190  double sx_cell[depos_order + 1] = {0.};
191  int j_node = 0;
192  int j_cell = 0;
193  if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) {
194  j_node = compute_shape_factor(sx_node, xmid);
195  }
196  if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
197  j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
198  }
199 
200  amrex::Real sx_jx[depos_order + 1] = {0._rt};
201  amrex::Real sx_jy[depos_order + 1] = {0._rt};
202  amrex::Real sx_jz[depos_order + 1] = {0._rt};
203  for (int ix=0; ix<=depos_order; ix++)
204  {
205  sx_jx[ix] = ((jx_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
206  sx_jy[ix] = ((jy_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
207  sx_jz[ix] = ((jz_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
208  }
209 
210  int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell);
211  int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell);
212  int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell);
213 #endif //AMREX_SPACEDIM >= 2
214 
215 #if defined(WARPX_DIM_3D)
216  // y direction
217  // Keep these double to avoid bug in single precision
218  const double ymid = ((yp - ymin) + relative_time*vy)*dyi;
219  double sy_node[depos_order + 1] = {0.};
220  double sy_cell[depos_order + 1] = {0.};
221  int k_node = 0;
222  int k_cell = 0;
223  if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) {
224  k_node = compute_shape_factor(sy_node, ymid);
225  }
226  if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
227  k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
228  }
229  amrex::Real sy_jx[depos_order + 1] = {0._rt};
230  amrex::Real sy_jy[depos_order + 1] = {0._rt};
231  amrex::Real sy_jz[depos_order + 1] = {0._rt};
232  for (int iy=0; iy<=depos_order; iy++)
233  {
234  sy_jx[iy] = ((jx_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
235  sy_jy[iy] = ((jy_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
236  sy_jz[iy] = ((jz_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
237  }
238  int const k_jx = ((jx_type[1] == NODE) ? k_node : k_cell);
239  int const k_jy = ((jy_type[1] == NODE) ? k_node : k_cell);
240  int const k_jz = ((jz_type[1] == NODE) ? k_node : k_cell);
241 #endif
242 
243  // z direction
244  // Keep these double to avoid bug in single precision
245  const double zmid = ((zp - zmin) + relative_time*vz)*dzi;
246  double sz_node[depos_order + 1] = {0.};
247  double sz_cell[depos_order + 1] = {0.};
248  int l_node = 0;
249  int l_cell = 0;
250  if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) {
251  l_node = compute_shape_factor(sz_node, zmid);
252  }
253  if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
254  l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
255  }
256  amrex::Real sz_jx[depos_order + 1] = {0._rt};
257  amrex::Real sz_jy[depos_order + 1] = {0._rt};
258  amrex::Real sz_jz[depos_order + 1] = {0._rt};
259  for (int iz=0; iz<=depos_order; iz++)
260  {
261  sz_jx[iz] = ((jx_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
262  sz_jy[iz] = ((jy_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
263  sz_jz[iz] = ((jz_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
264  }
265  int const l_jx = ((jx_type[zdir] == NODE) ? l_node : l_cell);
266  int const l_jy = ((jy_type[zdir] == NODE) ? l_node : l_cell);
267  int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell);
268 
269  // Deposit current into jx_arr, jy_arr and jz_arr
270 #if defined(WARPX_DIM_1D_Z)
271  for (int iz=0; iz<=depos_order; iz++){
272  amrex::Gpu::Atomic::AddNoRet(
273  &jx_arr(lo.x+l_jx+iz, 0, 0, 0),
274  sz_jx[iz]*wqx);
275  amrex::Gpu::Atomic::AddNoRet(
276  &jy_arr(lo.x+l_jy+iz, 0, 0, 0),
277  sz_jy[iz]*wqy);
278  amrex::Gpu::Atomic::AddNoRet(
279  &jz_arr(lo.x+l_jz+iz, 0, 0, 0),
280  sz_jz[iz]*wqz);
281  }
282 #endif
283 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
284  for (int iz=0; iz<=depos_order; iz++){
285  for (int ix=0; ix<=depos_order; ix++){
286  amrex::Gpu::Atomic::AddNoRet(
287  &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
288  sx_jx[ix]*sz_jx[iz]*wqx);
289  amrex::Gpu::Atomic::AddNoRet(
290  &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
291  sx_jy[ix]*sz_jy[iz]*wqy);
292  amrex::Gpu::Atomic::AddNoRet(
293  &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
294  sx_jz[ix]*sz_jz[iz]*wqz);
295 #if defined(WARPX_DIM_RZ)
296  Complex xy = xy0; // Note that xy is equal to e^{i m theta}
297  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
298  // The factor 2 on the weighting comes from the normalization of the modes
299  amrex::Gpu::Atomic::AddNoRet( &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 2*imode-1), 2._rt*sx_jx[ix]*sz_jx[iz]*wqx*xy.real());
300  amrex::Gpu::Atomic::AddNoRet( &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 2*imode ), 2._rt*sx_jx[ix]*sz_jx[iz]*wqx*xy.imag());
301  amrex::Gpu::Atomic::AddNoRet( &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 2*imode-1), 2._rt*sx_jy[ix]*sz_jy[iz]*wqy*xy.real());
302  amrex::Gpu::Atomic::AddNoRet( &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 2*imode ), 2._rt*sx_jy[ix]*sz_jy[iz]*wqy*xy.imag());
303  amrex::Gpu::Atomic::AddNoRet( &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 2*imode-1), 2._rt*sx_jz[ix]*sz_jz[iz]*wqz*xy.real());
304  amrex::Gpu::Atomic::AddNoRet( &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 2*imode ), 2._rt*sx_jz[ix]*sz_jz[iz]*wqz*xy.imag());
305  xy = xy*xy0;
306  }
307 #endif
308  }
309  }
310 #elif defined(WARPX_DIM_3D)
311  for (int iz=0; iz<=depos_order; iz++){
312  for (int iy=0; iy<=depos_order; iy++){
313  for (int ix=0; ix<=depos_order; ix++){
314  amrex::Gpu::Atomic::AddNoRet(
315  &jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz),
316  sx_jx[ix]*sy_jx[iy]*sz_jx[iz]*wqx);
317  amrex::Gpu::Atomic::AddNoRet(
318  &jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz),
319  sx_jy[ix]*sy_jy[iy]*sz_jy[iz]*wqy);
320  amrex::Gpu::Atomic::AddNoRet(
321  &jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz),
322  sx_jz[ix]*sy_jz[iy]*sz_jz[iz]*wqz);
323  }
324  }
325  }
326 #endif
327  }
328  );
329 #if defined(WARPX_USE_GPUCLOCK)
330  if( load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::GpuClock) {
331  amrex::Gpu::streamSynchronize();
332  *cost += *cost_real;
333  amrex::The_Managed_Arena()->free(cost_real);
334  }
335 #endif
336 }
337 
364 template <int depos_order>
365 void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition,
366  const amrex::ParticleReal * const wp,
367  const amrex::ParticleReal * const uxp,
368  const amrex::ParticleReal * const uyp,
369  const amrex::ParticleReal * const uzp,
370  const int * const ion_lev,
371  const amrex::Array4<amrex::Real>& Jx_arr,
372  const amrex::Array4<amrex::Real>& Jy_arr,
373  const amrex::Array4<amrex::Real>& Jz_arr,
374  const long np_to_depose,
375  const amrex::Real dt,
376  const amrex::Real relative_time,
377  const std::array<amrex::Real,3>& dx,
378  const std::array<amrex::Real, 3> xyzmin,
379  const amrex::Dim3 lo,
380  const amrex::Real q,
381  const int n_rz_azimuthal_modes,
382  amrex::Real * const cost,
383  const long load_balance_costs_update_algo)
384 {
385  using namespace amrex;
386 #if !defined(WARPX_DIM_RZ)
387  ignore_unused(n_rz_azimuthal_modes);
388 #endif
389 
390 #if !defined(AMREX_USE_GPU)
391  amrex::ignore_unused(cost, load_balance_costs_update_algo);
392 #endif
393 
394  // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
395  // (do_ionization=1)
396  bool const do_ionization = ion_lev;
397 #if !defined(WARPX_DIM_1D_Z)
398  Real const dxi = 1.0_rt / dx[0];
399 #endif
400 #if !defined(WARPX_DIM_1D_Z)
401  Real const xmin = xyzmin[0];
402 #endif
403 #if defined(WARPX_DIM_3D)
404  Real const dyi = 1.0_rt / dx[1];
405  Real const ymin = xyzmin[1];
406 #endif
407  Real const dzi = 1.0_rt / dx[2];
408  Real const zmin = xyzmin[2];
409 
410 #if defined(WARPX_DIM_3D)
411  Real const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]);
412  Real const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]);
413  Real const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]);
414 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
415  Real const invdtdx = 1.0_rt / (dt*dx[2]);
416  Real const invdtdz = 1.0_rt / (dt*dx[0]);
417  Real const invvol = 1.0_rt / (dx[0]*dx[2]);
418 #elif defined(WARPX_DIM_1D_Z)
419  Real const invdtdz = 1.0_rt / (dt*dx[0]);
420  Real const invvol = 1.0_rt / (dx[2]);
421 #endif
422 
423 #if defined(WARPX_DIM_RZ)
424  Complex const I = Complex{0._rt, 1._rt};
425 #endif
426 
427  Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c );
428 #if !defined(WARPX_DIM_1D_Z)
429  Real constexpr one_third = 1.0_rt / 3.0_rt;
430  Real constexpr one_sixth = 1.0_rt / 6.0_rt;
431 #endif
432 
433  // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
434 #if defined(WARPX_USE_GPUCLOCK)
435  amrex::Real* cost_real = nullptr;
436  if( load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::GpuClock) {
437  cost_real = (amrex::Real *) amrex::The_Managed_Arena()->alloc(sizeof(amrex::Real));
438  *cost_real = 0._rt;
439  }
440 #endif
441  amrex::ParallelFor(
442  np_to_depose,
443  [=] AMREX_GPU_DEVICE (long const ip) {
444 #if defined(WARPX_USE_GPUCLOCK)
445  KernelTimer kernelTimer(cost && load_balance_costs_update_algo
446  == LoadBalanceCostsUpdateAlgo::GpuClock, cost_real);
447 #endif
448 
449  // --- Get particle quantities
450  Real const gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
451  + uyp[ip]*uyp[ip]*clightsq
452  + uzp[ip]*uzp[ip]*clightsq);
453 
454  // wqx, wqy wqz are particle current in each direction
455  Real wq = q*wp[ip];
456  if (do_ionization){
457  wq *= ion_lev[ip];
458  }
459 
460  ParticleReal xp, yp, zp;
461  GetPosition(ip, xp, yp, zp);
462 
463 #if !defined(WARPX_DIM_1D_Z)
464  Real const wqx = wq*invdtdx;
465 #endif
466 #if defined(WARPX_DIM_3D)
467  Real const wqy = wq*invdtdy;
468 #endif
469  Real const wqz = wq*invdtdz;
470 
471  // computes current and old position in grid units
472 #if defined(WARPX_DIM_RZ)
473  Real const xp_new = xp + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv;
474  Real const yp_new = yp + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv;
475  Real const xp_mid = xp_new - 0.5_rt*dt*uxp[ip]*gaminv;
476  Real const yp_mid = yp_new - 0.5_rt*dt*uyp[ip]*gaminv;
477  Real const xp_old = xp_new - dt*uxp[ip]*gaminv;
478  Real const yp_old = yp_new - dt*uyp[ip]*gaminv;
479  Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
480  Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
481  Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
482  Real costheta_new, sintheta_new;
483  if (rp_new > 0._rt) {
484  costheta_new = xp_new/rp_new;
485  sintheta_new = yp_new/rp_new;
486  } else {
487  costheta_new = 1._rt;
488  sintheta_new = 0._rt;
489  }
490  amrex::Real costheta_mid, sintheta_mid;
491  if (rp_mid > 0._rt) {
492  costheta_mid = xp_mid/rp_mid;
493  sintheta_mid = yp_mid/rp_mid;
494  } else {
495  costheta_mid = 1._rt;
496  sintheta_mid = 0._rt;
497  }
498  amrex::Real costheta_old, sintheta_old;
499  if (rp_old > 0._rt) {
500  costheta_old = xp_old/rp_old;
501  sintheta_old = yp_old/rp_old;
502  } else {
503  costheta_old = 1._rt;
504  sintheta_old = 0._rt;
505  }
506  const Complex xy_new0 = Complex{costheta_new, sintheta_new};
507  const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
508  const Complex xy_old0 = Complex{costheta_old, sintheta_old};
509  // Keep these double to avoid bug in single precision
510  double const x_new = (rp_new - xmin)*dxi;
511  double const x_old = (rp_old - xmin)*dxi;
512 #else
513 #if !defined(WARPX_DIM_1D_Z)
514  // Keep these double to avoid bug in single precision
515  double const x_new = (xp - xmin + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv)*dxi;
516  double const x_old = x_new - dt*dxi*uxp[ip]*gaminv;
517 #endif
518 #endif
519 #if defined(WARPX_DIM_3D)
520  // Keep these double to avoid bug in single precision
521  double const y_new = (yp - ymin + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv)*dyi;
522  double const y_old = y_new - dt*dyi*uyp[ip]*gaminv;
523 #endif
524  // Keep these double to avoid bug in single precision
525  double const z_new = (zp - zmin + (relative_time + 0.5_rt*dt)*uzp[ip]*gaminv)*dzi;
526  double const z_old = z_new - dt*dzi*uzp[ip]*gaminv;
527 
528 #if defined(WARPX_DIM_RZ)
529  Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
530 #elif defined(WARPX_DIM_XZ)
531  Real const vy = uyp[ip]*gaminv;
532 #elif defined(WARPX_DIM_1D_Z)
533  Real const vx = uxp[ip]*gaminv;
534  Real const vy = uyp[ip]*gaminv;
535 #endif
536 
537  // Shape factor arrays
538  // Note that there are extra values above and below
539  // to possibly hold the factor for the old particle
540  // which can be at a different grid location.
541  // Keep these double to avoid bug in single precision
542 #if !defined(WARPX_DIM_1D_Z)
543  double sx_new[depos_order + 3] = {0.};
544  double sx_old[depos_order + 3] = {0.};
545 #endif
546 #if defined(WARPX_DIM_3D)
547  // Keep these double to avoid bug in single precision
548  double sy_new[depos_order + 3] = {0.};
549  double sy_old[depos_order + 3] = {0.};
550 #endif
551  // Keep these double to avoid bug in single precision
552  double sz_new[depos_order + 3] = {0.};
553  double sz_old[depos_order + 3] = {0.};
554 
555  // --- Compute shape factors
556  // Compute shape factors for position as they are now and at old positions
557  // [ijk]_new: leftmost grid point that the particle touches
558  Compute_shape_factor< depos_order > compute_shape_factor;
559  Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor;
560 
561 #if !defined(WARPX_DIM_1D_Z)
562  const int i_new = compute_shape_factor(sx_new+1, x_new);
563  const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
564 #endif
565 #if defined(WARPX_DIM_3D)
566  const int j_new = compute_shape_factor(sy_new+1, y_new);
567  const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
568 #endif
569  const int k_new = compute_shape_factor(sz_new+1, z_new);
570  const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new);
571 
572  // computes min/max positions of current contributions
573 #if !defined(WARPX_DIM_1D_Z)
574  int dil = 1, diu = 1;
575  if (i_old < i_new) dil = 0;
576  if (i_old > i_new) diu = 0;
577 #endif
578 #if defined(WARPX_DIM_3D)
579  int djl = 1, dju = 1;
580  if (j_old < j_new) djl = 0;
581  if (j_old > j_new) dju = 0;
582 #endif
583  int dkl = 1, dku = 1;
584  if (k_old < k_new) dkl = 0;
585  if (k_old > k_new) dku = 0;
586 
587 #if defined(WARPX_DIM_3D)
588 
589  for (int k=dkl; k<=depos_order+2-dku; k++) {
590  for (int j=djl; j<=depos_order+2-dju; j++) {
591  amrex::Real sdxi = 0._rt;
592  for (int i=dil; i<=depos_order+1-diu; i++) {
593  sdxi += wqx*(sx_old[i] - sx_new[i])*(
594  one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k])
595  +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k]));
596  amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdxi);
597  }
598  }
599  }
600  for (int k=dkl; k<=depos_order+2-dku; k++) {
601  for (int i=dil; i<=depos_order+2-diu; i++) {
602  amrex::Real sdyj = 0._rt;
603  for (int j=djl; j<=depos_order+1-dju; j++) {
604  sdyj += wqy*(sy_old[j] - sy_new[j])*(
605  one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k])
606  +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k]));
607  amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdyj);
608  }
609  }
610  }
611  for (int j=djl; j<=depos_order+2-dju; j++) {
612  for (int i=dil; i<=depos_order+2-diu; i++) {
613  amrex::Real sdzk = 0._rt;
614  for (int k=dkl; k<=depos_order+1-dku; k++) {
615  sdzk += wqz*(sz_old[k] - sz_new[k])*(
616  one_third*(sx_new[i]*sy_new[j] + sx_old[i]*sy_old[j])
617  +one_sixth*(sx_new[i]*sy_old[j] + sx_old[i]*sy_new[j]));
618  amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdzk);
619  }
620  }
621  }
622 
623 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
624 
625  for (int k=dkl; k<=depos_order+2-dku; k++) {
626  amrex::Real sdxi = 0._rt;
627  for (int i=dil; i<=depos_order+1-diu; i++) {
628  sdxi += wqx*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]);
629  amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi);
630 #if defined(WARPX_DIM_RZ)
631  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
632  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
633  // The factor 2 comes from the normalization of the modes
634  const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
635  amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real());
636  amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
637  xy_mid = xy_mid*xy_mid0;
638  }
639 #endif
640  }
641  }
642  for (int k=dkl; k<=depos_order+2-dku; k++) {
643  for (int i=dil; i<=depos_order+2-diu; i++) {
644  Real const sdyj = wq*vy*invvol*(
645  one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k])
646  +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k]));
647  amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
648 #if defined(WARPX_DIM_RZ)
649  Complex xy_new = xy_new0;
650  Complex xy_mid = xy_mid0;
651  Complex xy_old = xy_old0;
652  // Throughout the following loop, xy_ takes the value e^{i m theta_}
653  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
654  // The factor 2 comes from the normalization of the modes
655  // The minus sign comes from the different convention with respect to Davidson et al.
656  const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode
657  *(Complex(sx_new[i]*sz_new[k], 0._rt)*(xy_new - xy_mid)
658  + Complex(sx_old[i]*sz_old[k], 0._rt)*(xy_mid - xy_old));
659  amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real());
660  amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
661  xy_new = xy_new*xy_new0;
662  xy_mid = xy_mid*xy_mid0;
663  xy_old = xy_old*xy_old0;
664  }
665 #endif
666  }
667  }
668  for (int i=dil; i<=depos_order+2-diu; i++) {
669  Real sdzk = 0._rt;
670  for (int k=dkl; k<=depos_order+1-dku; k++) {
671  sdzk += wqz*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]);
672  amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
673 #if defined(WARPX_DIM_RZ)
674  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
675  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
676  // The factor 2 comes from the normalization of the modes
677  const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
678  amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real());
679  amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
680  xy_mid = xy_mid*xy_mid0;
681  }
682 #endif
683  }
684  }
685 #elif defined(WARPX_DIM_1D_Z)
686 
687  for (int k=dkl; k<=depos_order+2-dku; k++) {
688  amrex::Real sdxi = 0._rt;
689  sdxi += wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
690  amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+k_new-1+k, 0, 0, 0), sdxi);
691  }
692  for (int k=dkl; k<=depos_order+2-dku; k++) {
693  amrex::Real sdyj = 0._rt;
694  sdyj += wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
695  amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+k_new-1+k, 0, 0, 0), sdyj);
696  }
697  for (int k=dkl; k<=depos_order+1-dku; k++) {
698  amrex::Real sdzk = 0._rt;
699  sdzk += wqz*(sz_old[k] - sz_new[k]);
700  amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk);
701  }
702 #endif
703  }
704  );
705 #if defined(WARPX_USE_GPUCLOCK)
706  if( load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::GpuClock) {
707  amrex::Gpu::streamSynchronize();
708  *cost += *cost_real;
709  amrex::The_Managed_Arena()->free(cost_real);
710  }
711 #endif
712 }
713 
743 template <int depos_order>
744 void doVayDepositionShapeN (const GetParticlePosition& GetPosition,
745  const amrex::ParticleReal* const wp,
746  const amrex::ParticleReal* const uxp,
747  const amrex::ParticleReal* const uyp,
748  const amrex::ParticleReal* const uzp,
749  const int* const ion_lev,
750  amrex::FArrayBox& jx_fab,
751  amrex::FArrayBox& jy_fab,
752  amrex::FArrayBox& jz_fab,
753  const long np_to_depose,
754  const amrex::Real dt,
755  const amrex::Real relative_time,
756  const std::array<amrex::Real,3>& dx,
757  const std::array<amrex::Real,3>& xyzmin,
758  const amrex::Dim3 lo,
759  const amrex::Real q,
760  const int n_rz_azimuthal_modes,
761  amrex::Real* cost,
762  const long load_balance_costs_update_algo)
763 {
764 #if defined(WARPX_DIM_RZ)
765  amrex::ignore_unused(GetPosition,
766  wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab,
767  np_to_depose, dt, relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes);
768  amrex::Abort("Vay deposition not implemented in RZ geometry");
769 #endif
770 
771 #if defined(WARPX_DIM_1D_Z)
772  amrex::ignore_unused(GetPosition,
773  wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab,
774  np_to_depose, dt, relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes);
775  amrex::Abort("Vay deposition not implemented in cartesian 1D geometry");
776 #endif
777 
778 #if !defined(AMREX_USE_GPU)
779  amrex::ignore_unused(cost, load_balance_costs_update_algo);
780 #endif
781 
782 #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z)
783  amrex::ignore_unused(n_rz_azimuthal_modes);
784 
785  // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1
786  const bool do_ionization = ion_lev;
787 
788  // Inverse cell volume in each direction
789  const amrex::Real dxi = 1._rt / dx[0];
790  const amrex::Real dzi = 1._rt / dx[2];
791 #if defined(WARPX_DIM_3D)
792  const amrex::Real dyi = 1._rt / dx[1];
793 #endif
794 
795  // Inverse of time step
796  const amrex::Real invdt = 1._rt / dt;
797 
798  // Total inverse cell volume
799 #if defined(WARPX_DIM_XZ)
800  const amrex::Real invvol = dxi * dzi;
801 #elif defined(WARPX_DIM_3D)
802  const amrex::Real invvol = dxi * dyi * dzi;
803 #endif
804 
805  // Lower bound of physical domain in each direction
806  const amrex::Real xmin = xyzmin[0];
807  const amrex::Real zmin = xyzmin[2];
808 #if defined(WARPX_DIM_3D)
809  const amrex::Real ymin = xyzmin[1];
810 #endif
811 
812  // Auxiliary constants
813 #if defined(WARPX_DIM_3D)
814  const amrex::Real onethird = 1._rt / 3._rt;
815  const amrex::Real onesixth = 1._rt / 6._rt;
816 #endif
817 
818  // Inverse of light speed squared
819  const amrex::Real invcsq = 1._rt / (PhysConst::c * PhysConst::c);
820 
821  // Arrays where D will be stored
822  amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
823  amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
824  amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
825 
826  // Loop over particles and deposit (Dx,Dy,Dz) into jx_fab, jy_fab and jz_fab
827 #if defined(WARPX_USE_GPUCLOCK)
828  amrex::Real* cost_real = nullptr;
829  if( load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::GpuClock) {
830  cost_real = (amrex::Real *) amrex::The_Managed_Arena()->alloc(sizeof(amrex::Real));
831  *cost_real = 0._rt;
832  }
833 #endif
834  amrex::ParallelFor(np_to_depose, [=] AMREX_GPU_DEVICE (long ip)
835  {
836 #if defined(WARPX_USE_GPUCLOCK)
837  KernelTimer kernelTimer(cost && load_balance_costs_update_algo
838  == LoadBalanceCostsUpdateAlgo::GpuClock, cost_real);
839 #endif
840 
841  // Inverse of Lorentz factor gamma
842  const amrex::Real invgam = 1._rt / std::sqrt(1._rt + uxp[ip] * uxp[ip] * invcsq
843  + uyp[ip] * uyp[ip] * invcsq
844  + uzp[ip] * uzp[ip] * invcsq);
845  // Product of particle charges and weights
846  amrex::Real wq = q * wp[ip];
847  if (do_ionization) wq *= ion_lev[ip];
848 
849  // Current particle positions (in physical units)
850  amrex::ParticleReal xp, yp, zp;
851  GetPosition(ip, xp, yp, zp);
852 
853  // Particle velocities
854  const amrex::Real vx = uxp[ip] * invgam;
855  const amrex::Real vy = uyp[ip] * invgam;
856  const amrex::Real vz = uzp[ip] * invgam;
857 
858  // Modify the particle position to match the time of the deposition
859  xp += relative_time * vx;
860  yp += relative_time * vy;
861  zp += relative_time * vz;
862 
863  // Particle current densities
864 #if defined(WARPX_DIM_XZ)
865  const amrex::Real wqy = wq * vy * invvol;
866 #endif
867 
868  // Current and old particle positions in grid units
869  // Keep these double to avoid bug in single precision.
870  double const x_new = (xp - xmin + 0.5_rt*dt*vx) * dxi;
871  double const x_old = (xp - xmin - 0.5_rt*dt*vx) * dxi;
872 #if defined(WARPX_DIM_3D)
873  // Keep these double to avoid bug in single precision.
874  double const y_new = (yp - ymin + 0.5_rt*dt*vy) * dyi;
875  double const y_old = (yp - ymin - 0.5_rt*dt*vy) * dyi;
876 #endif
877  // Keep these double to avoid bug in single precision.
878  double const z_new = (zp - zmin + 0.5_rt*dt*vz) * dzi;
879  double const z_old = (zp - zmin - 0.5_rt*dt*vz) * dzi;
880 
881  // Shape factor arrays for current and old positions (nodal)
882  // Keep these double to avoid bug in single precision.
883  double sx_new[depos_order+1] = {0.};
884  double sx_old[depos_order+1] = {0.};
885 #if defined(WARPX_DIM_3D)
886  // Keep these double to avoid bug in single precision.
887  double sy_new[depos_order+1] = {0.};
888  double sy_old[depos_order+1] = {0.};
889 #endif
890  // Keep these double to avoid bug in single precision.
891  double sz_new[depos_order+1] = {0.};
892  double sz_old[depos_order+1] = {0.};
893 
894  // Compute shape factors for current positions
895 
896  // i_new leftmost grid point in x that the particle touches
897  // sx_new shape factor along x for the centering of each current
898  Compute_shape_factor< depos_order > const compute_shape_factor;
899  const int i_new = compute_shape_factor(sx_new, x_new);
900 #if defined(WARPX_DIM_3D)
901  // j_new leftmost grid point in y that the particle touches
902  // sy_new shape factor along y for the centering of each current
903  const int j_new = compute_shape_factor(sy_new, y_new);
904 #endif
905  // k_new leftmost grid point in z that the particle touches
906  // sz_new shape factor along z for the centering of each current
907  const int k_new = compute_shape_factor(sz_new, z_new);
908 
909  // Compute shape factors for old positions
910 
911  // i_old leftmost grid point in x that the particle touches
912  // sx_old shape factor along x for the centering of each current
913  const int i_old = compute_shape_factor(sx_old, x_old);
914 #if defined(WARPX_DIM_3D)
915  // j_old leftmost grid point in y that the particle touches
916  // sy_old shape factor along y for the centering of each current
917  const int j_old = compute_shape_factor(sy_old, y_old);
918 #endif
919  // k_old leftmost grid point in z that the particle touches
920  // sz_old shape factor along z for the centering of each current
921  const int k_old = compute_shape_factor(sz_old, z_old);
922 
923  // Deposit current into jx_arr, jy_arr and jz_arr
924 #if defined(WARPX_DIM_XZ)
925 
926  for (int k=0; k<=depos_order; k++) {
927  for (int i=0; i<=depos_order; i++) {
928 
929  // Re-casting sx_new and sz_new from double to amrex::Real so that
930  // Atomic::Add has consistent types in its argument
931  auto const sxn_szn = static_cast<amrex::Real>(sx_new[i] * sz_new[k]);
932  auto const sxo_szn = static_cast<amrex::Real>(sx_old[i] * sz_new[k]);
933  auto const sxn_szo = static_cast<amrex::Real>(sx_new[i] * sz_old[k]);
934  auto const sxo_szo = static_cast<amrex::Real>(sx_old[i] * sz_old[k]);
935 
936  // Jx
937  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
938  wq * invvol * invdt * 0.5_rt * sxn_szn);
939 
940  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0),
941  - wq * invvol * invdt * 0.5_rt * sxo_szn);
942 
943  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0),
944  wq * invvol * invdt * 0.5_rt * sxn_szo);
945 
946  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
947  - wq * invvol * invdt * 0.5_rt * sxo_szo);
948 
949  // Jy
950  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
951  wqy * 0.25_rt * sxn_szn);
952 
953  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0),
954  wqy * 0.25_rt * sxn_szo);
955 
956  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0),
957  wqy * 0.25_rt * sxo_szn);
958 
959  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
960  wqy * 0.25_rt * sxo_szo);
961 
962  // Jz
963  amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
964  wq * invvol * invdt * 0.5_rt * sxn_szn);
965 
966  amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x+i_new+i,lo.y+k_old+k,0,0),
967  - wq * invvol * invdt * 0.5_rt * sxn_szo);
968 
969  amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x+i_old+i,lo.y+k_new+k,0,0),
970  wq * invvol * invdt * 0.5_rt * sxo_szn);
971 
972  amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
973  - wq * invvol * invdt * 0.5_rt * sxo_szo);
974  }
975  }
976 
977 #elif defined(WARPX_DIM_3D)
978 
979  for (int k=0; k<=depos_order; k++) {
980  for (int j=0; j<=depos_order; j++) {
981 
982  auto const syn_szn = static_cast<amrex::Real>(sy_new[j] * sz_new[k]);
983  auto const syo_szn = static_cast<amrex::Real>(sy_old[j] * sz_new[k]);
984  auto const syn_szo = static_cast<amrex::Real>(sy_new[j] * sz_old[k]);
985  auto const syo_szo = static_cast<amrex::Real>(sy_old[j] * sz_old[k]);
986 
987  for (int i=0; i<=depos_order; i++) {
988 
989  auto const sxn_syn_szn = static_cast<amrex::Real>(sx_new[i]) * syn_szn;
990  auto const sxo_syn_szn = static_cast<amrex::Real>(sx_old[i]) * syn_szn;
991  auto const sxn_syo_szn = static_cast<amrex::Real>(sx_new[i]) * syo_szn;
992  auto const sxo_syo_szn = static_cast<amrex::Real>(sx_old[i]) * syo_szn;
993  auto const sxn_syn_szo = static_cast<amrex::Real>(sx_new[i]) * syn_szo;
994  auto const sxo_syn_szo = static_cast<amrex::Real>(sx_old[i]) * syn_szo;
995  auto const sxn_syo_szo = static_cast<amrex::Real>(sx_new[i]) * syo_szo;
996  auto const sxo_syo_szo = static_cast<amrex::Real>(sx_old[i]) * syo_szo;
997 
998  // Jx
999  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k),
1000  wq * invvol * invdt * onethird * sxn_syn_szn);
1001 
1002  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k),
1003  - wq * invvol * invdt * onethird * sxo_syn_szn);
1004 
1005  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k),
1006  wq * invvol * invdt * onesixth * sxn_syo_szn);
1007 
1008  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_old + i, lo.y + j_old + j,lo.z + k_new + k),
1009  - wq * invvol * invdt * onesixth * sxo_syo_szn);
1010 
1011  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k),
1012  wq * invvol * invdt * onesixth * sxn_syn_szo);
1013 
1014  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k),
1015  - wq * invvol * invdt * onesixth * sxo_syn_szo);
1016 
1017  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k),
1018  wq * invvol * invdt * onethird * sxn_syo_szo);
1019 
1020  amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k),
1021  - wq * invvol * invdt * onethird * sxo_syo_szo);
1022 
1023  // Jy
1024  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k),
1025  wq * invvol * invdt * onethird * sxn_syn_szn);
1026 
1027  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k),
1028  - wq * invvol * invdt * onethird * sxn_syo_szn);
1029 
1030  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k),
1031  wq * invvol * invdt * onesixth * sxo_syn_szn);
1032 
1033  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k),
1034  - wq * invvol * invdt * onesixth * sxo_syo_szn);
1035 
1036  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k),
1037  wq * invvol * invdt * onesixth * sxn_syn_szo);
1038 
1039  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k),
1040  - wq * invvol * invdt * onesixth * sxn_syo_szo);
1041 
1042  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k),
1043  wq * invvol * invdt * onethird * sxo_syn_szo);
1044 
1045  amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k),
1046  - wq * invvol * invdt * onethird * sxo_syo_szo);
1047 
1048  // Jz
1049  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k),
1050  wq * invvol * invdt * onethird * sxn_syn_szn);
1051 
1052  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k),
1053  - wq * invvol * invdt * onethird * sxn_syn_szo);
1054 
1055  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k),
1056  wq * invvol * invdt * onesixth * sxo_syn_szn);
1057 
1058  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k),
1059  - wq * invvol * invdt * onesixth * sxo_syn_szo);
1060 
1061  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k),
1062  wq * invvol * invdt * onesixth * sxn_syo_szn);
1063 
1064  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k),
1065  - wq * invvol * invdt * onesixth * sxn_syo_szo);
1066 
1067  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k),
1068  wq * invvol * invdt * onethird * sxo_syo_szn);
1069 
1070  amrex::Gpu::Atomic::AddNoRet(&jz_arr( lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k),
1071  - wq * invvol * invdt * onethird * sxo_syo_szo);
1072  }
1073  }
1074  }
1075 #endif
1076  } );
1077 # if defined(WARPX_USE_GPUCLOCK)
1078  if( load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::GpuClock) {
1079  amrex::Gpu::streamSynchronize();
1080  *cost += *cost_real;
1081  amrex::The_Managed_Arena()->free(cost_real);
1082  }
1083 # endif
1084 #endif // #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z)
1085 }
1086 #endif // CURRENTDEPOSITION_H_
amrex::BaseFab::box
const Box & box() const noexcept
Complex
amrex::GpuComplex< amrex::Real > Complex
Definition: WarpX_Complex.H:22
amrex::Gpu::Atomic::AddNoRet
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddNoRet(T *const sum, T const value) noexcept
amrex::BaseFab::array
AMREX_FORCE_INLINE Array4< T const > array() const noexcept
GetAndSetPosition.H
WarpXConst.H
LoadBalanceCostsUpdateAlgo::GpuClock
Definition: WarpXAlgorithmSelection.H:107
ablastr::constant::SI::c
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
compute_domain.dx
int dx
Definition: compute_domain.py:35
Stencil.dt
int dt
Definition: Stencil.py:468
doDepositionShapeN
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 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, amrex::Real *cost, const long load_balance_costs_update_algo)
Current Deposition for thread thread_num.
Definition: CurrentDeposition.H:53
doEsirkepovDepositionShapeN
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 *const 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 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, amrex::Real *const cost, const long load_balance_costs_update_algo)
Esirkepov Current Deposition for thread thread_num.
Definition: CurrentDeposition.H:365
Compute_shifted_shape_factor
Definition: ShapeFactors.H:81
amrex::Arena::free
virtual void free(void *pt)=0
doVayDepositionShapeN
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 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, amrex::Real *cost, const long load_balance_costs_update_algo)
Vay current deposition (Vay et al, 2013) for thread thread_num: deposit D in real space and store the...
Definition: CurrentDeposition.H:744
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMReX_REAL.H
Compute_shape_factor
Definition: ShapeFactors.H:26
amrex::GpuComplex::imag
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T imag() const noexcept
amrex::ParallelFor
std::enable_if_t< std::is_integral< T >::value > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
amrex::Abort
void Abort(const std::string &msg)
WarpXAlgorithmSelection.H
check_interp_points_and_weights.i
i
Definition: check_interp_points_and_weights.py:173
KernelTimer.H
amrex::Arena::alloc
virtual void * alloc(std::size_t sz)=0
ShapeFactors.H
amrex::GpuComplex::real
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T real() const noexcept
WarpX_Complex.H
AMReX_Arena.H
GetParticlePosition
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel...
Definition: GetAndSetPosition.H:52
amrex::Gpu::streamSynchronize
void streamSynchronize() noexcept
KernelTimer
Defines a timer object to be used on GPU; measures summed thread cycles.
Definition: KernelTimer.H:23
amrex::Box::type
AMREX_GPU_HOST_DEVICE IntVect type() const noexcept
amrex::The_Managed_Arena
Arena * The_Managed_Arena()
AMReX_Array4.H
NODE
NODE
Generated by   doxygen 1.8.13