WarpX
Loading...
Searching...
No Matches
FieldGather.H
Go to the documentation of this file.
1/* Copyright 2019 Axel Huebl, David Grote, Maxence Thevenet
2 * Revathi Jambunathan, Weiqun Zhang
3 *
4 * This file is part of WarpX.
5 *
6 * License: BSD-3-Clause-LBNL
7 */
8#ifndef WARPX_FIELDGATHER_H_
9#define WARPX_FIELDGATHER_H_
10
11#include "Particles/Gather/GetExternalFields.H"
12#include "Particles/Pusher/GetAndSetPosition.H"
13#include "Particles/ShapeFactors.H"
14#include "Utils/WarpX_Complex.H"
15#include "Utils/ParticleUtils.H"
16
17#include <AMReX.H>
18
34template <int depos_order_perp, int depos_order_para>
35AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
36void doDirectGatherVectorField (
37 [[maybe_unused]] const amrex::ParticleReal xp,
38 [[maybe_unused]] const amrex::ParticleReal yp,
39 [[maybe_unused]] const amrex::ParticleReal zp,
40 amrex::ParticleReal& Fxp,
41 amrex::ParticleReal& Fyp,
42 amrex::ParticleReal& Fzp,
43 amrex::Array4<amrex::Real const> const& Fx_arr,
44 amrex::Array4<amrex::Real const> const& Fy_arr,
45 amrex::Array4<amrex::Real const> const& Fz_arr,
46 const amrex::IndexType Fx_type,
47 const amrex::IndexType Fy_type,
48 const amrex::IndexType Fz_type,
49 const amrex::XDim3 & dinv,
50 const amrex::XDim3 & xyzmin,
51 const amrex::Dim3& lo,
52 [[maybe_unused]] const int n_rz_azimuthal_modes)
53{
54 using namespace amrex;
55 constexpr int NODE = amrex::IndexType::NODE;
56 constexpr int CELL = amrex::IndexType::CELL;
57
58 // Compute particle position in grid units
59#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
60 const double rp = std::sqrt(xp*xp + yp*yp);
61 const double x_bar = (rp - xyzmin.x)*dinv.x;
62#elif defined(WARPX_DIM_RSPHERE)
63 const double rp = std::sqrt(xp*xp + yp*yp + zp*zp);
64 const double x_bar = (rp - xyzmin.x)*dinv.x;
65#elif !defined(WARPX_DIM_1D_Z)
66 const double x_bar = (xp - xyzmin.x)*dinv.x;
67#endif
68#if defined(WARPX_DIM_3D)
69 const double y_bar = (yp - xyzmin.y)*dinv.y;
70#endif
71#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
72 const double z_bar = (zp - xyzmin.z)*dinv.z;
73#endif
74
75 const Compute_shape_factor< depos_order_perp > shape_factor_perp;
76 const Compute_shape_factor< depos_order_para > shape_factor_para;
77
78#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
79 constexpr int zdir = WARPX_ZINDEX;
80
81 // Compute shape factors for z-direction
82 int kp_node = 0;
83 int kp_cell = 0;
84 double sz_node[depos_order_perp+1] = {0.};
85 double sz_cell[depos_order_perp+1] = {0.};
86 if (Fx_type[zdir] == NODE || Fy_type[zdir] == NODE) {
87 kp_node = shape_factor_perp(sz_node, z_bar);
88 }
89 if (Fx_type[zdir] == CELL || Fy_type[zdir] == CELL) {
90 kp_cell = shape_factor_perp(sz_cell, z_bar - 0.5);
91 }
92
93 const int kp_Fx = ((Fx_type[zdir] == NODE) ? kp_node : kp_cell);
94 const double (&sz_Fx)[depos_order_perp+1] = ((Fx_type[zdir] == NODE) ? sz_node : sz_cell);
95
96 const int kp_Fy = ((Fy_type[zdir] == NODE) ? kp_node : kp_cell);
97 const double (&sz_Fy)[depos_order_perp+1] = ((Fy_type[zdir] == NODE) ? sz_node : sz_cell);
98
99 int kp_Fz = 0;
100 double sz_Fz[depos_order_para+1] = {0.};
101 if (Fz_type[zdir] == NODE) { kp_Fz = shape_factor_para(sz_Fz, z_bar); }
102 else { kp_Fz = shape_factor_para(sz_Fz, z_bar - 0.5); }
103#endif
104
105#if !defined(WARPX_DIM_1D_Z)
106 // Compute shape factors for x-direction
107 int ip_node = 0;
108 int ip_cell = 0;
109 double sx_node[depos_order_perp+1] = {0.};
110 double sx_cell[depos_order_perp+1] = {0.};
111 if (Fy_type[0] == NODE || Fz_type[0] == NODE) {
112 ip_node = shape_factor_perp(sx_node, x_bar);
113 }
114 if (Fy_type[0] == CELL || Fz_type[0] == CELL) {
115 ip_cell = shape_factor_perp(sx_cell, x_bar - 0.5);
116 }
117
118 const int ip_Fy = ((Fy_type[0] == NODE) ? ip_node : ip_cell);
119 const double (&sx_Fy)[depos_order_perp+1] = ((Fy_type[0] == NODE) ? sx_node : sx_cell);
120
121 const int ip_Fz = ((Fz_type[0] == NODE) ? ip_node : ip_cell);
122 const double (&sx_Fz)[depos_order_perp+1] = ((Fz_type[0] == NODE) ? sx_node : sx_cell);
123
124 int ip_Fx = 0;
125 double sx_Fx[depos_order_para + 1] = {0.};
126 if (Fx_type[0] == NODE) { ip_Fx = shape_factor_para(sx_Fx, x_bar); }
127 else { ip_Fx = shape_factor_para(sx_Fx, x_bar - 0.5); }
128#endif
129
130#if defined(WARPX_DIM_3D)
131 // Compute shape factors for y-direction
132 int jp_node = 0;
133 int jp_cell = 0;
134 double sy_node[depos_order_perp+1] = {0.};
135 double sy_cell[depos_order_perp+1] = {0.};
136 if (Fx_type[1] == NODE || Fz_type[1] == NODE) {
137 jp_node = shape_factor_perp(sy_node, y_bar);
138 }
139 if (Fx_type[1] == CELL || Fz_type[1] == CELL) {
140 jp_cell = shape_factor_perp(sy_cell, y_bar - 0.5);
141 }
142
143 const int jp_Fx = ((Fx_type[1] == NODE) ? jp_node : jp_cell);
144 const double (&sy_Fx)[depos_order_perp+1] = ((Fx_type[1] == NODE) ? sy_node : sy_cell);
145
146 const int jp_Fz = ((Fz_type[1] == NODE) ? jp_node : jp_cell);
147 const double (&sy_Fz)[depos_order_perp+1] = ((Fz_type[1] == NODE) ? sy_node : sy_cell);
148
149 int jp_Fy = 0;
150 double sy_Fy[depos_order_para + 1] = {0.};
151 if (Fy_type[1] == NODE) { jp_Fy = shape_factor_para(sy_Fy, y_bar); }
152 else { jp_Fy = shape_factor_para(sy_Fy, y_bar - 0.5); }
153
154 // Gather vector field F
155 amrex::Real weight;
156 for (int k=0; k<=depos_order_perp; k++) {
157 for (int j=0; j<=depos_order_perp; j++) {
158 for (int i=0; i<=depos_order_para; i++) {
159 weight = static_cast<amrex::Real>(sx_Fx[i]*sy_Fx[j]*sz_Fx[k]);
160 Fxp += Fx_arr(lo.x+ip_Fx+i, lo.y+jp_Fx+j, lo.z+kp_Fx+k)*weight;
161 }
162 }
163 }
164 for (int k=0; k<=depos_order_perp; k++) {
165 for (int j=0; j<=depos_order_para; j++) {
166 for (int i=0; i<=depos_order_perp; i++) {
167 weight = static_cast<amrex::Real>(sx_Fy[i]*sy_Fy[j]*sz_Fy[k]);
168 Fyp += Fy_arr(lo.x+ip_Fy+i, lo.y+jp_Fy+j, lo.z+kp_Fy+k)*weight;
169 }
170 }
171 }
172 for (int k=0; k<=depos_order_para; k++) {
173 for (int j=0; j<=depos_order_perp; j++) {
174 for (int i=0; i<=depos_order_perp; i++) {
175 weight = static_cast<amrex::Real>(sx_Fz[i]*sy_Fz[j]*sz_Fz[k]);
176 Fzp += Fz_arr(lo.x+ip_Fz+i, lo.y+jp_Fz+j, lo.z+kp_Fz+k)*weight;
177 }
178 }
179 }
180
181#elif defined(WARPX_DIM_XZ)
182
183 // Gather vector field F
184 for (int k=0; k<=depos_order_perp; k++) {
185 for (int i=0; i<=depos_order_para; i++) {
186 const auto weight_Fx = static_cast<amrex::Real>(sx_Fx[i]*sz_Fx[k]);
187 Fxp += Fx_arr(lo.x+ip_Fx+i, lo.y+kp_Fx+k, 0, 0)*weight_Fx;
188 }
189 }
190 for (int k=0; k<=depos_order_perp; k++) {
191 for (int i=0; i<=depos_order_perp; i++) {
192 const auto weight_Fy = static_cast<amrex::Real>(sx_Fy[i]*sz_Fy[k]);
193 Fyp += Fy_arr(lo.x+ip_Fy+i, lo.y+kp_Fy+k, 0, 0)*weight_Fy;
194 }
195 }
196 for (int k=0; k<=depos_order_para; k++) {
197 for (int i=0; i<=depos_order_perp; i++) {
198 const auto weight_Fz = static_cast<amrex::Real>(sx_Fz[i]*sz_Fz[k]);
199 Fzp += Fz_arr(lo.x+ip_Fz+i, lo.y+kp_Fz+k, 0, 0)*weight_Fz;
200 }
201 }
202
203#elif defined(WARPX_DIM_RZ)
204
205 amrex::ParticleReal Frp = 0.;
206 amrex::ParticleReal Fthp = 0.;
207
208 // Gather vector field F for azimuthal mode = 0
209 for (int k=0; k<=depos_order_perp; k++) {
210 for (int i=0; i<=depos_order_para; i++) {
211 const auto weight_Fr = static_cast<amrex::Real>(sx_Fx[i]*sz_Fx[k]);
212 Frp += Fx_arr(lo.x+ip_Fx+i, lo.y+kp_Fx+k, 0, 0)*weight_Fr;
213 }
214 }
215 for (int k=0; k<=depos_order_perp; k++) {
216 for (int i=0; i<=depos_order_perp; i++) {
217 const auto weight_Fth = static_cast<amrex::Real>(sx_Fy[i]*sz_Fy[k]);
218 Fthp += Fy_arr(lo.x+ip_Fy+i, lo.y+kp_Fy+k, 0, 0)*weight_Fth;
219 }
220 }
221 for (int k=0; k<=depos_order_para; k++) {
222 for (int i=0; i<=depos_order_perp; i++) {
223 const auto weight_Fz = static_cast<amrex::Real>(sx_Fz[i]*sz_Fz[k]);
224 Fzp += Fz_arr(lo.x+ip_Fz+i, lo.y+kp_Fz+k, 0, 0)*weight_Fz;
225 }
226 }
227
228 const auto costheta = static_cast<amrex::Real>(rp > 0._rt ? xp/rp: 1._rt);
229 const auto sintheta = static_cast<amrex::Real>(rp > 0._rt ? yp/rp: 0._rt);
230 const Complex xy0 = Complex{costheta, -sintheta};
231 Complex xy = xy0; // Throughout the following loop, xy takes the value e^{-i m theta}
232
233 // Gather vector field F for azimuthal modes > 0
234 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
235
236 for (int k=0; k<=depos_order_perp; k++) {
237 for (int i=0; i<=depos_order_para; i++) {
238 const auto weight_Fr = static_cast<amrex::Real>(sx_Fx[i]*sz_Fx[k]);
239 const auto dFr = (+ Fx_arr(lo.x+ip_Fx+i, lo.y+kp_Fx+k, 0, 2*imode-1)*xy.real()
240 - Fx_arr(lo.x+ip_Fx+i, lo.y+kp_Fx+k, 0, 2*imode)*xy.imag());
241 Frp += weight_Fr*dFr;
242 }
243 }
244 for (int k=0; k<=depos_order_perp; k++) {
245 for (int i=0; i<=depos_order_perp; i++) {
246 const auto weight_Fth = static_cast<amrex::Real>(sx_Fy[i]*sz_Fy[k]);
247 const auto dFth = (+ Fy_arr(lo.x+ip_Fy+i, lo.y+kp_Fy+k, 0, 2*imode-1)*xy.real()
248 - Fy_arr(lo.x+ip_Fy+i, lo.y+kp_Fy+k, 0, 2*imode)*xy.imag());
249 Fthp += weight_Fth*dFth;
250 }
251 }
252 for (int k=0; k<=depos_order_para; k++) {
253 for (int i=0; i<=depos_order_perp; i++) {
254 const auto weight_Fz = static_cast<amrex::Real>(sx_Fz[i]*sz_Fz[k]);
255 const auto dFz = (+ Fz_arr(lo.x+ip_Fz+i, lo.y+kp_Fz+k, 0, 2*imode-1)*xy.real()
256 - Fz_arr(lo.x+ip_Fz+i, lo.y+kp_Fz+k, 0, 2*imode)*xy.imag());
257 Fzp += weight_Fz*dFz;
258 }
259 }
260 xy = xy*xy0;
261
262 }
263
264 // Convert Frp and Fthp to Fxp and Fyp
265 Fxp += costheta*Frp - sintheta*Fthp;
266 Fyp += costheta*Fthp + sintheta*Frp;
267
268#elif defined(WARPX_DIM_1D_Z)
269
270 // Gather vector field F
271 for (int k=0; k<=depos_order_perp; k++) {
272 const auto weight_Fx = static_cast<amrex::Real>(sz_Fx[k]);
273 Fxp += Fx_arr(lo.x+kp_Fx+k, 0, 0)*weight_Fx;
274 //
275 const auto weight_Fy = static_cast<amrex::Real>(sz_Fy[k]);
276 Fyp += Fy_arr(lo.x+kp_Fy+k, 0, 0)*weight_Fy;
277 }
278 for (int k=0; k<=depos_order_para; k++) {
279 const auto weight_Fz = static_cast<amrex::Real>(sz_Fz[k]);
280 Fzp += Fz_arr(lo.x+kp_Fz+k, 0, 0)*weight_Fz;
281 }
282
283#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
284
285 amrex::ParticleReal Frp = 0.;
286 amrex::ParticleReal Fthetap = 0.;
287 // Z for CYLINDER, phi for SPHERE
288 amrex::ParticleReal F3p = 0.;
289
290 // Gather vector field F
291 for (int i=0; i<=depos_order_para; i++) {
292 const auto weight_Fx = static_cast<amrex::Real>(sx_Fx[i]);
293 Frp += Fx_arr(lo.x+ip_Fx+i, 0, 0)*weight_Fx;
294 }
295 for (int i=0; i<=depos_order_perp; i++) {
296 const auto weight_Fy = static_cast<amrex::Real>(sx_Fy[i]);
297 Fthetap += Fy_arr(lo.x+ip_Fy+i, 0, 0)*weight_Fy;
298 //
299 const auto weight_Fz = static_cast<amrex::Real>(sx_Fz[i]);
300 F3p += Fz_arr(lo.x+ip_Fz+i, 0, 0)*weight_Fz;
301 }
302
303#if defined(WARPX_DIM_RCYLINDER)
304
305 amrex::Real const costheta = (rp > 0. ? xp/rp: 1._rt);
306 amrex::Real const sintheta = (rp > 0. ? yp/rp: 0.);
307
308 // Convert Frp and Fthetap to Fx and Fy
309 Fxp += costheta*Frp - sintheta*Fthetap;
310 Fyp += costheta*Fthetap + sintheta*Frp;
311 Fzp += F3p;
312
313#elif defined(WARPX_DIM_RSPHERE)
314
315 amrex::Real const rpxy = std::sqrt(xp*xp + yp*yp);
316 amrex::Real const costheta = (rpxy > 0. ? xp/rpxy : 1.);
317 amrex::Real const sintheta = (rpxy > 0. ? yp/rpxy : 0.);
318 amrex::Real const cosphi = (rp > 0. ? rpxy/rp : 1.);
319 amrex::Real const sinphi = (rp > 0. ? zp/rp : 0.);
320
321 // Convert Frp, Fthetap, and Fphi to Fx, Fy, and Fz
322 Fxp += costheta*cosphi*Frp - sintheta*Fthetap - costheta*sinphi*F3p;
323 Fyp += sintheta*cosphi*Frp + costheta*Fthetap - sintheta*sinphi*F3p;
324 Fzp += sinphi*Frp + cosphi*F3p;
325
326#endif
327#endif
328}
329
348template <int depos_order, int galerkin_interpolation>
349AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
350void doGatherShapeN ([[maybe_unused]] const amrex::ParticleReal xp,
351 [[maybe_unused]] const amrex::ParticleReal yp,
352 [[maybe_unused]] const amrex::ParticleReal zp,
353 amrex::ParticleReal& Exp,
354 amrex::ParticleReal& Eyp,
355 amrex::ParticleReal& Ezp,
356 amrex::ParticleReal& Bxp,
357 amrex::ParticleReal& Byp,
358 amrex::ParticleReal& Bzp,
359 amrex::Array4<amrex::Real const> const& ex_arr,
360 amrex::Array4<amrex::Real const> const& ey_arr,
361 amrex::Array4<amrex::Real const> const& ez_arr,
362 amrex::Array4<amrex::Real const> const& bx_arr,
363 amrex::Array4<amrex::Real const> const& by_arr,
364 amrex::Array4<amrex::Real const> const& bz_arr,
365 const amrex::IndexType ex_type,
366 const amrex::IndexType ey_type,
367 const amrex::IndexType ez_type,
368 const amrex::IndexType bx_type,
369 const amrex::IndexType by_type,
370 const amrex::IndexType bz_type,
371 const amrex::XDim3 & dinv,
372 const amrex::XDim3 & xyzmin,
373 const amrex::Dim3& lo,
374 [[maybe_unused]] const int n_rz_azimuthal_modes)
375{
376 using namespace amrex;
377
378#ifdef WARPX_ZINDEX
379 constexpr int zdir = WARPX_ZINDEX;
380#endif
381
382 constexpr int NODE = amrex::IndexType::NODE;
383 constexpr int CELL = amrex::IndexType::CELL;
384
385 // --- Compute shape factors
386
387 Compute_shape_factor< depos_order > const compute_shape_factor;
388 Compute_shape_factor<depos_order - galerkin_interpolation > const compute_shape_factor_galerkin;
389
390#if !defined(WARPX_DIM_1D_Z)
391 // x direction
392 // Get particle position
393#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
394 const amrex::Real rp = std::sqrt(xp*xp + yp*yp);
395 const amrex::Real x = (rp - xyzmin.x)*dinv.x;
396#elif defined(WARPX_DIM_RSPHERE)
397 const amrex::Real rp = std::sqrt(xp*xp + yp*yp + zp*zp);
398 const amrex::Real x = (rp - xyzmin.x)*dinv.x;
399#else
400 const amrex::Real x = (xp - xyzmin.x)*dinv.x;
401#endif
402
403 // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current
404 // sx_[eb][xyz] shape factor along x for the centering of each current
405 // There are only two possible centerings, node or cell centered, so at most only two shape factor
406 // arrays will be needed.
407 amrex::Real sx_node[depos_order + 1] = {0._rt};
408 amrex::Real sx_cell[depos_order + 1] = {0._rt};
409 amrex::Real sx_node_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt};
410 amrex::Real sx_cell_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt};
411
412 int j_node = 0;
413 int j_cell = 0;
414 int j_node_v = 0;
415 int j_cell_v = 0;
416 if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) {
417 j_node = compute_shape_factor(sx_node, x);
418 }
419 if ((ey_type[0] == CELL) || (ez_type[0] == CELL) || (bx_type[0] == CELL)) {
420 j_cell = compute_shape_factor(sx_cell, x - 0.5_rt);
421 }
422 if ((ex_type[0] == NODE) || (by_type[0] == NODE) || (bz_type[0] == NODE)) {
423 j_node_v = compute_shape_factor_galerkin(sx_node_galerkin, x);
424 }
425 if ((ex_type[0] == CELL) || (by_type[0] == CELL) || (bz_type[0] == CELL)) {
426 j_cell_v = compute_shape_factor_galerkin(sx_cell_galerkin, x - 0.5_rt);
427 }
428 const amrex::Real (&sx_ex)[depos_order + 1 - galerkin_interpolation] = ((ex_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin);
429 const amrex::Real (&sx_ey)[depos_order + 1 ] = ((ey_type[0] == NODE) ? sx_node : sx_cell );
430 const amrex::Real (&sx_ez)[depos_order + 1 ] = ((ez_type[0] == NODE) ? sx_node : sx_cell );
431 const amrex::Real (&sx_bx)[depos_order + 1 ] = ((bx_type[0] == NODE) ? sx_node : sx_cell );
432 const amrex::Real (&sx_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin);
433 const amrex::Real (&sx_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin);
434 int const j_ex = ((ex_type[0] == NODE) ? j_node_v : j_cell_v);
435 int const j_ey = ((ey_type[0] == NODE) ? j_node : j_cell );
436 int const j_ez = ((ez_type[0] == NODE) ? j_node : j_cell );
437 int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell );
438 int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v);
439 int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v);
440#endif
441
442#if defined(WARPX_DIM_3D)
443 // y direction
444 const amrex::Real y = (yp-xyzmin.y)*dinv.y;
445 amrex::Real sy_node[depos_order + 1] = {0._rt};
446 amrex::Real sy_cell[depos_order + 1] = {0._rt};
447 amrex::Real sy_node_v[depos_order + 1 - galerkin_interpolation] = {0._rt};
448 amrex::Real sy_cell_v[depos_order + 1 - galerkin_interpolation] = {0._rt};
449 int k_node = 0;
450 int k_cell = 0;
451 int k_node_v = 0;
452 int k_cell_v = 0;
453 if ((ex_type[1] == NODE) || (ez_type[1] == NODE) || (by_type[1] == NODE)) {
454 k_node = compute_shape_factor(sy_node, y);
455 }
456 if ((ex_type[1] == CELL) || (ez_type[1] == CELL) || (by_type[1] == CELL)) {
457 k_cell = compute_shape_factor(sy_cell, y - 0.5_rt);
458 }
459 if ((ey_type[1] == NODE) || (bx_type[1] == NODE) || (bz_type[1] == NODE)) {
460 k_node_v = compute_shape_factor_galerkin(sy_node_v, y);
461 }
462 if ((ey_type[1] == CELL) || (bx_type[1] == CELL) || (bz_type[1] == CELL)) {
463 k_cell_v = compute_shape_factor_galerkin(sy_cell_v, y - 0.5_rt);
464 }
465 const amrex::Real (&sy_ex)[depos_order + 1 ] = ((ex_type[1] == NODE) ? sy_node : sy_cell );
466 const amrex::Real (&sy_ey)[depos_order + 1 - galerkin_interpolation] = ((ey_type[1] == NODE) ? sy_node_v : sy_cell_v);
467 const amrex::Real (&sy_ez)[depos_order + 1 ] = ((ez_type[1] == NODE) ? sy_node : sy_cell );
468 const amrex::Real (&sy_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[1] == NODE) ? sy_node_v : sy_cell_v);
469 const amrex::Real (&sy_by)[depos_order + 1 ] = ((by_type[1] == NODE) ? sy_node : sy_cell );
470 const amrex::Real (&sy_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[1] == NODE) ? sy_node_v : sy_cell_v);
471 int const k_ex = ((ex_type[1] == NODE) ? k_node : k_cell );
472 int const k_ey = ((ey_type[1] == NODE) ? k_node_v : k_cell_v);
473 int const k_ez = ((ez_type[1] == NODE) ? k_node : k_cell );
474 int const k_bx = ((bx_type[1] == NODE) ? k_node_v : k_cell_v);
475 int const k_by = ((by_type[1] == NODE) ? k_node : k_cell );
476 int const k_bz = ((bz_type[1] == NODE) ? k_node_v : k_cell_v);
477
478#endif
479
480#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
481 // z direction
482 const amrex::Real z = (zp-xyzmin.z)*dinv.z;
483 amrex::Real sz_node[depos_order + 1] = {0._rt};
484 amrex::Real sz_cell[depos_order + 1] = {0._rt};
485 amrex::Real sz_node_v[depos_order + 1 - galerkin_interpolation] = {0._rt};
486 amrex::Real sz_cell_v[depos_order + 1 - galerkin_interpolation] = {0._rt};
487 int l_node = 0;
488 int l_cell = 0;
489 int l_node_v = 0;
490 int l_cell_v = 0;
491 if ((ex_type[zdir] == NODE) || (ey_type[zdir] == NODE) || (bz_type[zdir] == NODE)) {
492 l_node = compute_shape_factor(sz_node, z);
493 }
494 if ((ex_type[zdir] == CELL) || (ey_type[zdir] == CELL) || (bz_type[zdir] == CELL)) {
495 l_cell = compute_shape_factor(sz_cell, z - 0.5_rt);
496 }
497 if ((ez_type[zdir] == NODE) || (bx_type[zdir] == NODE) || (by_type[zdir] == NODE)) {
498 l_node_v = compute_shape_factor_galerkin(sz_node_v, z);
499 }
500 if ((ez_type[zdir] == CELL) || (bx_type[zdir] == CELL) || (by_type[zdir] == CELL)) {
501 l_cell_v = compute_shape_factor_galerkin(sz_cell_v, z - 0.5_rt);
502 }
503 const amrex::Real (&sz_ex)[depos_order + 1 ] = ((ex_type[zdir] == NODE) ? sz_node : sz_cell );
504 const amrex::Real (&sz_ey)[depos_order + 1 ] = ((ey_type[zdir] == NODE) ? sz_node : sz_cell );
505 const amrex::Real (&sz_ez)[depos_order + 1 - galerkin_interpolation] = ((ez_type[zdir] == NODE) ? sz_node_v : sz_cell_v);
506 const amrex::Real (&sz_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[zdir] == NODE) ? sz_node_v : sz_cell_v);
507 const amrex::Real (&sz_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[zdir] == NODE) ? sz_node_v : sz_cell_v);
508 const amrex::Real (&sz_bz)[depos_order + 1 ] = ((bz_type[zdir] == NODE) ? sz_node : sz_cell );
509 int const l_ex = ((ex_type[zdir] == NODE) ? l_node : l_cell );
510 int const l_ey = ((ey_type[zdir] == NODE) ? l_node : l_cell );
511 int const l_ez = ((ez_type[zdir] == NODE) ? l_node_v : l_cell_v);
512 int const l_bx = ((bx_type[zdir] == NODE) ? l_node_v : l_cell_v);
513 int const l_by = ((by_type[zdir] == NODE) ? l_node_v : l_cell_v);
514 int const l_bz = ((bz_type[zdir] == NODE) ? l_node : l_cell );
515#endif
516
517 // Each field is gathered in a separate block of
518 // AMREX_SPACEDIM nested loops because the deposition
519 // order can differ for each component of each field
520 // when galerkin_interpolation is set to 1
521
522#if defined(WARPX_DIM_1D_Z)
523 // Gather field on particle Eyp from field on grid ey_arr
524 // Gather field on particle Exp from field on grid ex_arr
525 // Gather field on particle Bzp from field on grid bz_arr
526 for (int iz=0; iz<=depos_order; iz++){
527 Eyp += sz_ey[iz]*
528 ey_arr(lo.x+l_ey+iz, 0, 0, 0);
529 Exp += sz_ex[iz]*
530 ex_arr(lo.x+l_ex+iz, 0, 0, 0);
531 Bzp += sz_bz[iz]*
532 bz_arr(lo.x+l_bz+iz, 0, 0, 0);
533 }
534
535 // Gather field on particle Byp from field on grid by_arr
536 // Gather field on particle Ezp from field on grid ez_arr
537 // Gather field on particle Bxp from field on grid bx_arr
538 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
539 Ezp += sz_ez[iz]*
540 ez_arr(lo.x+l_ez+iz, 0, 0, 0);
541 Bxp += sz_bx[iz]*
542 bx_arr(lo.x+l_bx+iz, 0, 0, 0);
543 Byp += sz_by[iz]*
544 by_arr(lo.x+l_by+iz, 0, 0, 0);
545 }
546
547#elif defined(WARPX_DIM_XZ)
548 // Gather field on particle Eyp from field on grid ey_arr
549 for (int iz=0; iz<=depos_order; iz++){
550 for (int ix=0; ix<=depos_order; ix++){
551 Eyp += sx_ey[ix]*sz_ey[iz]*
552 ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 0);
553 }
554 }
555 // Gather field on particle Exp from field on grid ex_arr
556 // Gather field on particle Bzp from field on grid bz_arr
557 for (int iz=0; iz<=depos_order; iz++){
558 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
559 Exp += sx_ex[ix]*sz_ex[iz]*
560 ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 0);
561 Bzp += sx_bz[ix]*sz_bz[iz]*
562 bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 0);
563 }
564 }
565 // Gather field on particle Ezp from field on grid ez_arr
566 // Gather field on particle Bxp from field on grid bx_arr
567 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
568 for (int ix=0; ix<=depos_order; ix++){
569 Ezp += sx_ez[ix]*sz_ez[iz]*
570 ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 0);
571 Bxp += sx_bx[ix]*sz_bx[iz]*
572 bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 0);
573 }
574 }
575 // Gather field on particle Byp from field on grid by_arr
576 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
577 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
578 Byp += sx_by[ix]*sz_by[iz]*
579 by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 0);
580 }
581 }
582
583#elif defined(WARPX_DIM_RZ)
584
585 amrex::ParticleReal Erp = 0.;
586 amrex::ParticleReal Ethetap = 0.;
587 amrex::ParticleReal Brp = 0.;
588 amrex::ParticleReal Bthetap = 0.;
589
590 // Gather field on particle Ethetap from field on grid ey_arr
591 for (int iz=0; iz<=depos_order; iz++){
592 for (int ix=0; ix<=depos_order; ix++){
593 Ethetap += sx_ey[ix]*sz_ey[iz]*
594 ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 0);
595 }
596 }
597 // Gather field on particle Erp from field on grid ex_arr
598 // Gather field on particle Bzp from field on grid bz_arr
599 for (int iz=0; iz<=depos_order; iz++){
600 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
601 Erp += sx_ex[ix]*sz_ex[iz]*
602 ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 0);
603 Bzp += sx_bz[ix]*sz_bz[iz]*
604 bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 0);
605 }
606 }
607 // Gather field on particle Ezp from field on grid ez_arr
608 // Gather field on particle Brp from field on grid bx_arr
609 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
610 for (int ix=0; ix<=depos_order; ix++){
611 Ezp += sx_ez[ix]*sz_ez[iz]*
612 ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 0);
613 Brp += sx_bx[ix]*sz_bx[iz]*
614 bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 0);
615 }
616 }
617 // Gather field on particle Bthetap from field on grid by_arr
618 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
619 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
620 Bthetap += sx_by[ix]*sz_by[iz]*
621 by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 0);
622 }
623 }
624
625 amrex::Real costheta;
626 amrex::Real sintheta;
627 if (rp > 0.) {
628 costheta = xp/rp;
629 sintheta = yp/rp;
630 } else {
631 costheta = 1.;
632 sintheta = 0.;
633 }
634 const Complex xy0 = Complex{costheta, -sintheta};
635 Complex xy = xy0;
636
637 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
638
639 // Gather field on particle Ethetap from field on grid ey_arr
640 for (int iz=0; iz<=depos_order; iz++){
641 for (int ix=0; ix<=depos_order; ix++){
642 const amrex::Real dEy = (+ ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode-1)*xy.real()
643 - ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode)*xy.imag());
644 Ethetap += sx_ey[ix]*sz_ey[iz]*dEy;
645 }
646 }
647 // Gather field on particle Erp from field on grid ex_arr
648 // Gather field on particle Bzp from field on grid bz_arr
649 for (int iz=0; iz<=depos_order; iz++){
650 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
651 const amrex::Real dEx = (+ ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode-1)*xy.real()
652 - ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode)*xy.imag());
653 Erp += sx_ex[ix]*sz_ex[iz]*dEx;
654 const amrex::Real dBz = (+ bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode-1)*xy.real()
655 - bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode)*xy.imag());
656 Bzp += sx_bz[ix]*sz_bz[iz]*dBz;
657 }
658 }
659 // Gather field on particle Ezp from field on grid ez_arr
660 // Gather field on particle Brp from field on grid bx_arr
661 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
662 for (int ix=0; ix<=depos_order; ix++){
663 const amrex::Real dEz = (+ ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode-1)*xy.real()
664 - ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode)*xy.imag());
665 Ezp += sx_ez[ix]*sz_ez[iz]*dEz;
666 const amrex::Real dBx = (+ bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode-1)*xy.real()
667 - bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode)*xy.imag());
668 Brp += sx_bx[ix]*sz_bx[iz]*dBx;
669 }
670 }
671 // Gather field on particle Bthetap from field on grid by_arr
672 for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){
673 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
674 const amrex::Real dBy = (+ by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode-1)*xy.real()
675 - by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode)*xy.imag());
676 Bthetap += sx_by[ix]*sz_by[iz]*dBy;
677 }
678 }
679 xy = xy*xy0;
680 }
681
682 // Convert Erp, Ethetap, Brp, Bthetap to Exp, Eyp, Bxp, Byp
683 Exp += costheta*Erp - sintheta*Ethetap;
684 Eyp += costheta*Ethetap + sintheta*Erp;
685 Bxp += costheta*Brp - sintheta*Bthetap;
686 Byp += costheta*Bthetap + sintheta*Brp;
687
688#elif defined(WARPX_DIM_RCYLINDER)
689
690 amrex::ParticleReal Erp = 0.;
691 amrex::ParticleReal Ethetap = 0.;
692 amrex::ParticleReal Brp = 0.;
693 amrex::ParticleReal Bthetap = 0.;
694
695 // Gather field on particle Ethetap from field on grid ey_arr
696 for (int ix=0; ix<=depos_order; ix++){
697 Ethetap += sx_ey[ix]*ey_arr(lo.x+j_ey+ix, 0, 0, 0);
698 }
699 // Gather field on particle Erp from field on grid ex_arr
700 // Gather field on particle Bzp from field on grid bz_arr
701 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
702 Erp += sx_ex[ix]*ex_arr(lo.x+j_ex+ix, 0, 0, 0);
703 Bzp += sx_bz[ix]*bz_arr(lo.x+j_bz+ix, 0, 0, 0);
704 }
705 // Gather field on particle Ezp from field on grid ez_arr
706 // Gather field on particle Brp from field on grid bx_arr
707 for (int ix=0; ix<=depos_order; ix++){
708 Ezp += sx_ez[ix]*ez_arr(lo.x+j_ez+ix, 0, 0, 0);
709 Brp += sx_bx[ix]*bx_arr(lo.x+j_bx+ix, 0, 0, 0);
710 }
711 // Gather field on particle Bthetap from field on grid by_arr
712 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
713 Bthetap += sx_by[ix]*by_arr(lo.x+j_by+ix, 0, 0, 0);
714 }
715
716 amrex::Real const costheta = (rp > 0. ? xp/rp : 1.);
717 amrex::Real const sintheta = (rp > 0. ? yp/rp : 0.);
718
719 // Convert Erp, Ethetap, Brp, Bthetap to Exp, Eyp, Bxp, Byp
720 Exp += costheta*Erp - sintheta*Ethetap;
721 Eyp += costheta*Ethetap + sintheta*Erp;
722 Bxp += costheta*Brp - sintheta*Bthetap;
723 Byp += costheta*Bthetap + sintheta*Brp;
724
725#elif defined(WARPX_DIM_RSPHERE)
726
727 amrex::ParticleReal Erp = 0.;
728 amrex::ParticleReal Ethetap = 0.;
729 amrex::ParticleReal Ephip = 0.;
730 amrex::ParticleReal Brp = 0.;
731 amrex::ParticleReal Bthetap = 0.;
732 amrex::ParticleReal Bphip = 0.;
733
734 // Gather field on particle Ethetap from field on grid ey_arr
735 for (int ix=0; ix<=depos_order; ix++){
736 Ethetap += sx_ey[ix]*ey_arr(lo.x+j_ey+ix, 0, 0, 0);
737 }
738 // Gather field on particle Erp from field on grid ex_arr
739 // Gather field on particle Bphip from field on grid bz_arr
740 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
741 Erp += sx_ex[ix]*ex_arr(lo.x+j_ex+ix, 0, 0, 0);
742 Bphip += sx_bz[ix]*bz_arr(lo.x+j_bz+ix, 0, 0, 0);
743 }
744 // Gather field on particle Ephip from field on grid ez_arr
745 // Gather field on particle Brp from field on grid bx_arr
746 for (int ix=0; ix<=depos_order; ix++){
747 Ephip += sx_ez[ix]*ez_arr(lo.x+j_ez+ix, 0, 0, 0);
748 Brp += sx_bx[ix]*bx_arr(lo.x+j_bx+ix, 0, 0, 0);
749 }
750 // Gather field on particle Bthetap from field on grid by_arr
751 for (int ix=0; ix<=depos_order-galerkin_interpolation; ix++){
752 Bthetap += sx_by[ix]*by_arr(lo.x+j_by+ix, 0, 0, 0);
753 }
754
755 const amrex::Real rpxy = std::sqrt(xp*xp + yp*yp);
756 amrex::Real const costheta = (rpxy > 0. ? xp/rpxy : 1.);
757 amrex::Real const sintheta = (rpxy > 0. ? yp/rpxy : 0.);
758 amrex::Real const cosphi = (rp > 0. ? rpxy/rp : 1.);
759 amrex::Real const sinphi = (rp > 0. ? zp/rp : 0.);
760
761 // Convert Erp, Ethetap, Ephip, Brp, Bthetap, Bphip to Exp, Eyp, Ezp, Bxp, Byp, Bzp
762 Exp += costheta*cosphi*Erp - sintheta*Ethetap - costheta*sinphi*Ephip;
763 Eyp += sintheta*cosphi*Erp + costheta*Ethetap - sintheta*sinphi*Ephip;
764 Ezp += sinphi*Erp + cosphi*Ephip;
765 Bxp += costheta*cosphi*Brp - sintheta*Bthetap - costheta*sinphi*Bphip;
766 Byp += sintheta*cosphi*Brp + costheta*Bthetap - sintheta*sinphi*Bphip;
767 Bzp += sinphi*Brp + cosphi*Bphip;
768
769#else // defined(WARPX_DIM_3D)
770 // Gather field on particle Exp from field on grid ex_arr
771 for (int iz=0; iz<=depos_order; iz++){
772 for (int iy=0; iy<=depos_order; iy++){
773 for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
774 Exp += sx_ex[ix]*sy_ex[iy]*sz_ex[iz]*
775 ex_arr(lo.x+j_ex+ix, lo.y+k_ex+iy, lo.z+l_ex+iz);
776 }
777 }
778 }
779 // Gather field on particle Eyp from field on grid ey_arr
780 for (int iz=0; iz<=depos_order; iz++){
781 for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
782 for (int ix=0; ix<=depos_order; ix++){
783 Eyp += sx_ey[ix]*sy_ey[iy]*sz_ey[iz]*
784 ey_arr(lo.x+j_ey+ix, lo.y+k_ey+iy, lo.z+l_ey+iz);
785 }
786 }
787 }
788 // Gather field on particle Ezp from field on grid ez_arr
789 for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
790 for (int iy=0; iy<=depos_order; iy++){
791 for (int ix=0; ix<=depos_order; ix++){
792 Ezp += sx_ez[ix]*sy_ez[iy]*sz_ez[iz]*
793 ez_arr(lo.x+j_ez+ix, lo.y+k_ez+iy, lo.z+l_ez+iz);
794 }
795 }
796 }
797 // Gather field on particle Bzp from field on grid bz_arr
798 for (int iz=0; iz<=depos_order; iz++){
799 for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
800 for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
801 Bzp += sx_bz[ix]*sy_bz[iy]*sz_bz[iz]*
802 bz_arr(lo.x+j_bz+ix, lo.y+k_bz+iy, lo.z+l_bz+iz);
803 }
804 }
805 }
806 // Gather field on particle Byp from field on grid by_arr
807 for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
808 for (int iy=0; iy<=depos_order; iy++){
809 for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){
810 Byp += sx_by[ix]*sy_by[iy]*sz_by[iz]*
811 by_arr(lo.x+j_by+ix, lo.y+k_by+iy, lo.z+l_by+iz);
812 }
813 }
814 }
815 // Gather field on particle Bxp from field on grid bx_arr
816 for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){
817 for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){
818 for (int ix=0; ix<=depos_order; ix++){
819 Bxp += sx_bx[ix]*sy_bx[iy]*sz_bx[iz]*
820 bx_arr(lo.x+j_bx+ix, lo.y+k_bx+iy, lo.z+l_bx+iz);
821 }
822 }
823 }
824#endif
825}
826
847template <int depos_order>
848AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
849void doGatherShapeNEsirkepovStencilImplicit (
850 [[maybe_unused]] const amrex::ParticleReal xp_n,
851 [[maybe_unused]] const amrex::ParticleReal yp_n,
852 [[maybe_unused]] const amrex::ParticleReal zp_n,
853 [[maybe_unused]] const amrex::ParticleReal xp_nph,
854 [[maybe_unused]] const amrex::ParticleReal yp_nph,
855 [[maybe_unused]] const amrex::ParticleReal zp_nph,
856 amrex::ParticleReal& Exp,
857 amrex::ParticleReal& Eyp,
858 amrex::ParticleReal& Ezp,
859 amrex::ParticleReal& Bxp,
860 amrex::ParticleReal& Byp,
861 amrex::ParticleReal& Bzp,
862 amrex::Array4<amrex::Real const> const& Ex_arr,
863 amrex::Array4<amrex::Real const> const& Ey_arr,
864 amrex::Array4<amrex::Real const> const& Ez_arr,
865 amrex::Array4<amrex::Real const> const& Bx_arr,
866 amrex::Array4<amrex::Real const> const& By_arr,
867 amrex::Array4<amrex::Real const> const& Bz_arr,
868 [[maybe_unused]] const amrex::IndexType Ex_type,
869 [[maybe_unused]] const amrex::IndexType Ey_type,
870 [[maybe_unused]] const amrex::IndexType Ez_type,
871 [[maybe_unused]] const amrex::IndexType Bx_type,
872 [[maybe_unused]] const amrex::IndexType By_type,
873 [[maybe_unused]] const amrex::IndexType Bz_type,
874 const amrex::XDim3 & dinv,
875 const amrex::XDim3 & xyzmin,
876 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
877 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
878 const amrex::Dim3& lo,
879 const int n_rz_azimuthal_modes)
880{
881 using namespace amrex;
882#if !defined(WARPX_DIM_RZ)
883 ignore_unused(n_rz_azimuthal_modes);
884#endif
885
886#if (AMREX_SPACEDIM > 1)
887 Real constexpr one_third = 1.0_rt / 3.0_rt;
888 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
889#endif
890
891#if !defined(WARPX_DIM_1D_Z)
892 const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
893#endif
894#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
895 const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
896#endif
897#if !defined(WARPX_DIM_RCYLINDER)
898 const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
899#endif
900
901 // computes current and old position in grid units
902#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
903 amrex::Real const xp_new = xp_np1;
904 amrex::Real const yp_new = yp_np1;
905 amrex::Real const xp_old = xp_n;
906 amrex::Real const yp_old = yp_n;
907 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
908 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
909 // Keep these double to avoid bug in single precision
910 double x_new = (rp_new - xyzmin.x)*dinv.x;
911 double const x_old = (rp_old - xyzmin.x)*dinv.x;
912#elif defined(WARPX_DIM_RSPHERE)
913 amrex::Real const xp_new = xp_np1;
914 amrex::Real const yp_new = yp_np1;
915 amrex::Real const zp_new = zp_np1;
916 amrex::Real const xp_old = xp_n;
917 amrex::Real const yp_old = yp_n;
918 amrex::Real const zp_old = zp_n;
919 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
920 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
921
922 // Keep these double to avoid bug in single precision
923 double x_new = (rp_new - xyzmin.x)*dinv.x;
924 double const x_old = (rp_old - xyzmin.x)*dinv.x;
925#else
926#if !defined(WARPX_DIM_1D_Z)
927 // Keep these double to avoid bug in single precision
928 double x_new = (xp_np1 - xyzmin.x)*dinv.x;
929 double const x_old = (xp_n - xyzmin.x)*dinv.x;
930#endif
931#endif
932#if defined(WARPX_DIM_3D)
933 // Keep these double to avoid bug in single precision
934 double y_new = (yp_np1 - xyzmin.y)*dinv.y;
935 double const y_old = (yp_n - xyzmin.y)*dinv.y;
936#endif
937#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
938 // Keep these double to avoid bug in single precision
939 double z_new = (zp_np1 - xyzmin.z)*dinv.z;
940 double const z_old = (zp_n - xyzmin.z)*dinv.z;
941#endif
942
943 // Crop the particles at the boundaries if it is absorbing.
944 // Calculate the change in the length of the trajectory
945 // so that the time step size can be appropriately scaled.
946#if defined(WARPX_DIM_3D)
947 ParticleUtils::crop_at_boundary(x_old, y_old, z_old, x_new, y_new, z_new,
948 domain_double[0][0], domain_double[0][1], do_cropping[0]);
949 ParticleUtils::crop_at_boundary(y_old, z_old, x_old, y_new, z_new, x_new,
950 domain_double[1][0], domain_double[1][1], do_cropping[1]);
951 ParticleUtils::crop_at_boundary(z_old, x_old, y_old, z_new, x_new, y_new,
952 domain_double[2][0], domain_double[2][1], do_cropping[2]);
953#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
954 ParticleUtils::crop_at_boundary(x_old, z_old, x_new, z_new,
955 domain_double[0][0], domain_double[0][1], do_cropping[0]);
956 ParticleUtils::crop_at_boundary(z_old, x_old, z_new, x_new,
957 domain_double[1][0], domain_double[1][1], do_cropping[1]);
958#elif defined(WARPX_DIM_1D_Z)
959 ParticleUtils::crop_at_boundary(z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
960#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
961 ParticleUtils::crop_at_boundary(x_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
962#endif
963
964 // Shape factor arrays
965 // Note that there are extra values above and below
966 // to possibly hold the factor for the old particle
967 // which can be at a different grid location.
968 // Keep these double to avoid bug in single precision
969#if !defined(WARPX_DIM_1D_Z)
970 double sx_E_new[depos_order + 3] = {0.};
971 double sx_E_old[depos_order + 3] = {0.};
972#endif
973#if defined(WARPX_DIM_3D)
974 // Keep these double to avoid bug in single precision
975 double sy_E_new[depos_order + 3] = {0.};
976 double sy_E_old[depos_order + 3] = {0.};
977#endif
978#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
979 // Keep these double to avoid bug in single precision
980 double sz_E_new[depos_order + 3] = {0.};
981 double sz_E_old[depos_order + 3] = {0.};
982#endif
983
984#if defined(WARPX_DIM_3D)
985 double sx_B_new[depos_order + 3] = {0.};
986 double sx_B_old[depos_order + 3] = {0.};
987 double sy_B_new[depos_order + 3] = {0.};
988 double sy_B_old[depos_order + 3] = {0.};
989 double sz_B_new[depos_order + 3] = {0.};
990 double sz_B_old[depos_order + 3] = {0.};
991#endif
992
993#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
994 // Special shape functions are needed for By which is cell
995 // centered in both x and z. One lower order shape function is used.
996 double sx_By_new[depos_order + 2] = {0.};
997 double sx_By_old[depos_order + 2] = {0.};
998 double sz_By_new[depos_order + 2] = {0.};
999 double sz_By_old[depos_order + 2] = {0.};
1000#endif
1001
1002 // --- Compute shape factors
1003 // Compute shape factors for position as they are now and at old positions
1004 // [ijk]_new: leftmost grid point that the particle touches
1005 const Compute_shape_factor< depos_order > compute_shape_factor;
1006 const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor;
1007
1008#if !defined(WARPX_DIM_1D_Z)
1009 const int i_E_new = compute_shape_factor(sx_E_new+1, x_new);
1010 const int i_E_old = compute_shifted_shape_factor(sx_E_old, x_old, i_E_new);
1011#endif
1012#if defined(WARPX_DIM_3D)
1013 const int j_E_new = compute_shape_factor(sy_E_new+1, y_new);
1014 const int j_E_old = compute_shifted_shape_factor(sy_E_old, y_old, j_E_new);
1015#endif
1016#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1017 const int k_E_new = compute_shape_factor(sz_E_new+1, z_new);
1018 const int k_E_old = compute_shifted_shape_factor(sz_E_old, z_old, k_E_new);
1019#endif
1020
1021#if defined(WARPX_DIM_3D)
1022 const int i_B_new = compute_shape_factor(sx_B_new+1, x_new - 0.5_rt);
1023 const int i_B_old = compute_shifted_shape_factor(sx_B_old, x_old - 0.5_rt, i_B_new);
1024 const int j_B_new = compute_shape_factor(sy_B_new+1, y_new - 0.5_rt);
1025 const int j_B_old = compute_shifted_shape_factor(sy_B_old, y_old - 0.5_rt, j_B_new);
1026 const int k_B_new = compute_shape_factor(sz_B_new+1, z_new - 0.5_rt);
1027 const int k_B_old = compute_shifted_shape_factor(sz_B_old, z_old - 0.5_rt, k_B_new);
1028#endif
1029
1030 // computes min/max positions of current contributions
1031#if !defined(WARPX_DIM_1D_Z)
1032 int dil_E = 1, diu_E = 1;
1033 if (i_E_old < i_E_new) { dil_E = 0; }
1034 if (i_E_old > i_E_new) { diu_E = 0; }
1035#endif
1036#if defined(WARPX_DIM_3D)
1037 int djl_E = 1, dju_E = 1;
1038 if (j_E_old < j_E_new) { djl_E = 0; }
1039 if (j_E_old > j_E_new) { dju_E = 0; }
1040#endif
1041#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1042 int dkl_E = 1, dku_E = 1;
1043 if (k_E_old < k_E_new) { dkl_E = 0; }
1044 if (k_E_old > k_E_new) { dku_E = 0; }
1045#endif
1046
1047#if defined(WARPX_DIM_3D)
1048 int dil_B = 1, diu_B = 1;
1049 if (i_B_old < i_B_new) { dil_B = 0; }
1050 if (i_B_old > i_B_new) { diu_B = 0; }
1051 int djl_B = 1, dju_B = 1;
1052 if (j_B_old < j_B_new) { djl_B = 0; }
1053 if (j_B_old > j_B_new) { dju_B = 0; }
1054 int dkl_B = 1, dku_B = 1;
1055 if (k_B_old < k_B_new) { dkl_B = 0; }
1056 if (k_B_old > k_B_new) { dku_B = 0; }
1057#endif
1058
1059#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1060 const Compute_shape_factor< depos_order-1 > compute_shape_factor_By;
1061 const Compute_shifted_shape_factor< depos_order-1 > compute_shifted_shape_factor_By;
1062 const int i_By_new = compute_shape_factor_By(sx_By_new+1, x_new - 0.5_rt);
1063 const int i_By_old = compute_shifted_shape_factor_By(sx_By_old, x_old - 0.5_rt, i_By_new);
1064 const int k_By_new = compute_shape_factor_By(sz_By_new+1, z_new - 0.5_rt);
1065 const int k_By_old = compute_shifted_shape_factor_By(sz_By_old, z_old - 0.5_rt, k_By_new);
1066 int dil_By = 1, diu_By = 1;
1067 if (i_By_old < i_By_new) { dil_By = 0; }
1068 if (i_By_old > i_By_new) { diu_By = 0; }
1069 int dkl_By = 1, dku_By = 1;
1070 if (k_By_old < k_By_new) { dkl_By = 0; }
1071 if (k_By_old > k_By_new) { dku_By = 0; }
1072#endif
1073
1074#if defined(WARPX_DIM_3D)
1075
1076 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1077 for (int j=djl_E; j<=depos_order+2-dju_E; j++) {
1078 const amrex::Real sdzjk = one_third*(sy_E_new[j]*sz_E_new[k] + sy_E_old[j]*sz_E_old[k])
1079 +one_sixth*(sy_E_new[j]*sz_E_old[k] + sy_E_old[j]*sz_E_new[k]);
1080 amrex::Real sdxi = 0._rt;
1081 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1082 sdxi += (sx_E_old[i] - sx_E_new[i]);
1083 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1084 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;
1085 }
1086 }
1087 }
1088 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1089 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1090 const amrex::Real sdyik = one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
1091 +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]);
1092 amrex::Real sdyj = 0._rt;
1093 for (int j=djl_E; j<=depos_order+1-dju_E; j++) {
1094 sdyj += (sy_E_old[j] - sy_E_new[j]);
1095 auto sdyjov = static_cast<amrex::Real>((y_new - y_old) == 0. ? 1. : sdyj/(y_new - y_old));
1096 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;
1097 }
1098 }
1099 }
1100 for (int j=djl_E; j<=depos_order+2-dju_E; j++) {
1101 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1102 const amrex::Real sdzij = one_third*(sx_E_new[i]*sy_E_new[j] + sx_E_old[i]*sy_E_old[j])
1103 +one_sixth*(sx_E_new[i]*sy_E_old[j] + sx_E_old[i]*sy_E_new[j]);
1104 amrex::Real sdzk = 0._rt;
1105 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1106 sdzk += (sz_E_old[k] - sz_E_new[k]);
1107 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1108 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;
1109 }
1110 }
1111 }
1112 for (int k=dkl_B; k<=depos_order+2-dku_B; k++) {
1113 for (int j=djl_B; j<=depos_order+2-dju_B; j++) {
1114 const amrex::Real sdzjk = one_third*(sy_B_new[j]*sz_B_new[k] + sy_B_old[j]*sz_B_old[k])
1115 +one_sixth*(sy_B_new[j]*sz_B_old[k] + sy_B_old[j]*sz_B_new[k]);
1116 amrex::Real sdxi = 0._rt;
1117 for (int i=dil_B; i<=depos_order+1-diu_B; i++) {
1118 sdxi += (sx_B_old[i] - sx_B_new[i]);
1119 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1120 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;
1121 }
1122 }
1123 }
1124 for (int k=dkl_B; k<=depos_order+2-dku_B; k++) {
1125 for (int i=dil_B; i<=depos_order+2-diu_B; i++) {
1126 const amrex::Real sdyik = one_third*(sx_B_new[i]*sz_B_new[k] + sx_B_old[i]*sz_B_old[k])
1127 +one_sixth*(sx_B_new[i]*sz_B_old[k] + sx_B_old[i]*sz_B_new[k]);
1128 amrex::Real sdyj = 0._rt;
1129 for (int j=djl_B; j<=depos_order+1-dju_B; j++) {
1130 sdyj += (sy_B_old[j] - sy_B_new[j]);
1131 auto sdyjov = static_cast<amrex::Real>((y_new - y_old) == 0. ? 1. : sdyj/(y_new - y_old));
1132 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;
1133 }
1134 }
1135 }
1136 for (int j=djl_B; j<=depos_order+2-dju_B; j++) {
1137 for (int i=dil_B; i<=depos_order+2-diu_B; i++) {
1138 const amrex::Real sdzij = one_third*(sx_B_new[i]*sy_B_new[j] + sx_B_old[i]*sy_B_old[j])
1139 +one_sixth*(sx_B_new[i]*sy_B_old[j] + sx_B_old[i]*sy_B_new[j]);
1140 amrex::Real sdzk = 0._rt;
1141 for (int k=dkl_B; k<=depos_order+1-dku_B; k++) {
1142 sdzk += (sz_B_old[k] - sz_B_new[k]);
1143 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1144 Bzp += Bz_arr(lo.x+i_B_new-1+i, lo.y+j_B_new-1+j, lo.z+k_B_new-1+k)*sdzkov*sdzij;
1145 }
1146 }
1147 }
1148
1149#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1150
1151 amrex::ParticleReal E1p = 0._prt;
1152 amrex::ParticleReal E2p = 0._prt;
1153 amrex::ParticleReal B1p = 0._prt;
1154 amrex::ParticleReal B2p = 0._prt;
1155
1156 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1157 const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
1158 amrex::Real sdxi = 0._rt;
1159 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1160 sdxi += (sx_E_old[i] - sx_E_new[i]);
1161 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1162 E1p += Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
1163 Bzp += Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
1164 }
1165 }
1166 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1167 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1168 Real const sdyj = (
1169 one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
1170 +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]));
1171 E2p += Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdyj;
1172 }
1173 }
1174 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1175 const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]);
1176 amrex::Real sdzk = 0._rt;
1177 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1178 sdzk += (sz_E_old[k] - sz_E_new[k]);
1179 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1180 Ezp += Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
1181 B1p += Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
1182 }
1183 }
1184 for (int k=dkl_By; k<=depos_order+1-dku_By; k++) {
1185 for (int i=dil_By; i<=depos_order+1-diu_By; i++) {
1186 Real const sdyj = (
1187 one_third*(sx_By_new[i]*sz_By_new[k] + sx_By_old[i]*sz_By_old[k])
1188 +one_sixth*(sx_By_new[i]*sz_By_old[k] + sx_By_old[i]*sz_By_new[k]));
1189 B2p += By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 0)*sdyj;
1190 }
1191 }
1192
1193#ifdef WARPX_DIM_RZ
1194 const amrex::Real xp_mid = xp_nph;
1195 const amrex::Real yp_mid = yp_nph;
1196 const amrex::Real rp_mid = (rp_new + rp_old)/2._rt;
1197 const amrex::Real costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid: 1._rt);
1198 const amrex::Real sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid: 0._rt);
1199 const Complex xy_mid0 = Complex{costheta_mid, -sintheta_mid};
1200 Complex xy_mid = xy_mid0;
1201
1202 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1203
1204 // Gather field on particle E1p from field on grid ex_arr
1205 // Gather field on particle Bzp from field on grid bz_arr
1206 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1207 const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
1208 amrex::Real sdxi = 0._rt;
1209 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1210 sdxi += (sx_E_old[i] - sx_E_new[i]);
1211 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1212 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()
1213 - Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1214 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()
1215 - Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1216 E1p += dEx*sdxiov*sdzk;
1217 Bzp += dBz*sdxiov*sdzk;
1218 }
1219 }
1220 // Gather field on particle E2p from field on grid ey_arr
1221 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1222 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1223 Real const sdyj = (
1224 one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
1225 +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]));
1226 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()
1227 - Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1228 E2p += dEy*sdyj;
1229 }
1230 }
1231 // Gather field on particle Ezp from field on grid ez_arr
1232 // Gather field on particle B1p from field on grid bx_arr
1233 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1234 const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]);
1235 amrex::Real sdzk = 0._rt;
1236 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1237 sdzk += (sz_E_old[k] - sz_E_new[k]);
1238 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1239 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()
1240 - Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1241 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()
1242 - Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1243 Ezp += dEz*sdzkov*sdxi;
1244 B1p += dBx*sdzkov*sdxi;
1245 }
1246 }
1247 // Gather field on particle B2p from field on grid by_arr
1248 for (int k=dkl_By; k<=depos_order+1-dku_By; k++) {
1249 for (int i=dil_By; i<=depos_order+1-diu_By; i++) {
1250 Real const sdyj = (
1251 one_third*(sx_By_new[i]*sz_By_new[k] + sx_By_old[i]*sz_By_old[k])
1252 +one_sixth*(sx_By_new[i]*sz_By_old[k] + sx_By_old[i]*sz_By_new[k]));
1253 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()
1254 - By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 2*imode)*xy_mid.imag());
1255 B2p += dBy*sdyj;
1256 }
1257 }
1258 xy_mid = xy_mid*xy_mid0;
1259 }
1260
1261 // Convert E1p and E2p (which are actually Erp and Ethp) to Exp and Eyp
1262 Exp += costheta_mid*E1p - sintheta_mid*E2p;
1263 Eyp += costheta_mid*E2p + sintheta_mid*E1p;
1264
1265 // Convert B1p and B2p (which are actually Brp and Bthp) to Bxp and Byp
1266 Bxp += costheta_mid*B1p - sintheta_mid*B2p;
1267 Byp += costheta_mid*B2p + sintheta_mid*B1p;
1268#else
1269 Exp += E1p;
1270 Eyp += E2p;
1271
1272 Bxp += B1p;
1273 Byp += B2p;
1274#endif
1275
1276#elif defined(WARPX_DIM_1D_Z)
1277
1278 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1279 amrex::Real const sdzk = 0.5_rt*(sz_E_old[k] + sz_E_new[k]);
1280 Exp += Ex_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
1281 Eyp += Ey_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
1282 Bzp += Bz_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
1283 }
1284 amrex::Real sdzk = 0._rt;
1285 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1286 sdzk += (sz_E_old[k] - sz_E_new[k]);
1287 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1288 Bxp += Bx_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
1289 Byp += By_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
1290 Ezp += Ez_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
1291 }
1292
1293#elif defined(WARPX_DIM_RCYLINDER)
1294
1295 const amrex::Real xp_mid = xp_nph;
1296 const amrex::Real yp_mid = yp_nph;
1297 const amrex::Real rp_mid = (rp_new + rp_old)*0.5_rt;
1298 const amrex::Real costheta_mid = (rp_mid > 0. ? xp_mid/rp_mid: 1._rt);
1299 const amrex::Real sintheta_mid = (rp_mid > 0. ? yp_mid/rp_mid: 0.);
1300
1301 amrex::ParticleReal Erp = 0._prt;
1302 amrex::ParticleReal Ethetap = 0._prt;
1303 amrex::ParticleReal Brp = 0._prt;
1304 amrex::ParticleReal Bthetap = 0._prt;
1305
1306 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1307 amrex::Real const sdxi = 0.5_rt*(sx_E_old[i] + sx_E_new[i]);
1308 Brp += Bx_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1309 Ethetap += Ey_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1310 Ezp += Ez_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1311 }
1312 amrex::Real sdxi = 0._rt;
1313 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1314 sdxi += (sx_E_old[i] - sx_E_new[i]);
1315 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1316 Erp += Ex_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1317 Bthetap += By_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1318 Bzp += Bz_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1319 }
1320
1321 // Convert Erp, Ethetap, Brp, Bthetap to Exp, Eyp, Bxp, Byp
1322 Exp += costheta_mid*Erp - sintheta_mid*Ethetap;
1323 Eyp += costheta_mid*Ethetap + sintheta_mid*Erp;
1324 Bxp += costheta_mid*Brp - sintheta_mid*Bthetap;
1325 Byp += costheta_mid*Bthetap + sintheta_mid*Brp;
1326
1327#elif defined(WARPX_DIM_RSPHERE)
1328
1329 amrex::Real const xp_mid = xp_nph;
1330 amrex::Real const yp_mid = yp_nph;
1331 amrex::Real const zp_mid = zp_nph;
1332 amrex::Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
1333 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1334 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
1335 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
1336 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
1337 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
1338
1339 amrex::ParticleReal Erp = 0._prt;
1340 amrex::ParticleReal Ethetap = 0._prt;
1341 amrex::ParticleReal Ephip = 0._prt;
1342 amrex::ParticleReal Brp = 0._prt;
1343 amrex::ParticleReal Bthetap = 0._prt;
1344 amrex::ParticleReal Bphip = 0._prt;
1345
1346 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1347 amrex::Real const sdxi = 0.5_rt*(sx_E_old[i] + sx_E_new[i]);
1348 Brp += Bx_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1349 Ethetap += Ey_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1350 Ephip += Ez_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1351 }
1352 amrex::Real sdxi = 0._rt;
1353 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1354 sdxi += (sx_E_old[i] - sx_E_new[i]);
1355 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1356 Erp += Ex_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1357 Bthetap += By_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1358 Bphip += Bz_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1359 }
1360
1361 // Convert Erp, Ethetap, Ephip, Brp, Bthetap, Bphip to Exp, Eyp, Ezp, Bxp, Byp, Bzp
1362 Exp += costheta_mid*cosphi_mid*Erp - sintheta_mid*Ethetap - costheta_mid*sinphi_mid*Ephip;
1363 Eyp += sintheta_mid*cosphi_mid*Erp + costheta_mid*Ethetap - sintheta_mid*sinphi_mid*Ephip;
1364 Ezp += sinphi_mid*Erp + cosphi_mid*Ephip;
1365 Bxp += costheta_mid*cosphi_mid*Brp - sintheta_mid*Bthetap - costheta_mid*sinphi_mid*Bphip;
1366 Byp += sintheta_mid*cosphi_mid*Brp + costheta_mid*Bthetap - sintheta_mid*sinphi_mid*Bphip;
1367 Bzp += sinphi_mid*Brp + cosphi_mid*Bphip;
1368
1369#endif
1370}
1371
1393template <int depos_order>
1394AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1395void doGatherPicnicShapeN (
1396 [[maybe_unused]] const amrex::ParticleReal xp_n,
1397 [[maybe_unused]] const amrex::ParticleReal yp_n,
1398 [[maybe_unused]] const amrex::ParticleReal zp_n,
1399 [[maybe_unused]] const amrex::ParticleReal xp_nph,
1400 [[maybe_unused]] const amrex::ParticleReal yp_nph,
1401 [[maybe_unused]] const amrex::ParticleReal zp_nph,
1402 amrex::ParticleReal& Exp,
1403 amrex::ParticleReal& Eyp,
1404 amrex::ParticleReal& Ezp,
1405 amrex::ParticleReal& Bxp,
1406 amrex::ParticleReal& Byp,
1407 amrex::ParticleReal& Bzp,
1408 amrex::Array4<amrex::Real const> const& Ex_arr,
1409 amrex::Array4<amrex::Real const> const& Ey_arr,
1410 amrex::Array4<amrex::Real const> const& Ez_arr,
1411 amrex::Array4<amrex::Real const> const& Bx_arr,
1412 amrex::Array4<amrex::Real const> const& By_arr,
1413 amrex::Array4<amrex::Real const> const& Bz_arr,
1414 [[maybe_unused]] const amrex::IndexType Ex_type,
1415 [[maybe_unused]] const amrex::IndexType Ey_type,
1416 [[maybe_unused]] const amrex::IndexType Ez_type,
1417 const amrex::IndexType Bx_type,
1418 const amrex::IndexType By_type,
1419 const amrex::IndexType Bz_type,
1420 const amrex::XDim3 & dinv,
1421 const amrex::XDim3 & xyzmin,
1422 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
1423 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
1424 const amrex::Dim3& lo,
1425 const int n_rz_azimuthal_modes)
1426{
1427 using namespace amrex;
1428#if !defined(WARPX_DIM_RZ)
1429 ignore_unused(n_rz_azimuthal_modes);
1430#endif
1431
1432#if !defined(WARPX_DIM_1D_Z)
1433 const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
1434#endif
1435#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1436 const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
1437#endif
1438#if !defined(WARPX_DIM_RCYLINDER)
1439 const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
1440#endif
1441
1442#if (AMREX_SPACEDIM > 1)
1443 Real constexpr one_third = 1.0_rt / 3.0_rt;
1444 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
1445#endif
1446
1447 // computes current and old position in grid units
1448#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
1449 amrex::Real const xp_new = xp_np1;
1450 amrex::Real const yp_new = yp_np1;
1451 amrex::Real const xp_old = xp_n;
1452 amrex::Real const yp_old = yp_n;
1453 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
1454 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
1455 amrex::Real const xp_mid = xp_nph;
1456 amrex::Real const yp_mid = yp_nph;
1457 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1458 amrex::Real const costheta_mid = (rp_mid > 0. ? xp_mid/rp_mid: 1._rt);
1459 amrex::Real const sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid: 0._rt);
1460#if defined(WARPX_DIM_RZ)
1461 const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
1462#endif
1463
1464 // Keep these double to avoid bug in single precision
1465 double x_new = (rp_new - xyzmin.x)*dinv.x;
1466 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1467#elif defined(WARPX_DIM_RSPHERE)
1468 amrex::Real const xp_new = xp_np1;
1469 amrex::Real const yp_new = yp_np1;
1470 amrex::Real const zp_new = zp_np1;
1471 amrex::Real const xp_mid = xp_nph;
1472 amrex::Real const yp_mid = yp_nph;
1473 amrex::Real const zp_mid = zp_nph;
1474 amrex::Real const xp_old = xp_n;
1475 amrex::Real const yp_old = yp_n;
1476 amrex::Real const zp_old = zp_n;
1477 amrex::Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
1478 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
1479 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
1480 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1481
1482 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
1483 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
1484 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
1485 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
1486
1487 // Keep these double to avoid bug in single precision
1488 double x_new = (rp_new - xyzmin.x)*dinv.x;
1489 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1490#elif !defined(WARPX_DIM_1D_Z)
1491 // Keep these double to avoid bug in single precision
1492 double x_new = (xp_np1 - xyzmin.x)*dinv.x;
1493 double const x_old = (xp_n - xyzmin.x)*dinv.x;
1494#endif
1495#if defined(WARPX_DIM_3D)
1496 // Keep these double to avoid bug in single precision
1497 double y_new = (yp_np1 - xyzmin.y)*dinv.y;
1498 double const y_old = (yp_n - xyzmin.y)*dinv.y;
1499#endif
1500#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1501 // Keep these double to avoid bug in single precision
1502 double z_new = (zp_np1 - xyzmin.z)*dinv.z;
1503 double const z_old = (zp_n - xyzmin.z)*dinv.z;
1504#endif
1505
1506 // compute total change in particle position and the initial cell
1507 // locations in each direction used to find the position at cell crossings.
1508#if !defined(WARPX_DIM_1D_Z)
1509 const double dxp = x_new - x_old;
1510#endif
1511#if defined(WARPX_DIM_3D)
1512 const double dyp = y_new - y_old;
1513#endif
1514#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1515 const double dzp = z_new - z_old;
1516#endif
1517
1518#if defined(WARPX_DIM_3D)
1519 ParticleUtils::crop_at_boundary(x_old, y_old, z_old, x_new, y_new, z_new,
1520 domain_double[0][0], domain_double[0][1], do_cropping[0]);
1521 ParticleUtils::crop_at_boundary(y_old, z_old, x_old, y_new, z_new, x_new,
1522 domain_double[1][0], domain_double[1][1], do_cropping[1]);
1523 ParticleUtils::crop_at_boundary(z_old, x_old, y_old, z_new, x_new, y_new,
1524 domain_double[2][0], domain_double[2][1], do_cropping[2]);
1525#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1526 ParticleUtils::crop_at_boundary(x_old, z_old, x_new, z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
1527 ParticleUtils::crop_at_boundary(z_old, x_old, z_new, x_new, domain_double[1][0], domain_double[1][1], do_cropping[1]);
1528#elif defined(WARPX_DIM_1D_Z)
1529 ParticleUtils::crop_at_boundary(z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
1530#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1531 ParticleUtils::crop_at_boundary(x_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
1532#endif
1533
1534 // 1) Determine the number of segments.
1535 // 2) Loop over segments and gather electric field.
1536 // 3) Gather magnetic field.
1537
1538 // cell crossings are defined at cell edges if depos_order is odd
1539 // cell crossings are defined at cell centers if depos_order is even
1540
1541 int num_segments = 1;
1542 double shift = 0.0;
1543 if ( (depos_order % 2) == 0 ) { shift = 0.5; }
1544
1545#if defined(WARPX_DIM_3D)
1546
1547 // compute cell crossings in X-direction
1548 const auto i_old = static_cast<int>(x_old-shift);
1549 const auto i_new = static_cast<int>(x_new-shift);
1550 const int cell_crossings_x = std::abs(i_new-i_old);
1551 num_segments += cell_crossings_x;
1552
1553 // compute cell crossings in Y-direction
1554 const auto j_old = static_cast<int>(y_old-shift);
1555 const auto j_new = static_cast<int>(y_new-shift);
1556 const int cell_crossings_y = std::abs(j_new-j_old);
1557 num_segments += cell_crossings_y;
1558
1559 // compute cell crossings in Z-direction
1560 const auto k_old = static_cast<int>(z_old-shift);
1561 const auto k_new = static_cast<int>(z_new-shift);
1562 const int cell_crossings_z = std::abs(k_new-k_old);
1563 num_segments += cell_crossings_z;
1564
1565 // need to assert that the number of cell crossings in each direction
1566 // is within the range permitted by the number of guard cells
1567 // e.g., if (num_segments > 7) ...
1568
1569 // compute the initial cell location used to find the cell crossings.
1570 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1571 const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
1572 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1573 double Xcell = 0., Ycell = 0., Zcell = 0.;
1574 if (num_segments > 1) {
1575 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1576 Ycell = static_cast<double>(j_old) + shift + 0.5*(1.-dirY_sign);
1577 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1578 }
1579
1580 // loop over the number of segments and deposit
1581 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1582 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1583 double dxp_seg, dyp_seg, dzp_seg;
1584 double x0_new, y0_new, z0_new;
1585 double x0_old = x_old;
1586 double y0_old = y_old;
1587 double z0_old = z_old;
1588
1589 for (int ns=0; ns<num_segments; ns++) {
1590
1591 if (ns == num_segments-1) { // final segment
1592
1593 x0_new = x_new;
1594 y0_new = y_new;
1595 z0_new = z_new;
1596 dxp_seg = x0_new - x0_old;
1597 dyp_seg = y0_new - y0_old;
1598 dzp_seg = z0_new - z0_old;
1599
1600 }
1601 else {
1602
1603 x0_new = Xcell + dirX_sign;
1604 y0_new = Ycell + dirY_sign;
1605 z0_new = Zcell + dirZ_sign;
1606 dxp_seg = x0_new - x0_old;
1607 dyp_seg = y0_new - y0_old;
1608 dzp_seg = z0_new - z0_old;
1609
1610 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
1611 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
1612 Xcell = x0_new;
1613 dyp_seg = dyp/dxp*dxp_seg;
1614 dzp_seg = dzp/dxp*dxp_seg;
1615 y0_new = y0_old + dyp_seg;
1616 z0_new = z0_old + dzp_seg;
1617 }
1618 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
1619 Ycell = y0_new;
1620 dxp_seg = dxp/dyp*dyp_seg;
1621 dzp_seg = dzp/dyp*dyp_seg;
1622 x0_new = x0_old + dxp_seg;
1623 z0_new = z0_old + dzp_seg;
1624 }
1625 else {
1626 Zcell = z0_new;
1627 dxp_seg = dxp/dzp*dzp_seg;
1628 dyp_seg = dyp/dzp*dzp_seg;
1629 x0_new = x0_old + dxp_seg;
1630 y0_new = y0_old + dyp_seg;
1631 }
1632
1633 }
1634
1635 // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, dyp, or dzp)
1636 const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1637 const auto seg_factor_y = static_cast<double>(dyp == 0. ? 1. : dyp_seg/dyp);
1638 const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1639
1640 // compute cell-based weights using the average segment position
1641 double sx_cell[depos_order] = {0.};
1642 double sy_cell[depos_order] = {0.};
1643 double sz_cell[depos_order] = {0.};
1644 double const x0_bar = (x0_new + x0_old)/2.0;
1645 double const y0_bar = (y0_new + y0_old)/2.0;
1646 double const z0_bar = (z0_new + z0_old)/2.0;
1647 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1648 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
1649 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1650
1651 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1652 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1653 double sx_old_cell[depos_order] = {0.};
1654 double sx_new_cell[depos_order] = {0.};
1655 double sy_old_cell[depos_order] = {0.};
1656 double sy_new_cell[depos_order] = {0.};
1657 double sz_old_cell[depos_order] = {0.};
1658 double sz_new_cell[depos_order] = {0.};
1659 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1660 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
1661 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1662 ignore_unused(i0_cell_2, j0_cell_2, k0_cell_2);
1663 for (int m=0; m<depos_order; m++) {
1664 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1665 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
1666 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1667 }
1668 }
1669
1670 // compute node-based weights using the old and new segment positions
1671 double sx_old_node[depos_order+1] = {0.};
1672 double sx_new_node[depos_order+1] = {0.};
1673 double sy_old_node[depos_order+1] = {0.};
1674 double sy_new_node[depos_order+1] = {0.};
1675 double sz_old_node[depos_order+1] = {0.};
1676 double sz_new_node[depos_order+1] = {0.};
1677 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1678 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
1679 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1680
1681 // gather Ex for this segment
1682 amrex::Real weight;
1683 for (int i=0; i<=depos_order-1; i++) {
1684 for (int j=0; j<=depos_order; j++) {
1685 for (int k=0; k<=depos_order; k++) {
1686 weight = sx_cell[i]*( sy_old_node[j]*sz_old_node[k]*one_third
1687 + sy_old_node[j]*sz_new_node[k]*one_sixth
1688 + sy_new_node[j]*sz_old_node[k]*one_sixth
1689 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
1690 Exp += Ex_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k)*weight;
1691 }
1692 }
1693 }
1694
1695 // gather Ey for this segment
1696 for (int i=0; i<=depos_order; i++) {
1697 for (int j=0; j<=depos_order-1; j++) {
1698 for (int k=0; k<=depos_order; k++) {
1699 weight = sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
1700 + sx_old_node[i]*sz_new_node[k]*one_sixth
1701 + sx_new_node[i]*sz_old_node[k]*one_sixth
1702 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1703 Eyp += Ey_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k)*weight;
1704 }
1705 }
1706 }
1707
1708 // gather Ez for this segment
1709 for (int i=0; i<=depos_order; i++) {
1710 for (int j=0; j<=depos_order; j++) {
1711 for (int k=0; k<=depos_order-1; k++) {
1712 weight = sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
1713 + sx_old_node[i]*sy_new_node[j]*one_sixth
1714 + sx_new_node[i]*sy_old_node[j]*one_sixth
1715 + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
1716 Ezp += Ez_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k)*weight;
1717 }
1718 }
1719 }
1720
1721 // update old segment values
1722 if (ns < num_segments-1) {
1723 x0_old = x0_new;
1724 y0_old = y0_new;
1725 z0_old = z0_new;
1726 }
1727
1728 } // end loop over segments
1729
1730#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1731
1732 // XZ: E1p = Exp, E2p = Eyp
1733 // RZ: E1p = Erp, E2p = Ethp
1734 amrex::ParticleReal E1p = 0.0_prt;
1735 amrex::ParticleReal E2p = 0.0_prt;
1736
1737 // compute cell crossings in X-direction
1738 const auto i_old = static_cast<int>(x_old-shift);
1739 const auto i_new = static_cast<int>(x_new-shift);
1740 const int cell_crossings_x = std::abs(i_new-i_old);
1741 num_segments += cell_crossings_x;
1742
1743 // compute cell crossings in Z-direction
1744 const auto k_old = static_cast<int>(z_old-shift);
1745 const auto k_new = static_cast<int>(z_new-shift);
1746 const int cell_crossings_z = std::abs(k_new-k_old);
1747 num_segments += cell_crossings_z;
1748
1749 // need to assert that the number of cell crossings in each direction
1750 // is within the range permitted by the number of guard cells
1751 // e.g., if (num_segments > 5) ...
1752
1753 // compute the initial cell location used to find the cell crossings.
1754 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1755 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1756 double Xcell = 0., Zcell = 0.;
1757 if (num_segments > 1) {
1758 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1759 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1760 }
1761
1762 // loop over the number of segments and deposit
1763 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1764 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1765 double dxp_seg, dzp_seg;
1766 double x0_new, z0_new;
1767 double x0_old = x_old;
1768 double z0_old = z_old;
1769
1770 for (int ns=0; ns<num_segments; ns++) {
1771
1772 if (ns == num_segments-1) { // final segment
1773
1774 x0_new = x_new;
1775 z0_new = z_new;
1776 dxp_seg = x0_new - x0_old;
1777 dzp_seg = z0_new - z0_old;
1778
1779 }
1780 else {
1781
1782 x0_new = Xcell + dirX_sign;
1783 z0_new = Zcell + dirZ_sign;
1784 dxp_seg = x0_new - x0_old;
1785 dzp_seg = z0_new - z0_old;
1786
1787 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1788 Xcell = x0_new;
1789 dzp_seg = dzp/dxp*dxp_seg;
1790 z0_new = z0_old + dzp_seg;
1791 }
1792 else {
1793 Zcell = z0_new;
1794 dxp_seg = dxp/dzp*dzp_seg;
1795 x0_new = x0_old + dxp_seg;
1796 }
1797
1798 }
1799
1800 // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1801 const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1802 const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1803
1804 // compute cell-based weights using the average segment position
1805 double sx_cell[depos_order] = {0.};
1806 double sz_cell[depos_order] = {0.};
1807 double const x0_bar = (x0_new + x0_old)/2.0;
1808 double const z0_bar = (z0_new + z0_old)/2.0;
1809 const int i0_cell = compute_shape_factor_cell(sx_cell, x0_bar-0.5);
1810 const int k0_cell = compute_shape_factor_cell(sz_cell, z0_bar-0.5);
1811
1812 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1813 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1814 double sx_old_cell[depos_order] = {0.};
1815 double sx_new_cell[depos_order] = {0.};
1816 double sz_old_cell[depos_order] = {0.};
1817 double sz_new_cell[depos_order] = {0.};
1818 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1819 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1820 ignore_unused(i0_cell_2, k0_cell_2);
1821 for (int m=0; m<depos_order; m++) {
1822 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1823 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1824 }
1825 }
1826
1827 // compute node-based weights using the old and new segment positions
1828 double sx_old_node[depos_order+1] = {0.};
1829 double sx_new_node[depos_order+1] = {0.};
1830 double sz_old_node[depos_order+1] = {0.};
1831 double sz_new_node[depos_order+1] = {0.};
1832 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1833 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1834
1835 // gather Ex for this segment
1836 amrex::Real weight;
1837 for (int i=0; i<=depos_order-1; i++) {
1838 for (int k=0; k<=depos_order; k++) {
1839 weight = sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1840 E1p += Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 0)*weight;
1841#if defined(WARPX_DIM_RZ)
1842 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{-i m theta}
1843 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1844 const auto dEx = (+ Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode-1)*xy_mid.real()
1845 - Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode)*xy_mid.imag());
1846 E1p += weight*dEx;
1847 xy_mid = xy_mid*xy_mid0;
1848 }
1849#endif
1850 }
1851 }
1852
1853 // gather out-of-plane Ey for this segment
1854 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1855 for (int i=0; i<=depos_order; i++) {
1856 for (int k=0; k<=depos_order; k++) {
1857 weight = ( sx_old_node[i]*sz_old_node[k]*one_third
1858 + sx_old_node[i]*sz_new_node[k]*one_sixth
1859 + sx_new_node[i]*sz_old_node[k]*one_sixth
1860 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1861 E2p += Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 0)*weight;
1862#if defined(WARPX_DIM_RZ)
1863 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{-i m theta}
1864 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1865 const auto dEy = (+ Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode-1)*xy_mid.real()
1866 - Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode)*xy_mid.imag());
1867 E2p += weight*dEy;
1868 xy_mid = xy_mid*xy_mid0;
1869 }
1870#endif
1871 }
1872 }
1873
1874 // gather Ez for this segment
1875 for (int i=0; i<=depos_order; i++) {
1876 for (int k=0; k<=depos_order-1; k++) {
1877 weight = sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1878 Ezp += Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 0)*weight;
1879#if defined(WARPX_DIM_RZ)
1880 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{-i m theta}
1881 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1882 const auto dEz = (+ Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode-1)*xy_mid.real()
1883 - Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode)*xy_mid.imag());
1884 Ezp += weight*dEz;
1885 xy_mid = xy_mid*xy_mid0;
1886 }
1887#endif
1888 }
1889 }
1890
1891 // update old segment values
1892 if (ns < num_segments-1) {
1893 x0_old = x0_new;
1894 z0_old = z0_new;
1895 }
1896
1897 }
1898
1899#if defined(WARPX_DIM_RZ)
1900 // Convert E1p = Erp and E2p = Ethp to Exp and Eyp
1901 Exp += costheta_mid*E1p - sintheta_mid*E2p;
1902 Eyp += costheta_mid*E2p + sintheta_mid*E1p;
1903#else
1904 Exp += E1p;
1905 Eyp += E2p;
1906#endif
1907
1908#elif defined(WARPX_DIM_1D_Z)
1909
1910 // compute cell crossings in Z-direction
1911 const auto k_old = static_cast<int>(z_old-shift);
1912 const auto k_new = static_cast<int>(z_new-shift);
1913 const int cell_crossings_z = std::abs(k_new-k_old);
1914 num_segments += cell_crossings_z;
1915
1916 // need to assert that the number of cell crossings in each direction
1917 // is within the range permitted by the number of guard cells
1918 // e.g., if (num_segments > 3) ...
1919
1920 // compute the initial cell location used to find the cell crossings.
1921 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1922 double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1923
1924 // loop over the number of segments and deposit
1925 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1926 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1927 double dzp_seg;
1928 double z0_new;
1929 double z0_old = z_old;
1930
1931 for (int ns=0; ns<num_segments; ns++) {
1932
1933 if (ns == num_segments-1) { // final segment
1934 z0_new = z_new;
1935 dzp_seg = z0_new - z0_old;
1936 }
1937 else {
1938 Zcell = Zcell + dirZ_sign;
1939 z0_new = Zcell;
1940 dzp_seg = z0_new - z0_old;
1941 }
1942
1943 // compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1944 const auto seg_factor = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1945
1946 // compute cell-based weights using the average segment position
1947 double sz_cell[depos_order] = {0.};
1948 double const z0_bar = (z0_new + z0_old)/2.0;
1949 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1950
1951 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1952 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1953 double sz_old_cell[depos_order] = {0.};
1954 double sz_new_cell[depos_order] = {0.};
1955 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1956 ignore_unused(k0_cell_2);
1957 for (int m=0; m<depos_order; m++) {
1958 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1959 }
1960 }
1961
1962 // compute node-based weights using the old and new segment positions
1963 double sz_old_node[depos_order+1] = {0.};
1964 double sz_new_node[depos_order+1] = {0.};
1965 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1966
1967 // gather out-of-plane Ex and Ey for this segment
1968 for (int k=0; k<=depos_order; k++) {
1969 auto weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1970 Exp += Ex_arr(lo.x+k0_node+k, 0, 0)*weight;
1971 Eyp += Ey_arr(lo.x+k0_node+k, 0, 0)*weight;
1972 }
1973
1974 // gather Ez for this segment
1975 for (int k=0; k<=depos_order-1; k++) {
1976 auto weight = sz_cell[k]*seg_factor;
1977 Ezp += Ez_arr(lo.x+k0_cell+k, 0, 0)*weight;
1978 }
1979
1980 // update old segment values
1981 if (ns < num_segments-1) {
1982 z0_old = z0_new;
1983 }
1984
1985 }
1986
1987#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1988
1989 // RCYLINDER: E1p = Erp, E2p = Ethp, E3p = Ezp
1990 // RSPHERE: E1p = Erp, E2p = Ethp, E3p = Ephp
1991 amrex::ParticleReal E1p = 0.0_prt;
1992 amrex::ParticleReal E2p = 0.0_prt;
1993 amrex::ParticleReal E3p = 0.0_prt;
1994
1995 // compute cell crossings in X-direction
1996 const auto i_old = static_cast<int>(x_old-shift);
1997 const auto i_new = static_cast<int>(x_new-shift);
1998 const int cell_crossings_x = std::abs(i_new-i_old);
1999 num_segments += cell_crossings_x;
2000
2001 // need to assert that the number of cell crossings in each direction
2002 // is within the range permitted by the number of guard cells
2003 // e.g., if (num_segments > 3) ...
2004
2005 // compute the initial cell location used to find the cell crossings.
2006 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
2007 double Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
2008
2009 // loop over the number of segments and deposit
2010 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
2011 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
2012 double dxp_seg;
2013 double x0_new;
2014 double x0_old = x_old;
2015
2016 for (int ns=0; ns<num_segments; ns++) {
2017
2018 if (ns == num_segments-1) { // final segment
2019 x0_new = x_new;
2020 dxp_seg = x0_new - x0_old;
2021 }
2022 else {
2023 Xcell = Xcell + dirX_sign;
2024 x0_new = Xcell;
2025 dxp_seg = x0_new - x0_old;
2026 }
2027
2028 // compute the segment factor (equal to dt_seg/dt for nonzero dxp)
2029 const auto seg_factor = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
2030
2031 // compute cell-based weights using the average segment position
2032 double sx_cell[depos_order] = {0.};
2033 double const x0_bar = (x0_new + x0_old)/2.0;
2034 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
2035
2036 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
2037 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
2038 double sx_old_cell[depos_order] = {0.};
2039 double sx_new_cell[depos_order] = {0.};
2040 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
2041 ignore_unused(i0_cell_2);
2042 for (int m=0; m<depos_order; m++) {
2043 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
2044 }
2045 }
2046
2047 // compute node-based weights using the old and new segment positions
2048 double sx_old_node[depos_order+1] = {0.};
2049 double sx_new_node[depos_order+1] = {0.};
2050 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
2051
2052 // gather out-of-plane Ey and Ez for this segment
2053 for (int i=0; i<=depos_order; i++) {
2054 auto weight = 0.5_rt*(sx_old_node[i] + sx_new_node[i])*seg_factor;
2055 E2p += Ey_arr(lo.x+i0_node+i, 0, 0)*weight;
2056 E3p += Ez_arr(lo.x+i0_node+i, 0, 0)*weight;
2057 }
2058
2059 // gather Ex for this segment
2060 for (int i=0; i<=depos_order-1; i++) {
2061 auto weight = sx_cell[i]*seg_factor;
2062 E1p += Ex_arr(lo.x+i0_cell+i, 0, 0)*weight;
2063 }
2064
2065 // update old segment values
2066 if (ns < num_segments-1) {
2067 x0_old = x0_new;
2068 }
2069
2070 }
2071
2072#if defined(WARPX_DIM_RCYLINDER)
2073 // Convert E1p = Erp and E2p = Ethp to Exp and Eyp
2074 Exp += costheta_mid*E1p - sintheta_mid*E2p;
2075 Eyp += costheta_mid*E2p + sintheta_mid*E1p;
2076 Ezp += E3p;
2077#elif defined(WARPX_DIM_RSPHERE)
2078 // Convert E1p = Erp, E2p = Ethp, and E3p = Ephp to Exp, Eyp, and Ezp
2079 Exp += costheta_mid*cosphi_mid*E1p - sintheta_mid*E2p - costheta_mid*sinphi_mid*E3p;
2080 Eyp += sintheta_mid*cosphi_mid*E1p + costheta_mid*E2p - sintheta_mid*sinphi_mid*E3p;
2081 Ezp += sinphi_mid*E1p + cosphi_mid*E3p;
2082#endif
2083
2084#endif
2085
2086 // Gather magnetic field using direct gather with shape = 1.
2087 // This passes in the uncropped values of the half step position.
2088 // A future fix could pass in the average of the cropped values so
2089 // that the position would be within the domain.
2090 const int depos_order_perp = 1;
2091 const int depos_order_para = 1;
2092 amrex::ParticleReal B1p = 0._prt;
2093 amrex::ParticleReal B2p = 0._prt;
2094 amrex::ParticleReal B3p = 0._prt;
2095 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
2096#if defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
2097 rp_mid, 0._prt, 0._prt,
2098#elif defined(WARPX_DIM_RZ)
2099 rp_mid, 0._prt, zp_nph,
2100#else
2101 xp_nph, yp_nph, zp_nph,
2102#endif
2103 B1p, B2p, B3p,
2104 Bx_arr, By_arr, Bz_arr,
2105 Bx_type, By_type, Bz_type,
2106 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2107
2108 // Because we pass rp_mid and 0. instead of xp_nph and yp_nph above for
2109 // axisymmetric geometries, the returned fields on the particle are in mapped space.
2110#if defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RZ)
2111 // Convert B1p = Brp and B2p = Bthp to Bxp and Byp
2112 Bxp += costheta_mid*B1p - sintheta_mid*B2p;
2113 Byp += costheta_mid*B2p + sintheta_mid*B1p;
2114 Bzp += B3p;
2115#elif defined(WARPX_DIM_RSPHERE)
2116 // Convert B1p = Brp, B2p = Bthp, and B3p = Bphip to Bxp, Byp, and Bzp
2117 Bxp += costheta_mid*cosphi_mid*B1p - sintheta_mid*B2p - costheta_mid*sinphi_mid*B3p;
2118 Byp += sintheta_mid*cosphi_mid*B1p + costheta_mid*B2p - sintheta_mid*sinphi_mid*B3p;
2119 Bzp += sinphi_mid*B1p + cosphi_mid*B3p;
2120#else
2121 Bxp += B1p;
2122 Byp += B2p;
2123 Bzp += B3p;
2124#endif
2125
2126}
2127
2145AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2146void doGatherShapeN (const amrex::ParticleReal xp,
2147 const amrex::ParticleReal yp,
2148 const amrex::ParticleReal zp,
2149 amrex::ParticleReal& Exp,
2150 amrex::ParticleReal& Eyp,
2151 amrex::ParticleReal& Ezp,
2152 amrex::ParticleReal& Bxp,
2153 amrex::ParticleReal& Byp,
2154 amrex::ParticleReal& Bzp,
2155 amrex::Array4<amrex::Real const> const& ex_arr,
2156 amrex::Array4<amrex::Real const> const& ey_arr,
2157 amrex::Array4<amrex::Real const> const& ez_arr,
2158 amrex::Array4<amrex::Real const> const& bx_arr,
2159 amrex::Array4<amrex::Real const> const& by_arr,
2160 amrex::Array4<amrex::Real const> const& bz_arr,
2161 const amrex::IndexType ex_type,
2162 const amrex::IndexType ey_type,
2163 const amrex::IndexType ez_type,
2164 const amrex::IndexType bx_type,
2165 const amrex::IndexType by_type,
2166 const amrex::IndexType bz_type,
2167 const amrex::XDim3 & dinv,
2168 const amrex::XDim3 & xyzmin,
2169 const amrex::Dim3& lo,
2170 const int n_rz_azimuthal_modes,
2171 const int nox,
2172 const bool galerkin_interpolation)
2173{
2174 if (galerkin_interpolation) {
2175 if (nox == 1) {
2176 doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2177 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2178 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2179 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2180 } else if (nox == 2) {
2181 doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2182 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2183 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2184 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2185 } else if (nox == 3) {
2186 doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2187 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2188 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2189 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2190 } else if (nox == 4) {
2191 doGatherShapeN<4,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2192 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2193 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2194 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2195 }
2196 } else {
2197 if (nox == 1) {
2198 doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2199 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2200 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2201 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2202 } else if (nox == 2) {
2203 doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2204 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2205 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2206 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2207 } else if (nox == 3) {
2208 doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2209 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2210 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2211 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2212 } else if (nox == 4) {
2213 doGatherShapeN<4,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2214 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2215 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2216 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2217 }
2218 }
2219}
2220
2221
2242AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2243void doGatherShapeNImplicit (
2244 const amrex::ParticleReal xp_n,
2245 const amrex::ParticleReal yp_n,
2246 const amrex::ParticleReal zp_n,
2247 const amrex::ParticleReal xp_nph,
2248 const amrex::ParticleReal yp_nph,
2249 const amrex::ParticleReal zp_nph,
2250 amrex::ParticleReal& Exp,
2251 amrex::ParticleReal& Eyp,
2252 amrex::ParticleReal& Ezp,
2253 amrex::ParticleReal& Bxp,
2254 amrex::ParticleReal& Byp,
2255 amrex::ParticleReal& Bzp,
2256 amrex::Array4<amrex::Real const> const& ex_arr,
2257 amrex::Array4<amrex::Real const> const& ey_arr,
2258 amrex::Array4<amrex::Real const> const& ez_arr,
2259 amrex::Array4<amrex::Real const> const& bx_arr,
2260 amrex::Array4<amrex::Real const> const& by_arr,
2261 amrex::Array4<amrex::Real const> const& bz_arr,
2262 const amrex::IndexType ex_type,
2263 const amrex::IndexType ey_type,
2264 const amrex::IndexType ez_type,
2265 const amrex::IndexType bx_type,
2266 const amrex::IndexType by_type,
2267 const amrex::IndexType bz_type,
2268 const amrex::XDim3 & dinv,
2269 const amrex::XDim3 & xyzmin,
2270 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
2271 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
2272 const amrex::Dim3& lo,
2273 const int n_rz_azimuthal_modes,
2274 const int nox,
2275 const CurrentDepositionAlgo depos_type )
2276{
2277 if (depos_type == CurrentDepositionAlgo::Esirkepov) {
2278 if (nox == 1) {
2279 doGatherShapeNEsirkepovStencilImplicit<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2280 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2281 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2282 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2283 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2284 } else if (nox == 2) {
2285 doGatherShapeNEsirkepovStencilImplicit<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2286 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2287 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2288 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2289 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2290 } else if (nox == 3) {
2291 doGatherShapeNEsirkepovStencilImplicit<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2292 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2293 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2294 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2295 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2296 } else if (nox == 4) {
2297 doGatherShapeNEsirkepovStencilImplicit<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2298 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2299 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2300 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2301 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2302 }
2303 }
2304 else if (depos_type == CurrentDepositionAlgo::Villasenor) {
2305 if (nox == 1) {
2306 doGatherPicnicShapeN<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2307 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2308 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2309 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2310 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2311 } else if (nox == 2) {
2312 doGatherPicnicShapeN<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2313 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2314 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2315 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2316 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2317 } else if (nox == 3) {
2318 doGatherPicnicShapeN<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2319 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2320 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2321 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2322 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2323 } else if (nox == 4) {
2324 doGatherPicnicShapeN<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2325 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2326 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2327 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2328 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2329 }
2330 }
2331 else if (depos_type == CurrentDepositionAlgo::Direct) {
2332 if (nox == 1) {
2333 doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2334 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2335 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2336 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2337 } else if (nox == 2) {
2338 doGatherShapeN<2,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2339 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2340 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2341 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2342 } else if (nox == 3) {
2343 doGatherShapeN<3,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2344 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2345 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2346 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2347 } else if (nox == 4) {
2348 doGatherShapeN<4,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2349 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2350 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2351 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2352 }
2353 }
2354}
2355
2356#endif // WARPX_FIELDGATHER_H_
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
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::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Field gather for a single particle.
Definition FieldGather.H:350
doGatherPicnicShapeN
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherPicnicShapeN(const amrex::ParticleReal xp_n, const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, const amrex::ParticleReal xp_nph, const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &Ex_arr, amrex::Array4< amrex::Real const > const &Ey_arr, amrex::Array4< amrex::Real const > const &Ez_arr, amrex::Array4< amrex::Real const > const &Bx_arr, amrex::Array4< amrex::Real const > const &By_arr, amrex::Array4< amrex::Real const > const &Bz_arr, const amrex::IndexType Ex_type, const amrex::IndexType Ey_type, const amrex::IndexType Ez_type, const amrex::IndexType Bx_type, const amrex::IndexType By_type, const amrex::IndexType Bz_type, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::GpuArray< amrex::GpuArray< double, 2 >, 3 > &domain_double, const amrex::GpuArray< amrex::GpuArray< bool, 2 >, 3 > &do_cropping, 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:1395
doGatherShapeNEsirkepovStencilImplicit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeNEsirkepovStencilImplicit(const amrex::ParticleReal xp_n, const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, const amrex::ParticleReal xp_nph, const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &Ex_arr, amrex::Array4< amrex::Real const > const &Ey_arr, amrex::Array4< amrex::Real const > const &Ez_arr, amrex::Array4< amrex::Real const > const &Bx_arr, amrex::Array4< amrex::Real const > const &By_arr, amrex::Array4< amrex::Real const > const &Bz_arr, const amrex::IndexType Ex_type, const amrex::IndexType Ey_type, const amrex::IndexType Ez_type, const amrex::IndexType Bx_type, const amrex::IndexType By_type, const amrex::IndexType Bz_type, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::GpuArray< amrex::GpuArray< double, 2 >, 3 > &domain_double, const amrex::GpuArray< amrex::GpuArray< bool, 2 >, 3 > &do_cropping, 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:849
doGatherShapeNImplicit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeNImplicit(const amrex::ParticleReal xp_n, const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, const amrex::ParticleReal xp_nph, const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &ex_arr, amrex::Array4< amrex::Real const > const &ey_arr, amrex::Array4< amrex::Real const > const &ez_arr, amrex::Array4< amrex::Real const > const &bx_arr, amrex::Array4< amrex::Real const > const &by_arr, amrex::Array4< amrex::Real const > const &bz_arr, const amrex::IndexType ex_type, const amrex::IndexType ey_type, const amrex::IndexType ez_type, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::GpuArray< amrex::GpuArray< double, 2 >, 3 > &domain_double, const amrex::GpuArray< amrex::GpuArray< bool, 2 >, 3 > &do_cropping, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes, const int nox, const CurrentDepositionAlgo depos_type)
Field gather for a single particle.
Definition FieldGather.H:2243
doDirectGatherVectorField
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doDirectGatherVectorField(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::ParticleReal &Fxp, amrex::ParticleReal &Fyp, amrex::ParticleReal &Fzp, amrex::Array4< amrex::Real const > const &Fx_arr, amrex::Array4< amrex::Real const > const &Fy_arr, amrex::Array4< amrex::Real const > const &Fz_arr, const amrex::IndexType Fx_type, const amrex::IndexType Fy_type, const amrex::IndexType Fz_type, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Gather vector field F for a single particle.
Definition FieldGather.H:36
WarpX_Complex.H
Complex
amrex::GpuComplex< amrex::Real > Complex
Definition WarpX_Complex.H:22
CurrentDepositionAlgo
CurrentDepositionAlgo
Definition WarpXAlgorithmSelection.H:86
CurrentDepositionAlgo::Esirkepov
@ Esirkepov
Definition WarpXAlgorithmSelection.H:86
CurrentDepositionAlgo::Villasenor
@ Villasenor
Definition WarpXAlgorithmSelection.H:86
CurrentDepositionAlgo::Direct
@ Direct
Definition WarpXAlgorithmSelection.H:86
amrex::Real
amrex_real Real
amrex::ParticleReal
amrex_particle_real ParticleReal
amrex::Array4
ArrayND< T, 4, true > Array4
ParticleUtils::crop_at_boundary
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void crop_at_boundary(double &x1, double xmin, double xmax, amrex::GpuArray< bool, 2 > const &do_cropping)
Definition ParticleUtils.H:256
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
amrex::shift
__host__ __device__ BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
amrex::IndexType
IndexTypeND< 3 > IndexType
Compute_shape_factor_pair
Definition ShapeFactors.H:168
Compute_shape_factor
Definition ShapeFactors.H:29
Compute_shifted_shape_factor
Definition ShapeFactors.H:95
amrex::GpuComplex::real
__host__ __device__ constexpr T real() const noexcept
amrex::GpuComplex::imag
__host__ __device__ constexpr T imag() const noexcept
Generated by doxygen 1.13.2