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