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 
471 template <int depos_order, int lower_in_v>
472 void doGatherShapeN(const GetParticlePosition<PIdx>& getPosition,
473  const GetExternalEBField& getExternalEB,
474  amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp,
475  amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp,
476  amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp,
477  amrex::FArrayBox const * const exfab,
478  amrex::FArrayBox const * const eyfab,
479  amrex::FArrayBox const * const ezfab,
480  amrex::FArrayBox const * const bxfab,
481  amrex::FArrayBox const * const byfab,
482  amrex::FArrayBox const * const bzfab,
483  const long np_to_gather,
484  const std::array<amrex::Real, 3>& dx,
485  const std::array<amrex::Real, 3> xyzmin,
486  const amrex::Dim3 lo,
487  const int n_rz_azimuthal_modes)
488 {
489 
490  amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
491  amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
492 
493  amrex::Array4<const amrex::Real> const& ex_arr = exfab->array();
494  amrex::Array4<const amrex::Real> const& ey_arr = eyfab->array();
495  amrex::Array4<const amrex::Real> const& ez_arr = ezfab->array();
496  amrex::Array4<const amrex::Real> const& bx_arr = bxfab->array();
497  amrex::Array4<const amrex::Real> const& by_arr = byfab->array();
498  amrex::Array4<const amrex::Real> const& bz_arr = bzfab->array();
499 
500  amrex::IndexType const ex_type = exfab->box().ixType();
501  amrex::IndexType const ey_type = eyfab->box().ixType();
502  amrex::IndexType const ez_type = ezfab->box().ixType();
503  amrex::IndexType const bx_type = bxfab->box().ixType();
504  amrex::IndexType const by_type = byfab->box().ixType();
505  amrex::IndexType const bz_type = bzfab->box().ixType();
506 
507  // Loop over particles and gather fields from
508  // {e,b}{x,y,z}_arr to {E,B}{xyz}p.
509  amrex::ParallelFor(
510  np_to_gather,
511  [=] AMREX_GPU_DEVICE (long ip) {
512 
513  amrex::ParticleReal xp, yp, zp;
514  getPosition(ip, xp, yp, zp);
515  getExternalEB(ip, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip]);
516 
517  doGatherShapeN<depos_order, lower_in_v>(
518  xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip],
519  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
520  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
521  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
522  }
523  );
524 }
525 
543 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
544 void doGatherShapeN (const amrex::ParticleReal xp,
545  const amrex::ParticleReal yp,
546  const amrex::ParticleReal zp,
547  amrex::ParticleReal& Exp,
548  amrex::ParticleReal& Eyp,
549  amrex::ParticleReal& Ezp,
550  amrex::ParticleReal& Bxp,
551  amrex::ParticleReal& Byp,
552  amrex::ParticleReal& Bzp,
553  amrex::Array4<amrex::Real const> const& ex_arr,
554  amrex::Array4<amrex::Real const> const& ey_arr,
555  amrex::Array4<amrex::Real const> const& ez_arr,
556  amrex::Array4<amrex::Real const> const& bx_arr,
557  amrex::Array4<amrex::Real const> const& by_arr,
558  amrex::Array4<amrex::Real const> const& bz_arr,
559  const amrex::IndexType ex_type,
560  const amrex::IndexType ey_type,
561  const amrex::IndexType ez_type,
562  const amrex::IndexType bx_type,
563  const amrex::IndexType by_type,
564  const amrex::IndexType bz_type,
565  const amrex::GpuArray<amrex::Real, 3>& dx_arr,
566  const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
567  const amrex::Dim3& lo,
568  const int n_rz_azimuthal_modes,
569  const int nox,
570  const bool galerkin_interpolation)
571 {
572  if (galerkin_interpolation) {
573  if (nox == 1) {
574  doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
575  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
576  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
577  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
578  } else if (nox == 2) {
579  doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
580  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
581  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
582  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
583  } else if (nox == 3) {
584  doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
585  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
586  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
587  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
588  }
589  } else {
590  if (nox == 1) {
591  doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
592  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
593  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
594  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
595  } else if (nox == 2) {
596  doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
597  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
598  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
599  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
600  } else if (nox == 3) {
601  doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
602  ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
603  ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
604  dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
605  }
606  }
607 }
608 
609 #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
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 &...)
stencil.dx
tuple dx
lab frame
Definition: stencil.py:429
stencil.xmin
xmin
Definition: stencil.py:424
Compute_shape_factor
Definition: ShapeFactors.H:29
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