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 FIELDGATHER_H_
9 #define 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 (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,
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::GpuArray<amrex::Real, 3>& dx,
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)
67  amrex::ignore_unused(yp);
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>
473 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
474 void doGatherShapeNEsirkepovStencilImplicit (
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,
499  const amrex::GpuArray<amrex::Real, 3>& dx,
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  ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
529 #endif
530 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
531  ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
532 #endif
533  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  Compute_shape_factor< depos_order > compute_shape_factor;
614  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  Compute_shape_factor< depos_order-1 > compute_shape_factor_By;
665  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  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  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  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  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  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  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  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  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  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  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 
900 template <int depos_order, int lower_in_v>
901 void doGatherShapeN(const GetParticlePosition<PIdx>& getPosition,
902  const GetExternalEBField& getExternalEB,
903  amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp,
904  amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp,
905  amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp,
906  amrex::FArrayBox const * const exfab,
907  amrex::FArrayBox const * const eyfab,
908  amrex::FArrayBox const * const ezfab,
909  amrex::FArrayBox const * const bxfab,
910  amrex::FArrayBox const * const byfab,
911  amrex::FArrayBox const * const bzfab,
912  const long np_to_gather,
913  const std::array<amrex::Real, 3>& dx,
914  const std::array<amrex::Real, 3> xyzmin,
915  const amrex::Dim3 lo,
916  const int n_rz_azimuthal_modes)
917 {
918 
919  amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
920  amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
921 
922  amrex::Array4<const amrex::Real> const& ex_arr = exfab->array();
923  amrex::Array4<const amrex::Real> const& ey_arr = eyfab->array();
924  amrex::Array4<const amrex::Real> const& ez_arr = ezfab->array();
925  amrex::Array4<const amrex::Real> const& bx_arr = bxfab->array();
926  amrex::Array4<const amrex::Real> const& by_arr = byfab->array();
927  amrex::Array4<const amrex::Real> const& bz_arr = bzfab->array();
928 
929  amrex::IndexType const ex_type = exfab->box().ixType();
930  amrex::IndexType const ey_type = eyfab->box().ixType();
931  amrex::IndexType const ez_type = ezfab->box().ixType();
932  amrex::IndexType const bx_type = bxfab->box().ixType();
933  amrex::IndexType const by_type = byfab->box().ixType();
934  amrex::IndexType const bz_type = bzfab->box().ixType();
935 
936  // Loop over particles and gather fields from
937  // {e,b}{x,y,z}_arr to {E,B}{xyz}p.
938  amrex::ParallelFor(
939  np_to_gather,
940  [=] AMREX_GPU_DEVICE (long ip) {
941 
942  amrex::ParticleReal xp, yp, zp;
943  getPosition(ip, xp, yp, zp);
944  getExternalEB(ip, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip]);
945 
946  doGatherShapeN<depos_order, lower_in_v>(
947  xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip],
948  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
949  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
950  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
951  }
952  );
953 }
954 
972 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
973 void doGatherShapeN (const amrex::ParticleReal xp,
974  const amrex::ParticleReal yp,
975  const amrex::ParticleReal zp,
976  amrex::ParticleReal& Exp,
977  amrex::ParticleReal& Eyp,
978  amrex::ParticleReal& Ezp,
979  amrex::ParticleReal& Bxp,
980  amrex::ParticleReal& Byp,
981  amrex::ParticleReal& Bzp,
982  amrex::Array4<amrex::Real const> const& ex_arr,
983  amrex::Array4<amrex::Real const> const& ey_arr,
984  amrex::Array4<amrex::Real const> const& ez_arr,
985  amrex::Array4<amrex::Real const> const& bx_arr,
986  amrex::Array4<amrex::Real const> const& by_arr,
987  amrex::Array4<amrex::Real const> const& bz_arr,
988  const amrex::IndexType ex_type,
989  const amrex::IndexType ey_type,
990  const amrex::IndexType ez_type,
991  const amrex::IndexType bx_type,
992  const amrex::IndexType by_type,
993  const amrex::IndexType bz_type,
994  const amrex::GpuArray<amrex::Real, 3>& dx_arr,
995  const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
996  const amrex::Dim3& lo,
997  const int n_rz_azimuthal_modes,
998  const int nox,
999  const bool galerkin_interpolation)
1000 {
1001  if (galerkin_interpolation) {
1002  if (nox == 1) {
1003  doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1004  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1005  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1006  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1007  } else if (nox == 2) {
1008  doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1009  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1010  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1011  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1012  } else if (nox == 3) {
1013  doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1014  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1015  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1016  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1017  }
1018  } else {
1019  if (nox == 1) {
1020  doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1021  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1022  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1023  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1024  } else if (nox == 2) {
1025  doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1026  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1027  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1028  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1029  } else if (nox == 3) {
1030  doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1031  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1032  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1033  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1034  }
1035  }
1036 }
1037 
1038 
1057 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1058 void doGatherShapeNImplicit (
1059  const amrex::ParticleReal xp_n,
1060  const amrex::ParticleReal yp_n,
1061  const amrex::ParticleReal zp_n,
1062  const amrex::ParticleReal xp_nph,
1063  const amrex::ParticleReal yp_nph,
1064  const amrex::ParticleReal zp_nph,
1065  amrex::ParticleReal& Exp,
1066  amrex::ParticleReal& Eyp,
1067  amrex::ParticleReal& Ezp,
1068  amrex::ParticleReal& Bxp,
1069  amrex::ParticleReal& Byp,
1070  amrex::ParticleReal& Bzp,
1071  amrex::Array4<amrex::Real const> const& ex_arr,
1072  amrex::Array4<amrex::Real const> const& ey_arr,
1073  amrex::Array4<amrex::Real const> const& ez_arr,
1074  amrex::Array4<amrex::Real const> const& bx_arr,
1075  amrex::Array4<amrex::Real const> const& by_arr,
1076  amrex::Array4<amrex::Real const> const& bz_arr,
1077  const amrex::IndexType ex_type,
1078  const amrex::IndexType ey_type,
1079  const amrex::IndexType ez_type,
1080  const amrex::IndexType bx_type,
1081  const amrex::IndexType by_type,
1082  const amrex::IndexType bz_type,
1083  const amrex::GpuArray<amrex::Real, 3>& dx_arr,
1084  const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
1085  const amrex::Dim3& lo,
1086  const int n_rz_azimuthal_modes,
1087  const int nox,
1088  const bool galerkin_interpolation)
1089 {
1090  if (galerkin_interpolation) {
1091  if (nox == 1) {
1092  doGatherShapeNEsirkepovStencilImplicit<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1093  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1094  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1095  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1096  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1097  } else if (nox == 2) {
1098  doGatherShapeNEsirkepovStencilImplicit<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1099  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1100  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1101  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1102  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1103  } else if (nox == 3) {
1104  doGatherShapeNEsirkepovStencilImplicit<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1105  Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1106  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1107  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1108  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1109  }
1110  } else {
1111  if (nox == 1) {
1112  doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1113  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1114  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1115  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1116  } else if (nox == 2) {
1117  doGatherShapeN<2,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1118  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1119  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1120  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1121  } else if (nox == 3) {
1122  doGatherShapeN<3,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1123  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1124  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1125  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1126  }
1127  }
1128 }
1129 
1130 #endif // 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
doGatherShapeN
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
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::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 bool galerkin_interpolation)
Field gather for a single particle.
Definition: FieldGather.H:1058
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::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
WarpX_Complex.H
NODE
NODE
BaseFab< Real >::array
AMREX_FORCE_INLINE Array4< Real const > array() const noexcept
BaseFab< Real >::box
const Box & box() const noexcept
amrex::Box::ixType
AMREX_GPU_HOST_DEVICE IndexType ixType() const noexcept
amrex::ParallelFor
std::enable_if_t< std::is_integral< T >::value > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
check_interp_points_and_weights.i
i
Definition: check_interp_points_and_weights.py:174
stencil.dx
tuple dx
lab frame
Definition: stencil.py:429
stencil.xmin
xmin
Definition: stencil.py:424
Compute_shape_factor
Definition: ShapeFactors.H:29
Compute_shifted_shape_factor
Definition: ShapeFactors.H:84
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