WarpX
Loading...
Searching...
No Matches
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 WARPX_CURRENTDEPOSITION_H_
9#define WARPX_CURRENTDEPOSITION_H_
10
11#include "Particles/Deposition/SharedDepositionUtils.H"
12#include "Particles/Pusher/GetAndSetPosition.H"
13#include "Particles/Pusher/UpdatePosition.H"
14#include "Particles/ShapeFactors.H"
15#include "Utils/TextMsg.H"
16#include "Utils/WarpXAlgorithmSelection.H"
17#include "Utils/WarpXConst.H"
18#ifdef WARPX_DIM_RZ
19# include "Utils/WarpX_Complex.H"
20#endif
21
22#include <AMReX.H>
23#include <AMReX_Arena.H>
24#include <AMReX_Array4.H>
25#include <AMReX_Dim3.H>
26#include <AMReX_REAL.H>
27
46template <int depos_order>
47AMREX_GPU_HOST_DEVICE AMREX_INLINE
48void doDepositionShapeNKernel([[maybe_unused]] const amrex::ParticleReal xp,
49 [[maybe_unused]] const amrex::ParticleReal yp,
50 [[maybe_unused]] const amrex::ParticleReal zp,
51 const amrex::ParticleReal wq,
52 const amrex::ParticleReal vx,
53 const amrex::ParticleReal vy,
54 const amrex::ParticleReal vz,
55 amrex::Array4<amrex::Real> const& jx_arr,
56 amrex::Array4<amrex::Real> const& jy_arr,
57 amrex::Array4<amrex::Real> const& jz_arr,
58 amrex::IntVect const& jx_type,
59 amrex::IntVect const& jy_type,
60 amrex::IntVect const& jz_type,
61 const amrex::Real relative_time,
62 const amrex::XDim3 & dinv,
63 const amrex::XDim3 & xyzmin,
64 const amrex::Real invvol,
65 const amrex::Dim3 lo,
66 [[maybe_unused]] const int n_rz_azimuthal_modes)
67{
68 using namespace amrex::literals;
69
70 constexpr int NODE = amrex::IndexType::NODE;
71 constexpr int CELL = amrex::IndexType::CELL;
72
73 // wqx, wqy wqz are particle current in each direction
74#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
75 // In RZ and RCYLINDER, wqx is actually wqr, and wqy is wqtheta
76 // Convert to cylindrical at the mid point
77 const amrex::Real xpmid = xp + relative_time*vx;
78 const amrex::Real ypmid = yp + relative_time*vy;
79 const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
80 const amrex::Real costheta = (rpmid > 0._rt ? xpmid/rpmid : 1._rt);
81 const amrex::Real sintheta = (rpmid > 0._rt ? ypmid/rpmid : 0._rt);
82 const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
83 const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
84 const amrex::Real wqz = wq*invvol*vz;
85 #if defined(WARPX_DIM_RZ)
86 const Complex xy0 = Complex{costheta, sintheta};
87 #endif
88#elif defined(WARPX_DIM_RSPHERE)
89 // Convert to cylindrical at the mid point
90 const amrex::Real xpmid = xp + relative_time*vx;
91 const amrex::Real ypmid = yp + relative_time*vy;
92 const amrex::Real zpmid = zp + relative_time*vz;
93 const amrex::Real rpxymid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
94 const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid + zpmid*zpmid);
95 const amrex::Real costheta = (rpxymid > 0._rt ? xpmid/rpxymid : 1._rt);
96 const amrex::Real sintheta = (rpxymid > 0._rt ? ypmid/rpxymid : 0._rt);
97 const amrex::Real cosphi = (rpmid > 0._rt ? rpxymid/rpmid : 1._rt);
98 const amrex::Real sinphi = (rpmid > 0._rt ? zpmid/rpmid : 0._rt);
99 // convert from Cartesian to spherical
100 const amrex::Real wqx = wq*invvol*(+vx*costheta*cosphi + vy*sintheta*cosphi + vz*sinphi);
101 const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
102 const amrex::Real wqz = wq*invvol*(-vx*costheta*sinphi - vy*sintheta*sinphi + vz*cosphi);
103#else
104 const amrex::Real wqx = wq*invvol*vx;
105 const amrex::Real wqy = wq*invvol*vy;
106 const amrex::Real wqz = wq*invvol*vz;
107#endif
108
109 // --- Compute shape factors
110 Compute_shape_factor< depos_order > const compute_shape_factor;
111#if !defined(WARPX_DIM_1D_Z)
112 // x direction
113 // Get particle position after 1/2 push back in position
114 // Keep these double to avoid bug in single precision
115
116#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
117 const double xmid = (rpmid - xyzmin.x)*dinv.x;
118#else
119 const double xmid = ((xp - xyzmin.x) + relative_time*vx)*dinv.x;
120#endif
121
122 // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
123 // sx_j[xyz] shape factor along x for the centering of each current
124 // There are only two possible centerings, node or cell centered, so at most only two shape factor
125 // arrays will be needed.
126 // Keep these double to avoid bug in single precision
127 double sx_node[depos_order + 1] = {0.};
128 double sx_cell[depos_order + 1] = {0.};
129 int j_node = 0;
130 int j_cell = 0;
131 if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) {
132 j_node = compute_shape_factor(sx_node, xmid);
133 }
134 if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
135 j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
136 }
137
138 amrex::Real sx_jx[depos_order + 1] = {0._rt};
139 amrex::Real sx_jy[depos_order + 1] = {0._rt};
140 amrex::Real sx_jz[depos_order + 1] = {0._rt};
141 for (int ix=0; ix<=depos_order; ix++)
142 {
143 sx_jx[ix] = ((jx_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
144 sx_jy[ix] = ((jy_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
145 sx_jz[ix] = ((jz_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
146 }
147
148 int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell);
149 int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell);
150 int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell);
151#endif
152
153#if defined(WARPX_DIM_3D)
154 // y direction
155 // Keep these double to avoid bug in single precision
156 const double ymid = ((yp - xyzmin.y) + relative_time*vy)*dinv.y;
157 double sy_node[depos_order + 1] = {0.};
158 double sy_cell[depos_order + 1] = {0.};
159 int k_node = 0;
160 int k_cell = 0;
161 if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) {
162 k_node = compute_shape_factor(sy_node, ymid);
163 }
164 if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
165 k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
166 }
167 amrex::Real sy_jx[depos_order + 1] = {0._rt};
168 amrex::Real sy_jy[depos_order + 1] = {0._rt};
169 amrex::Real sy_jz[depos_order + 1] = {0._rt};
170 for (int iy=0; iy<=depos_order; iy++)
171 {
172 sy_jx[iy] = ((jx_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
173 sy_jy[iy] = ((jy_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
174 sy_jz[iy] = ((jz_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
175 }
176 int const k_jx = ((jx_type[1] == NODE) ? k_node : k_cell);
177 int const k_jy = ((jy_type[1] == NODE) ? k_node : k_cell);
178 int const k_jz = ((jz_type[1] == NODE) ? k_node : k_cell);
179#endif
180
181#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
182 // z direction
183 // Keep these double to avoid bug in single precision
184 constexpr int zdir = WARPX_ZINDEX;
185 const double zmid = ((zp - xyzmin.z) + relative_time*vz)*dinv.z;
186 double sz_node[depos_order + 1] = {0.};
187 double sz_cell[depos_order + 1] = {0.};
188 int l_node = 0;
189 int l_cell = 0;
190 if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) {
191 l_node = compute_shape_factor(sz_node, zmid);
192 }
193 if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
194 l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
195 }
196 amrex::Real sz_jx[depos_order + 1] = {0._rt};
197 amrex::Real sz_jy[depos_order + 1] = {0._rt};
198 amrex::Real sz_jz[depos_order + 1] = {0._rt};
199 for (int iz=0; iz<=depos_order; iz++)
200 {
201 sz_jx[iz] = ((jx_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
202 sz_jy[iz] = ((jy_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
203 sz_jz[iz] = ((jz_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
204 }
205 int const l_jx = ((jx_type[zdir] == NODE) ? l_node : l_cell);
206 int const l_jy = ((jy_type[zdir] == NODE) ? l_node : l_cell);
207 int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell);
208#endif
209
210 // Deposit current into jx_arr, jy_arr and jz_arr
211#if defined(WARPX_DIM_1D_Z)
212 for (int iz=0; iz<=depos_order; iz++){
213 amrex::Gpu::Atomic::AddNoRet(
214 &jx_arr(lo.x+l_jx+iz, 0, 0, 0),
215 sz_jx[iz]*wqx);
216 amrex::Gpu::Atomic::AddNoRet(
217 &jy_arr(lo.x+l_jy+iz, 0, 0, 0),
218 sz_jy[iz]*wqy);
219 amrex::Gpu::Atomic::AddNoRet(
220 &jz_arr(lo.x+l_jz+iz, 0, 0, 0),
221 sz_jz[iz]*wqz);
222 }
223#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
224 for (int ix=0; ix<=depos_order; ix++){
225 amrex::Gpu::Atomic::AddNoRet( &jx_arr(lo.x+j_jx+ix, 0, 0, 0), sx_jx[ix]*wqx);
226 amrex::Gpu::Atomic::AddNoRet( &jy_arr(lo.x+j_jy+ix, 0, 0, 0), sx_jy[ix]*wqy);
227 amrex::Gpu::Atomic::AddNoRet( &jz_arr(lo.x+j_jz+ix, 0, 0, 0), sx_jz[ix]*wqz);
228 }
229#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
230 for (int iz=0; iz<=depos_order; iz++){
231 for (int ix=0; ix<=depos_order; ix++){
232 amrex::Gpu::Atomic::AddNoRet(
233 &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
234 sx_jx[ix]*sz_jx[iz]*wqx);
235 amrex::Gpu::Atomic::AddNoRet(
236 &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
237 sx_jy[ix]*sz_jy[iz]*wqy);
238 amrex::Gpu::Atomic::AddNoRet(
239 &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
240 sx_jz[ix]*sz_jz[iz]*wqz);
241#if defined(WARPX_DIM_RZ)
242 Complex xy = xy0; // Note that xy is equal to e^{i m theta}
243 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
244 // The factor 2 on the weighting comes from the normalization of the modes
245 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());
246 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());
247 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());
248 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());
249 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());
250 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());
251 xy = xy*xy0;
252 }
253#endif
254 }
255 }
256#elif defined(WARPX_DIM_3D)
257 for (int iz=0; iz<=depos_order; iz++){
258 for (int iy=0; iy<=depos_order; iy++){
259 for (int ix=0; ix<=depos_order; ix++){
260 amrex::Gpu::Atomic::AddNoRet(
261 &jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz),
262 sx_jx[ix]*sy_jx[iy]*sz_jx[iz]*wqx);
263 amrex::Gpu::Atomic::AddNoRet(
264 &jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz),
265 sx_jy[ix]*sy_jy[iy]*sz_jy[iz]*wqy);
266 amrex::Gpu::Atomic::AddNoRet(
267 &jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz),
268 sx_jz[ix]*sy_jz[iy]*sz_jz[iz]*wqz);
269 }
270 }
271 }
272#endif
273}
274
297template <int depos_order>
298void doDepositionShapeN (const GetParticlePosition<PIdx>& GetPosition,
299 const amrex::ParticleReal * const wp,
300 const amrex::ParticleReal * const uxp,
301 const amrex::ParticleReal * const uyp,
302 const amrex::ParticleReal * const uzp,
303 const int* ion_lev,
304 amrex::FArrayBox& jx_fab,
305 amrex::FArrayBox& jy_fab,
306 amrex::FArrayBox& jz_fab,
307 long np_to_deposit,
308 amrex::Real relative_time,
309 const amrex::XDim3 & dinv,
310 const amrex::XDim3 & xyzmin,
311 amrex::Dim3 lo,
312 amrex::Real q,
313 [[maybe_unused]]int n_rz_azimuthal_modes)
314{
315 using namespace amrex::literals;
316
317 // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
318 // (do_ionization=1)
319 const bool do_ionization = ion_lev;
320
321 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
322
323 const amrex::Real clightsq = 1.0_rt/PhysConst::c/PhysConst::c;
324
325 amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
326 amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
327 amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
328 amrex::IntVect const jx_type = jx_fab.box().type();
329 amrex::IntVect const jy_type = jy_fab.box().type();
330 amrex::IntVect const jz_type = jz_fab.box().type();
331
332 // Loop over particles and deposit into jx_fab, jy_fab and jz_fab
333 amrex::ParallelFor(
334 np_to_deposit,
335 [=] AMREX_GPU_DEVICE (long ip) {
336 amrex::ParticleReal xp, yp, zp;
337 GetPosition(ip, xp, yp, zp);
338
339 // --- Get particle quantities
340 const amrex::Real gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
341 + uyp[ip]*uyp[ip]*clightsq
342 + uzp[ip]*uzp[ip]*clightsq);
343 const amrex::Real vx = uxp[ip]*gaminv;
344 const amrex::Real vy = uyp[ip]*gaminv;
345 const amrex::Real vz = uzp[ip]*gaminv;
346
347 amrex::Real wq = q*wp[ip];
348 if (do_ionization){
349 wq *= ion_lev[ip];
350 }
351
352 doDepositionShapeNKernel<depos_order>(xp, yp, zp, wq, vx, vy, vz, jx_arr, jy_arr, jz_arr,
353 jx_type, jy_type, jz_type,
354 relative_time, dinv, xyzmin,
355 invvol, lo, n_rz_azimuthal_modes);
356
357 }
358 );
359}
360
382template <int depos_order>
383void doDepositionShapeNImplicit(const GetParticlePosition<PIdx>& GetPosition,
384 const amrex::ParticleReal * const wp,
385 const amrex::ParticleReal * const uxp_n,
386 const amrex::ParticleReal * const uyp_n,
387 const amrex::ParticleReal * const uzp_n,
388 const amrex::ParticleReal * const uxp_nph,
389 const amrex::ParticleReal * const uyp_nph,
390 const amrex::ParticleReal * const uzp_nph,
391 const int * const ion_lev,
392 amrex::FArrayBox& jx_fab,
393 amrex::FArrayBox& jy_fab,
394 amrex::FArrayBox& jz_fab,
395 const long np_to_deposit,
396 const amrex::XDim3 & dinv,
397 const amrex::XDim3 & xyzmin,
398 const amrex::Dim3 lo,
399 const amrex::Real q,
400 [[maybe_unused]]const int n_rz_azimuthal_modes)
401{
402 using namespace amrex::literals;
403
404 // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
405 // (do_ionization=1)
406 const bool do_ionization = ion_lev;
407
408 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
409
410 amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
411 amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
412 amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
413 amrex::IntVect const jx_type = jx_fab.box().type();
414 amrex::IntVect const jy_type = jy_fab.box().type();
415 amrex::IntVect const jz_type = jz_fab.box().type();
416
417 // Loop over particles and deposit into jx_fab, jy_fab and jz_fab
418 amrex::ParallelFor(
419 np_to_deposit,
420 [=] AMREX_GPU_DEVICE (long ip) {
421 amrex::ParticleReal xp, yp, zp;
422 GetPosition(ip, xp, yp, zp);
423
424 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
425 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
426 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
427
428 const amrex::Real vx = uxp_nph[ip]*gaminv;
429 const amrex::Real vy = uyp_nph[ip]*gaminv;
430 const amrex::Real vz = uzp_nph[ip]*gaminv;
431
432 amrex::Real wq = q*wp[ip];
433 if (do_ionization){
434 wq *= ion_lev[ip];
435 }
436
437 const amrex::Real relative_time = 0._rt;
438 doDepositionShapeNKernel<depos_order>(xp, yp, zp, wq, vx, vy, vz, jx_arr, jy_arr, jz_arr,
439 jx_type, jy_type, jz_type,
440 relative_time, dinv, xyzmin,
441 invvol, lo, n_rz_azimuthal_modes);
442
443 }
444 );
445}
446
475template <int depos_order>
476void doDepositionSharedShapeN (const GetParticlePosition<PIdx>& GetPosition,
477 const amrex::ParticleReal * const wp,
478 const amrex::ParticleReal * const uxp,
479 const amrex::ParticleReal * const uyp,
480 const amrex::ParticleReal * const uzp,
481 const int* ion_lev,
482 amrex::FArrayBox& jx_fab,
483 amrex::FArrayBox& jy_fab,
484 amrex::FArrayBox& jz_fab,
485 long np_to_deposit,
486 const amrex::Real relative_time,
487 const amrex::XDim3 & dinv,
488 const amrex::XDim3 & xyzmin,
489 amrex::Dim3 lo,
490 amrex::Real q,
491 int n_rz_azimuthal_modes,
492 const amrex::DenseBins<WarpXParticleContainer::ParticleTileType::ParticleTileDataType>& a_bins,
493 const amrex::Box& box,
494 const amrex::Geometry& geom,
495 const amrex::IntVect& a_tbox_max_size,
496 const int threads_per_block,
497 const amrex::IntVect bin_size)
498{
499 using namespace amrex::literals;
500
501#if (defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)) && \
502 !(defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE))
503 using namespace amrex;
504
505 auto permutation = a_bins.permutationPtr();
506
507 amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
508 amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
509 amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
510 amrex::IntVect const jx_type = jx_fab.box().type();
511 amrex::IntVect const jy_type = jy_fab.box().type();
512 amrex::IntVect const jz_type = jz_fab.box().type();
513
514 constexpr int zdir = WARPX_ZINDEX;
515 constexpr int NODE = amrex::IndexType::NODE;
516 constexpr int CELL = amrex::IndexType::CELL;
517
518 // Loop over particles and deposit into jx_fab, jy_fab and jz_fab
519 const auto dxiarr = geom.InvCellSizeArray();
520 const auto plo = geom.ProbLoArray();
521 const auto domain = geom.Domain();
522
523 amrex::Box sample_tbox(IntVect(AMREX_D_DECL(0,0,0)), a_tbox_max_size - 1);
524 sample_tbox.grow(depos_order);
525
526 amrex::Box sample_tbox_x = convert(sample_tbox, jx_type);
527 amrex::Box sample_tbox_y = convert(sample_tbox, jy_type);
528 amrex::Box sample_tbox_z = convert(sample_tbox, jz_type);
529
530 const auto npts = amrex::max(sample_tbox_x.numPts(), sample_tbox_y.numPts(), sample_tbox_z.numPts());
531
532 const int nblocks = a_bins.numBins();
533 const auto offsets_ptr = a_bins.offsetsPtr();
534
535 const std::size_t shared_mem_bytes = npts*sizeof(amrex::Real);
536 const std::size_t max_shared_mem_bytes = amrex::Gpu::Device::sharedMemPerBlock();
537 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(shared_mem_bytes <= max_shared_mem_bytes,
538 "Tile size too big for GPU shared memory current deposition");
539
540 amrex::ignore_unused(np_to_deposit);
541 // Launch one thread-block per bin
542 amrex::launch(
543 nblocks, threads_per_block, shared_mem_bytes, amrex::Gpu::gpuStream(),
544 [=] AMREX_GPU_DEVICE () noexcept {
545 const int bin_id = blockIdx.x;
546 const unsigned int bin_start = offsets_ptr[bin_id];
547 const unsigned int bin_stop = offsets_ptr[bin_id+1];
548
549 if (bin_start == bin_stop) { return; /*this bin has no particles*/ }
550
551 // These boxes define the index space for the shared memory buffers
552 amrex::Box buffer_box;
553 {
554 ParticleReal xp, yp, zp;
555 GetPosition(permutation[bin_start], xp, yp, zp);
556#if defined(WARPX_DIM_3D)
557 IntVect iv = IntVect(int( amrex::Math::floor((xp-plo[0]) * dxiarr[0]) ),
558 int( amrex::Math::floor((yp-plo[1]) * dxiarr[1]) ),
559 int( amrex::Math::floor((zp-plo[2]) * dxiarr[2]) ));
560#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
561 IntVect iv = IntVect(int( amrex::Math::floor((xp-plo[0]) * dxiarr[0]) ),
562 int( amrex::Math::floor((zp-plo[1]) * dxiarr[1]) ));
563#elif defined(WARPX_DIM_1D_Z)
564 IntVect iv = IntVect(int( amrex::Math::floor((zp-plo[0]) * dxiarr[0]) ));
565#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
566 IntVect iv = IntVect(int( amrex::Math::floor((xp-plo[0]) * dxiarr[0]) ));
567#endif
568 iv += domain.smallEnd();
569 getTileIndex(iv, box, true, bin_size, buffer_box);
570 }
571
572 buffer_box.grow(depos_order);
573 Box tbox_x = convert(buffer_box, jx_type);
574 Box tbox_y = convert(buffer_box, jy_type);
575 Box tbox_z = convert(buffer_box, jz_type);
576
577 Gpu::SharedMemory<amrex::Real> gsm;
578 amrex::Real* const shared = gsm.dataPtr();
579
580 amrex::Array4<amrex::Real> const jx_buff(shared,
581 amrex::begin(tbox_x), amrex::end(tbox_x), 1);
582 amrex::Array4<amrex::Real> const jy_buff(shared,
583 amrex::begin(tbox_y), amrex::end(tbox_y), 1);
584 amrex::Array4<amrex::Real> const jz_buff(shared,
585 amrex::begin(tbox_z), amrex::end(tbox_z), 1);
586
587 // Zero-initialize the temporary array in shared memory
588 volatile amrex::Real* vs = shared;
589 for (int i = threadIdx.x; i < npts; i += blockDim.x){
590 vs[i] = 0.0;
591 }
592 __syncthreads();
593 for (unsigned int ip_orig = bin_start+threadIdx.x; ip_orig<bin_stop; ip_orig += blockDim.x)
594 {
595 const unsigned int ip = permutation[ip_orig];
596 depositComponent<depos_order>(GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_buff, jx_type,
597 relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes,
598 ip, zdir, NODE, CELL, 0);
599 }
600
601 __syncthreads();
602 addLocalToGlobal(tbox_x, jx_arr, jx_buff);
603 for (int i = threadIdx.x; i < npts; i += blockDim.x){
604 vs[i] = 0.0;
605 }
606
607 __syncthreads();
608 for (unsigned int ip_orig = bin_start+threadIdx.x; ip_orig<bin_stop; ip_orig += blockDim.x)
609 {
610 const unsigned int ip = permutation[ip_orig];
611 depositComponent<depos_order>(GetPosition, wp, uxp, uyp, uzp, ion_lev, jy_buff, jy_type,
612 relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes,
613 ip, zdir, NODE, CELL, 1);
614 }
615
616 __syncthreads();
617 addLocalToGlobal(tbox_y, jy_arr, jy_buff);
618 for (int i = threadIdx.x; i < npts; i += blockDim.x){
619 vs[i] = 0.0;
620 }
621
622 __syncthreads();
623 for (unsigned int ip_orig = bin_start+threadIdx.x; ip_orig<bin_stop; ip_orig += blockDim.x)
624 {
625 const unsigned int ip = permutation[ip_orig];
626 depositComponent<depos_order>(GetPosition, wp, uxp, uyp, uzp, ion_lev, jz_buff, jz_type,
627 relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes,
628 ip, zdir, NODE, CELL, 2);
629 }
630
631 __syncthreads();
632 addLocalToGlobal(tbox_z, jz_arr, jz_buff);
633 });
634#else // not using hip/cuda
635 // Note, you should never reach this part of the code. This funcion cannot be called unless
636 // using HIP/CUDA, and those things are checked prior
637 //don't use any args
638 ignore_unused(
639 GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab,
640 np_to_deposit, relative_time, dinv, xyzmin, lo, q,
641 n_rz_azimuthal_modes, a_bins, box, geom, a_tbox_max_size,
642 threads_per_block, bin_size);
643 WARPX_ABORT_WITH_MESSAGE("Shared memory only implemented for HIP/CUDA");
644#endif
645}
646
674template <int depos_order>
675void doEsirkepovDepositionShapeN (const GetParticlePosition<PIdx>& GetPosition,
676 const amrex::ParticleReal * const wp,
677 const amrex::ParticleReal * const uxp,
678 const amrex::ParticleReal * const uyp,
679 const amrex::ParticleReal * const uzp,
680 const int* ion_lev,
681 const amrex::Array4<amrex::Real>& Jx_arr,
682 const amrex::Array4<amrex::Real>& Jy_arr,
683 const amrex::Array4<amrex::Real>& Jz_arr,
684 long np_to_deposit,
685 amrex::Real dt,
686 amrex::Real relative_time,
687 const amrex::XDim3 & dinv,
688 const amrex::XDim3 & xyzmin,
689 amrex::Dim3 lo,
690 amrex::Real q,
691 [[maybe_unused]] int n_rz_azimuthal_modes,
692 const amrex::Array4<const int>& reduced_particle_shape_mask,
693 bool enable_reduced_shape
694 )
695{
696 using namespace amrex;
697 using namespace amrex::literals;
698
699 // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
700 // (do_ionization=1)
701 bool const do_ionization = ion_lev;
702#if !defined(WARPX_DIM_3D)
703 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
704#endif
705
706 amrex::XDim3 const invdtd = amrex::XDim3{(1.0_rt/dt)*dinv.y*dinv.z,
707 (1.0_rt/dt)*dinv.x*dinv.z,
708 (1.0_rt/dt)*dinv.x*dinv.y};
709
710 Real constexpr clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c );
711
712#if (AMREX_SPACEDIM > 1)
713 Real constexpr one_third = 1.0_rt / 3.0_rt;
714 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
715#endif
716
717 // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
718
719 // (Compile 2 versions of the kernel: with and without reduced shape)
720 enum eb_flags : int { has_reduced_shape, no_reduced_shape };
721 const int reduce_shape_runtime_flag = (enable_reduced_shape && (depos_order>1))? has_reduced_shape : no_reduced_shape;
722
723 amrex::ParallelFor( TypeList<CompileTimeOptions<has_reduced_shape,no_reduced_shape>>{},
724 {reduce_shape_runtime_flag},
725 np_to_deposit, [=] AMREX_GPU_DEVICE (long ip, auto reduce_shape_control) {
726 // --- Get particle quantities
727 Real const gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
728 + uyp[ip]*uyp[ip]*clightsq
729 + uzp[ip]*uzp[ip]*clightsq);
730
731 Real wq = q*wp[ip];
732 if (do_ionization){
733 wq *= ion_lev[ip];
734 }
735
736 ParticleReal xp, yp, zp;
737 GetPosition(ip, xp, yp, zp);
738
739 // computes current and old position in grid units
740#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
741 Real const xp_new = xp + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv;
742 Real const yp_new = yp + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv;
743 Real const xp_mid = xp_new - 0.5_rt*dt*uxp[ip]*gaminv;
744 Real const yp_mid = yp_new - 0.5_rt*dt*uyp[ip]*gaminv;
745 Real const xp_old = xp_new - dt*uxp[ip]*gaminv;
746 Real const yp_old = yp_new - dt*uyp[ip]*gaminv;
747 Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
748 Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
749 Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
750 const amrex::Real costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt);
751 const amrex::Real sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt);
752 // Keep these double to avoid bug in single precision
753 double const x_new = (rp_new - xyzmin.x)*dinv.x;
754 double const x_old = (rp_old - xyzmin.x)*dinv.x;
755#if defined(WARPX_DIM_RZ)
756 const amrex::Real costheta_new = (rp_new > 0._rt ? xp_new/rp_new : 1._rt);
757 const amrex::Real sintheta_new = (rp_new > 0._rt ? yp_new/rp_new : 0._rt);
758 const amrex::Real costheta_old = (rp_old > 0._rt ? xp_old/rp_old : 1._rt);
759 const amrex::Real sintheta_old = (rp_old > 0._rt ? yp_old/rp_old : 0._rt);
760 const Complex xy_new0 = Complex{costheta_new, sintheta_new};
761 const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
762 const Complex xy_old0 = Complex{costheta_old, sintheta_old};
763#endif
764#elif defined(WARPX_DIM_RSPHERE)
765 Real const xp_new = xp + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv;
766 Real const yp_new = yp + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv;
767 Real const zp_new = zp + (relative_time + 0.5_rt*dt)*uzp[ip]*gaminv;
768 Real const xp_mid = xp_new - 0.5_rt*dt*uxp[ip]*gaminv;
769 Real const yp_mid = yp_new - 0.5_rt*dt*uyp[ip]*gaminv;
770 Real const zp_mid = zp_new - 0.5_rt*dt*uzp[ip]*gaminv;
771 Real const xp_old = xp_new - dt*uxp[ip]*gaminv;
772 Real const yp_old = yp_new - dt*uyp[ip]*gaminv;
773 Real const zp_old = zp_new - dt*uzp[ip]*gaminv;
774 Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
775 Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
776 Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
777 Real const rp_mid = (rp_new + rp_old)*0.5_rt;
778
779 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
780 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
781 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
782 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
783
784 // Keep these double to avoid bug in single precision
785 double const x_new = (rp_new - xyzmin.x)*dinv.x;
786 double const x_old = (rp_old - xyzmin.x)*dinv.x;
787#else
788#if !defined(WARPX_DIM_1D_Z)
789 // Keep these double to avoid bug in single precision
790 double const x_new = (xp - xyzmin.x + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv)*dinv.x;
791 double const x_old = x_new - dt*dinv.x*uxp[ip]*gaminv;
792#endif
793#endif
794#if defined(WARPX_DIM_3D)
795 // Keep these double to avoid bug in single precision
796 double const y_new = (yp - xyzmin.y + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv)*dinv.y;
797 double const y_old = y_new - dt*dinv.y*uyp[ip]*gaminv;
798#endif
799#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
800 // Keep these double to avoid bug in single precision
801 double const z_new = (zp - xyzmin.z + (relative_time + 0.5_rt*dt)*uzp[ip]*gaminv)*dinv.z;
802 double const z_old = z_new - dt*dinv.z*uzp[ip]*gaminv;
803#endif
804
805 // Check whether the particle is close to the EB at the old and new position
806 bool reduce_shape_old, reduce_shape_new;
807#ifdef AMREX_USE_CUDA
808 amrex::ignore_unused(reduced_particle_shape_mask, lo); // Needed to avoid compilation error with nvcc
809#endif
810 if constexpr (reduce_shape_control == has_reduced_shape) {
811#if defined(WARPX_DIM_3D)
812 reduce_shape_old = reduced_particle_shape_mask(
813 lo.x + int(amrex::Math::floor(x_old)),
814 lo.y + int(amrex::Math::floor(y_old)),
815 lo.z + int(amrex::Math::floor(z_old)));
816 reduce_shape_new = reduced_particle_shape_mask(
817 lo.x + int(amrex::Math::floor(x_new)),
818 lo.y + int(amrex::Math::floor(y_new)),
819 lo.z + int(amrex::Math::floor(z_new)));
820#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
821 reduce_shape_old = reduced_particle_shape_mask(
822 lo.x + int(amrex::Math::floor(x_old)),
823 lo.y + int(amrex::Math::floor(z_old)),
824 0);
825 reduce_shape_new = reduced_particle_shape_mask(
826 lo.x + int(amrex::Math::floor(x_new)),
827 lo.y + int(amrex::Math::floor(z_new)),
828 0);
829#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
830 reduce_shape_old = reduced_particle_shape_mask(
831 lo.x + int(amrex::Math::floor(x_old)),
832 0, 0);
833 reduce_shape_new = reduced_particle_shape_mask(
834 lo.x + int(amrex::Math::floor(x_new)),
835 0, 0);
836#elif defined(WARPX_DIM_1D_Z)
837 reduce_shape_old = reduced_particle_shape_mask(
838 lo.x + int(amrex::Math::floor(z_old)),
839 0, 0);
840 reduce_shape_new = reduced_particle_shape_mask(
841 lo.x + int(amrex::Math::floor(z_new)),
842 0, 0);
843#endif
844 } else {
845 reduce_shape_old = false;
846 reduce_shape_new = false;
847 }
848
849#if defined(WARPX_DIM_RZ)
850 Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
851#elif defined(WARPX_DIM_XZ)
852 Real const vy = uyp[ip]*gaminv;
853#elif defined(WARPX_DIM_1D_Z)
854 Real const vx = uxp[ip]*gaminv;
855 Real const vy = uyp[ip]*gaminv;
856#elif defined(WARPX_DIM_RCYLINDER)
857 Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
858 Real const vz = uzp[ip]*gaminv;
859#elif defined(WARPX_DIM_RSPHERE)
860 // convert from Cartesian to spherical
861 Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
862 Real const vz = (-uxp[ip]*costheta_mid*sinphi_mid - uyp[ip]*sintheta_mid*sinphi_mid + uzp[ip]*cosphi_mid)*gaminv;
863#endif
864
865 // --- Compute shape factors
866 // Compute shape factors for position as they are now and at old positions
867 // [ijk]_new: leftmost grid point that the particle touches
868 const Compute_shape_factor< depos_order > compute_shape_factor;
869 const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor;
870 // In cells marked by reduced_particle_shape_mask, we need order 1 deposition
871 const Compute_shifted_shape_factor< 1 > compute_shifted_shape_factor_order1;
872 amrex::ignore_unused(compute_shifted_shape_factor_order1); // unused for `no_reduced_shape`
873
874 // Shape factor arrays
875 // Note that there are extra values above and below
876 // to possibly hold the factor for the old particle
877 // which can be at a different grid location.
878 // Keep these double to avoid bug in single precision
879#if !defined(WARPX_DIM_1D_Z)
880 double sx_new[depos_order + 3] = {0.};
881 double sx_old[depos_order + 3] = {0.};
882 const int i_new = compute_shape_factor(sx_new+1, x_new );
883 const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
884 // If particle is close to the embedded boundary, recompute deposition with order 1 shape
885 if constexpr (reduce_shape_control == has_reduced_shape) {
886 if (reduce_shape_new) {
887 for (int i=0; i<depos_order+3; i++) {sx_new[i] = 0.;} // Erase previous deposition
888 compute_shifted_shape_factor_order1( sx_new+depos_order/2, x_new, i_new+depos_order/2 ); // Redeposit with order 1
889 }
890 if (reduce_shape_old) {
891 for (int i=0; i<depos_order+3; i++) {sx_old[i] = 0.;} // Erase previous deposition
892 compute_shifted_shape_factor_order1( sx_old+depos_order/2, x_old, i_new+depos_order/2 ); // Redeposit with order 1
893 }
894 // Note: depos_order/2 in the above code corresponds to the shift between the index of the lowest point
895 // to which the particle can deposit, with shape of order `depos_order` vs with shape of order 1
896 }
897#endif
898#if defined(WARPX_DIM_3D)
899 double sy_new[depos_order + 3] = {0.};
900 double sy_old[depos_order + 3] = {0.};
901 const int j_new = compute_shape_factor(sy_new+1, y_new);
902 const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
903 // If particle is close to the embedded boundary, recompute deposition with order 1 shape
904 if constexpr (reduce_shape_control == has_reduced_shape) {
905 if (reduce_shape_new) {
906 for (int j=0; j<depos_order+3; j++) {sy_new[j] = 0.;} // Erase previous deposition
907 compute_shifted_shape_factor_order1( sy_new+depos_order/2, y_new, j_new+depos_order/2 ); // Redeposit with order 1
908 }
909 if (reduce_shape_old) {
910 for (int j=0; j<depos_order+3; j++) {sy_old[j] = 0.;} // Erase previous deposition
911 compute_shifted_shape_factor_order1( sy_old+depos_order/2, y_old, j_new+depos_order/2 ); // Redeposit with order 1
912 }
913 // Note: depos_order/2 in the above code corresponds to the shift between the index of the lowest point
914 // to which the particle can deposit, with shape of order `depos_order` vs with shape of order 1
915 }
916#endif
917#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
918 double sz_new[depos_order + 3] = {0.};
919 double sz_old[depos_order + 3] = {0.};
920 const int k_new = compute_shape_factor(sz_new+1, z_new );
921 const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new );
922 // If particle is close to the embedded boundary, recompute deposition with order 1 shape
923 if constexpr (reduce_shape_control == has_reduced_shape) {
924 if (reduce_shape_new) {
925 for (int k=0; k<depos_order+3; k++) {sz_new[k] = 0.;} // Erase previous deposition
926 compute_shifted_shape_factor_order1( sz_new+depos_order/2, z_new, k_new+depos_order/2 ); // Redeposit with order 1
927 }
928 if (reduce_shape_old) {
929 for (int k=0; k<depos_order+3; k++) {sz_old[k] = 0.;} // Erase previous deposition
930 compute_shifted_shape_factor_order1( sz_old+depos_order/2, z_old, k_new+depos_order/2 ); // Redeposit with order 1
931 }
932 // Note: depos_order/2 in the above code corresponds to the shift between the index of the lowest point
933 // to which the particle can deposit, with shape of order `depos_order` vs with shape of order 1
934 }
935#endif
936
937 // computes min/max positions of current contributions
938#if !defined(WARPX_DIM_1D_Z)
939 int dil = 1, diu = 1;
940 if (i_old < i_new) { dil = 0; }
941 if (i_old > i_new) { diu = 0; }
942#endif
943#if defined(WARPX_DIM_3D)
944 int djl = 1, dju = 1;
945 if (j_old < j_new) { djl = 0; }
946 if (j_old > j_new) { dju = 0; }
947#endif
948#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
949 int dkl = 1, dku = 1;
950 if (k_old < k_new) { dkl = 0; }
951 if (k_old > k_new) { dku = 0; }
952#endif
953
954#if defined(WARPX_DIM_3D)
955
956 for (int k=dkl; k<=depos_order+2-dku; k++) {
957 for (int j=djl; j<=depos_order+2-dju; j++) {
958 amrex::Real sdxi = 0._rt;
959 for (int i=dil; i<=depos_order+1-diu; i++) {
960 sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*(
961 one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k])
962 +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k]));
963 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);
964 }
965 }
966 }
967 for (int k=dkl; k<=depos_order+2-dku; k++) {
968 for (int i=dil; i<=depos_order+2-diu; i++) {
969 amrex::Real sdyj = 0._rt;
970 for (int j=djl; j<=depos_order+1-dju; j++) {
971 sdyj += wq*invdtd.y*(sy_old[j] - sy_new[j])*(
972 one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k])
973 +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k]));
974 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);
975 }
976 }
977 }
978 for (int j=djl; j<=depos_order+2-dju; j++) {
979 for (int i=dil; i<=depos_order+2-diu; i++) {
980 amrex::Real sdzk = 0._rt;
981 for (int k=dkl; k<=depos_order+1-dku; k++) {
982 sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*(
983 one_third*(sx_new[i]*sy_new[j] + sx_old[i]*sy_old[j])
984 +one_sixth*(sx_new[i]*sy_old[j] + sx_old[i]*sy_new[j]));
985 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);
986 }
987 }
988 }
989
990#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
991
992 for (int k=dkl; k<=depos_order+2-dku; k++) {
993 amrex::Real sdxi = 0._rt;
994 for (int i=dil; i<=depos_order+1-diu; i++) {
995 sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]);
996 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi);
997#if defined(WARPX_DIM_RZ)
998 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
999 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1000 // The factor 2 comes from the normalization of the modes
1001 const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
1002 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());
1003 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
1004 xy_mid = xy_mid*xy_mid0;
1005 }
1006#endif
1007 }
1008 }
1009 for (int k=dkl; k<=depos_order+2-dku; k++) {
1010 for (int i=dil; i<=depos_order+2-diu; i++) {
1011 Real const sdyj = wq*vy*invvol*(
1012 one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k])
1013 +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k]));
1014 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
1015#if defined(WARPX_DIM_RZ)
1016 Complex const I = Complex{0._rt, 1._rt};
1017 Complex xy_new = xy_new0;
1018 Complex xy_mid = xy_mid0;
1019 Complex xy_old = xy_old0;
1020 // Throughout the following loop, xy_ takes the value e^{i m theta_}
1021 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1022 // The factor 2 comes from the normalization of the modes
1023 // The minus sign comes from the different convention with respect to Davidson et al.
1024 const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xyzmin.x*dinv.x)*wq*invdtd.x/(amrex::Real)imode
1025 *(Complex(sx_new[i]*sz_new[k], 0._rt)*(xy_new - xy_mid)
1026 + Complex(sx_old[i]*sz_old[k], 0._rt)*(xy_mid - xy_old));
1027 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());
1028 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
1029 xy_new = xy_new*xy_new0;
1030 xy_mid = xy_mid*xy_mid0;
1031 xy_old = xy_old*xy_old0;
1032 }
1033#endif
1034 }
1035 }
1036 for (int i=dil; i<=depos_order+2-diu; i++) {
1037 Real sdzk = 0._rt;
1038 for (int k=dkl; k<=depos_order+1-dku; k++) {
1039 sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]);
1040 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
1041#if defined(WARPX_DIM_RZ)
1042 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1043 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1044 // The factor 2 comes from the normalization of the modes
1045 const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
1046 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());
1047 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
1048 xy_mid = xy_mid*xy_mid0;
1049 }
1050#endif
1051 }
1052 }
1053#elif defined(WARPX_DIM_1D_Z)
1054
1055 for (int k=dkl; k<=depos_order+2-dku; k++) {
1056 amrex::Real const sdxi = wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1057 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+k_new-1+k, 0, 0, 0), sdxi);
1058 }
1059 for (int k=dkl; k<=depos_order+2-dku; k++) {
1060 amrex::Real const sdyj = wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1061 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+k_new-1+k, 0, 0, 0), sdyj);
1062 }
1063 amrex::Real sdzk = 0._rt;
1064 for (int k=dkl; k<=depos_order+1-dku; k++) {
1065 sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k]);
1066 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk);
1067 }
1068
1069#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1070
1071 amrex::Real sdri = 0._rt;
1072 for (int i=dil; i<=depos_order+1-diu; i++) {
1073 sdri += wq*invdtd.x*(sx_old[i] - sx_new[i]);
1074 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, 0, 0, 0), sdri);
1075 }
1076 for (int i=dil; i<=depos_order+2-diu; i++) {
1077 amrex::Real const sdyj = wq*vy*invvol*0.5_rt*(sx_old[i] + sx_new[i]);
1078 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, 0, 0, 0), sdyj);
1079 }
1080 for (int i=dil; i<=depos_order+2-diu; i++) {
1081 amrex::Real const sdzi = wq*vz*invvol*0.5_rt*(sx_old[i] + sx_new[i]);
1082 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, 0, 0, 0), sdzi);
1083 }
1084#endif
1085 }
1086 );
1087}
1088
1113template <int depos_order>
1114void doChargeConservingDepositionShapeNImplicit ([[maybe_unused]]const amrex::ParticleReal * const xp_n,
1115 [[maybe_unused]]const amrex::ParticleReal * const yp_n,
1116 [[maybe_unused]]const amrex::ParticleReal * const zp_n,
1117 const GetParticlePosition<PIdx>& GetPosition,
1118 const amrex::ParticleReal * const wp,
1119 [[maybe_unused]]const amrex::ParticleReal * const uxp_n,
1120 [[maybe_unused]]const amrex::ParticleReal * const uyp_n,
1121 [[maybe_unused]]const amrex::ParticleReal * const uzp_n,
1122 [[maybe_unused]]const amrex::ParticleReal * const uxp_nph,
1123 [[maybe_unused]]const amrex::ParticleReal * const uyp_nph,
1124 [[maybe_unused]]const amrex::ParticleReal * const uzp_nph,
1125 const int * const ion_lev,
1126 const amrex::Array4<amrex::Real>& Jx_arr,
1127 const amrex::Array4<amrex::Real>& Jy_arr,
1128 const amrex::Array4<amrex::Real>& Jz_arr,
1129 const long np_to_deposit,
1130 const amrex::Real dt,
1131 const amrex::XDim3 & dinv,
1132 const amrex::XDim3 & xyzmin,
1133 const amrex::Dim3 lo,
1134 const amrex::Real q,
1135 [[maybe_unused]] const int n_rz_azimuthal_modes)
1136{
1137 using namespace amrex;
1138
1139 // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
1140 // (do_ionization=1)
1141 bool const do_ionization = ion_lev;
1142
1143#if !defined(WARPX_DIM_3D)
1144 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
1145#endif
1146
1147 amrex::XDim3 const invdtd = amrex::XDim3{(1.0_rt/dt)*dinv.y*dinv.z,
1148 (1.0_rt/dt)*dinv.x*dinv.z,
1149 (1.0_rt/dt)*dinv.x*dinv.y};
1150
1151#if (AMREX_SPACEDIM > 1)
1152 Real constexpr one_third = 1.0_rt / 3.0_rt;
1153 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
1154#endif
1155
1156 // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
1157 amrex::ParallelFor(
1158 np_to_deposit,
1159 [=] AMREX_GPU_DEVICE (long const ip) {
1160
1161#if !defined(WARPX_DIM_3D)
1162 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
1163 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
1164 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
1165#endif
1166
1167 Real wq = q*wp[ip];
1168 if (do_ionization){
1169 wq *= ion_lev[ip];
1170 }
1171
1172 ParticleReal xp_nph, yp_nph, zp_nph;
1173 GetPosition(ip, xp_nph, yp_nph, zp_nph);
1174
1175#if !defined(WARPX_DIM_1D_Z)
1176 ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n[ip];
1177#endif
1178#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1179 ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n[ip];
1180#endif
1181#if !defined(WARPX_DIM_RCYLINDER)
1182 ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n[ip];
1183#endif
1184
1185 // computes current and old position in grid units
1186#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
1187 amrex::Real const xp_new = xp_np1;
1188 amrex::Real const yp_new = yp_np1;
1189 amrex::Real const xp_mid = xp_nph;
1190 amrex::Real const yp_mid = yp_nph;
1191 amrex::Real const xp_old = xp_n[ip];
1192 amrex::Real const yp_old = yp_n[ip];
1193 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
1194 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
1195 amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
1196 const amrex::Real costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt);
1197 const amrex::Real sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt);
1198 // Keep these double to avoid bug in single precision
1199 double const x_new = (rp_new - xyzmin.x)*dinv.x;
1200 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1201#if defined(WARPX_DIM_RZ)
1202 const amrex::Real costheta_new = (rp_new > 0._rt ? xp_new/rp_new : 1._rt);
1203 const amrex::Real sintheta_new = (rp_new > 0._rt ? yp_new/rp_new : 0._rt);
1204 const amrex::Real costheta_old = (rp_old > 0._rt ? xp_old/rp_old : 1._rt);
1205 const amrex::Real sintheta_old = (rp_old > 0._rt ? yp_old/rp_old : 0._rt);
1206 const Complex xy_new0 = Complex{costheta_new, sintheta_new};
1207 const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
1208 const Complex xy_old0 = Complex{costheta_old, sintheta_old};
1209#endif
1210#elif defined(WARPX_DIM_RSPHERE)
1211 amrex::Real const xp_new = xp_np1;
1212 amrex::Real const yp_new = yp_np1;
1213 amrex::Real const zp_new = zp_np1;
1214 amrex::Real const xp_mid = xp_nph;
1215 amrex::Real const yp_mid = yp_nph;
1216 amrex::Real const zp_mid = zp_nph;
1217 amrex::Real const xp_old = xp_n[ip];
1218 amrex::Real const yp_old = yp_n[ip];
1219 amrex::Real const zp_old = zp_n[ip];
1220 amrex::Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
1221 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
1222 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
1223 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1224
1225 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
1226 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
1227 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
1228 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
1229
1230 // Keep these double to avoid bug in single precision
1231 double const x_new = (rp_new - xyzmin.x)*dinv.x;
1232 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1233#else
1234#if !defined(WARPX_DIM_1D_Z)
1235 // Keep these double to avoid bug in single precision
1236 double const x_new = (xp_np1 - xyzmin.x)*dinv.x;
1237 double const x_old = (xp_n[ip] - xyzmin.x)*dinv.x;
1238#endif
1239#endif
1240#if defined(WARPX_DIM_3D)
1241 // Keep these double to avoid bug in single precision
1242 double const y_new = (yp_np1 - xyzmin.y)*dinv.y;
1243 double const y_old = (yp_n[ip] - xyzmin.y)*dinv.y;
1244#endif
1245#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1246 // Keep these double to avoid bug in single precision
1247 double const z_new = (zp_np1 - xyzmin.z)*dinv.z;
1248 double const z_old = (zp_n[ip] - xyzmin.z)*dinv.z;
1249#endif
1250
1251#if defined(WARPX_DIM_RZ)
1252 amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
1253#elif defined(WARPX_DIM_XZ)
1254 amrex::Real const vy = uyp_nph[ip]*gaminv;
1255#elif defined(WARPX_DIM_1D_Z)
1256 amrex::Real const vx = uxp_nph[ip]*gaminv;
1257 amrex::Real const vy = uyp_nph[ip]*gaminv;
1258#elif defined(WARPX_DIM_RCYLINDER)
1259 amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
1260 amrex::Real const vz = uzp_nph[ip]*gaminv;
1261#elif defined(WARPX_DIM_RSPHERE)
1262 // convert from Cartesian to spherical
1263 amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
1264 amrex::Real const vz = (-uxp_nph[ip]*costheta_mid*sinphi_mid - uyp_nph[ip]*sintheta_mid*sinphi_mid + uzp_nph[ip]*cosphi_mid)*gaminv;
1265#endif
1266
1267 // --- Compute shape factors
1268 // Compute shape factors for position as they are now and at old positions
1269 // [ijk]_new: leftmost grid point that the particle touches
1270 const Compute_shape_factor< depos_order > compute_shape_factor;
1271 const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor;
1272
1273 // Shape factor arrays
1274 // Note that there are extra values above and below
1275 // to possibly hold the factor for the old particle
1276 // which can be at a different grid location.
1277 // Keep these double to avoid bug in single precision
1278#if !defined(WARPX_DIM_1D_Z)
1279 double sx_new[depos_order + 3] = {0.};
1280 double sx_old[depos_order + 3] = {0.};
1281 const int i_new = compute_shape_factor(sx_new+1, x_new);
1282 const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
1283#endif
1284#if defined(WARPX_DIM_3D)
1285 double sy_new[depos_order + 3] = {0.};
1286 double sy_old[depos_order + 3] = {0.};
1287 const int j_new = compute_shape_factor(sy_new+1, y_new);
1288 const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
1289#endif
1290#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1291 double sz_new[depos_order + 3] = {0.};
1292 double sz_old[depos_order + 3] = {0.};
1293 const int k_new = compute_shape_factor(sz_new+1, z_new);
1294 const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new);
1295#endif
1296
1297 // computes min/max positions of current contributions
1298#if !defined(WARPX_DIM_1D_Z)
1299 int dil = 1, diu = 1;
1300 if (i_old < i_new) { dil = 0; }
1301 if (i_old > i_new) { diu = 0; }
1302#endif
1303#if defined(WARPX_DIM_3D)
1304 int djl = 1, dju = 1;
1305 if (j_old < j_new) { djl = 0; }
1306 if (j_old > j_new) { dju = 0; }
1307#endif
1308#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1309 int dkl = 1, dku = 1;
1310 if (k_old < k_new) { dkl = 0; }
1311 if (k_old > k_new) { dku = 0; }
1312#endif
1313
1314#if defined(WARPX_DIM_3D)
1315
1316 for (int k=dkl; k<=depos_order+2-dku; k++) {
1317 for (int j=djl; j<=depos_order+2-dju; j++) {
1318 amrex::Real sdxi = 0._rt;
1319 for (int i=dil; i<=depos_order+1-diu; i++) {
1320 sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*(
1321 one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k])
1322 +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k]));
1323 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);
1324 }
1325 }
1326 }
1327 for (int k=dkl; k<=depos_order+2-dku; k++) {
1328 for (int i=dil; i<=depos_order+2-diu; i++) {
1329 amrex::Real sdyj = 0._rt;
1330 for (int j=djl; j<=depos_order+1-dju; j++) {
1331 sdyj += wq*invdtd.y*(sy_old[j] - sy_new[j])*(
1332 one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k])
1333 +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k]));
1334 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);
1335 }
1336 }
1337 }
1338 for (int j=djl; j<=depos_order+2-dju; j++) {
1339 for (int i=dil; i<=depos_order+2-diu; i++) {
1340 amrex::Real sdzk = 0._rt;
1341 for (int k=dkl; k<=depos_order+1-dku; k++) {
1342 sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*(
1343 one_third*(sx_new[i]*sy_new[j] + sx_old[i]*sy_old[j])
1344 +one_sixth*(sx_new[i]*sy_old[j] + sx_old[i]*sy_new[j]));
1345 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);
1346 }
1347 }
1348 }
1349
1350#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1351
1352 for (int k=dkl; k<=depos_order+2-dku; k++) {
1353 amrex::Real sdxi = 0._rt;
1354 for (int i=dil; i<=depos_order+1-diu; i++) {
1355 sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]);
1356 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi);
1357#if defined(WARPX_DIM_RZ)
1358 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1359 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1360 // The factor 2 comes from the normalization of the modes
1361 const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
1362 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());
1363 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
1364 xy_mid = xy_mid*xy_mid0;
1365 }
1366#endif
1367 }
1368 }
1369 for (int k=dkl; k<=depos_order+2-dku; k++) {
1370 for (int i=dil; i<=depos_order+2-diu; i++) {
1371 Real const sdyj = wq*vy*invvol*(
1372 one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k])
1373 +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k]));
1374 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
1375#if defined(WARPX_DIM_RZ)
1376 Complex const I = Complex{0._rt, 1._rt};
1377 Complex xy_new = xy_new0;
1378 Complex xy_mid = xy_mid0;
1379 Complex xy_old = xy_old0;
1380 // Throughout the following loop, xy_ takes the value e^{i m theta_}
1381 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1382 // The factor 2 comes from the normalization of the modes
1383 // The minus sign comes from the different convention with respect to Davidson et al.
1384 const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xyzmin.x*dinv.x)*wq*invdtd.x/(amrex::Real)imode
1385 *(Complex(sx_new[i]*sz_new[k], 0._rt)*(xy_new - xy_mid)
1386 + Complex(sx_old[i]*sz_old[k], 0._rt)*(xy_mid - xy_old));
1387 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());
1388 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
1389 xy_new = xy_new*xy_new0;
1390 xy_mid = xy_mid*xy_mid0;
1391 xy_old = xy_old*xy_old0;
1392 }
1393#endif
1394 }
1395 }
1396 for (int i=dil; i<=depos_order+2-diu; i++) {
1397 Real sdzk = 0._rt;
1398 for (int k=dkl; k<=depos_order+1-dku; k++) {
1399 sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]);
1400 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
1401#if defined(WARPX_DIM_RZ)
1402 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1403 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1404 // The factor 2 comes from the normalization of the modes
1405 const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
1406 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());
1407 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
1408 xy_mid = xy_mid*xy_mid0;
1409 }
1410#endif
1411 }
1412 }
1413#elif defined(WARPX_DIM_1D_Z)
1414
1415 for (int k=dkl; k<=depos_order+2-dku; k++) {
1416 amrex::Real const sdxi = wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1417 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+k_new-1+k, 0, 0, 0), sdxi);
1418 }
1419 for (int k=dkl; k<=depos_order+2-dku; k++) {
1420 amrex::Real const sdyj = wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1421 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+k_new-1+k, 0, 0, 0), sdyj);
1422 }
1423 amrex::Real sdzk = 0._rt;
1424 for (int k=dkl; k<=depos_order+1-dku; k++) {
1425 sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k]);
1426 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk);
1427 }
1428
1429#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1430
1431 amrex::Real sdri = 0._rt;
1432 for (int i=dil; i<=depos_order+1-diu; i++) {
1433 sdri += wq*invdtd.x*(sx_old[i] - sx_new[i]);
1434 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, 0, 0, 0), sdri);
1435 }
1436 for (int i=dil; i<=depos_order+2-diu; i++) {
1437 amrex::Real const sdyj = wq*vy*invvol*0.5_rt*(sx_old[i] + sx_new[i]);
1438 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, 0, 0, 0), sdyj);
1439 }
1440 for (int i=dil; i<=depos_order+2-diu; i++) {
1441 amrex::Real const sdzk = wq*vz*invvol*0.5_rt*(sx_old[i] + sx_new[i]);
1442 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, 0, 0, 0), sdzk);
1443 }
1444#endif
1445 }
1446 );
1447}
1448
1469template <int depos_order>
1470AMREX_GPU_HOST_DEVICE AMREX_INLINE
1471void VillasenorDepositionShapeNKernel ([[maybe_unused]]amrex::ParticleReal const xp_old,
1472 [[maybe_unused]]amrex::ParticleReal const yp_old,
1473 [[maybe_unused]]amrex::ParticleReal const zp_old,
1474 [[maybe_unused]]amrex::ParticleReal const xp_new,
1475 [[maybe_unused]]amrex::ParticleReal const yp_new,
1476 [[maybe_unused]]amrex::ParticleReal const zp_new,
1477 amrex::ParticleReal const wq,
1478 [[maybe_unused]]amrex::ParticleReal const uxp_mid,
1479 [[maybe_unused]]amrex::ParticleReal const uyp_mid,
1480 [[maybe_unused]]amrex::ParticleReal const uzp_mid,
1481 [[maybe_unused]]amrex::ParticleReal const gaminv,
1482 amrex::Array4<amrex::Real>const & Jx_arr,
1483 amrex::Array4<amrex::Real>const & Jy_arr,
1484 amrex::Array4<amrex::Real>const & Jz_arr,
1485 amrex::Real const dt,
1486 amrex::XDim3 const & dinv,
1487 amrex::XDim3 const & xyzmin,
1488 amrex::Dim3 const lo,
1489 amrex::Real const invvol,
1490 [[maybe_unused]] int const n_rz_azimuthal_modes)
1491{
1492
1493 using namespace amrex::literals;
1494
1495#if (AMREX_SPACEDIM > 1)
1496 amrex::Real constexpr one_third = 1.0_rt / 3.0_rt;
1497 amrex::Real constexpr one_sixth = 1.0_rt / 6.0_rt;
1498#endif
1499
1500 // computes current and old position in grid units
1501#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
1502 amrex::Real const xp_mid = (xp_new + xp_old)*0.5_rt;
1503 amrex::Real const yp_mid = (yp_new + yp_old)*0.5_rt;
1504 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
1505 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
1506 amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
1507 amrex::Real const costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt);
1508 amrex::Real const sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt);
1509#if defined(WARPX_DIM_RZ)
1510 Complex const xy_mid0 = Complex{costheta_mid, sintheta_mid};
1511#endif
1512
1513 // Keep these double to avoid bug in single precision
1514 double const x_new = (rp_new - xyzmin.x)*dinv.x;
1515 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1516 amrex::Real const vx = (rp_new - rp_old)/dt;
1517 amrex::Real const vy = (-uxp_mid*sintheta_mid + uyp_mid*costheta_mid)*gaminv;
1518#if defined(WARPX_DIM_RCYLINDER)
1519 amrex::Real const vz = uzp_mid*gaminv;
1520#endif
1521#elif defined(WARPX_DIM_RSPHERE)
1522 amrex::Real const xp_mid = (xp_new + xp_old)*0.5_rt;
1523 amrex::Real const yp_mid = (yp_new + yp_old)*0.5_rt;
1524 amrex::Real const zp_mid = (zp_new + zp_old)*0.5_rt;
1525 amrex::Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
1526 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
1527 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
1528 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1529
1530 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
1531 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
1532 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
1533 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
1534
1535 // Keep these double to avoid bug in single precision
1536 double const x_new = (rp_new - xyzmin.x)*dinv.x;
1537 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1538 amrex::Real const vx = (rp_new - rp_old)/dt;
1539 // convert from Cartesian to spherical
1540 amrex::Real const vy = (-uxp_mid*sintheta_mid + uyp_mid*costheta_mid)*gaminv;
1541 amrex::Real const vz = (-uxp_mid*costheta_mid*sinphi_mid - uyp_mid*sintheta_mid*sinphi_mid + uzp_mid*cosphi_mid)*gaminv;
1542#elif defined(WARPX_DIM_XZ)
1543 // Keep these double to avoid bug in single precision
1544 double const x_new = (xp_new - xyzmin.x)*dinv.x;
1545 double const x_old = (xp_old - xyzmin.x)*dinv.x;
1546 amrex::Real const vx = (xp_new - xp_old)/dt;
1547 amrex::Real const vy = uyp_mid*gaminv;
1548#elif defined(WARPX_DIM_1D_Z)
1549 amrex::Real const vx = uxp_mid*gaminv;
1550 amrex::Real const vy = uyp_mid*gaminv;
1551#elif defined(WARPX_DIM_3D)
1552 // Keep these double to avoid bug in single precision
1553 double const x_new = (xp_new - xyzmin.x)*dinv.x;
1554 double const x_old = (xp_old - xyzmin.x)*dinv.x;
1555 double const y_new = (yp_new - xyzmin.y)*dinv.y;
1556 double const y_old = (yp_old - xyzmin.y)*dinv.y;
1557 amrex::Real const vx = (xp_new - xp_old)/dt;
1558 amrex::Real const vy = (yp_new - yp_old)/dt;
1559#endif
1560
1561#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1562 // Keep these double to avoid bug in single precision
1563 double const z_new = (zp_new - xyzmin.z)*dinv.z;
1564 double const z_old = (zp_old - xyzmin.z)*dinv.z;
1565 amrex::Real const vz = (zp_new - zp_old)/dt;
1566#endif
1567
1568 // Define velocity kernals to deposit
1569 amrex::Real const wqx = wq*vx*invvol;
1570 amrex::Real const wqy = wq*vy*invvol;
1571 amrex::Real const wqz = wq*vz*invvol;
1572
1573 // 1) Determine the number of segments.
1574 // 2) Loop over segments and deposit current.
1575
1576 // cell crossings are defined at cell edges if depos_order is odd
1577 // cell crossings are defined at cell centers if depos_order is even
1578
1579 int num_segments = 1;
1580 double shift = 0.0;
1581 if ( (depos_order % 2) == 0 ) { shift = 0.5; }
1582
1583#if defined(WARPX_DIM_3D)
1584
1585 // compute cell crossings in X-direction
1586 const auto i_old = static_cast<int>(x_old-shift);
1587 const auto i_new = static_cast<int>(x_new-shift);
1588 const int cell_crossings_x = std::abs(i_new-i_old);
1589 num_segments += cell_crossings_x;
1590
1591 // compute cell crossings in Y-direction
1592 const auto j_old = static_cast<int>(y_old-shift);
1593 const auto j_new = static_cast<int>(y_new-shift);
1594 const int cell_crossings_y = std::abs(j_new-j_old);
1595 num_segments += cell_crossings_y;
1596
1597 // compute cell crossings in Z-direction
1598 const auto k_old = static_cast<int>(z_old-shift);
1599 const auto k_new = static_cast<int>(z_new-shift);
1600 const int cell_crossings_z = std::abs(k_new-k_old);
1601 num_segments += cell_crossings_z;
1602
1603 // need to assert that the number of cell crossings in each direction
1604 // is within the range permitted by the number of guard cells
1605 // e.g., if (num_segments > 7) ...
1606
1607 // compute total change in particle position and the initial cell
1608 // locations in each direction used to find the position at cell crossings.
1609 const double dxp = x_new - x_old;
1610 const double dyp = y_new - y_old;
1611 const double dzp = z_new - z_old;
1612 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1613 const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
1614 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1615 double Xcell = 0., Ycell = 0., Zcell = 0.;
1616 if (num_segments > 1) {
1617 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1618 Ycell = static_cast<double>(j_old) + shift + 0.5*(1.-dirY_sign);
1619 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1620 }
1621
1622 // loop over the number of segments and deposit
1623 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1624 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1625 double dxp_seg, dyp_seg, dzp_seg;
1626 double x0_new, y0_new, z0_new;
1627 double x0_old = x_old;
1628 double y0_old = y_old;
1629 double z0_old = z_old;
1630
1631 for (int ns=0; ns<num_segments; ns++) {
1632
1633 if (ns == num_segments-1) { // final segment
1634
1635 x0_new = x_new;
1636 y0_new = y_new;
1637 z0_new = z_new;
1638 dxp_seg = x0_new - x0_old;
1639 dyp_seg = y0_new - y0_old;
1640 dzp_seg = z0_new - z0_old;
1641
1642 }
1643 else {
1644
1645 x0_new = Xcell + dirX_sign;
1646 y0_new = Ycell + dirY_sign;
1647 z0_new = Zcell + dirZ_sign;
1648 dxp_seg = x0_new - x0_old;
1649 dyp_seg = y0_new - y0_old;
1650 dzp_seg = z0_new - z0_old;
1651
1652 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
1653 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
1654 Xcell = x0_new;
1655 dyp_seg = dyp/dxp*dxp_seg;
1656 dzp_seg = dzp/dxp*dxp_seg;
1657 y0_new = y0_old + dyp_seg;
1658 z0_new = z0_old + dzp_seg;
1659 }
1660 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
1661 Ycell = y0_new;
1662 dxp_seg = dxp/dyp*dyp_seg;
1663 dzp_seg = dzp/dyp*dyp_seg;
1664 x0_new = x0_old + dxp_seg;
1665 z0_new = z0_old + dzp_seg;
1666 }
1667 else {
1668 Zcell = z0_new;
1669 dxp_seg = dxp/dzp*dzp_seg;
1670 dyp_seg = dyp/dzp*dzp_seg;
1671 x0_new = x0_old + dxp_seg;
1672 y0_new = y0_old + dyp_seg;
1673 }
1674
1675 }
1676
1677 // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, dyp, or dzp)
1678 const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1679 const auto seg_factor_y = static_cast<double>(dyp == 0. ? 1. : dyp_seg/dyp);
1680 const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1681
1682 // compute cell-based weights using the average segment position
1683 double sx_cell[depos_order] = {0.};
1684 double sy_cell[depos_order] = {0.};
1685 double sz_cell[depos_order] = {0.};
1686 double const x0_bar = (x0_new + x0_old)/2.0;
1687 double const y0_bar = (y0_new + y0_old)/2.0;
1688 double const z0_bar = (z0_new + z0_old)/2.0;
1689 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1690 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
1691 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1692
1693 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1694 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1695 double sx_old_cell[depos_order] = {0.};
1696 double sx_new_cell[depos_order] = {0.};
1697 double sy_old_cell[depos_order] = {0.};
1698 double sy_new_cell[depos_order] = {0.};
1699 double sz_old_cell[depos_order] = {0.};
1700 double sz_new_cell[depos_order] = {0.};
1701 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1702 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
1703 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1704 amrex::ignore_unused(i0_cell_2, j0_cell_2, k0_cell_2);
1705 for (int m=0; m<depos_order; m++) {
1706 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1707 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
1708 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1709 }
1710 }
1711
1712 // compute node-based weights using the old and new segment positions
1713 double sx_old_node[depos_order+1] = {0.};
1714 double sx_new_node[depos_order+1] = {0.};
1715 double sy_old_node[depos_order+1] = {0.};
1716 double sy_new_node[depos_order+1] = {0.};
1717 double sz_old_node[depos_order+1] = {0.};
1718 double sz_new_node[depos_order+1] = {0.};
1719 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1720 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
1721 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1722
1723 // deposit Jx for this segment
1724 amrex::Real this_Jx;
1725 for (int i=0; i<=depos_order-1; i++) {
1726 for (int j=0; j<=depos_order; j++) {
1727 for (int k=0; k<=depos_order; k++) {
1728 this_Jx = wqx*sx_cell[i]*( sy_old_node[j]*sz_old_node[k]*one_third
1729 + sy_old_node[j]*sz_new_node[k]*one_sixth
1730 + sy_new_node[j]*sz_old_node[k]*one_sixth
1731 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
1732 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k), this_Jx);
1733 }
1734 }
1735 }
1736
1737 // deposit Jy for this segment
1738 amrex::Real this_Jy;
1739 for (int i=0; i<=depos_order; i++) {
1740 for (int j=0; j<=depos_order-1; j++) {
1741 for (int k=0; k<=depos_order; k++) {
1742 this_Jy = wqy*sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
1743 + sx_old_node[i]*sz_new_node[k]*one_sixth
1744 + sx_new_node[i]*sz_old_node[k]*one_sixth
1745 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1746 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k), this_Jy);
1747 }
1748 }
1749 }
1750
1751 // deposit Jz for this segment
1752 amrex::Real this_Jz;
1753 for (int i=0; i<=depos_order; i++) {
1754 for (int j=0; j<=depos_order; j++) {
1755 for (int k=0; k<=depos_order-1; k++) {
1756 this_Jz = wqz*sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
1757 + sx_old_node[i]*sy_new_node[j]*one_sixth
1758 + sx_new_node[i]*sy_old_node[j]*one_sixth
1759 + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
1760 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k), this_Jz);
1761 }
1762 }
1763 }
1764
1765 // update old segment values
1766 if (ns < num_segments-1) {
1767 x0_old = x0_new;
1768 y0_old = y0_new;
1769 z0_old = z0_new;
1770 }
1771
1772 } // end loop over segments
1773
1774#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1775
1776 // compute cell crossings in X-direction
1777 const auto i_old = static_cast<int>(x_old-shift);
1778 const auto i_new = static_cast<int>(x_new-shift);
1779 const int cell_crossings_x = std::abs(i_new-i_old);
1780 num_segments += cell_crossings_x;
1781
1782 // compute cell crossings in Z-direction
1783 const auto k_old = static_cast<int>(z_old-shift);
1784 const auto k_new = static_cast<int>(z_new-shift);
1785 const int cell_crossings_z = std::abs(k_new-k_old);
1786 num_segments += cell_crossings_z;
1787
1788 // need to assert that the number of cell crossings in each direction
1789 // is within the range permitted by the number of guard cells
1790 // e.g., if (num_segments > 5) ...
1791
1792 // compute total change in particle position and the initial cell
1793 // locations in each direction used to find the position at cell crossings.
1794 const double dxp = x_new - x_old;
1795 const double dzp = z_new - z_old;
1796 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1797 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1798 double Xcell = 0., Zcell = 0.;
1799 if (num_segments > 1) {
1800 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1801 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1802 }
1803
1804 // loop over the number of segments and deposit
1805 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1806 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1807 double dxp_seg, dzp_seg;
1808 double x0_new, z0_new;
1809 double x0_old = x_old;
1810 double z0_old = z_old;
1811
1812 for (int ns=0; ns<num_segments; ns++) {
1813
1814 if (ns == num_segments-1) { // final segment
1815
1816 x0_new = x_new;
1817 z0_new = z_new;
1818 dxp_seg = x0_new - x0_old;
1819 dzp_seg = z0_new - z0_old;
1820
1821 }
1822 else {
1823
1824 x0_new = Xcell + dirX_sign;
1825 z0_new = Zcell + dirZ_sign;
1826 dxp_seg = x0_new - x0_old;
1827 dzp_seg = z0_new - z0_old;
1828
1829 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1830 Xcell = x0_new;
1831 dzp_seg = dzp/dxp*dxp_seg;
1832 z0_new = z0_old + dzp_seg;
1833 }
1834 else {
1835 Zcell = z0_new;
1836 dxp_seg = dxp/dzp*dzp_seg;
1837 x0_new = x0_old + dxp_seg;
1838 }
1839
1840 }
1841
1842 // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1843 const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1844 const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1845
1846 // compute cell-based weights using the average segment position
1847 double sx_cell[depos_order] = {0.};
1848 double sz_cell[depos_order] = {0.};
1849 double const x0_bar = (x0_new + x0_old)/2.0;
1850 double const z0_bar = (z0_new + z0_old)/2.0;
1851 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1852 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1853
1854 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1855 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1856 double sx_old_cell[depos_order] = {0.};
1857 double sx_new_cell[depos_order] = {0.};
1858 double sz_old_cell[depos_order] = {0.};
1859 double sz_new_cell[depos_order] = {0.};
1860 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1861 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1862 amrex::ignore_unused(i0_cell_2, k0_cell_2);
1863 for (int m=0; m<depos_order; m++) {
1864 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1865 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1866 }
1867 }
1868
1869 // compute node-based weights using the old and new segment positions
1870 double sx_old_node[depos_order+1] = {0.};
1871 double sx_new_node[depos_order+1] = {0.};
1872 double sz_old_node[depos_order+1] = {0.};
1873 double sz_new_node[depos_order+1] = {0.};
1874 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1875 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1876
1877 // deposit Jx for this segment
1878 amrex::Real this_Jx;
1879 for (int i=0; i<=depos_order-1; i++) {
1880 for (int k=0; k<=depos_order; k++) {
1881 this_Jx = wqx*sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1882 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 0), this_Jx);
1883#if defined(WARPX_DIM_RZ)
1884 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1885 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1886 // The factor 2 comes from the normalization of the modes
1887 const Complex djr_cmplx = 2._rt*this_Jx*xy_mid;
1888 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode-1), djr_cmplx.real());
1889 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode), djr_cmplx.imag());
1890 xy_mid = xy_mid*xy_mid0;
1891 }
1892#endif
1893 }
1894 }
1895
1896 // deposit out-of-plane Jy for this segment
1897 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1898 amrex::Real this_Jy;
1899 for (int i=0; i<=depos_order; i++) {
1900 for (int k=0; k<=depos_order; k++) {
1901 this_Jy = wqy*( sx_old_node[i]*sz_old_node[k]*one_third
1902 + sx_old_node[i]*sz_new_node[k]*one_sixth
1903 + sx_new_node[i]*sz_old_node[k]*one_sixth
1904 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1905 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 0), this_Jy);
1906#if defined(WARPX_DIM_RZ)
1907 Complex xy_mid = xy_mid0;
1908 // Throughout the following loop, xy_ takes the value e^{i m theta_}
1909 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1910 // The factor 2 comes from the normalization of the modes
1911 const Complex djy_cmplx = 2._rt*this_Jy*xy_mid;
1912 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode-1), djy_cmplx.real());
1913 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode), djy_cmplx.imag());
1914 xy_mid = xy_mid*xy_mid0;
1915 }
1916#endif
1917 }
1918 }
1919
1920 // deposit Jz for this segment
1921 amrex::Real this_Jz;
1922 for (int i=0; i<=depos_order; i++) {
1923 for (int k=0; k<=depos_order-1; k++) {
1924 this_Jz = wqz*sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1925 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 0), this_Jz);
1926#if defined(WARPX_DIM_RZ)
1927 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1928 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1929 // The factor 2 comes from the normalization of the modes
1930 const Complex djz_cmplx = 2._rt*this_Jz*xy_mid;
1931 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode-1), djz_cmplx.real());
1932 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode), djz_cmplx.imag());
1933 xy_mid = xy_mid*xy_mid0;
1934 }
1935#endif
1936 }
1937 }
1938
1939 // update old segment values
1940 if (ns < num_segments-1) {
1941 x0_old = x0_new;
1942 z0_old = z0_new;
1943 }
1944
1945 } // end loop over segments
1946
1947#elif defined(WARPX_DIM_1D_Z)
1948
1949 // compute cell crossings in Z-direction
1950 const auto k_old = static_cast<int>(z_old-shift);
1951 const auto k_new = static_cast<int>(z_new-shift);
1952 const int cell_crossings_z = std::abs(k_new-k_old);
1953 num_segments += cell_crossings_z;
1954
1955 // need to assert that the number of cell crossings in each direction
1956 // is within the range permitted by the number of guard cells
1957 // e.g., if (num_segments > 3) ...
1958
1959 // compute dzp and the initial cell location used to find the cell crossings.
1960 double const dzp = z_new - z_old;
1961 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1962 double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1963
1964 // loop over the number of segments and deposit
1965 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1966 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1967 double dzp_seg;
1968 double z0_new;
1969 double z0_old = z_old;
1970
1971 for (int ns=0; ns<num_segments; ns++) {
1972
1973 if (ns == num_segments-1) { // final segment
1974 z0_new = z_new;
1975 dzp_seg = z0_new - z0_old;
1976 }
1977 else {
1978 Zcell = Zcell + dirZ_sign;
1979 z0_new = Zcell;
1980 dzp_seg = z0_new - z0_old;
1981 }
1982
1983 // compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1984 const auto seg_factor = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1985
1986 // compute cell-based weights using the average segment position
1987 double sz_cell[depos_order] = {0.};
1988 double const z0_bar = (z0_new + z0_old)/2.0;
1989 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1990
1991 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1992 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1993 double sz_old_cell[depos_order] = {0.};
1994 double sz_new_cell[depos_order] = {0.};
1995 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1996 amrex::ignore_unused(k0_cell_2);
1997 for (int m=0; m<depos_order; m++) {
1998 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1999 }
2000 }
2001
2002 // compute node-based weights using the old and new segment positions
2003 double sz_old_node[depos_order+1] = {0.};
2004 double sz_new_node[depos_order+1] = {0.};
2005 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
2006
2007 // deposit out-of-plane Jx and Jy for this segment
2008 for (int k=0; k<=depos_order; k++) {
2009 const amrex::Real weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
2010 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+k0_node+k, 0, 0), wqx*weight);
2011 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+k0_node+k, 0, 0), wqy*weight);
2012 }
2013
2014 // deposit Jz for this segment
2015 for (int k=0; k<=depos_order-1; k++) {
2016 const amrex::Real this_Jz = wqz*sz_cell[k]*seg_factor;
2017 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k0_cell+k, 0, 0), this_Jz);
2018 }
2019
2020 // update old segment values
2021 if (ns < num_segments-1) {
2022 z0_old = z0_new;
2023 }
2024
2025 }
2026
2027#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
2028
2029 // compute cell crossings in X-direction
2030 const auto i_old = static_cast<int>(x_old-shift);
2031 const auto i_new = static_cast<int>(x_new-shift);
2032 const int cell_crossings_x = std::abs(i_new-i_old);
2033 num_segments += cell_crossings_x;
2034
2035 // need to assert that the number of cell crossings in each direction
2036 // is within the range permitted by the number of guard cells
2037 // e.g., if (num_segments > 3) ...
2038
2039 // compute dxp and the initial cell location used to find the cell crossings.
2040 double const dxp = x_new - x_old;
2041 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
2042 double Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
2043
2044 // loop over the number of segments and deposit
2045 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
2046 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
2047 double dxp_seg;
2048 double x0_new;
2049 double x0_old = x_old;
2050
2051 for (int ns=0; ns<num_segments; ns++) {
2052
2053 if (ns == num_segments-1) { // final segment
2054 x0_new = x_new;
2055 dxp_seg = x0_new - x0_old;
2056 }
2057 else {
2058 Xcell = Xcell + dirX_sign;
2059 x0_new = Xcell;
2060 dxp_seg = x0_new - x0_old;
2061 }
2062
2063 // compute the segment factor (equal to dt_seg/dt for nonzero dxp)
2064 const auto seg_factor = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
2065
2066 // compute cell-based weights using the average segment position
2067 double sx_cell[depos_order] = {0.};
2068 double const x0_bar = (x0_new + x0_old)/2.0;
2069 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
2070
2071 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
2072 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
2073 double sx_old_cell[depos_order] = {0.};
2074 double sx_new_cell[depos_order] = {0.};
2075 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
2076 amrex::ignore_unused(i0_cell_2);
2077 for (int m=0; m<depos_order; m++) {
2078 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
2079 }
2080 }
2081
2082 // compute node-based weights using the old and new segment positions
2083 double sx_old_node[depos_order+1] = {0.};
2084 double sx_new_node[depos_order+1] = {0.};
2085 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
2086
2087 // deposit out-of-plane Jy and Jz for this segment
2088 for (int i=0; i<=depos_order; i++) {
2089 const amrex::Real weight = 0.5_rt*(sx_old_node[i] + sx_new_node[i])*seg_factor;
2090 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, 0, 0), wqy*weight);
2091 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, 0, 0), wqz*weight);
2092 }
2093
2094 // deposit Jx for this segment
2095 for (int i=0; i<=depos_order-1; i++) {
2096 const amrex::Real this_Jx = wqx*sx_cell[i]*seg_factor;
2097 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, 0, 0), this_Jx);
2098 }
2099
2100 // update old segment values
2101 if (ns < num_segments-1) {
2102 x0_old = x0_new;
2103 }
2104
2105 }
2106
2107#endif
2108}
2109
2135template <int depos_order>
2136void doVillasenorDepositionShapeNExplicit (const GetParticlePosition<PIdx>& GetPosition,
2137 const amrex::ParticleReal * const wp,
2138 [[maybe_unused]]const amrex::ParticleReal * const uxp,
2139 [[maybe_unused]]const amrex::ParticleReal * const uyp,
2140 [[maybe_unused]]const amrex::ParticleReal * const uzp,
2141 const int * const ion_lev,
2142 const amrex::Array4<amrex::Real>& Jx_arr,
2143 const amrex::Array4<amrex::Real>& Jy_arr,
2144 const amrex::Array4<amrex::Real>& Jz_arr,
2145 const long np_to_deposit,
2146 const amrex::Real dt,
2147 const amrex::Real relative_time,
2148 const amrex::XDim3 & dinv,
2149 const amrex::XDim3 & xyzmin,
2150 const amrex::Dim3 lo,
2151 const amrex::Real q,
2152 [[maybe_unused]] const int n_rz_azimuthal_modes)
2153{
2154 using namespace amrex::literals;
2155
2156 // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
2157 // (do_ionization=1)
2158 bool const do_ionization = ion_lev;
2159
2160 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
2161
2162 // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
2163 amrex::ParallelFor(
2164 np_to_deposit,
2165 [=] AMREX_GPU_DEVICE (long const ip) {
2166
2167 constexpr amrex::ParticleReal inv_c2 = 1._prt/(PhysConst::c*PhysConst::c);
2168 const amrex::Real gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*inv_c2
2169 + uyp[ip]*uyp[ip]*inv_c2
2170 + uzp[ip]*uzp[ip]*inv_c2);
2171
2172 amrex::Real wq = q*wp[ip];
2173 if (do_ionization){
2174 wq *= ion_lev[ip];
2175 }
2176
2177 amrex::ParticleReal xp, yp, zp;
2178 GetPosition(ip, xp, yp, zp);
2179
2180 // computes current and old position
2181 amrex::Real const xp_new = xp + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv;
2182 amrex::Real const xp_old = xp_new - dt*uxp[ip]*gaminv;
2183 amrex::Real const yp_new = yp + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv;
2184 amrex::Real const yp_old = yp_new - dt*uyp[ip]*gaminv;
2185 amrex::Real const zp_new = zp + (relative_time + 0.5_rt*dt)*uzp[ip]*gaminv;
2186 amrex::Real const zp_old = zp_new - dt*uzp[ip]*gaminv;
2187
2188 VillasenorDepositionShapeNKernel<depos_order>(xp_old, yp_old, zp_old, xp_new, yp_new, zp_new, wq,
2189 uxp[ip], uyp[ip], uzp[ip], gaminv,
2190 Jx_arr, Jy_arr, Jz_arr,
2191 dt, dinv, xyzmin, lo, invvol, n_rz_azimuthal_modes);
2192
2193 });
2194}
2195
2222template <int depos_order>
2223void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::ParticleReal * const xp_n_data,
2224 [[maybe_unused]]const amrex::ParticleReal * const yp_n_data,
2225 [[maybe_unused]]const amrex::ParticleReal * const zp_n_data,
2226 const GetParticlePosition<PIdx>& GetPosition,
2227 const amrex::ParticleReal * const wp,
2228 [[maybe_unused]]const amrex::ParticleReal * const uxp_n,
2229 [[maybe_unused]]const amrex::ParticleReal * const uyp_n,
2230 [[maybe_unused]]const amrex::ParticleReal * const uzp_n,
2231 [[maybe_unused]]const amrex::ParticleReal * const uxp_nph,
2232 [[maybe_unused]]const amrex::ParticleReal * const uyp_nph,
2233 [[maybe_unused]]const amrex::ParticleReal * const uzp_nph,
2234 const int * const ion_lev,
2235 const amrex::Array4<amrex::Real>& Jx_arr,
2236 const amrex::Array4<amrex::Real>& Jy_arr,
2237 const amrex::Array4<amrex::Real>& Jz_arr,
2238 const long np_to_deposit,
2239 const amrex::Real dt,
2240 const amrex::XDim3 & dinv,
2241 const amrex::XDim3 & xyzmin,
2242 const amrex::Dim3 lo,
2243 const amrex::Real q,
2244 [[maybe_unused]] const int n_rz_azimuthal_modes)
2245{
2246 using namespace amrex::literals;
2247
2248 // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
2249 // (do_ionization=1)
2250 bool const do_ionization = ion_lev;
2251
2252 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
2253
2254 // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
2255 amrex::ParallelFor(
2256 np_to_deposit,
2257 [=] AMREX_GPU_DEVICE (long const ip) {
2258
2259 // Skip particles with zero weight.
2260 // This should only be the case for particles that will be suborbited.
2261 if (wp[ip] == 0.) { return; }
2262
2263#if !defined(WARPX_DIM_3D)
2264
2265 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
2266 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
2267 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
2268#else
2269 // gaminv is unused in 3D
2270 const amrex::ParticleReal gaminv = 1.;
2271#endif
2272
2273 amrex::Real wq = q*wp[ip];
2274 if (do_ionization){
2275 wq *= ion_lev[ip];
2276 }
2277
2278 amrex::ParticleReal xp_nph, yp_nph, zp_nph;
2279 GetPosition(ip, xp_nph, yp_nph, zp_nph);
2280
2281 amrex::ParticleReal const xp_n = (xp_n_data ? xp_n_data[ip] : 0._prt);
2282 amrex::ParticleReal const yp_n = (yp_n_data ? yp_n_data[ip] : 0._prt);
2283 amrex::ParticleReal const zp_n = (zp_n_data ? zp_n_data[ip] : 0._prt);
2284
2285 amrex::ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n;
2286 amrex::ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n;
2287 amrex::ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n;
2288
2289 VillasenorDepositionShapeNKernel<depos_order>(xp_n, yp_n, zp_n, xp_np1, yp_np1, zp_np1, wq,
2290 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip], gaminv,
2291 Jx_arr, Jy_arr, Jz_arr,
2292 dt, dinv, xyzmin, lo, invvol, n_rz_azimuthal_modes);
2293
2294 });
2295}
2296
2323template <int depos_order>
2324void doVayDepositionShapeN (const GetParticlePosition<PIdx>& GetPosition,
2325 const amrex::ParticleReal* const wp,
2326 const amrex::ParticleReal* const uxp,
2327 const amrex::ParticleReal* const uyp,
2328 const amrex::ParticleReal* const uzp,
2329 const int* const ion_lev,
2330 amrex::FArrayBox& Dx_fab,
2331 amrex::FArrayBox& Dy_fab,
2332 amrex::FArrayBox& Dz_fab,
2333 long np_to_deposit,
2334 amrex::Real dt,
2335 amrex::Real relative_time,
2336 const amrex::XDim3 & dinv,
2337 const amrex::XDim3 & xyzmin,
2338 amrex::Dim3 lo,
2339 amrex::Real q,
2340 [[maybe_unused]]int n_rz_azimuthal_modes)
2341{
2342 using namespace amrex::literals;
2343
2344#if defined(WARPX_DIM_RZ)
2345 amrex::ignore_unused(GetPosition,
2346 wp, uxp, uyp, uzp, ion_lev, Dx_fab, Dy_fab, Dz_fab,
2347 np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q);
2348 WARPX_ABORT_WITH_MESSAGE("Vay deposition not implemented in RZ geometry");
2349#endif
2350
2351#if defined(WARPX_DIM_1D_Z) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
2352 amrex::ignore_unused(GetPosition,
2353 wp, uxp, uyp, uzp, ion_lev, Dx_fab, Dy_fab, Dz_fab,
2354 np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q);
2355 WARPX_ABORT_WITH_MESSAGE("Vay deposition not implemented in 1D geometry");
2356#endif
2357
2358#if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z || defined WARPX_DIM_RCYLINDER || defined WARPX_DIM_RSPHERE)
2359
2360 // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1
2361 const bool do_ionization = ion_lev;
2362
2363 // Inverse of time step
2364 const amrex::Real invdt = 1._rt / dt;
2365
2366 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
2367
2368 // Allocate temporary arrays
2369#if defined(WARPX_DIM_3D)
2370 AMREX_ALWAYS_ASSERT(Dx_fab.box() == Dy_fab.box() && Dx_fab.box() == Dz_fab.box());
2371 amrex::FArrayBox temp_fab{Dx_fab.box(), 4};
2372#elif defined(WARPX_DIM_XZ)
2373 AMREX_ALWAYS_ASSERT(Dx_fab.box() == Dz_fab.box());
2374 amrex::FArrayBox temp_fab{Dx_fab.box(), 2};
2375#endif
2376 temp_fab.setVal<amrex::RunOn::Device>(0._rt);
2377 amrex::Array4<amrex::Real> const& temp_arr = temp_fab.array();
2378
2379 // Arrays where D will be stored
2380 amrex::Array4<amrex::Real> const& Dx_arr = Dx_fab.array();
2381 amrex::Array4<amrex::Real> const& Dy_arr = Dy_fab.array();
2382 amrex::Array4<amrex::Real> const& Dz_arr = Dz_fab.array();
2383
2384 // Loop over particles and deposit (Dx,Dy,Dz) into Dx_fab, Dy_fab and Dz_fab
2385 amrex::ParallelFor(np_to_deposit, [=] AMREX_GPU_DEVICE (long ip)
2386 {
2387 // Inverse of Lorentz factor gamma
2388 const amrex::Real invgam = 1._rt / std::sqrt(1._rt + uxp[ip] * uxp[ip] * PhysConst::inv_c2
2389 + uyp[ip] * uyp[ip] * PhysConst::inv_c2
2390 + uzp[ip] * uzp[ip] * PhysConst::inv_c2);
2391 // Product of particle charges and weights
2392 amrex::Real wq = q * wp[ip];
2393 if (do_ionization) { wq *= ion_lev[ip]; }
2394
2395 // Current particle positions (in physical units)
2396 amrex::ParticleReal xp, yp, zp;
2397 GetPosition(ip, xp, yp, zp);
2398
2399 // Particle velocities
2400 const amrex::Real vx = uxp[ip] * invgam;
2401 const amrex::Real vy = uyp[ip] * invgam;
2402 const amrex::Real vz = uzp[ip] * invgam;
2403
2404 // Modify the particle position to match the time of the deposition
2405 xp += relative_time * vx;
2406 yp += relative_time * vy;
2407 zp += relative_time * vz;
2408
2409 // Current and old particle positions in grid units
2410 // Keep these double to avoid bug in single precision.
2411 double const x_new = (xp - xyzmin.x + 0.5_rt*dt*vx) * dinv.x;
2412 double const x_old = (xp - xyzmin.x - 0.5_rt*dt*vx) * dinv.x;
2413#if defined(WARPX_DIM_3D)
2414 // Keep these double to avoid bug in single precision.
2415 double const y_new = (yp - xyzmin.y + 0.5_rt*dt*vy) * dinv.y;
2416 double const y_old = (yp - xyzmin.y - 0.5_rt*dt*vy) * dinv.y;
2417#endif
2418 // Keep these double to avoid bug in single precision.
2419 double const z_new = (zp - xyzmin.z + 0.5_rt*dt*vz) * dinv.z;
2420 double const z_old = (zp - xyzmin.z - 0.5_rt*dt*vz) * dinv.z;
2421
2422 // Shape factor arrays for current and old positions (nodal)
2423 // Keep these double to avoid bug in single precision.
2424 double sx_new[depos_order+1] = {0.};
2425 double sx_old[depos_order+1] = {0.};
2426#if defined(WARPX_DIM_3D)
2427 // Keep these double to avoid bug in single precision.
2428 double sy_new[depos_order+1] = {0.};
2429 double sy_old[depos_order+1] = {0.};
2430#endif
2431 // Keep these double to avoid bug in single precision.
2432 double sz_new[depos_order+1] = {0.};
2433 double sz_old[depos_order+1] = {0.};
2434
2435 // Compute shape factors for current positions
2436
2437 // i_new leftmost grid point in x that the particle touches
2438 // sx_new shape factor along x for the centering of each current
2439 Compute_shape_factor< depos_order > const compute_shape_factor;
2440 const int i_new = compute_shape_factor(sx_new, x_new);
2441#if defined(WARPX_DIM_3D)
2442 // j_new leftmost grid point in y that the particle touches
2443 // sy_new shape factor along y for the centering of each current
2444 const int j_new = compute_shape_factor(sy_new, y_new);
2445#endif
2446 // k_new leftmost grid point in z that the particle touches
2447 // sz_new shape factor along z for the centering of each current
2448 const int k_new = compute_shape_factor(sz_new, z_new);
2449
2450 // Compute shape factors for old positions
2451
2452 // i_old leftmost grid point in x that the particle touches
2453 // sx_old shape factor along x for the centering of each current
2454 const int i_old = compute_shape_factor(sx_old, x_old);
2455#if defined(WARPX_DIM_3D)
2456 // j_old leftmost grid point in y that the particle touches
2457 // sy_old shape factor along y for the centering of each current
2458 const int j_old = compute_shape_factor(sy_old, y_old);
2459#endif
2460 // k_old leftmost grid point in z that the particle touches
2461 // sz_old shape factor along z for the centering of each current
2462 const int k_old = compute_shape_factor(sz_old, z_old);
2463
2464 // Deposit current into Dx_arr, Dy_arr and Dz_arr
2465#if defined(WARPX_DIM_XZ)
2466
2467 const amrex::Real wqy = wq * vy * invvol;
2468 for (int k=0; k<=depos_order; k++) {
2469 for (int i=0; i<=depos_order; i++) {
2470
2471 // Re-casting sx_new and sz_new from double to amrex::Real so that
2472 // Atomic::Add has consistent types in its argument
2473 auto const sxn_szn = static_cast<amrex::Real>(sx_new[i] * sz_new[k]);
2474 auto const sxo_szn = static_cast<amrex::Real>(sx_old[i] * sz_new[k]);
2475 auto const sxn_szo = static_cast<amrex::Real>(sx_new[i] * sz_old[k]);
2476 auto const sxo_szo = static_cast<amrex::Real>(sx_old[i] * sz_old[k]);
2477
2478 if (i_new == i_old && k_new == k_old) {
2479 // temp arrays for Dx and Dz
2480 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
2481 wq * invvol * invdt * (sxn_szn - sxo_szo));
2482
2483 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 1),
2484 wq * invvol * invdt * (sxn_szo - sxo_szn));
2485
2486 // Dy
2487 amrex::Gpu::Atomic::AddNoRet(&Dy_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
2488 wqy * 0.25_rt * (sxn_szn + sxn_szo + sxo_szn + sxo_szo));
2489 } else {
2490 // temp arrays for Dx and Dz
2491 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
2492 wq * invvol * invdt * sxn_szn);
2493
2494 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
2495 - wq * invvol * invdt * sxo_szo);
2496
2497 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 1),
2498 wq * invvol * invdt * sxn_szo);
2499
2500 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 1),
2501 - wq * invvol * invdt * sxo_szn);
2502
2503 // Dy
2504 amrex::Gpu::Atomic::AddNoRet(&Dy_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0),
2505 wqy * 0.25_rt * sxn_szn);
2506
2507 amrex::Gpu::Atomic::AddNoRet(&Dy_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0),
2508 wqy * 0.25_rt * sxn_szo);
2509
2510 amrex::Gpu::Atomic::AddNoRet(&Dy_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0),
2511 wqy * 0.25_rt * sxo_szn);
2512
2513 amrex::Gpu::Atomic::AddNoRet(&Dy_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0),
2514 wqy * 0.25_rt * sxo_szo);
2515 }
2516
2517 }
2518 }
2519
2520#elif defined(WARPX_DIM_3D)
2521
2522 for (int k=0; k<=depos_order; k++) {
2523 for (int j=0; j<=depos_order; j++) {
2524
2525 auto const syn_szn = static_cast<amrex::Real>(sy_new[j] * sz_new[k]);
2526 auto const syo_szn = static_cast<amrex::Real>(sy_old[j] * sz_new[k]);
2527 auto const syn_szo = static_cast<amrex::Real>(sy_new[j] * sz_old[k]);
2528 auto const syo_szo = static_cast<amrex::Real>(sy_old[j] * sz_old[k]);
2529
2530 for (int i=0; i<=depos_order; i++) {
2531
2532 auto const sxn_syn_szn = static_cast<amrex::Real>(sx_new[i]) * syn_szn;
2533 auto const sxo_syn_szn = static_cast<amrex::Real>(sx_old[i]) * syn_szn;
2534 auto const sxn_syo_szn = static_cast<amrex::Real>(sx_new[i]) * syo_szn;
2535 auto const sxo_syo_szn = static_cast<amrex::Real>(sx_old[i]) * syo_szn;
2536 auto const sxn_syn_szo = static_cast<amrex::Real>(sx_new[i]) * syn_szo;
2537 auto const sxo_syn_szo = static_cast<amrex::Real>(sx_old[i]) * syn_szo;
2538 auto const sxn_syo_szo = static_cast<amrex::Real>(sx_new[i]) * syo_szo;
2539 auto const sxo_syo_szo = static_cast<amrex::Real>(sx_old[i]) * syo_szo;
2540
2541 if (i_new == i_old && j_new == j_old && k_new == k_old) {
2542 // temp arrays for Dx, Dy and Dz
2543 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k, 0),
2544 wq * invvol * invdt * (sxn_syn_szn - sxo_syo_szo));
2545
2546 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k, 1),
2547 wq * invvol * invdt * (sxn_syn_szo - sxo_syo_szn));
2548
2549 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k, 2),
2550 wq * invvol * invdt * (sxn_syo_szn - sxo_syn_szo));
2551
2552 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k, 3),
2553 wq * invvol * invdt * (sxo_syn_szn - sxn_syo_szo));
2554 } else {
2555 // temp arrays for Dx, Dy and Dz
2556 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k, 0),
2557 wq * invvol * invdt * sxn_syn_szn);
2558
2559 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k, 0),
2560 - wq * invvol * invdt * sxo_syo_szo);
2561
2562 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k, 1),
2563 wq * invvol * invdt * sxn_syn_szo);
2564
2565 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k, 1),
2566 - wq * invvol * invdt * sxo_syo_szn);
2567
2568 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k, 2),
2569 wq * invvol * invdt * sxn_syo_szn);
2570
2571 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k, 2),
2572 - wq * invvol * invdt * sxo_syn_szo);
2573
2574 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k, 3),
2575 wq * invvol * invdt * sxo_syn_szn);
2576
2577 amrex::Gpu::Atomic::AddNoRet(&temp_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k, 3),
2578 - wq * invvol * invdt * sxn_syo_szo);
2579 }
2580 }
2581 }
2582 }
2583#endif
2584 } );
2585
2586#if defined(WARPX_DIM_3D)
2587 amrex::ParallelFor(Dx_fab.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
2588 {
2589 const amrex::Real t_a = temp_arr(i,j,k,0);
2590 const amrex::Real t_b = temp_arr(i,j,k,1);
2591 const amrex::Real t_c = temp_arr(i,j,k,2);
2592 const amrex::Real t_d = temp_arr(i,j,k,3);
2593 Dx_arr(i,j,k) += (1._rt/6._rt)*(2_rt*t_a + t_b + t_c - 2._rt*t_d);
2594 Dy_arr(i,j,k) += (1._rt/6._rt)*(2_rt*t_a + t_b - 2._rt*t_c + t_d);
2595 Dz_arr(i,j,k) += (1._rt/6._rt)*(2_rt*t_a - 2._rt*t_b + t_c + t_d);
2596 });
2597#elif defined(WARPX_DIM_XZ)
2598 amrex::ParallelFor(Dx_fab.box(), [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
2599 {
2600 const amrex::Real t_a = temp_arr(i,j,0,0);
2601 const amrex::Real t_b = temp_arr(i,j,0,1);
2602 Dx_arr(i,j,0) += (0.5_rt)*(t_a + t_b);
2603 Dz_arr(i,j,0) += (0.5_rt)*(t_a - t_b);
2604 });
2605#endif
2606 // Synchronize so that temp_fab can be safely deallocated in its destructor
2607 amrex::Gpu::streamSynchronize();
2608
2609#endif // #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z)
2610}
2611#endif // WARPX_CURRENTDEPOSITION_H_
AMReX_Arena.H
AMReX_Array4.H
AMREX_ALWAYS_ASSERT
#define AMREX_ALWAYS_ASSERT(EX)
AMReX_Dim3.H
AMREX_INLINE
#define AMREX_INLINE
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
AMReX_REAL.H
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
VillasenorDepositionShapeNKernel
AMREX_GPU_HOST_DEVICE AMREX_INLINE void VillasenorDepositionShapeNKernel(amrex::ParticleReal const xp_old, amrex::ParticleReal const yp_old, amrex::ParticleReal const zp_old, amrex::ParticleReal const xp_new, amrex::ParticleReal const yp_new, amrex::ParticleReal const zp_new, amrex::ParticleReal const wq, amrex::ParticleReal const uxp_mid, amrex::ParticleReal const uyp_mid, amrex::ParticleReal const uzp_mid, amrex::ParticleReal const gaminv, amrex::Array4< amrex::Real >const &Jx_arr, amrex::Array4< amrex::Real >const &Jy_arr, amrex::Array4< amrex::Real >const &Jz_arr, amrex::Real const dt, amrex::XDim3 const &dinv, amrex::XDim3 const &xyzmin, amrex::Dim3 const lo, amrex::Real const invvol, int const n_rz_azimuthal_modes)
Villasenor and Buneman Current Deposition for thread thread_num kernal. This is a charge-conserving d...
Definition CurrentDeposition.H:1471
doEsirkepovDepositionShapeN
void doEsirkepovDepositionShapeN(const GetParticlePosition< PIdx > &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, long np_to_deposit, amrex::Real dt, amrex::Real relative_time, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes, const amrex::Array4< const int > &reduced_particle_shape_mask, bool enable_reduced_shape)
Esirkepov Current Deposition for thread thread_num.
Definition CurrentDeposition.H:675
doDepositionSharedShapeN
void doDepositionSharedShapeN(const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp, const amrex::ParticleReal *const uyp, const amrex::ParticleReal *const uzp, const int *ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, long np_to_deposit, const amrex::Real relative_time, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes, const amrex::DenseBins< WarpXParticleContainer::ParticleTileType::ParticleTileDataType > &a_bins, const amrex::Box &box, const amrex::Geometry &geom, const amrex::IntVect &a_tbox_max_size, const int threads_per_block, const amrex::IntVect bin_size)
Current Deposition for thread thread_num using shared memory.
Definition CurrentDeposition.H:476
doVillasenorDepositionShapeNExplicit
void doVillasenorDepositionShapeNExplicit(const GetParticlePosition< PIdx > &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_deposit, const amrex::Real dt, const amrex::Real relative_time, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const int n_rz_azimuthal_modes)
Villasenor and Buneman Current Deposition for thread thread_num for explicit scheme....
Definition CurrentDeposition.H:2136
doVayDepositionShapeN
void doVayDepositionShapeN(const GetParticlePosition< PIdx > &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 &Dx_fab, amrex::FArrayBox &Dy_fab, amrex::FArrayBox &Dz_fab, long np_to_deposit, amrex::Real dt, amrex::Real relative_time, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, amrex::Dim3 lo, amrex::Real q, int 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:2324
doChargeConservingDepositionShapeNImplicit
void doChargeConservingDepositionShapeNImplicit(const amrex::ParticleReal *const xp_n, const amrex::ParticleReal *const yp_n, const amrex::ParticleReal *const zp_n, const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp_n, const amrex::ParticleReal *const uyp_n, const amrex::ParticleReal *const uzp_n, const amrex::ParticleReal *const uxp_nph, const amrex::ParticleReal *const uyp_nph, const amrex::ParticleReal *const uzp_nph, 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_deposit, const amrex::Real dt, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const int n_rz_azimuthal_modes)
Esirkepov Current Deposition for thread thread_num for implicit scheme The difference from doEsirkepo...
Definition CurrentDeposition.H:1114
doDepositionShapeNImplicit
void doDepositionShapeNImplicit(const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp_n, const amrex::ParticleReal *const uyp_n, const amrex::ParticleReal *const uzp_n, const amrex::ParticleReal *const uxp_nph, const amrex::ParticleReal *const uyp_nph, const amrex::ParticleReal *const uzp_nph, const int *const ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, const long np_to_deposit, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const int n_rz_azimuthal_modes)
Direct current deposition for thread thread_num for the implicit scheme The only difference from doDe...
Definition CurrentDeposition.H:383
doDepositionShapeNKernel
AMREX_GPU_HOST_DEVICE AMREX_INLINE void doDepositionShapeNKernel(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, const amrex::ParticleReal wq, const amrex::ParticleReal vx, const amrex::ParticleReal vy, const amrex::ParticleReal vz, amrex::Array4< amrex::Real > const &jx_arr, amrex::Array4< amrex::Real > const &jy_arr, amrex::Array4< amrex::Real > const &jz_arr, amrex::IntVect const &jx_type, amrex::IntVect const &jy_type, amrex::IntVect const &jz_type, const amrex::Real relative_time, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Real invvol, const amrex::Dim3 lo, const int n_rz_azimuthal_modes)
Kernel for the direct current deposition for thread thread_num.
Definition CurrentDeposition.H:48
doVillasenorDepositionShapeNImplicit
void doVillasenorDepositionShapeNImplicit(const amrex::ParticleReal *const xp_n_data, const amrex::ParticleReal *const yp_n_data, const amrex::ParticleReal *const zp_n_data, const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp_n, const amrex::ParticleReal *const uyp_n, const amrex::ParticleReal *const uzp_n, const amrex::ParticleReal *const uxp_nph, const amrex::ParticleReal *const uyp_nph, const amrex::ParticleReal *const uzp_nph, 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_deposit, const amrex::Real dt, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const int n_rz_azimuthal_modes)
Villasenor and Buneman Current Deposition for thread thread_num for implicit scheme....
Definition CurrentDeposition.H:2223
doDepositionShapeN
void doDepositionShapeN(const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp, const amrex::ParticleReal *const uyp, const amrex::ParticleReal *const uzp, const int *ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, long np_to_deposit, amrex::Real relative_time, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes)
Current Deposition for thread thread_num.
Definition CurrentDeposition.H:298
RigidAdvanceMode::vz
@ vz
Definition RigidInjectedParticleContainer.H:27
GetImplicitGammaInverse
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal GetImplicitGammaInverse(const amrex::ParticleReal uxp_n, const amrex::ParticleReal uyp_n, const amrex::ParticleReal uzp_n, const amrex::ParticleReal uxp_nph, const amrex::ParticleReal uyp_nph, const amrex::ParticleReal uzp_nph) noexcept
Compute the inverse Lorentz factor for the position update in the implicit methods,...
Definition UpdatePosition.H:77
TextMsg.H
WARPX_ABORT_WITH_MESSAGE
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition TextMsg.H:15
WARPX_ALWAYS_ASSERT_WITH_MESSAGE
#define WARPX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition TextMsg.H:13
WarpX_Complex.H
Complex
amrex::GpuComplex< amrex::Real > Complex
Definition WarpX_Complex.H:22
NODE
NODE
amrex::BaseFab::box
const Box & box() const noexcept
amrex::BaseFab::array
Array4< T const > array() const noexcept
amrex::BaseFab::setVal
void setVal(T const &x, const Box &bx, int dcomp, int ncomp) noexcept
amrex::BoxND::grow
__host__ __device__ BoxND & grow(int i) noexcept
amrex::BoxND::numPts
__host__ __device__ Long numPts() const noexcept
amrex::BoxND::type
__host__ __device__ IntVectND< dim > type() const noexcept
amrex::CoordSys::InvCellSizeArray
GpuArray< Real, 3 > InvCellSizeArray() const noexcept
amrex::DenseBins::offsetsPtr
index_type * offsetsPtr() noexcept
amrex::DenseBins::permutationPtr
index_type * permutationPtr() noexcept
amrex::DenseBins::numBins
Long numBins() const noexcept
amrex::Geometry::Domain
const Box & Domain() const noexcept
amrex::Geometry::ProbLoArray
GpuArray< Real, 3 > ProbLoArray() const noexcept
amrex::Gpu::Device::sharedMemPerBlock
static std::size_t sharedMemPerBlock() noexcept
ablastr::constant::SI::c
constexpr auto c
vacuum speed of light [m/s]
Definition constant.H:149
ablastr::constant::SI::inv_c2
constexpr auto inv_c2
inverse of the square of the vacuum speed of light [s^2/m^2]
Definition constant.H:216
amrex::Gpu::Atomic::AddNoRet
__host__ __device__ AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
amrex::Gpu::streamSynchronize
void streamSynchronize() noexcept
amrex::Gpu::gpuStream
gpuStream_t gpuStream() noexcept
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
amrex::convert
__host__ __device__ BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
amrex::getTileIndex
__host__ __device__ int getTileIndex(const IntVect &iv, const Box &box, const bool a_do_tiling, const IntVect &a_tile_size, Box &tbx)
amrex::ParallelFor
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::Box
BoxND< 3 > Box
amrex::launch
void launch(T const &n, L &&f) noexcept
amrex::begin
__host__ __device__ Dim3 begin(BoxND< dim > const &box) noexcept
amrex::shift
__host__ __device__ BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
amrex::IntVect
IntVectND< 3 > IntVect
amrex::max
__host__ __device__ constexpr const T & max(const T &a, const T &b) noexcept
amrex::end
__host__ __device__ Dim3 end(BoxND< dim > const &box) noexcept
Compute_shape_factor_pair
Definition ShapeFactors.H:168
Compute_shape_factor
Definition ShapeFactors.H:29
Compute_shifted_shape_factor
Definition ShapeFactors.H:95
GetParticlePosition
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel.
Definition GetAndSetPosition.H:75
amrex::GpuComplex::real
__host__ __device__ constexpr T real() const noexcept
amrex::GpuComplex::imag
__host__ __device__ constexpr T imag() const noexcept
amrex::Gpu::SharedMemory::dataPtr
__device__ T * dataPtr() noexcept
Generated by doxygen 1.13.2