WarpX
FieldGather.H
Go to the documentation of this file.
1 /* Copyright 2019 Axel Huebl, David Grote, Maxence Thevenet
2  * Revathi Jambunathan, Weiqun Zhang
3  *
4  * This file is part of WarpX.
5  *
6  * License: BSD-3-Clause-LBNL
7  */
8 #ifndef WARPX_FIELDGATHER_H_
9 #define WARPX_FIELDGATHER_H_
10 
13 #include "Particles/ShapeFactors.H"
14 #include "Utils/WarpX_Complex.H"
15 
16 #include <AMReX.H>
17 
36 template <int depos_order, int galerkin_interpolation>
38 void doGatherShapeN (const amrex::ParticleReal xp,
39  const amrex::ParticleReal yp,
40  const amrex::ParticleReal zp,
41  amrex::ParticleReal& Exp,
42  amrex::ParticleReal& Eyp,
43  amrex::ParticleReal& Ezp,
44  amrex::ParticleReal& Bxp,
45  amrex::ParticleReal& Byp,
46  amrex::ParticleReal& Bzp,
53  const amrex::IndexType ex_type,
54  const amrex::IndexType ey_type,
55  const amrex::IndexType ez_type,
56  const amrex::IndexType bx_type,
57  const amrex::IndexType by_type,
58  const amrex::IndexType bz_type,
60  const amrex::GpuArray<amrex::Real, 3>& xyzmin,
61  const amrex::Dim3& lo,
62  const int n_rz_azimuthal_modes)
63 {
64  using namespace amrex;
65 
66 #if defined(WARPX_DIM_XZ)
68 #endif
69 
70 #if defined(WARPX_DIM_1D_Z)
71  amrex::ignore_unused(xp,yp);
72 #endif
73 
74 #ifndef WARPX_DIM_RZ
75  amrex::ignore_unused(n_rz_azimuthal_modes);
76 #endif
77 
78 #if (AMREX_SPACEDIM >= 2)
79  const amrex::Real dxi = 1.0_rt/dx[0];
80 #endif
81  const amrex::Real dzi = 1.0_rt/dx[2];
82 #if defined(WARPX_DIM_3D)
83  const amrex::Real dyi = 1.0_rt/dx[1];
84 #endif
85 
86 #if (AMREX_SPACEDIM >= 2)
87  const amrex::Real xmin = xyzmin[0];
88 #endif
89 #if defined(WARPX_DIM_3D)
90  const amrex::Real ymin = xyzmin[1];
91 #endif
92  const amrex::Real zmin = xyzmin[2];
93 
94  constexpr int zdir = WARPX_ZINDEX;
95  constexpr int NODE = amrex::IndexType::NODE;
96  constexpr int CELL = amrex::IndexType::CELL;
97 
98  // --- Compute shape factors
99 
100  Compute_shape_factor< depos_order > const compute_shape_factor;
101  Compute_shape_factor<depos_order - galerkin_interpolation > const compute_shape_factor_galerkin;
102 
103 #if (AMREX_SPACEDIM >= 2)
104  // x direction
105  // Get particle position
106 #ifdef WARPX_DIM_RZ
107  const amrex::Real rp = std::sqrt(xp*xp + yp*yp);
108  const amrex::Real x = (rp - xmin)*dxi;
109 #else
110  const amrex::Real x = (xp-xmin)*dxi;
111 #endif
112 
113  // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current
114  // sx_[eb][xyz] shape factor along x for the centering of each current
115  // There are only two possible centerings, node or cell centered, so at most only two shape factor
116  // arrays will be needed.
117  amrex::Real sx_node[depos_order + 1];
118  amrex::Real sx_cell[depos_order + 1];
119  amrex::Real sx_node_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt};
120  amrex::Real sx_cell_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt};
121 
122  int j_node = 0;
123  int j_cell = 0;
124  int j_node_v = 0;
125  int j_cell_v = 0;
126  if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) {
127  j_node = compute_shape_factor(sx_node, x);
128  }
129  if ((ey_type[0] == CELL) || (ez_type[0] == CELL) || (bx_type[0] == CELL)) {
130  j_cell = compute_shape_factor(sx_cell, x - 0.5_rt);
131  }
132  if ((ex_type[0] == NODE) || (by_type[0] == NODE) || (bz_type[0] == NODE)) {
133  j_node_v = compute_shape_factor_galerkin(sx_node_galerkin, x);
134  }
135  if ((ex_type[0] == CELL) || (by_type[0] == CELL) || (bz_type[0] == CELL)) {
136  j_cell_v = compute_shape_factor_galerkin(sx_cell_galerkin, x - 0.5_rt);
137  }
138  const amrex::Real (&sx_ex)[depos_order + 1 - galerkin_interpolation] = ((ex_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin);
139  const amrex::Real (&sx_ey)[depos_order + 1 ] = ((ey_type[0] == NODE) ? sx_node : sx_cell );
140  const amrex::Real (&sx_ez)[depos_order + 1 ] = ((ez_type[0] == NODE) ? sx_node : sx_cell );
141  const amrex::Real (&sx_bx)[depos_order + 1 ] = ((bx_type[0] == NODE) ? sx_node : sx_cell );
142  const amrex::Real (&sx_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin);
143  const amrex::Real (&sx_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin);
144  int const j_ex = ((ex_type[0] == NODE) ? j_node_v : j_cell_v);
145  int const j_ey = ((ey_type[0] == NODE) ? j_node : j_cell );
146  int const j_ez = ((ez_type[0] == NODE) ? j_node : j_cell );
147  int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell );
148  int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v);
149  int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v);
150 #endif
151 
152 #if defined(WARPX_DIM_3D)
153  // y direction
154  const amrex::Real y = (yp-ymin)*dyi;
155  amrex::Real sy_node[depos_order + 1];
156  amrex::Real sy_cell[depos_order + 1];
157  amrex::Real sy_node_v[depos_order + 1 - galerkin_interpolation];
158  amrex::Real sy_cell_v[depos_order + 1 - galerkin_interpolation];
159  int k_node = 0;
160  int k_cell = 0;
161  int k_node_v = 0;
162  int k_cell_v = 0;
163  if ((ex_type[1] == NODE) || (ez_type[1] == NODE) || (by_type[1] == NODE)) {
164  k_node = compute_shape_factor(sy_node, y);
165  }
166  if ((ex_type[1] == CELL) || (ez_type[1] == CELL) || (by_type[1] == CELL)) {
167  k_cell = compute_shape_factor(sy_cell, y - 0.5_rt);
168  }
169  if ((ey_type[1] == NODE) || (bx_type[1] == NODE) || (bz_type[1] == NODE)) {
170  k_node_v = compute_shape_factor_galerkin(sy_node_v, y);
171  }
172  if ((ey_type[1] == CELL) || (bx_type[1] == CELL) || (bz_type[1] == CELL)) {
173  k_cell_v = compute_shape_factor_galerkin(sy_cell_v, y - 0.5_rt);
174  }
175  const amrex::Real (&sy_ex)[depos_order + 1 ] = ((ex_type[1] == NODE) ? sy_node : sy_cell );
176  const amrex::Real (&sy_ey)[depos_order + 1 - galerkin_interpolation] = ((ey_type[1] == NODE) ? sy_node_v : sy_cell_v);
177  const amrex::Real (&sy_ez)[depos_order + 1 ] = ((ez_type[1] == NODE) ? sy_node : sy_cell );
178  const amrex::Real (&sy_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[1] == NODE) ? sy_node_v : sy_cell_v);
179  const amrex::Real (&sy_by)[depos_order + 1 ] = ((by_type[1] == NODE) ? sy_node : sy_cell );
180  const amrex::Real (&sy_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[1] == NODE) ? sy_node_v : sy_cell_v);
181  int const k_ex = ((ex_type[1] == NODE) ? k_node : k_cell );
182  int const k_ey = ((ey_type[1] == NODE) ? k_node_v : k_cell_v);
183  int const k_ez = ((ez_type[1] == NODE) ? k_node : k_cell );
184  int const k_bx = ((bx_type[1] == NODE) ? k_node_v : k_cell_v);
185  int const k_by = ((by_type[1] == NODE) ? k_node : k_cell );
186  int const k_bz = ((bz_type[1] == NODE) ? k_node_v : k_cell_v);
187 
188 #endif
189  // z direction
190  const amrex::Real z = (zp-zmin)*dzi;
191  amrex::Real sz_node[depos_order + 1];
192  amrex::Real sz_cell[depos_order + 1];
193  amrex::Real sz_node_v[depos_order + 1 - galerkin_interpolation];
194  amrex::Real sz_cell_v[depos_order + 1 - galerkin_interpolation];
195  int l_node = 0;
196  int l_cell = 0;
197  int l_node_v = 0;
198  int l_cell_v = 0;
199  if ((ex_type[zdir] == NODE) || (ey_type[zdir] == NODE) || (bz_type[zdir] == NODE)) {
200  l_node = compute_shape_factor(sz_node, z);
201  }
202  if ((ex_type[zdir] == CELL) || (ey_type[zdir] == CELL) || (bz_type[zdir] == CELL)) {
203  l_cell = compute_shape_factor(sz_cell, z - 0.5_rt);
204  }
205  if ((ez_type[zdir] == NODE) || (bx_type[zdir] == NODE) || (by_type[zdir] == NODE)) {
206  l_node_v = compute_shape_factor_galerkin(sz_node_v, z);
207  }
208  if ((ez_type[zdir] == CELL) || (bx_type[zdir] == CELL) || (by_type[zdir] == CELL)) {
209  l_cell_v = compute_shape_factor_galerkin(sz_cell_v, z - 0.5_rt);
210  }
211  const amrex::Real (&sz_ex)[depos_order + 1 ] = ((ex_type[zdir] == NODE) ? sz_node : sz_cell );
212  const amrex::Real (&sz_ey)[depos_order + 1 ] = ((ey_type[zdir] == NODE) ? sz_node : sz_cell );
213  const amrex::Real (&sz_ez)[depos_order + 1 - galerkin_interpolation] = ((ez_type[zdir] == NODE) ? sz_node_v : sz_cell_v);
214  const amrex::Real (&sz_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[zdir] == NODE) ? sz_node_v : sz_cell_v);
215  const amrex::Real (&sz_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[zdir] == NODE) ? sz_node_v : sz_cell_v);
216  const amrex::Real (&sz_bz)[depos_order + 1 ] = ((bz_type[zdir] == NODE) ? sz_node : sz_cell );
217  int const l_ex = ((ex_type[zdir] == NODE) ? l_node : l_cell );
218  int const l_ey = ((ey_type[zdir] == NODE) ? l_node : l_cell );
219  int const l_ez = ((ez_type[zdir] == NODE) ? l_node_v : l_cell_v);
220  int const l_bx = ((bx_type[zdir] == NODE) ? l_node_v : l_cell_v);
221  int const l_by = ((by_type[zdir] == NODE) ? l_node_v : l_cell_v);
222  int const l_bz = ((bz_type[zdir] == NODE) ? l_node : l_cell );
223 
224 
225  // Each field is gathered in a separate block of
226  // AMREX_SPACEDIM nested loops because the deposition
227  // order can differ for each component of each field
228  // when galerkin_interpolation is set to 1
229 
230 #if defined(WARPX_DIM_1D_Z)
231  // Gather field on particle Eyp from field on grid ey_arr
232  // Gather field on particle Exp from field on grid ex_arr
233  // Gather field on particle Bzp from field on grid bz_arr
234  for (int iz=0; iz<=depos_order; iz++){
235  Eyp += sz_ey[iz]*
236  ey_arr(lo.x+l_ey+iz, 0, 0, 0);
237  Exp += sz_ex[iz]*
238  ex_arr(lo.x+l_ex+iz, 0, 0, 0);
239  Bzp += sz_bz[iz]*
240  bz_arr(lo.x+l_bz+iz, 0, 0, 0);
241  }
242 
243  // Gather field on particle Byp from field on grid by_arr
244  // Gather field on particle Ezp from field on grid ez_arr
245  // Gather field on particle Bxp from field on grid bx_arr
246  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
247  Ezp += sz_ez[iz]*
248  ez_arr(lo.x+l_ez+iz, 0, 0, 0);
249  Bxp += sz_bx[iz]*
250  bx_arr(lo.x+l_bx+iz, 0, 0, 0);
251  Byp += sz_by[iz]*
252  by_arr(lo.x+l_by+iz, 0, 0, 0);
253  }
254 
255 #elif defined(WARPX_DIM_XZ)
256  // Gather field on particle Eyp from field on grid ey_arr
257  for (int iz=0; iz<=depos_order; iz++){
258  for (int ix=0; ix<=depos_order; ix++){
259  Eyp += sx_ey[ix]*sz_ey[iz]*
260  ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 0);
261  }
262  }
263  // Gather field on particle Exp from field on grid ex_arr
264  // Gather field on particle Bzp from field on grid bz_arr
265  for (int iz=0; iz<=depos_order; iz++){
266  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
267  Exp += sx_ex[ix]*sz_ex[iz]*
268  ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 0);
269  Bzp += sx_bz[ix]*sz_bz[iz]*
270  bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 0);
271  }
272  }
273  // Gather field on particle Ezp from field on grid ez_arr
274  // Gather field on particle Bxp from field on grid bx_arr
275  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
276  for (int ix=0; ix<=depos_order; ix++){
277  Ezp += sx_ez[ix]*sz_ez[iz]*
278  ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 0);
279  Bxp += sx_bx[ix]*sz_bx[iz]*
280  bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 0);
281  }
282  }
283  // Gather field on particle Byp from field on grid by_arr
284  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
285  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
286  Byp += sx_by[ix]*sz_by[iz]*
287  by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 0);
288  }
289  }
290 
291 #elif defined(WARPX_DIM_RZ)
292 
293  amrex::ParticleReal Erp = 0.;
294  amrex::ParticleReal Ethetap = 0.;
295  amrex::ParticleReal Brp = 0.;
296  amrex::ParticleReal Bthetap = 0.;
297 
298  // Gather field on particle Ethetap from field on grid ey_arr
299  for (int iz=0; iz<=depos_order; iz++){
300  for (int ix=0; ix<=depos_order; ix++){
301  Ethetap += sx_ey[ix]*sz_ey[iz]*
302  ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 0);
303  }
304  }
305  // Gather field on particle Erp from field on grid ex_arr
306  // Gather field on particle Bzp from field on grid bz_arr
307  for (int iz=0; iz<=depos_order; iz++){
308  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
309  Erp += sx_ex[ix]*sz_ex[iz]*
310  ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 0);
311  Bzp += sx_bz[ix]*sz_bz[iz]*
312  bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 0);
313  }
314  }
315  // Gather field on particle Ezp from field on grid ez_arr
316  // Gather field on particle Brp from field on grid bx_arr
317  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
318  for (int ix=0; ix<=depos_order; ix++){
319  Ezp += sx_ez[ix]*sz_ez[iz]*
320  ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 0);
321  Brp += sx_bx[ix]*sz_bx[iz]*
322  bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 0);
323  }
324  }
325  // Gather field on particle Bthetap from field on grid by_arr
326  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
327  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
328  Bthetap += sx_by[ix]*sz_by[iz]*
329  by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 0);
330  }
331  }
332 
333  amrex::Real costheta;
334  amrex::Real sintheta;
335  if (rp > 0.) {
336  costheta = xp/rp;
337  sintheta = yp/rp;
338  } else {
339  costheta = 1.;
340  sintheta = 0.;
341  }
342  const Complex xy0 = Complex{costheta, -sintheta};
343  Complex xy = xy0;
344 
345  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
346 
347  // Gather field on particle Ethetap from field on grid ey_arr
348  for (int iz=0; iz<=depos_order; iz++){
349  for (int ix=0; ix<=depos_order; ix++){
350  const amrex::Real dEy = (+ ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode-1)*xy.real()
351  - ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode)*xy.imag());
352  Ethetap += sx_ey[ix]*sz_ey[iz]*dEy;
353  }
354  }
355  // Gather field on particle Erp from field on grid ex_arr
356  // Gather field on particle Bzp from field on grid bz_arr
357  for (int iz=0; iz<=depos_order; iz++){
358  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
359  const amrex::Real dEx = (+ ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode-1)*xy.real()
360  - ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode)*xy.imag());
361  Erp += sx_ex[ix]*sz_ex[iz]*dEx;
362  const amrex::Real dBz = (+ bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode-1)*xy.real()
363  - bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode)*xy.imag());
364  Bzp += sx_bz[ix]*sz_bz[iz]*dBz;
365  }
366  }
367  // Gather field on particle Ezp from field on grid ez_arr
368  // Gather field on particle Brp from field on grid bx_arr
369  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
370  for (int ix=0; ix<=depos_order; ix++){
371  const amrex::Real dEz = (+ ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode-1)*xy.real()
372  - ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode)*xy.imag());
373  Ezp += sx_ez[ix]*sz_ez[iz]*dEz;
374  const amrex::Real dBx = (+ bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode-1)*xy.real()
375  - bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode)*xy.imag());
376  Brp += sx_bx[ix]*sz_bx[iz]*dBx;
377  }
378  }
379  // Gather field on particle Bthetap from field on grid by_arr
380  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
381  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
382  const amrex::Real dBy = (+ by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode-1)*xy.real()
383  - by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode)*xy.imag());
384  Bthetap += sx_by[ix]*sz_by[iz]*dBy;
385  }
386  }
387  xy = xy*xy0;
388  }
389 
390  // Convert Erp and Ethetap to Ex and Ey
391  Exp += costheta*Erp - sintheta*Ethetap;
392  Eyp += costheta*Ethetap + sintheta*Erp;
393  Bxp += costheta*Brp - sintheta*Bthetap;
394  Byp += costheta*Bthetap + sintheta*Brp;
395 
396 #else // defined(WARPX_DIM_3D)
397  // Gather field on particle Exp from field on grid ex_arr
398  for (int iz=0; iz<=depos_order; iz++){
399  for (int iy=0; iy<=depos_order; iy++){
400  for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
401  Exp += sx_ex[ix]*sy_ex[iy]*sz_ex[iz]*
402  ex_arr(lo.x+j_ex+ix, lo.y+k_ex+iy, lo.z+l_ex+iz);
403  }
404  }
405  }
406  // Gather field on particle Eyp from field on grid ey_arr
407  for (int iz=0; iz<=depos_order; iz++){
408  for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
409  for (int ix=0; ix<=depos_order; ix++){
410  Eyp += sx_ey[ix]*sy_ey[iy]*sz_ey[iz]*
411  ey_arr(lo.x+j_ey+ix, lo.y+k_ey+iy, lo.z+l_ey+iz);
412  }
413  }
414  }
415  // Gather field on particle Ezp from field on grid ez_arr
416  for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
417  for (int iy=0; iy<=depos_order; iy++){
418  for (int ix=0; ix<=depos_order; ix++){
419  Ezp += sx_ez[ix]*sy_ez[iy]*sz_ez[iz]*
420  ez_arr(lo.x+j_ez+ix, lo.y+k_ez+iy, lo.z+l_ez+iz);
421  }
422  }
423  }
424  // Gather field on particle Bzp from field on grid bz_arr
425  for (int iz=0; iz<=depos_order; iz++){
426  for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
427  for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
428  Bzp += sx_bz[ix]*sy_bz[iy]*sz_bz[iz]*
429  bz_arr(lo.x+j_bz+ix, lo.y+k_bz+iy, lo.z+l_bz+iz);
430  }
431  }
432  }
433  // Gather field on particle Byp from field on grid by_arr
434  for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
435  for (int iy=0; iy<=depos_order; iy++){
436  for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
437  Byp += sx_by[ix]*sy_by[iy]*sz_by[iz]*
438  by_arr(lo.x+j_by+ix, lo.y+k_by+iy, lo.z+l_by+iz);
439  }
440  }
441  }
442  // Gather field on particle Bxp from field on grid bx_arr
443  for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
444  for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
445  for (int ix=0; ix<=depos_order; ix++){
446  Bxp += sx_bx[ix]*sy_bx[iy]*sz_bx[iz]*
447  bx_arr(lo.x+j_bx+ix, lo.y+k_bx+iy, lo.z+l_bx+iz);
448  }
449  }
450  }
451 #endif
452 }
453 
472 template <int depos_order>
475  [[maybe_unused]] const amrex::ParticleReal xp_n,
476  [[maybe_unused]] const amrex::ParticleReal yp_n,
477  const amrex::ParticleReal zp_n,
478  [[maybe_unused]] const amrex::ParticleReal xp_nph,
479  [[maybe_unused]] const amrex::ParticleReal yp_nph,
480  const amrex::ParticleReal zp_nph,
481  amrex::ParticleReal& Exp,
482  amrex::ParticleReal& Eyp,
483  amrex::ParticleReal& Ezp,
484  amrex::ParticleReal& Bxp,
485  amrex::ParticleReal& Byp,
486  amrex::ParticleReal& Bzp,
487  amrex::Array4<amrex::Real const> const& Ex_arr,
488  amrex::Array4<amrex::Real const> const& Ey_arr,
489  amrex::Array4<amrex::Real const> const& Ez_arr,
490  amrex::Array4<amrex::Real const> const& Bx_arr,
491  amrex::Array4<amrex::Real const> const& By_arr,
492  amrex::Array4<amrex::Real const> const& Bz_arr,
493  [[maybe_unused]] const amrex::IndexType Ex_type,
494  [[maybe_unused]] const amrex::IndexType Ey_type,
495  [[maybe_unused]] const amrex::IndexType Ez_type,
496  [[maybe_unused]] const amrex::IndexType Bx_type,
497  [[maybe_unused]] const amrex::IndexType By_type,
498  [[maybe_unused]] const amrex::IndexType Bz_type,
500  const amrex::GpuArray<amrex::Real, 3>& xyzmin,
501  const amrex::Dim3& lo,
502  const int n_rz_azimuthal_modes)
503 {
504  using namespace amrex;
505 #if !defined(WARPX_DIM_RZ)
506  ignore_unused(n_rz_azimuthal_modes);
507 #endif
508 
509 #if !defined(WARPX_DIM_1D_Z)
510  Real const dxi = 1.0_rt / dx[0];
511 #endif
512 #if !defined(WARPX_DIM_1D_Z)
513  Real const xmin = xyzmin[0];
514 #endif
515 #if defined(WARPX_DIM_3D)
516  Real const dyi = 1.0_rt / dx[1];
517  Real const ymin = xyzmin[1];
518 #endif
519  Real const dzi = 1.0_rt / dx[2];
520  Real const zmin = xyzmin[2];
521 
522 #if !defined(WARPX_DIM_1D_Z)
523  Real constexpr one_third = 1.0_rt / 3.0_rt;
524  Real constexpr one_sixth = 1.0_rt / 6.0_rt;
525 #endif
526 
527 #if !defined(WARPX_DIM_1D_Z)
528  const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
529 #endif
530 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
531  const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
532 #endif
533  const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
534 
535  // computes current and old position in grid units
536 #if defined(WARPX_DIM_RZ)
537  amrex::Real const xp_new = xp_np1;
538  amrex::Real const yp_new = yp_np1;
539  amrex::Real const xp_mid = xp_nph;
540  amrex::Real const yp_mid = yp_nph;
541  amrex::Real const xp_old = xp_n;
542  amrex::Real const yp_old = yp_n;
543  amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
544  amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
545  amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
546  amrex::Real costheta_mid, sintheta_mid;
547  if (rp_mid > 0._rt) {
548  costheta_mid = xp_mid/rp_mid;
549  sintheta_mid = yp_mid/rp_mid;
550  } else {
551  costheta_mid = 1._rt;
552  sintheta_mid = 0._rt;
553  }
554  const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
555  // Keep these double to avoid bug in single precision
556  double const x_new = (rp_new - xmin)*dxi;
557  double const x_old = (rp_old - xmin)*dxi;
558 #else
559 #if !defined(WARPX_DIM_1D_Z)
560  // Keep these double to avoid bug in single precision
561  double const x_new = (xp_np1 - xmin)*dxi;
562  double const x_old = (xp_n - xmin)*dxi;
563 #endif
564 #endif
565 #if defined(WARPX_DIM_3D)
566  // Keep these double to avoid bug in single precision
567  double const y_new = (yp_np1 - ymin)*dyi;
568  double const y_old = (yp_n - ymin)*dyi;
569 #endif
570  // Keep these double to avoid bug in single precision
571  double const z_new = (zp_np1 - zmin)*dzi;
572  double const z_old = (zp_n - zmin)*dzi;
573 
574  // Shape factor arrays
575  // Note that there are extra values above and below
576  // to possibly hold the factor for the old particle
577  // which can be at a different grid location.
578  // Keep these double to avoid bug in single precision
579 #if !defined(WARPX_DIM_1D_Z)
580  double sx_E_new[depos_order + 3] = {0.};
581  double sx_E_old[depos_order + 3] = {0.};
582 #endif
583 #if defined(WARPX_DIM_3D)
584  // Keep these double to avoid bug in single precision
585  double sy_E_new[depos_order + 3] = {0.};
586  double sy_E_old[depos_order + 3] = {0.};
587 #endif
588  // Keep these double to avoid bug in single precision
589  double sz_E_new[depos_order + 3] = {0.};
590  double sz_E_old[depos_order + 3] = {0.};
591 
592 #if defined(WARPX_DIM_3D)
593  double sx_B_new[depos_order + 3] = {0.};
594  double sx_B_old[depos_order + 3] = {0.};
595  double sy_B_new[depos_order + 3] = {0.};
596  double sy_B_old[depos_order + 3] = {0.};
597  double sz_B_new[depos_order + 3] = {0.};
598  double sz_B_old[depos_order + 3] = {0.};
599 #endif
600 
601 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
602  // Special shape functions are needed for By which is cell
603  // centered in both x and z. One lower order shape function is used.
604  double sx_By_new[depos_order + 2] = {0.};
605  double sx_By_old[depos_order + 2] = {0.};
606  double sz_By_new[depos_order + 2] = {0.};
607  double sz_By_old[depos_order + 2] = {0.};
608 #endif
609 
610  // --- Compute shape factors
611  // Compute shape factors for position as they are now and at old positions
612  // [ijk]_new: leftmost grid point that the particle touches
613  const Compute_shape_factor< depos_order > compute_shape_factor;
614  const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor;
615 
616 #if !defined(WARPX_DIM_1D_Z)
617  const int i_E_new = compute_shape_factor(sx_E_new+1, x_new);
618  const int i_E_old = compute_shifted_shape_factor(sx_E_old, x_old, i_E_new);
619 #endif
620 #if defined(WARPX_DIM_3D)
621  const int j_E_new = compute_shape_factor(sy_E_new+1, y_new);
622  const int j_E_old = compute_shifted_shape_factor(sy_E_old, y_old, j_E_new);
623 #endif
624  const int k_E_new = compute_shape_factor(sz_E_new+1, z_new);
625  const int k_E_old = compute_shifted_shape_factor(sz_E_old, z_old, k_E_new);
626 
627 #if defined(WARPX_DIM_3D)
628  const int i_B_new = compute_shape_factor(sx_B_new+1, x_new - 0.5_rt);
629  const int i_B_old = compute_shifted_shape_factor(sx_B_old, x_old - 0.5_rt, i_B_new);
630  const int j_B_new = compute_shape_factor(sy_B_new+1, y_new - 0.5_rt);
631  const int j_B_old = compute_shifted_shape_factor(sy_B_old, y_old - 0.5_rt, j_B_new);
632  const int k_B_new = compute_shape_factor(sz_B_new+1, z_new - 0.5_rt);
633  const int k_B_old = compute_shifted_shape_factor(sz_B_old, z_old - 0.5_rt, k_B_new);
634 #endif
635 
636  // computes min/max positions of current contributions
637 #if !defined(WARPX_DIM_1D_Z)
638  int dil_E = 1, diu_E = 1;
639  if (i_E_old < i_E_new) { dil_E = 0; }
640  if (i_E_old > i_E_new) { diu_E = 0; }
641 #endif
642 #if defined(WARPX_DIM_3D)
643  int djl_E = 1, dju_E = 1;
644  if (j_E_old < j_E_new) { djl_E = 0; }
645  if (j_E_old > j_E_new) { dju_E = 0; }
646 #endif
647  int dkl_E = 1, dku_E = 1;
648  if (k_E_old < k_E_new) { dkl_E = 0; }
649  if (k_E_old > k_E_new) { dku_E = 0; }
650 
651 #if defined(WARPX_DIM_3D)
652  int dil_B = 1, diu_B = 1;
653  if (i_B_old < i_B_new) { dil_B = 0; }
654  if (i_B_old > i_B_new) { diu_B = 0; }
655  int djl_B = 1, dju_B = 1;
656  if (j_B_old < j_B_new) { djl_B = 0; }
657  if (j_B_old > j_B_new) { dju_B = 0; }
658  int dkl_B = 1, dku_B = 1;
659  if (k_B_old < k_B_new) { dkl_B = 0; }
660  if (k_B_old > k_B_new) { dku_B = 0; }
661 #endif
662 
663 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
664  const Compute_shape_factor< depos_order-1 > compute_shape_factor_By;
665  const Compute_shifted_shape_factor< depos_order-1 > compute_shifted_shape_factor_By;
666  const int i_By_new = compute_shape_factor_By(sx_By_new+1, x_new - 0.5_rt);
667  const int i_By_old = compute_shifted_shape_factor_By(sx_By_old, x_old - 0.5_rt, i_By_new);
668  const int k_By_new = compute_shape_factor_By(sz_By_new+1, z_new - 0.5_rt);
669  const int k_By_old = compute_shifted_shape_factor_By(sz_By_old, z_old - 0.5_rt, k_By_new);
670  int dil_By = 1, diu_By = 1;
671  if (i_By_old < i_By_new) { dil_By = 0; }
672  if (i_By_old > i_By_new) { diu_By = 0; }
673  int dkl_By = 1, dku_By = 1;
674  if (k_By_old < k_By_new) { dkl_By = 0; }
675  if (k_By_old > k_By_new) { dku_By = 0; }
676 #endif
677 
678 #if defined(WARPX_DIM_3D)
679 
680  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
681  for (int j=djl_E; j<=depos_order+2-dju_E; j++) {
682  const amrex::Real sdzjk = one_third*(sy_E_new[j]*sz_E_new[k] + sy_E_old[j]*sz_E_old[k])
683  +one_sixth*(sy_E_new[j]*sz_E_old[k] + sy_E_old[j]*sz_E_new[k]);
684  amrex::Real sdxi = 0._rt;
685  for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
686  sdxi += (sx_E_old[i] - sx_E_new[i]);
687  auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
688  Exp += Ex_arr(lo.x+i_E_new-1+i, lo.y+j_E_new-1+j, lo.z+k_E_new-1+k)*sdxiov*sdzjk;
689  }
690  }
691  }
692  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
693  for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
694  const amrex::Real sdyik = one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
695  +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]);
696  amrex::Real sdyj = 0._rt;
697  for (int j=djl_E; j<=depos_order+1-dju_E; j++) {
698  sdyj += (sy_E_old[j] - sy_E_new[j]);
699  auto sdyjov = static_cast<amrex::Real>((y_new - y_old) == 0. ? 1. : sdyj/(y_new - y_old));
700  Eyp += Ey_arr(lo.x+i_E_new-1+i, lo.y+j_E_new-1+j, lo.z+k_E_new-1+k)*sdyjov*sdyik;
701  }
702  }
703  }
704  for (int j=djl_E; j<=depos_order+2-dju_E; j++) {
705  for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
706  const amrex::Real sdzij = one_third*(sx_E_new[i]*sy_E_new[j] + sx_E_old[i]*sy_E_old[j])
707  +one_sixth*(sx_E_new[i]*sy_E_old[j] + sx_E_old[i]*sy_E_new[j]);
708  amrex::Real sdzk = 0._rt;
709  for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
710  sdzk += (sz_E_old[k] - sz_E_new[k]);
711  auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
712  Ezp += Ez_arr(lo.x+i_E_new-1+i, lo.y+j_E_new-1+j, lo.z+k_E_new-1+k)*sdzkov*sdzij;
713  }
714  }
715  }
716  for (int k=dkl_B; k<=depos_order+2-dku_B; k++) {
717  for (int j=djl_B; j<=depos_order+2-dju_B; j++) {
718  const amrex::Real sdzjk = one_third*(sy_B_new[j]*sz_B_new[k] + sy_B_old[j]*sz_B_old[k])
719  +one_sixth*(sy_B_new[j]*sz_B_old[k] + sy_B_old[j]*sz_B_new[k]);
720  amrex::Real sdxi = 0._rt;
721  for (int i=dil_B; i<=depos_order+1-diu_B; i++) {
722  sdxi += (sx_B_old[i] - sx_B_new[i]);
723  auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
724  Bxp += Bx_arr(lo.x+i_B_new-1+i, lo.y+j_B_new-1+j, lo.z+k_B_new-1+k)*sdxiov*sdzjk;
725  }
726  }
727  }
728  for (int k=dkl_B; k<=depos_order+2-dku_B; k++) {
729  for (int i=dil_B; i<=depos_order+2-diu_B; i++) {
730  const amrex::Real sdyik = one_third*(sx_B_new[i]*sz_B_new[k] + sx_B_old[i]*sz_B_old[k])
731  +one_sixth*(sx_B_new[i]*sz_B_old[k] + sx_B_old[i]*sz_B_new[k]);
732  amrex::Real sdyj = 0._rt;
733  for (int j=djl_B; j<=depos_order+1-dju_B; j++) {
734  sdyj += (sy_B_old[j] - sy_B_new[j]);
735  auto sdyjov = static_cast<amrex::Real>((y_new - y_old) == 0. ? 1. : sdyj/(y_new - y_old));
736  Byp += By_arr(lo.x+i_B_new-1+i, lo.y+j_B_new-1+j, lo.z+k_B_new-1+k)*sdyjov*sdyik;
737  }
738  }
739  }
740  for (int j=djl_B; j<=depos_order+2-dju_B; j++) {
741  for (int i=dil_B; i<=depos_order+2-diu_B; i++) {
742  const amrex::Real sdzij = one_third*(sx_B_new[i]*sy_B_new[j] + sx_B_old[i]*sy_B_old[j])
743  +one_sixth*(sx_B_new[i]*sy_B_old[j] + sx_B_old[i]*sy_B_new[j]);
744  amrex::Real sdzk = 0._rt;
745  for (int k=dkl_B; k<=depos_order+1-dku_B; k++) {
746  sdzk += (sz_B_old[k] - sz_B_new[k]);
747  auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
748  Bzp += Bz_arr(lo.x+i_B_new-1+i, lo.y+j_B_new-1+j, lo.z+k_E_new-1+k)*sdzkov*sdzij;
749  }
750  }
751  }
752 
753 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
754 
755  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
756  const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
757  amrex::Real sdxi = 0._rt;
758  for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
759  sdxi += (sx_E_old[i] - sx_E_new[i]);
760  auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
761  Exp += Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
762  Bzp += Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
763  }
764  }
765  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
766  for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
767  Real const sdyj = (
768  one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
769  +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]));
770  Eyp += Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdyj;
771  }
772  }
773  for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
774  const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]);
775  amrex::Real sdzk = 0._rt;
776  for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
777  sdzk += (sz_E_old[k] - sz_E_new[k]);
778  auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
779  Ezp += Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
780  Bxp += Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
781  }
782  }
783  for (int k=dkl_By; k<=depos_order+1-dku_By; k++) {
784  for (int i=dil_By; i<=depos_order+1-diu_By; i++) {
785  Real const sdyj = (
786  one_third*(sx_By_new[i]*sz_By_new[k] + sx_By_old[i]*sz_By_old[k])
787  +one_sixth*(sx_By_new[i]*sz_By_old[k] + sx_By_old[i]*sz_By_new[k]));
788  Byp += By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 0)*sdyj;
789  }
790  }
791 
792 #ifdef WARPX_DIM_RZ
793  Complex xy_mid = xy_mid0;
794 
795  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
796 
797  // Gather field on particle Exp from field on grid ex_arr
798  // Gather field on particle Bzp from field on grid bz_arr
799  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
800  const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
801  amrex::Real sdxi = 0._rt;
802  for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
803  sdxi += (sx_E_old[i] - sx_E_new[i]);
804  auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
805  const amrex::Real dEx = (+ Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.real()
806  - Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
807  const amrex::Real dBz = (+ Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.real()
808  - Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
809  Exp += dEx*sdxiov*sdzk;
810  Bzp += dBz*sdxiov*sdzk;
811  }
812  }
813  // Gather field on particle Eyp from field on grid ey_arr
814  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
815  for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
816  Real const sdyj = (
817  one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
818  +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]));
819  const amrex::Real dEy = (+ Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.real()
820  - Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
821  Eyp += dEy*sdyj;
822  }
823  }
824  // Gather field on particle Ezp from field on grid ez_arr
825  // Gather field on particle Bxp from field on grid bx_arr
826  for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
827  const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]);
828  amrex::Real sdzk = 0._rt;
829  for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
830  sdzk += (sz_E_old[k] - sz_E_new[k]);
831  auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
832  const amrex::Real dEz = (+ Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.real()
833  - Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
834  const amrex::Real dBx = (+ Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.real()
835  - Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
836  Ezp += dEz*sdzkov*sdxi;
837  Bxp += dBx*sdzkov*sdxi;
838  }
839  }
840  // Gather field on particle Byp from field on grid by_arr
841  for (int k=dkl_By; k<=depos_order+1-dku_By; k++) {
842  for (int i=dil_By; i<=depos_order+1-diu_By; i++) {
843  Real const sdyj = (
844  one_third*(sx_By_new[i]*sz_By_new[k] + sx_By_old[i]*sz_By_old[k])
845  +one_sixth*(sx_By_new[i]*sz_By_old[k] + sx_By_old[i]*sz_By_new[k]));
846  const amrex::Real dBy = (+ By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 2*imode-1)*xy_mid.real()
847  - By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 2*imode)*xy_mid.imag());
848  Byp += dBy*sdyj;
849  }
850  }
851  xy_mid = xy_mid*xy_mid0;
852  }
853 
854  // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey
855  const amrex::Real Exp_save = Exp;
856  Exp = costheta_mid*Exp - sintheta_mid*Eyp;
857  Eyp = costheta_mid*Eyp + sintheta_mid*Exp_save;
858  const amrex::Real Bxp_save = Bxp;
859  Bxp = costheta_mid*Bxp - sintheta_mid*Byp;
860  Byp = costheta_mid*Byp + sintheta_mid*Bxp_save;
861 
862 #endif
863 
864 #elif defined(WARPX_DIM_1D_Z)
865 
866  for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
867  amrex::Real const sdzk = 0.5_rt*(sz_E_old[k] + sz_E_new[k]);
868  Exp += Ex_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
869  Eyp += Ey_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
870  Bzp += Bz_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
871  }
872  amrex::Real sdzk = 0._rt;
873  for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
874  sdzk += (sz_E_old[k] - sz_E_new[k]);
875  auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
876  Bxp += Bx_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
877  Byp += By_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
878  Ezp += Ez_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
879  }
880 #endif
881 }
882 
902 template <int depos_order>
905  [[maybe_unused]] const amrex::ParticleReal xp_n,
906  [[maybe_unused]] const amrex::ParticleReal yp_n,
907  const amrex::ParticleReal zp_n,
908  [[maybe_unused]] const amrex::ParticleReal xp_nph,
909  [[maybe_unused]] const amrex::ParticleReal yp_nph,
910  const amrex::ParticleReal zp_nph,
911  amrex::ParticleReal& Exp,
912  amrex::ParticleReal& Eyp,
913  amrex::ParticleReal& Ezp,
914  amrex::ParticleReal& Bxp,
915  amrex::ParticleReal& Byp,
916  amrex::ParticleReal& Bzp,
917  amrex::Array4<amrex::Real const> const& Ex_arr,
918  amrex::Array4<amrex::Real const> const& Ey_arr,
919  amrex::Array4<amrex::Real const> const& Ez_arr,
920  amrex::Array4<amrex::Real const> const& Bx_arr,
921  amrex::Array4<amrex::Real const> const& By_arr,
922  amrex::Array4<amrex::Real const> const& Bz_arr,
923  [[maybe_unused]] const amrex::IndexType Ex_type,
924  [[maybe_unused]] const amrex::IndexType Ey_type,
925  [[maybe_unused]] const amrex::IndexType Ez_type,
926  [[maybe_unused]] const amrex::IndexType Bx_type,
927  [[maybe_unused]] const amrex::IndexType By_type,
928  [[maybe_unused]] const amrex::IndexType Bz_type,
930  const amrex::GpuArray<amrex::Real, 3>& xyzmin,
931  const amrex::Dim3& lo,
932  const int n_rz_azimuthal_modes)
933 {
934  using namespace amrex;
935 #if !defined(WARPX_DIM_RZ)
936  ignore_unused(n_rz_azimuthal_modes);
937 #endif
938 
939 #if !defined(WARPX_DIM_1D_Z)
940  Real const dxi = 1.0_rt / dx[0];
941 #endif
942 #if !defined(WARPX_DIM_1D_Z)
943  Real const xmin = xyzmin[0];
944 #endif
945 #if defined(WARPX_DIM_3D)
946  Real const dyi = 1.0_rt / dx[1];
947  Real const ymin = xyzmin[1];
948 #endif
949  Real const dzi = 1.0_rt / dx[2];
950  Real const zmin = xyzmin[2];
951 
952 #if !defined(WARPX_DIM_1D_Z)
953  const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
954 #endif
955 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
956  const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
957 #endif
958  const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
959 
960 #if !defined(WARPX_DIM_1D_Z)
961  Real constexpr one_third = 1.0_rt / 3.0_rt;
962  Real constexpr one_sixth = 1.0_rt / 6.0_rt;
963 #endif
964 
965  // computes current and old position in grid units
966 #if defined(WARPX_DIM_RZ)
967  amrex::Real const xp_new = xp_np1;
968  amrex::Real const yp_new = yp_np1;
969  amrex::Real const xp_mid = xp_nph;
970  amrex::Real const yp_mid = yp_nph;
971  amrex::Real const xp_old = xp_n;
972  amrex::Real const yp_old = yp_n;
973  amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
974  amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
975  amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
976  amrex::Real costheta_mid, sintheta_mid;
977  if (rp_mid > 0._rt) {
978  costheta_mid = xp_mid/rp_mid;
979  sintheta_mid = yp_mid/rp_mid;
980  } else {
981  costheta_mid = 1._rt;
982  sintheta_mid = 0._rt;
983  }
984  const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
985  // Keep these double to avoid bug in single precision
986  double const x_new = (rp_new - xmin)*dxi;
987  double const x_old = (rp_old - xmin)*dxi;
988  double const x_bar = (rp_mid - xmin)*dxi;
989 #elif !defined(WARPX_DIM_1D_Z)
990  // Keep these double to avoid bug in single precision
991  double const x_new = (xp_np1 - xmin)*dxi;
992  double const x_old = (xp_n - xmin)*dxi;
993  double const x_bar = (xp_nph - xmin)*dxi;
994 #endif
995 #if defined(WARPX_DIM_3D)
996  // Keep these double to avoid bug in single precision
997  double const y_new = (yp_np1 - ymin)*dyi;
998  double const y_old = (yp_n - ymin)*dyi;
999  double const y_bar = (yp_nph - ymin)*dyi;
1000 #endif
1001  // Keep these double to avoid bug in single precision
1002  double const z_new = (zp_np1 - zmin)*dzi;
1003  double const z_old = (zp_n - zmin)*dzi;
1004  double const z_bar = (zp_nph - zmin)*dzi;
1005 
1006  // 1) Determine the number of segments.
1007  // 2) Loop over segments and gather electric field.
1008  // 3) Gather magnetic field.
1009 
1010  // cell crossings are defined at cell edges if depos_order is odd
1011  // cell crossings are defined at cell centers if depos_order is even
1012 
1013  int num_segments = 1;
1014  double shift = 0.0;
1015  if ( (depos_order % 2) == 0 ) { shift = 0.5; }
1016 
1017 #if defined(WARPX_DIM_3D)
1018 
1019  // compute cell crossings in X-direction
1020  const auto i_old = static_cast<int>(x_old-shift);
1021  const auto i_new = static_cast<int>(x_new-shift);
1022  const int cell_crossings_x = std::abs(i_new-i_old);
1023  num_segments += cell_crossings_x;
1024 
1025  // compute cell crossings in Y-direction
1026  const auto j_old = static_cast<int>(y_old-shift);
1027  const auto j_new = static_cast<int>(y_new-shift);
1028  const int cell_crossings_y = std::abs(j_new-j_old);
1029  num_segments += cell_crossings_y;
1030 
1031  // compute cell crossings in Z-direction
1032  const auto k_old = static_cast<int>(z_old-shift);
1033  const auto k_new = static_cast<int>(z_new-shift);
1034  const int cell_crossings_z = std::abs(k_new-k_old);
1035  num_segments += cell_crossings_z;
1036 
1037  // need to assert that the number of cell crossings in each direction
1038  // is within the range permitted by the number of guard cells
1039  // e.g., if (num_segments > 7) ...
1040 
1041  // compute total change in particle position and the initial cell
1042  // locations in each direction used to find the position at cell crossings.
1043  const double dxp = x_new - x_old;
1044  const double dyp = y_new - y_old;
1045  const double dzp = z_new - z_old;
1046  const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1047  const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
1048  const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1049  double Xcell = 0., Ycell = 0., Zcell = 0.;
1050  if (num_segments > 1) {
1051  Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1052  Ycell = static_cast<double>(j_old) + shift + 0.5*(1.-dirY_sign);
1053  Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1054  }
1055 
1056  // loop over the number of segments and deposit
1057  const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1058  const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1059  double dxp_seg, dyp_seg, dzp_seg;
1060  double x0_new, y0_new, z0_new;
1061  double x0_old = x_old;
1062  double y0_old = y_old;
1063  double z0_old = z_old;
1064 
1065  for (int ns=0; ns<num_segments; ns++) {
1066 
1067  if (ns == num_segments-1) { // final segment
1068 
1069  x0_new = x_new;
1070  y0_new = y_new;
1071  z0_new = z_new;
1072  dxp_seg = x0_new - x0_old;
1073  dyp_seg = y0_new - y0_old;
1074  dzp_seg = z0_new - z0_old;
1075 
1076  }
1077  else {
1078 
1079  x0_new = Xcell + dirX_sign;
1080  y0_new = Ycell + dirY_sign;
1081  z0_new = Zcell + dirZ_sign;
1082  dxp_seg = x0_new - x0_old;
1083  dyp_seg = y0_new - y0_old;
1084  dzp_seg = z0_new - z0_old;
1085 
1086  if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
1087  && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
1088  Xcell = x0_new;
1089  dyp_seg = dyp/dxp*dxp_seg;
1090  dzp_seg = dzp/dxp*dxp_seg;
1091  y0_new = y0_old + dyp_seg;
1092  z0_new = z0_old + dzp_seg;
1093  }
1094  else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
1095  Ycell = y0_new;
1096  dxp_seg = dxp/dyp*dyp_seg;
1097  dzp_seg = dzp/dyp*dyp_seg;
1098  x0_new = x0_old + dxp_seg;
1099  z0_new = z0_old + dzp_seg;
1100  }
1101  else {
1102  Zcell = z0_new;
1103  dxp_seg = dxp/dzp*dzp_seg;
1104  dyp_seg = dyp/dzp*dzp_seg;
1105  x0_new = x0_old + dxp_seg;
1106  y0_new = y0_old + dyp_seg;
1107  }
1108 
1109  }
1110 
1111  // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, dyp, or dzp)
1112  const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1113  const auto seg_factor_y = static_cast<double>(dyp == 0. ? 1. : dyp_seg/dyp);
1114  const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1115 
1116  // compute cell-based weights using the average segment position
1117  double sx_cell[depos_order] = {0.};
1118  double sy_cell[depos_order] = {0.};
1119  double sz_cell[depos_order] = {0.};
1120  double const x0_bar = (x0_new + x0_old)/2.0;
1121  double const y0_bar = (y0_new + y0_old)/2.0;
1122  double const z0_bar = (z0_new + z0_old)/2.0;
1123  const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1124  const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
1125  const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1126 
1127  if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1128  const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1129  double sx_old_cell[depos_order] = {0.};
1130  double sx_new_cell[depos_order] = {0.};
1131  double sy_old_cell[depos_order] = {0.};
1132  double sy_new_cell[depos_order] = {0.};
1133  double sz_old_cell[depos_order] = {0.};
1134  double sz_new_cell[depos_order] = {0.};
1135  const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1136  const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
1137  const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1138  ignore_unused(i0_cell_2, j0_cell_2, k0_cell_2);
1139  for (int m=0; m<depos_order; m++) {
1140  sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1141  sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
1142  sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1143  }
1144  }
1145 
1146  // compute node-based weights using the old and new segment positions
1147  double sx_old_node[depos_order+1] = {0.};
1148  double sx_new_node[depos_order+1] = {0.};
1149  double sy_old_node[depos_order+1] = {0.};
1150  double sy_new_node[depos_order+1] = {0.};
1151  double sz_old_node[depos_order+1] = {0.};
1152  double sz_new_node[depos_order+1] = {0.};
1153  const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1154  const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
1155  const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1156 
1157  // gather Ex for this segment
1158  amrex::Real weight;
1159  for (int i=0; i<=depos_order-1; i++) {
1160  for (int j=0; j<=depos_order; j++) {
1161  for (int k=0; k<=depos_order; k++) {
1162  weight = sx_cell[i]*( sy_old_node[j]*sz_old_node[k]*one_third
1163  + sy_old_node[j]*sz_new_node[k]*one_sixth
1164  + sy_new_node[j]*sz_old_node[k]*one_sixth
1165  + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
1166  Exp += Ex_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k)*weight;
1167  }
1168  }
1169  }
1170 
1171  // gather Ey for this segment
1172  for (int i=0; i<=depos_order; i++) {
1173  for (int j=0; j<=depos_order-1; j++) {
1174  for (int k=0; k<=depos_order; k++) {
1175  weight = sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
1176  + sx_old_node[i]*sz_new_node[k]*one_sixth
1177  + sx_new_node[i]*sz_old_node[k]*one_sixth
1178  + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1179  Eyp += Ey_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k)*weight;
1180  }
1181  }
1182  }
1183 
1184  // gather Ez for this segment
1185  for (int i=0; i<=depos_order; i++) {
1186  for (int j=0; j<=depos_order; j++) {
1187  for (int k=0; k<=depos_order-1; k++) {
1188  weight = sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
1189  + sx_old_node[i]*sy_new_node[j]*one_sixth
1190  + sx_new_node[i]*sy_old_node[j]*one_sixth
1191  + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
1192  Ezp += Ez_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k)*weight;
1193  }
1194  }
1195  }
1196 
1197  // update old segment values
1198  if (ns < num_segments-1) {
1199  x0_old = x0_new;
1200  y0_old = y0_new;
1201  z0_old = z0_new;
1202  }
1203 
1204  } // end loop over segments
1205 
1206  // gather magnetic field
1207  const int depos_order_B = 1; // second template parameter?
1208  const Compute_shape_factor< depos_order_B > compute_shape_factor_B;
1209  double sz_bar_node[depos_order_B+1] = {0.};
1210  double sz_bar_cell[depos_order_B+1] = {0.};
1211  const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar);
1212  const int k_bar_cell = compute_shape_factor_B(sz_bar_cell, z_bar-0.5);
1213  double sy_bar_node[depos_order_B+1] = {0.};
1214  double sy_bar_cell[depos_order_B+1] = {0.};
1215  const int j_bar_node = compute_shape_factor_B(sy_bar_node, y_bar);
1216  const int j_bar_cell = compute_shape_factor_B(sy_bar_cell, y_bar-0.5);
1217  double sx_bar_node[depos_order_B+1] = {0.};
1218  double sx_bar_cell[depos_order_B+1] = {0.};
1219  const int i_bar_node = compute_shape_factor_B(sx_bar_node, x_bar);
1220  const int i_bar_cell = compute_shape_factor_B(sx_bar_cell, x_bar-0.5);
1221 
1222  amrex::Real weight;
1223  for (int i=0; i<=depos_order_B; i++) {
1224  for (int j=0; j<=depos_order_B; j++) {
1225  for (int k=0; k<=depos_order_B; k++) {
1226  weight = static_cast<amrex::Real>(sx_bar_node[i]*sy_bar_cell[j]*sz_bar_cell[k]);
1227  Bxp += Bx_arr(lo.x+i_bar_node+i, lo.y+j_bar_cell+j, lo.z+k_bar_cell+k)*weight;
1228  //
1229  weight = static_cast<amrex::Real>(sx_bar_cell[i]*sy_bar_node[j]*sz_bar_cell[k]);
1230  Byp += By_arr(lo.x+i_bar_cell+i, lo.y+j_bar_node+j, lo.z+k_bar_cell+k)*weight;
1231  //
1232  weight = static_cast<amrex::Real>(sx_bar_cell[i]*sy_bar_cell[j]*sz_bar_node[k]);
1233  Bzp += Bz_arr(lo.x+i_bar_cell+i, lo.y+j_bar_cell+j, lo.z+k_bar_node+k)*weight;
1234  }
1235  }
1236  }
1237 
1238 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1239 
1240  // compute cell crossings in X-direction
1241  const auto i_old = static_cast<int>(x_old-shift);
1242  const auto i_new = static_cast<int>(x_new-shift);
1243  const int cell_crossings_x = std::abs(i_new-i_old);
1244  num_segments += cell_crossings_x;
1245 
1246  // compute cell crossings in Z-direction
1247  const auto k_old = static_cast<int>(z_old-shift);
1248  const auto k_new = static_cast<int>(z_new-shift);
1249  const int cell_crossings_z = std::abs(k_new-k_old);
1250  num_segments += cell_crossings_z;
1251 
1252  // need to assert that the number of cell crossings in each direction
1253  // is within the range permitted by the number of guard cells
1254  // e.g., if (num_segments > 5) ...
1255 
1256  // compute total change in particle position and the initial cell
1257  // locations in each direction used to find the position at cell crossings.
1258  const double dxp = x_new - x_old;
1259  const double dzp = z_new - z_old;
1260  const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1261  const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1262  double Xcell = 0., Zcell = 0.;
1263  if (num_segments > 1) {
1264  Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1265  Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1266  }
1267 
1268  // loop over the number of segments and deposit
1269  const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1270  const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1271  double dxp_seg, dzp_seg;
1272  double x0_new, z0_new;
1273  double x0_old = x_old;
1274  double z0_old = z_old;
1275 
1276  for (int ns=0; ns<num_segments; ns++) {
1277 
1278  if (ns == num_segments-1) { // final segment
1279 
1280  x0_new = x_new;
1281  z0_new = z_new;
1282  dxp_seg = x0_new - x0_old;
1283  dzp_seg = z0_new - z0_old;
1284 
1285  }
1286  else {
1287 
1288  x0_new = Xcell + dirX_sign;
1289  z0_new = Zcell + dirZ_sign;
1290  dxp_seg = x0_new - x0_old;
1291  dzp_seg = z0_new - z0_old;
1292 
1293  if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1294  Xcell = x0_new;
1295  dzp_seg = dzp/dxp*dxp_seg;
1296  z0_new = z0_old + dzp_seg;
1297  }
1298  else {
1299  Zcell = z0_new;
1300  dxp_seg = dxp/dzp*dzp_seg;
1301  x0_new = x0_old + dxp_seg;
1302  }
1303 
1304  }
1305 
1306  // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1307  const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1308  const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1309 
1310  // compute cell-based weights using the average segment position
1311  double sx_cell[depos_order] = {0.};
1312  double sz_cell[depos_order] = {0.};
1313  double const x0_bar = (x0_new + x0_old)/2.0;
1314  double const z0_bar = (z0_new + z0_old)/2.0;
1315  const int i0_cell = compute_shape_factor_cell(sx_cell, x0_bar-0.5);
1316  const int k0_cell = compute_shape_factor_cell(sz_cell, z0_bar-0.5);
1317 
1318  if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1319  const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1320  double sx_old_cell[depos_order] = {0.};
1321  double sx_new_cell[depos_order] = {0.};
1322  double sz_old_cell[depos_order] = {0.};
1323  double sz_new_cell[depos_order] = {0.};
1324  const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1325  const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1326  ignore_unused(i0_cell_2, k0_cell_2);
1327  for (int m=0; m<depos_order; m++) {
1328  sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1329  sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1330  }
1331  }
1332 
1333  // compute node-based weights using the old and new segment positions
1334  double sx_old_node[depos_order+1] = {0.};
1335  double sx_new_node[depos_order+1] = {0.};
1336  double sz_old_node[depos_order+1] = {0.};
1337  double sz_new_node[depos_order+1] = {0.};
1338  const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1339  const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1340 
1341  // gather Ex for this segment
1342  amrex::Real weight;
1343  for (int i=0; i<=depos_order-1; i++) {
1344  for (int k=0; k<=depos_order; k++) {
1345  weight = sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1346  Exp += Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 0)*weight;
1347 #if defined(WARPX_DIM_RZ)
1348  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1349  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1350  const auto dEx = (+ Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode-1)*xy_mid.real()
1351  - Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode)*xy_mid.imag());
1352  Exp += weight*dEx;
1353  xy_mid = xy_mid*xy_mid0;
1354  }
1355 #endif
1356  }
1357  }
1358 
1359  // gather out-of-plane Ey for this segment
1360  const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1361  for (int i=0; i<=depos_order; i++) {
1362  for (int k=0; k<=depos_order; k++) {
1363  weight = ( sx_old_node[i]*sz_old_node[k]*one_third
1364  + sx_old_node[i]*sz_new_node[k]*one_sixth
1365  + sx_new_node[i]*sz_old_node[k]*one_sixth
1366  + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1367  Eyp += Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 0)*weight;
1368 #if defined(WARPX_DIM_RZ)
1369  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1370  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1371  const auto dEy = (+ Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode-1)*xy_mid.real()
1372  - Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode)*xy_mid.imag());
1373  Eyp += weight*dEy;
1374  xy_mid = xy_mid*xy_mid0;
1375  }
1376 #endif
1377  }
1378  }
1379 
1380  // gather Ez for this segment
1381  for (int i=0; i<=depos_order; i++) {
1382  for (int k=0; k<=depos_order-1; k++) {
1383  weight = sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1384  Ezp += Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 0)*weight;
1385 #if defined(WARPX_DIM_RZ)
1386  Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
1387  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1388  const auto dEz = (+ Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode-1)*xy_mid.real()
1389  - Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode)*xy_mid.imag());
1390  Ezp += weight*dEz;
1391  xy_mid = xy_mid*xy_mid0;
1392  }
1393 #endif
1394  }
1395  }
1396 
1397  // update old segment values
1398  if (ns < num_segments-1) {
1399  x0_old = x0_new;
1400  z0_old = z0_new;
1401  }
1402 
1403  }
1404 
1405  // gather magnetic field and out-of-plane electric field
1406  const int depos_order_B = 1; // second template parameter?
1407  const Compute_shape_factor< depos_order_B > compute_shape_factor_B;
1408  double sz_bar_node[depos_order_B+1] = {0.};
1409  double sz_bar_cell[depos_order_B+1] = {0.};
1410  const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar);
1411  const int k_bar_cell = compute_shape_factor_B(sz_bar_cell, z_bar-0.5);
1412  double sx_bar_node[depos_order_B+1] = {0.};
1413  double sx_bar_cell[depos_order_B+1] = {0.};
1414  const int i_bar_node = compute_shape_factor_B(sx_bar_node, x_bar);
1415  const int i_bar_cell = compute_shape_factor_B(sx_bar_cell, x_bar-0.5);
1416 
1417  for (int i=0; i<=depos_order_B; i++) {
1418  for (int k=0; k<=depos_order_B; k++) {
1419  const auto weight_Bz = static_cast<amrex::Real>(sx_bar_cell[i]*sz_bar_node[k]);
1420  Bzp += Bz_arr(lo.x+i_bar_cell+i, lo.y+k_bar_node+k, 0, 0)*weight_Bz;
1421  //
1422  const auto weight_Bx = static_cast<amrex::Real>(sx_bar_node[i]*sz_bar_cell[k]);
1423  Bxp += Bx_arr(lo.x+i_bar_node+i, lo.y+k_bar_cell+k, 0, 0)*weight_Bx;
1424  //
1425  const auto weight_By = static_cast<amrex::Real>(sx_bar_cell[i]*sz_bar_cell[k]);
1426  Byp += By_arr(lo.x+i_bar_cell+i, lo.y+k_bar_cell+k, 0, 0)*weight_By;
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  const auto dBx = (+ Bx_arr(lo.x+i_bar_node+i, lo.y+k_bar_cell+k, 0, 2*imode-1)*xy_mid.real()
1431  - Bx_arr(lo.x+i_bar_node+i, lo.y+k_bar_cell+k, 0, 2*imode)*xy_mid.imag());
1432  const auto dBy = (+ By_arr(lo.x+i_bar_cell+i, lo.y+k_bar_cell+k, 0, 2*imode-1)*xy_mid.real()
1433  - By_arr(lo.x+i_bar_cell+i, lo.y+k_bar_cell+k, 0, 2*imode)*xy_mid.imag());
1434  const auto dBz = (+ Bz_arr(lo.x+i_bar_cell+i, lo.y+k_bar_node+k, 0, 2*imode-1)*xy_mid.real()
1435  - Bz_arr(lo.x+i_bar_cell+i, lo.y+k_bar_node+k, 0, 2*imode)*xy_mid.imag());
1436  Bxp += weight_Bx*dBx;
1437  Byp += weight_By*dBy;
1438  Bzp += weight_Bz*dBz;
1439  xy_mid = xy_mid*xy_mid0;
1440  }
1441 #endif
1442  }
1443  }
1444 
1445 #ifdef WARPX_DIM_RZ
1446 
1447  // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey
1448  const amrex::Real Exp_save = Exp;
1449  Exp = costheta_mid*Exp - sintheta_mid*Eyp;
1450  Eyp = costheta_mid*Eyp + sintheta_mid*Exp_save;
1451  const amrex::Real Bxp_save = Bxp;
1452  Bxp = costheta_mid*Bxp - sintheta_mid*Byp;
1453  Byp = costheta_mid*Byp + sintheta_mid*Bxp_save;
1454 
1455 #endif
1456 
1457 #elif defined(WARPX_DIM_1D_Z)
1458 
1459  // compute cell crossings in Z-direction
1460  const auto k_old = static_cast<int>(z_old-shift);
1461  const auto k_new = static_cast<int>(z_new-shift);
1462  const int cell_crossings_z = std::abs(k_new-k_old);
1463  num_segments += cell_crossings_z;
1464 
1465  // need to assert that the number of cell crossings in each direction
1466  // is within the range permitted by the number of guard cells
1467  // e.g., if (num_segments > 3) ...
1468 
1469  // compute dzp and the initial cell location used to find the cell crossings.
1470  double const dzp = z_new - z_old;
1471  const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1472  double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1473 
1474  // loop over the number of segments and deposit
1475  const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1476  const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1477  double dzp_seg;
1478  double z0_new;
1479  double z0_old = z_old;
1480 
1481  for (int ns=0; ns<num_segments; ns++) {
1482 
1483  if (ns == num_segments-1) { // final segment
1484  z0_new = z_new;
1485  dzp_seg = z0_new - z0_old;
1486  }
1487  else {
1488  Zcell = Zcell + dirZ_sign;
1489  z0_new = Zcell;
1490  dzp_seg = z0_new - z0_old;
1491  }
1492 
1493  // compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1494  const auto seg_factor = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1495 
1496  // compute cell-based weights using the average segment position
1497  double sz_cell[depos_order] = {0.};
1498  double const z0_bar = (z0_new + z0_old)/2.0;
1499  const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1500 
1501  if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1502  const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1503  double sz_old_cell[depos_order] = {0.};
1504  double sz_new_cell[depos_order] = {0.};
1505  const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1506  ignore_unused(k0_cell_2);
1507  for (int m=0; m<depos_order; m++) {
1508  sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1509  }
1510  }
1511 
1512  // compute node-based weights using the old and new segment positions
1513  double sz_old_node[depos_order+1] = {0.};
1514  double sz_new_node[depos_order+1] = {0.};
1515  const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1516 
1517  // gather out-of-plane Ex and Ey for this segment
1518  for (int k=0; k<=depos_order; k++) {
1519  auto weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1520  Exp += Ex_arr(lo.x+k0_node+k, 0, 0)*weight;
1521  Eyp += Ey_arr(lo.x+k0_node+k, 0, 0)*weight;
1522  }
1523 
1524  // gather Ez for this segment
1525  for (int k=0; k<=depos_order-1; k++) {
1526  auto weight = sz_cell[k]*seg_factor;
1527  Ezp += Ez_arr(lo.x+k0_cell+k, 0, 0)*weight;
1528  }
1529 
1530  // update old segment values
1531  if (ns < num_segments-1) {
1532  z0_old = z0_new;
1533  }
1534 
1535  }
1536 
1537  // gather magnetic field
1538  const int depos_order_B = 1; // second template parameter?
1539  const Compute_shape_factor< depos_order_B > compute_shape_factor_B;
1540  double sz_bar_node[depos_order_B+1] = {0.};
1541  double sz_bar_cell[depos_order_B+1] = {0.};
1542  const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar);
1543  const int k_bar_cell = compute_shape_factor_B(sz_bar_cell, z_bar-0.5_rt);
1544 
1545  amrex::Real weight;
1546  for (int k=0; k<=depos_order_B; k++) {
1547  weight = static_cast<amrex::Real>(sz_bar_node[k]);
1548  Bzp += Bz_arr(lo.x+k_bar_node+k, 0, 0)*weight;
1549  //
1550  weight = static_cast<amrex::Real>(sz_bar_cell[k]);
1551  Bxp += Bx_arr(lo.x+k_bar_cell+k, 0, 0)*weight;
1552  Byp += By_arr(lo.x+k_bar_cell+k, 0, 0)*weight;
1553  }
1554 
1555 #endif
1556 }
1557 
1575 template <int depos_order, int lower_in_v>
1577  const GetExternalEBField& getExternalEB,
1578  amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp,
1579  amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp,
1580  amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp,
1581  amrex::FArrayBox const * const exfab,
1582  amrex::FArrayBox const * const eyfab,
1583  amrex::FArrayBox const * const ezfab,
1584  amrex::FArrayBox const * const bxfab,
1585  amrex::FArrayBox const * const byfab,
1586  amrex::FArrayBox const * const bzfab,
1587  const long np_to_gather,
1588  const std::array<amrex::Real, 3>& dx,
1589  const std::array<amrex::Real, 3> xyzmin,
1590  const amrex::Dim3 lo,
1591  const int n_rz_azimuthal_modes)
1592 {
1593 
1594  const amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
1595  const amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
1596 
1597  amrex::Array4<const amrex::Real> const& ex_arr = exfab->array();
1598  amrex::Array4<const amrex::Real> const& ey_arr = eyfab->array();
1599  amrex::Array4<const amrex::Real> const& ez_arr = ezfab->array();
1600  amrex::Array4<const amrex::Real> const& bx_arr = bxfab->array();
1601  amrex::Array4<const amrex::Real> const& by_arr = byfab->array();
1602  amrex::Array4<const amrex::Real> const& bz_arr = bzfab->array();
1603 
1604  amrex::IndexType const ex_type = exfab->box().ixType();
1605  amrex::IndexType const ey_type = eyfab->box().ixType();
1606  amrex::IndexType const ez_type = ezfab->box().ixType();
1607  amrex::IndexType const bx_type = bxfab->box().ixType();
1608  amrex::IndexType const by_type = byfab->box().ixType();
1609  amrex::IndexType const bz_type = bzfab->box().ixType();
1610 
1611  // Loop over particles and gather fields from
1612  // {e,b}{x,y,z}_arr to {E,B}{xyz}p.
1614  np_to_gather,
1615  [=] AMREX_GPU_DEVICE (long ip) {
1616 
1617  amrex::ParticleReal xp, yp, zp;
1618  getPosition(ip, xp, yp, zp);
1619  getExternalEB(ip, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip]);
1620 
1621  doGatherShapeN<depos_order, lower_in_v>(
1622  xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip],
1623  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1624  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1625  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1626  }
1627  );
1628 }
1629 
1648 void doGatherShapeN (const amrex::ParticleReal xp,
1649  const amrex::ParticleReal yp,
1650  const amrex::ParticleReal zp,
1651  amrex::ParticleReal& Exp,
1652  amrex::ParticleReal& Eyp,
1653  amrex::ParticleReal& Ezp,
1654  amrex::ParticleReal& Bxp,
1655  amrex::ParticleReal& Byp,
1656  amrex::ParticleReal& Bzp,
1657  amrex::Array4<amrex::Real const> const& ex_arr,
1658  amrex::Array4<amrex::Real const> const& ey_arr,
1659  amrex::Array4<amrex::Real const> const& ez_arr,
1660  amrex::Array4<amrex::Real const> const& bx_arr,
1661  amrex::Array4<amrex::Real const> const& by_arr,
1662  amrex::Array4<amrex::Real const> const& bz_arr,
1663  const amrex::IndexType ex_type,
1664  const amrex::IndexType ey_type,
1665  const amrex::IndexType ez_type,
1666  const amrex::IndexType bx_type,
1667  const amrex::IndexType by_type,
1668  const amrex::IndexType bz_type,
1669  const amrex::GpuArray<amrex::Real, 3>& dx_arr,
1670  const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
1671  const amrex::Dim3& lo,
1672  const int n_rz_azimuthal_modes,
1673  const int nox,
1674  const bool galerkin_interpolation)
1675 {
1676  if (galerkin_interpolation) {
1677  if (nox == 1) {
1678  doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1679  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1680  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1681  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1682  } else if (nox == 2) {
1683  doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1684  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1685  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1686  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1687  } else if (nox == 3) {
1688  doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1689  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1690  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1691  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1692  } else if (nox == 4) {
1693  doGatherShapeN<4,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1694  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1695  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1696  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1697  }
1698  } else {
1699  if (nox == 1) {
1700  doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1701  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1702  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1703  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1704  } else if (nox == 2) {
1705  doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1706  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1707  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1708  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1709  } else if (nox == 3) {
1710  doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1711  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1712  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1713  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1714  } else if (nox == 4) {
1715  doGatherShapeN<4,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1716  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1717  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1718  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1719  }
1720  }
1721 }
1722 
1723 
1745  const amrex::ParticleReal xp_n,
1746  const amrex::ParticleReal yp_n,
1747  const amrex::ParticleReal zp_n,
1748  const amrex::ParticleReal xp_nph,
1749  const amrex::ParticleReal yp_nph,
1750  const amrex::ParticleReal zp_nph,
1751  amrex::ParticleReal& Exp,
1752  amrex::ParticleReal& Eyp,
1753  amrex::ParticleReal& Ezp,
1754  amrex::ParticleReal& Bxp,
1755  amrex::ParticleReal& Byp,
1756  amrex::ParticleReal& Bzp,
1757  amrex::Array4<amrex::Real const> const& ex_arr,
1758  amrex::Array4<amrex::Real const> const& ey_arr,
1759  amrex::Array4<amrex::Real const> const& ez_arr,
1760  amrex::Array4<amrex::Real const> const& bx_arr,
1761  amrex::Array4<amrex::Real const> const& by_arr,
1762  amrex::Array4<amrex::Real const> const& bz_arr,
1763  const amrex::IndexType ex_type,
1764  const amrex::IndexType ey_type,
1765  const amrex::IndexType ez_type,
1766  const amrex::IndexType bx_type,
1767  const amrex::IndexType by_type,
1768  const amrex::IndexType bz_type,
1769  const amrex::GpuArray<amrex::Real, 3>& dx_arr,
1770  const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
1771  const amrex::Dim3& lo,
1772  const int n_rz_azimuthal_modes,
1773  const int nox,
1774  const int depos_type )
1775 {
1776  if (depos_type==0) { // CurrentDepositionAlgo::Esirkepov
1777  if (nox == 1) {
1778  doGatherShapeNEsirkepovStencilImplicit<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1779  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1780  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1781  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1782  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1783  } else if (nox == 2) {
1784  doGatherShapeNEsirkepovStencilImplicit<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1785  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1786  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1787  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1788  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1789  } else if (nox == 3) {
1790  doGatherShapeNEsirkepovStencilImplicit<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1791  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1792  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1793  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1794  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1795  } else if (nox == 4) {
1796  doGatherShapeNEsirkepovStencilImplicit<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1797  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1798  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1799  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1800  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1801  }
1802  }
1803  else if (depos_type==3) { // CurrentDepositionAlgo::Villasenor
1804  if (nox == 1) {
1805  doGatherPicnicShapeN<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1806  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1807  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1808  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1809  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1810  } else if (nox == 2) {
1811  doGatherPicnicShapeN<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1812  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1813  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1814  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1815  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1816  } else if (nox == 3) {
1817  doGatherPicnicShapeN<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1818  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1819  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1820  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1821  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1822  } else if (nox == 4) {
1823  doGatherPicnicShapeN<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1824  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1825  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1826  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1827  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1828  }
1829  }
1830  else if (depos_type==1) { // CurrentDepositionAlgo::Direct
1831  if (nox == 1) {
1832  doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1833  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1834  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1835  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1836  } else if (nox == 2) {
1837  doGatherShapeN<2,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1838  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1839  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1840  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1841  } else if (nox == 3) {
1842  doGatherShapeN<3,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1843  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1844  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1845  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1846  } else if (nox == 4) {
1847  doGatherShapeN<4,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1848  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1849  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1850  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1851  }
1852  }
1853 }
1854 
1855 #endif // WARPX_FIELDGATHER_H_
#define AMREX_FORCE_INLINE
#define AMREX_GPU_DEVICE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeN(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &ex_arr, amrex::Array4< amrex::Real const > const &ey_arr, amrex::Array4< amrex::Real const > const &ez_arr, amrex::Array4< amrex::Real const > const &bx_arr, amrex::Array4< amrex::Real const > const &by_arr, amrex::Array4< amrex::Real const > const &bz_arr, const amrex::IndexType ex_type, const amrex::IndexType ey_type, const amrex::IndexType ez_type, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx, const amrex::GpuArray< amrex::Real, 3 > &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Field gather for a single particle.
Definition: FieldGather.H:38
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeNImplicit(const amrex::ParticleReal xp_n, const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, const amrex::ParticleReal xp_nph, const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &ex_arr, amrex::Array4< amrex::Real const > const &ey_arr, amrex::Array4< amrex::Real const > const &ez_arr, amrex::Array4< amrex::Real const > const &bx_arr, amrex::Array4< amrex::Real const > const &by_arr, amrex::Array4< amrex::Real const > const &bz_arr, const amrex::IndexType ex_type, const amrex::IndexType ey_type, const amrex::IndexType ez_type, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx_arr, const amrex::GpuArray< amrex::Real, 3 > &xyzmin_arr, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes, const int nox, const int depos_type)
Field gather for a single particle.
Definition: FieldGather.H:1744
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherPicnicShapeN([[maybe_unused]] const amrex::ParticleReal xp_n, [[maybe_unused]] const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, [[maybe_unused]] const amrex::ParticleReal xp_nph, [[maybe_unused]] const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &Ex_arr, amrex::Array4< amrex::Real const > const &Ey_arr, amrex::Array4< amrex::Real const > const &Ez_arr, amrex::Array4< amrex::Real const > const &Bx_arr, amrex::Array4< amrex::Real const > const &By_arr, amrex::Array4< amrex::Real const > const &Bz_arr, [[maybe_unused]] const amrex::IndexType Ex_type, [[maybe_unused]] const amrex::IndexType Ey_type, [[maybe_unused]] const amrex::IndexType Ez_type, [[maybe_unused]] const amrex::IndexType Bx_type, [[maybe_unused]] const amrex::IndexType By_type, [[maybe_unused]] const amrex::IndexType Bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx, const amrex::GpuArray< amrex::Real, 3 > &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Energy conserving field gather for thread thread_num for the implicit scheme This uses the same stenc...
Definition: FieldGather.H:904
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeNEsirkepovStencilImplicit([[maybe_unused]] const amrex::ParticleReal xp_n, [[maybe_unused]] const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, [[maybe_unused]] const amrex::ParticleReal xp_nph, [[maybe_unused]] const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &Ex_arr, amrex::Array4< amrex::Real const > const &Ey_arr, amrex::Array4< amrex::Real const > const &Ez_arr, amrex::Array4< amrex::Real const > const &Bx_arr, amrex::Array4< amrex::Real const > const &By_arr, amrex::Array4< amrex::Real const > const &Bz_arr, [[maybe_unused]] const amrex::IndexType Ex_type, [[maybe_unused]] const amrex::IndexType Ey_type, [[maybe_unused]] const amrex::IndexType Ez_type, [[maybe_unused]] const amrex::IndexType Bx_type, [[maybe_unused]] const amrex::IndexType By_type, [[maybe_unused]] const amrex::IndexType Bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx, const amrex::GpuArray< amrex::Real, 3 > &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Energy conserving field gather for thread thread_num for the implicit scheme This uses the same stenc...
Definition: FieldGather.H:474
NODE
AMREX_FORCE_INLINE Array4< Real const > array() const noexcept
const Box & box() const noexcept
AMREX_GPU_HOST_DEVICE IndexType ixType() const noexcept
constexpr int iz
constexpr int iy
constexpr int ix
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box shift(const Box &b, int dir, int nzones) noexcept
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
i
Definition: check_interp_points_and_weights.py:174
tuple dx
lab frame
Definition: stencil.py:429
xmin
Definition: stencil.py:424
Definition: ShapeFactors.H:168
Definition: ShapeFactors.H:29
Definition: ShapeFactors.H:95
Functor class that assigns external field values (E and B) to particles.
Definition: GetExternalFields.H:25
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T real() const noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T imag() const noexcept