ARTEMIS
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) || defined(WARPX_DIM_RZ)
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 #ifdef WARPX_DIM_RZ
292 
293  amrex::Real costheta;
294  amrex::Real sintheta;
295  if (rp > 0.) {
296  costheta = xp/rp;
297  sintheta = yp/rp;
298  } else {
299  costheta = 1.;
300  sintheta = 0.;
301  }
302  const Complex xy0 = Complex{costheta, -sintheta};
303  Complex xy = xy0;
304 
305  for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
306 
307  // Gather field on particle Eyp from field on grid ey_arr
308  for (int iz=0; iz<=depos_order; iz++){
309  for (int ix=0; ix<=depos_order; ix++){
310  const amrex::Real dEy = (+ ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode-1)*xy.real()
311  - ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode)*xy.imag());
312  Eyp += sx_ey[ix]*sz_ey[iz]*dEy;
313  }
314  }
315  // Gather field on particle Exp from field on grid ex_arr
316  // Gather field on particle Bzp from field on grid bz_arr
317  for (int iz=0; iz<=depos_order; iz++){
318  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
319  const amrex::Real dEx = (+ ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode-1)*xy.real()
320  - ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode)*xy.imag());
321  Exp += sx_ex[ix]*sz_ex[iz]*dEx;
322  const amrex::Real dBz = (+ bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode-1)*xy.real()
323  - bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode)*xy.imag());
324  Bzp += sx_bz[ix]*sz_bz[iz]*dBz;
325  }
326  }
327  // Gather field on particle Ezp from field on grid ez_arr
328  // Gather field on particle Bxp from field on grid bx_arr
329  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
330  for (int ix=0; ix<=depos_order; ix++){
331  const amrex::Real dEz = (+ ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode-1)*xy.real()
332  - ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode)*xy.imag());
333  Ezp += sx_ez[ix]*sz_ez[iz]*dEz;
334  const amrex::Real dBx = (+ bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode-1)*xy.real()
335  - bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode)*xy.imag());
336  Bxp += sx_bx[ix]*sz_bx[iz]*dBx;
337  }
338  }
339  // Gather field on particle Byp from field on grid by_arr
340  for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
341  for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
342  const amrex::Real dBy = (+ by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode-1)*xy.real()
343  - by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode)*xy.imag());
344  Byp += sx_by[ix]*sz_by[iz]*dBy;
345  }
346  }
347  xy = xy*xy0;
348  }
349 
350  // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey
351  const amrex::Real Exp_save = Exp;
352  Exp = costheta*Exp - sintheta*Eyp;
353  Eyp = costheta*Eyp + sintheta*Exp_save;
354  const amrex::Real Bxp_save = Bxp;
355  Bxp = costheta*Bxp - sintheta*Byp;
356  Byp = costheta*Byp + sintheta*Bxp_save;
357 #endif
358 
359 #else // defined(WARPX_DIM_3D)
360  // Gather field on particle Exp from field on grid ex_arr
361  for (int iz=0; iz<=depos_order; iz++){
362  for (int iy=0; iy<=depos_order; iy++){
363  for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
364  Exp += sx_ex[ix]*sy_ex[iy]*sz_ex[iz]*
365  ex_arr(lo.x+j_ex+ix, lo.y+k_ex+iy, lo.z+l_ex+iz);
366  }
367  }
368  }
369  // Gather field on particle Eyp from field on grid ey_arr
370  for (int iz=0; iz<=depos_order; iz++){
371  for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
372  for (int ix=0; ix<=depos_order; ix++){
373  Eyp += sx_ey[ix]*sy_ey[iy]*sz_ey[iz]*
374  ey_arr(lo.x+j_ey+ix, lo.y+k_ey+iy, lo.z+l_ey+iz);
375  }
376  }
377  }
378  // Gather field on particle Ezp from field on grid ez_arr
379  for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
380  for (int iy=0; iy<=depos_order; iy++){
381  for (int ix=0; ix<=depos_order; ix++){
382  Ezp += sx_ez[ix]*sy_ez[iy]*sz_ez[iz]*
383  ez_arr(lo.x+j_ez+ix, lo.y+k_ez+iy, lo.z+l_ez+iz);
384  }
385  }
386  }
387  // Gather field on particle Bzp from field on grid bz_arr
388  for (int iz=0; iz<=depos_order; iz++){
389  for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
390  for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
391  Bzp += sx_bz[ix]*sy_bz[iy]*sz_bz[iz]*
392  bz_arr(lo.x+j_bz+ix, lo.y+k_bz+iy, lo.z+l_bz+iz);
393  }
394  }
395  }
396  // Gather field on particle Byp from field on grid by_arr
397  for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
398  for (int iy=0; iy<=depos_order; iy++){
399  for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
400  Byp += sx_by[ix]*sy_by[iy]*sz_by[iz]*
401  by_arr(lo.x+j_by+ix, lo.y+k_by+iy, lo.z+l_by+iz);
402  }
403  }
404  }
405  // Gather field on particle Bxp from field on grid bx_arr
406  for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
407  for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
408  for (int ix=0; ix<=depos_order; ix++){
409  Bxp += sx_bx[ix]*sy_bx[iy]*sz_bx[iz]*
410  bx_arr(lo.x+j_bx+ix, lo.y+k_bx+iy, lo.z+l_bx+iz);
411  }
412  }
413  }
414 #endif
415 }
416 
434 template <int depos_order, int lower_in_v>
435 void doGatherShapeN(const GetParticlePosition& getPosition,
436  const GetExternalEBField& getExternalEB,
437  amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp,
438  amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp,
439  amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp,
440  amrex::FArrayBox const * const exfab,
441  amrex::FArrayBox const * const eyfab,
442  amrex::FArrayBox const * const ezfab,
443  amrex::FArrayBox const * const bxfab,
444  amrex::FArrayBox const * const byfab,
445  amrex::FArrayBox const * const bzfab,
446  const long np_to_gather,
447  const std::array<amrex::Real, 3>& dx,
448  const std::array<amrex::Real, 3> xyzmin,
449  const amrex::Dim3 lo,
450  const int n_rz_azimuthal_modes)
451 {
452 
453  amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
454  amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
455 
456  amrex::Array4<const amrex::Real> const& ex_arr = exfab->array();
457  amrex::Array4<const amrex::Real> const& ey_arr = eyfab->array();
458  amrex::Array4<const amrex::Real> const& ez_arr = ezfab->array();
459  amrex::Array4<const amrex::Real> const& bx_arr = bxfab->array();
460  amrex::Array4<const amrex::Real> const& by_arr = byfab->array();
461  amrex::Array4<const amrex::Real> const& bz_arr = bzfab->array();
462 
463  amrex::IndexType const ex_type = exfab->box().ixType();
464  amrex::IndexType const ey_type = eyfab->box().ixType();
465  amrex::IndexType const ez_type = ezfab->box().ixType();
466  amrex::IndexType const bx_type = bxfab->box().ixType();
467  amrex::IndexType const by_type = byfab->box().ixType();
468  amrex::IndexType const bz_type = bzfab->box().ixType();
469 
470  // Loop over particles and gather fields from
471  // {e,b}{x,y,z}_arr to {E,B}{xyz}p.
472  amrex::ParallelFor(
473  np_to_gather,
474  [=] AMREX_GPU_DEVICE (long ip) {
475 
476  amrex::ParticleReal xp, yp, zp;
477  getPosition(ip, xp, yp, zp);
478  getExternalEB(ip, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip]);
479 
480  doGatherShapeN<depos_order, lower_in_v>(
481  xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip],
482  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
483  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
484  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
485  }
486  );
487 }
488 
506 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
507 void doGatherShapeN (const amrex::ParticleReal xp,
508  const amrex::ParticleReal yp,
509  const amrex::ParticleReal zp,
510  amrex::ParticleReal& Exp,
511  amrex::ParticleReal& Eyp,
512  amrex::ParticleReal& Ezp,
513  amrex::ParticleReal& Bxp,
514  amrex::ParticleReal& Byp,
515  amrex::ParticleReal& Bzp,
516  amrex::Array4<amrex::Real const> const& ex_arr,
517  amrex::Array4<amrex::Real const> const& ey_arr,
518  amrex::Array4<amrex::Real const> const& ez_arr,
519  amrex::Array4<amrex::Real const> const& bx_arr,
520  amrex::Array4<amrex::Real const> const& by_arr,
521  amrex::Array4<amrex::Real const> const& bz_arr,
522  const amrex::IndexType ex_type,
523  const amrex::IndexType ey_type,
524  const amrex::IndexType ez_type,
525  const amrex::IndexType bx_type,
526  const amrex::IndexType by_type,
527  const amrex::IndexType bz_type,
528  const amrex::GpuArray<amrex::Real, 3>& dx_arr,
529  const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
530  const amrex::Dim3& lo,
531  const int n_rz_azimuthal_modes,
532  const int nox,
533  const bool galerkin_interpolation)
534 {
535  if (galerkin_interpolation) {
536  if (nox == 1) {
537  doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
538  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
539  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
540  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
541  } else if (nox == 2) {
542  doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
543  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
544  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
545  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
546  } else if (nox == 3) {
547  doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
548  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
549  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
550  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
551  }
552  } else {
553  if (nox == 1) {
554  doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
555  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
556  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
557  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
558  } else if (nox == 2) {
559  doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
560  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
561  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
562  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
563  } else if (nox == 3) {
564  doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
565  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
566  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
567  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
568  }
569  }
570 }
571 
572 #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
WarpX_Complex.H
NODE
NODE
amrex::BaseFab::box
const Box & box() const noexcept
amrex::BaseFab::array
AMREX_FORCE_INLINE Array4< T const > array() const noexcept
amrex::BoxND::ixType
AMREX_GPU_HOST_DEVICE IndexTypeND< dim > ixType() const noexcept
Excitation_Flag_Generator.y
def y
Definition: Excitation_Flag_Generator.py:76
Excitation_Flag_Generator.z
def z
Definition: Excitation_Flag_Generator.py:77
amrex::ParallelFor
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
stencil.nox
int nox
Definition: stencil.py:442
stencil.dx
int dx
Definition: stencil.py:436
Compute_shape_factor
Definition: ShapeFactors.H:27
GetExternalEBField
Functor class that assigns external field values (E and B) to particles.
Definition: GetExternalFields.H:25
GetParticlePosition
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel.
Definition: GetAndSetPosition.H:53
amrex::GpuComplex::real
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr T real() const noexcept
amrex::GpuComplex::imag
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr T imag() const noexcept
Generated by doxygen 1.9.1