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