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 amrex::Real costheta = (rp > 0. ? xp/rp: 1._rt);
229 const amrex::Real sintheta = (rp > 0. ? yp/rp: 0.);
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 and Ethetap to Ex and Ey
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 and Ethetap to Ex and Ey
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, and Ephi to Ex, Ey, and Ez
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_E_new-1+k)*sdzkov*sdzij;
1145 }
1146 }
1147 }
1148
1149#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1150
1151 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1152 const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
1153 amrex::Real sdxi = 0._rt;
1154 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1155 sdxi += (sx_E_old[i] - sx_E_new[i]);
1156 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1157 Exp += Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
1158 Bzp += Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
1159 }
1160 }
1161 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1162 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1163 Real const sdyj = (
1164 one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
1165 +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]));
1166 Eyp += Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdyj;
1167 }
1168 }
1169 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1170 const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]);
1171 amrex::Real sdzk = 0._rt;
1172 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1173 sdzk += (sz_E_old[k] - sz_E_new[k]);
1174 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1175 Ezp += Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
1176 Bxp += Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
1177 }
1178 }
1179 for (int k=dkl_By; k<=depos_order+1-dku_By; k++) {
1180 for (int i=dil_By; i<=depos_order+1-diu_By; i++) {
1181 Real const sdyj = (
1182 one_third*(sx_By_new[i]*sz_By_new[k] + sx_By_old[i]*sz_By_old[k])
1183 +one_sixth*(sx_By_new[i]*sz_By_old[k] + sx_By_old[i]*sz_By_new[k]));
1184 Byp += By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 0)*sdyj;
1185 }
1186 }
1187
1188#ifdef WARPX_DIM_RZ
1189 const amrex::Real xp_mid = xp_nph;
1190 const amrex::Real yp_mid = yp_nph;
1191 const amrex::Real rp_mid = (rp_new + rp_old)/2._rt;
1192 const amrex::Real costheta_mid = (rp_mid > 0. ? xp_mid/rp_mid: 1._rt);
1193 const amrex::Real sintheta_mid = (rp_mid > 0. ? yp_mid/rp_mid: 0.);
1194 const Complex xy_mid0 = Complex{costheta_mid, -sintheta_mid};
1195 Complex xy_mid = xy_mid0;
1196
1197 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1198
1199 // Gather field on particle Exp from field on grid ex_arr
1200 // Gather field on particle Bzp from field on grid bz_arr
1201 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1202 const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
1203 amrex::Real sdxi = 0._rt;
1204 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1205 sdxi += (sx_E_old[i] - sx_E_new[i]);
1206 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1207 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()
1208 - Ex_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1209 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()
1210 - Bz_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1211 Exp += dEx*sdxiov*sdzk;
1212 Bzp += dBz*sdxiov*sdzk;
1213 }
1214 }
1215 // Gather field on particle Eyp from field on grid ey_arr
1216 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1217 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1218 Real const sdyj = (
1219 one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k])
1220 +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]));
1221 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()
1222 - Ey_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1223 Eyp += dEy*sdyj;
1224 }
1225 }
1226 // Gather field on particle Ezp from field on grid ez_arr
1227 // Gather field on particle Bxp from field on grid bx_arr
1228 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1229 const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]);
1230 amrex::Real sdzk = 0._rt;
1231 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1232 sdzk += (sz_E_old[k] - sz_E_new[k]);
1233 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1234 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()
1235 - Ez_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1236 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()
1237 - Bx_arr(lo.x+i_E_new-1+i, lo.y+k_E_new-1+k, 0, 2*imode)*xy_mid.imag());
1238 Ezp += dEz*sdzkov*sdxi;
1239 Bxp += dBx*sdzkov*sdxi;
1240 }
1241 }
1242 // Gather field on particle Byp from field on grid by_arr
1243 for (int k=dkl_By; k<=depos_order+1-dku_By; k++) {
1244 for (int i=dil_By; i<=depos_order+1-diu_By; i++) {
1245 Real const sdyj = (
1246 one_third*(sx_By_new[i]*sz_By_new[k] + sx_By_old[i]*sz_By_old[k])
1247 +one_sixth*(sx_By_new[i]*sz_By_old[k] + sx_By_old[i]*sz_By_new[k]));
1248 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()
1249 - By_arr(lo.x+i_By_new-1+i, lo.y+k_By_new-1+k, 0, 2*imode)*xy_mid.imag());
1250 Byp += dBy*sdyj;
1251 }
1252 }
1253 xy_mid = xy_mid*xy_mid0;
1254 }
1255
1256 // Convert Exp and Eyp (which are actually Erp and Ethp) to Exp and Eyp
1257 const amrex::Real Erp = Exp;
1258 const amrex::Real Ethp = Eyp;
1259 Exp = costheta_mid*Erp - sintheta_mid*Ethp;
1260 Eyp = costheta_mid*Ethp + sintheta_mid*Erp;
1261 // Convert Bxp and Byp (which are actually Brp and Bthp) to Bxp and Byp
1262 const amrex::Real Brp = Bxp;
1263 const amrex::Real Bthp = Byp;
1264 Bxp = costheta_mid*Brp - sintheta_mid*Bthp;
1265 Byp = costheta_mid*Bthp + sintheta_mid*Brp;
1266
1267#endif
1268
1269#elif defined(WARPX_DIM_1D_Z)
1270
1271 for (int k=dkl_E; k<=depos_order+2-dku_E; k++) {
1272 amrex::Real const sdzk = 0.5_rt*(sz_E_old[k] + sz_E_new[k]);
1273 Exp += Ex_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
1274 Eyp += Ey_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
1275 Bzp += Bz_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzk;
1276 }
1277 amrex::Real sdzk = 0._rt;
1278 for (int k=dkl_E; k<=depos_order+1-dku_E; k++) {
1279 sdzk += (sz_E_old[k] - sz_E_new[k]);
1280 auto sdzkov = static_cast<amrex::Real>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
1281 Bxp += Bx_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
1282 Byp += By_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
1283 Ezp += Ez_arr(lo.x+k_E_new-1+k, 0, 0, 0)*sdzkov;
1284 }
1285
1286#elif defined(WARPX_DIM_RCYLINDER)
1287
1288 const amrex::Real xp_mid = xp_nph;
1289 const amrex::Real yp_mid = yp_nph;
1290 const amrex::Real rp_mid = (rp_new + rp_old)*0.5_rt;
1291 const amrex::Real costheta_mid = (rp_mid > 0. ? xp_mid/rp_mid: 1._rt);
1292 const amrex::Real sintheta_mid = (rp_mid > 0. ? yp_mid/rp_mid: 0.);
1293
1294 amrex::ParticleReal Erp = 0.;
1295 amrex::ParticleReal Ethetap = 0.;
1296 amrex::ParticleReal Brp = 0.;
1297 amrex::ParticleReal Bthetap = 0.;
1298
1299 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1300 amrex::Real const sdxi = 0.5_rt*(sx_E_old[i] + sx_E_new[i]);
1301 Brp += Bx_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1302 Ethetap += Ey_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1303 Ezp += Ez_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1304 }
1305 amrex::Real sdxi = 0._rt;
1306 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1307 sdxi += (sx_E_old[i] - sx_E_new[i]);
1308 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1309 Erp += Ex_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1310 Bthetap += By_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1311 Bzp += Bz_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1312 }
1313
1314 // Convert Erp and Ethetap to Ex and Ey
1315 Exp += costheta_mid*Erp - sintheta_mid*Ethetap;
1316 Eyp += costheta_mid*Ethetap + sintheta_mid*Erp;
1317 Bxp += costheta_mid*Brp - sintheta_mid*Bthetap;
1318 Byp += costheta_mid*Bthetap + sintheta_mid*Brp;
1319
1320#elif defined(WARPX_DIM_RSPHERE)
1321
1322 amrex::Real const xp_mid = xp_nph;
1323 amrex::Real const yp_mid = yp_nph;
1324 amrex::Real const zp_mid = zp_nph;
1325 amrex::Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
1326 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1327 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
1328 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
1329 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
1330 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
1331
1332 amrex::ParticleReal Erp = 0.;
1333 amrex::ParticleReal Ethetap = 0.;
1334 amrex::ParticleReal Ephip = 0.;
1335 amrex::ParticleReal Brp = 0.;
1336 amrex::ParticleReal Bthetap = 0.;
1337 amrex::ParticleReal Bphip = 0.;
1338
1339 for (int i=dil_E; i<=depos_order+2-diu_E; i++) {
1340 amrex::Real const sdxi = 0.5_rt*(sx_E_old[i] + sx_E_new[i]);
1341 Brp += Bx_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1342 Ethetap += Ey_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1343 Ephip += Ez_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxi;
1344 }
1345 amrex::Real sdxi = 0._rt;
1346 for (int i=dil_E; i<=depos_order+1-diu_E; i++) {
1347 sdxi += (sx_E_old[i] - sx_E_new[i]);
1348 auto sdxiov = static_cast<amrex::Real>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
1349 Erp += Ex_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1350 Bthetap += By_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1351 Bphip += Bz_arr(lo.x+i_E_new-1+i, 0, 0, 0)*sdxiov;
1352 }
1353
1354 // Convert Erp, Ethetap, and Ephi to Ex, Ey, and Ez
1355 Exp += costheta_mid*cosphi_mid*Erp - sintheta_mid*Ethetap - costheta_mid*sinphi_mid*Ephip;
1356 Eyp += sintheta_mid*cosphi_mid*Erp + costheta_mid*Ethetap - sintheta_mid*sinphi_mid*Ephip;
1357 Ezp += sinphi_mid*Erp + cosphi_mid*Ephip;
1358 Bxp += costheta_mid*cosphi_mid*Brp - sintheta_mid*Bthetap - costheta_mid*sinphi_mid*Bphip;
1359 Byp += sintheta_mid*cosphi_mid*Brp + costheta_mid*Bthetap - sintheta_mid*sinphi_mid*Bphip;
1360 Bzp += sinphi_mid*Brp + cosphi_mid*Bphip;
1361
1362#endif
1363}
1364
1386template <int depos_order>
1387AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1388void doGatherPicnicShapeN (
1389 [[maybe_unused]] const amrex::ParticleReal xp_n,
1390 [[maybe_unused]] const amrex::ParticleReal yp_n,
1391 [[maybe_unused]] const amrex::ParticleReal zp_n,
1392 [[maybe_unused]] const amrex::ParticleReal xp_nph,
1393 [[maybe_unused]] const amrex::ParticleReal yp_nph,
1394 [[maybe_unused]] const amrex::ParticleReal zp_nph,
1395 amrex::ParticleReal& Exp,
1396 amrex::ParticleReal& Eyp,
1397 amrex::ParticleReal& Ezp,
1398 amrex::ParticleReal& Bxp,
1399 amrex::ParticleReal& Byp,
1400 amrex::ParticleReal& Bzp,
1401 amrex::Array4<amrex::Real const> const& Ex_arr,
1402 amrex::Array4<amrex::Real const> const& Ey_arr,
1403 amrex::Array4<amrex::Real const> const& Ez_arr,
1404 amrex::Array4<amrex::Real const> const& Bx_arr,
1405 amrex::Array4<amrex::Real const> const& By_arr,
1406 amrex::Array4<amrex::Real const> const& Bz_arr,
1407 [[maybe_unused]] const amrex::IndexType Ex_type,
1408 [[maybe_unused]] const amrex::IndexType Ey_type,
1409 [[maybe_unused]] const amrex::IndexType Ez_type,
1410 const amrex::IndexType Bx_type,
1411 const amrex::IndexType By_type,
1412 const amrex::IndexType Bz_type,
1413 const amrex::XDim3 & dinv,
1414 const amrex::XDim3 & xyzmin,
1415 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
1416 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
1417 const amrex::Dim3& lo,
1418 const int n_rz_azimuthal_modes)
1419{
1420 using namespace amrex;
1421#if !defined(WARPX_DIM_RZ)
1422 ignore_unused(n_rz_azimuthal_modes);
1423#endif
1424
1425#if !defined(WARPX_DIM_1D_Z)
1426 const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
1427#endif
1428#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1429 const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
1430#endif
1431#if !defined(WARPX_DIM_RCYLINDER)
1432 const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
1433#endif
1434
1435#if (AMREX_SPACEDIM > 1)
1436 Real constexpr one_third = 1.0_rt / 3.0_rt;
1437 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
1438#endif
1439
1440 // computes current and old position in grid units
1441#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
1442 amrex::Real const xp_new = xp_np1;
1443 amrex::Real const yp_new = yp_np1;
1444 amrex::Real const xp_old = xp_n;
1445 amrex::Real const yp_old = yp_n;
1446 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
1447 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
1448 amrex::Real const xp_mid = xp_nph;
1449 amrex::Real const yp_mid = yp_nph;
1450 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1451 amrex::Real const costheta_mid = (rp_mid > 0. ? xp_mid/rp_mid: 1._rt);
1452 amrex::Real const sintheta_mid = (rp_mid > 0. ? yp_mid/rp_mid: 0.);
1453#if defined(WARPX_DIM_RZ)
1454 const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
1455#endif
1456
1457 // Keep these double to avoid bug in single precision
1458 double x_new = (rp_new - xyzmin.x)*dinv.x;
1459 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1460#elif defined(WARPX_DIM_RSPHERE)
1461 amrex::Real const xp_new = xp_np1;
1462 amrex::Real const yp_new = yp_np1;
1463 amrex::Real const zp_new = zp_np1;
1464 amrex::Real const xp_mid = xp_nph;
1465 amrex::Real const yp_mid = yp_nph;
1466 amrex::Real const zp_mid = zp_nph;
1467 amrex::Real const xp_old = xp_n;
1468 amrex::Real const yp_old = yp_n;
1469 amrex::Real const zp_old = zp_n;
1470 amrex::Real const rpxy_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
1471 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
1472 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
1473 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
1474
1475 amrex::Real const costheta_mid = (rpxy_mid > 0. ? xp_mid/rpxy_mid : 1._rt);
1476 amrex::Real const sintheta_mid = (rpxy_mid > 0. ? yp_mid/rpxy_mid : 0._rt);
1477 amrex::Real const cosphi_mid = (rp_mid > 0. ? rpxy_mid/rp_mid : 1._rt);
1478 amrex::Real const sinphi_mid = (rp_mid > 0. ? zp_mid/rp_mid : 0._rt);
1479
1480 // Keep these double to avoid bug in single precision
1481 double x_new = (rp_new - xyzmin.x)*dinv.x;
1482 double const x_old = (rp_old - xyzmin.x)*dinv.x;
1483#elif !defined(WARPX_DIM_1D_Z)
1484 // Keep these double to avoid bug in single precision
1485 double x_new = (xp_np1 - xyzmin.x)*dinv.x;
1486 double const x_old = (xp_n - xyzmin.x)*dinv.x;
1487#endif
1488#if defined(WARPX_DIM_3D)
1489 // Keep these double to avoid bug in single precision
1490 double y_new = (yp_np1 - xyzmin.y)*dinv.y;
1491 double const y_old = (yp_n - xyzmin.y)*dinv.y;
1492#endif
1493#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1494 // Keep these double to avoid bug in single precision
1495 double z_new = (zp_np1 - xyzmin.z)*dinv.z;
1496 double const z_old = (zp_n - xyzmin.z)*dinv.z;
1497#endif
1498
1499 // compute total change in particle position and the initial cell
1500 // locations in each direction used to find the position at cell crossings.
1501#if !defined(WARPX_DIM_1D_Z)
1502 const double dxp = x_new - x_old;
1503#endif
1504#if defined(WARPX_DIM_3D)
1505 const double dyp = y_new - y_old;
1506#endif
1507#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
1508 const double dzp = z_new - z_old;
1509#endif
1510
1511#if defined(WARPX_DIM_3D)
1512 ParticleUtils::crop_at_boundary(x_old, y_old, z_old, x_new, y_new, z_new,
1513 domain_double[0][0], domain_double[0][1], do_cropping[0]);
1514 ParticleUtils::crop_at_boundary(y_old, z_old, x_old, y_new, z_new, x_new,
1515 domain_double[1][0], domain_double[1][1], do_cropping[1]);
1516 ParticleUtils::crop_at_boundary(z_old, x_old, y_old, z_new, x_new, y_new,
1517 domain_double[2][0], domain_double[2][1], do_cropping[2]);
1518#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1519 ParticleUtils::crop_at_boundary(x_old, z_old, x_new, z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
1520 ParticleUtils::crop_at_boundary(z_old, x_old, z_new, x_new, domain_double[1][0], domain_double[1][1], do_cropping[1]);
1521#elif defined(WARPX_DIM_1D_Z)
1522 ParticleUtils::crop_at_boundary(z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
1523#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1524 ParticleUtils::crop_at_boundary(x_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
1525#endif
1526
1527 // 1) Determine the number of segments.
1528 // 2) Loop over segments and gather electric field.
1529 // 3) Gather magnetic field.
1530
1531 // cell crossings are defined at cell edges if depos_order is odd
1532 // cell crossings are defined at cell centers if depos_order is even
1533
1534 int num_segments = 1;
1535 double shift = 0.0;
1536 if ( (depos_order % 2) == 0 ) { shift = 0.5; }
1537
1538#if defined(WARPX_DIM_3D)
1539
1540 // compute cell crossings in X-direction
1541 const auto i_old = static_cast<int>(x_old-shift);
1542 const auto i_new = static_cast<int>(x_new-shift);
1543 const int cell_crossings_x = std::abs(i_new-i_old);
1544 num_segments += cell_crossings_x;
1545
1546 // compute cell crossings in Y-direction
1547 const auto j_old = static_cast<int>(y_old-shift);
1548 const auto j_new = static_cast<int>(y_new-shift);
1549 const int cell_crossings_y = std::abs(j_new-j_old);
1550 num_segments += cell_crossings_y;
1551
1552 // compute cell crossings in Z-direction
1553 const auto k_old = static_cast<int>(z_old-shift);
1554 const auto k_new = static_cast<int>(z_new-shift);
1555 const int cell_crossings_z = std::abs(k_new-k_old);
1556 num_segments += cell_crossings_z;
1557
1558 // need to assert that the number of cell crossings in each direction
1559 // is within the range permitted by the number of guard cells
1560 // e.g., if (num_segments > 7) ...
1561
1562 // compute the initial cell location used to find the cell crossings.
1563 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1564 const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
1565 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1566 double Xcell = 0., Ycell = 0., Zcell = 0.;
1567 if (num_segments > 1) {
1568 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1569 Ycell = static_cast<double>(j_old) + shift + 0.5*(1.-dirY_sign);
1570 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1571 }
1572
1573 // loop over the number of segments and deposit
1574 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1575 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1576 double dxp_seg, dyp_seg, dzp_seg;
1577 double x0_new, y0_new, z0_new;
1578 double x0_old = x_old;
1579 double y0_old = y_old;
1580 double z0_old = z_old;
1581
1582 for (int ns=0; ns<num_segments; ns++) {
1583
1584 if (ns == num_segments-1) { // final segment
1585
1586 x0_new = x_new;
1587 y0_new = y_new;
1588 z0_new = z_new;
1589 dxp_seg = x0_new - x0_old;
1590 dyp_seg = y0_new - y0_old;
1591 dzp_seg = z0_new - z0_old;
1592
1593 }
1594 else {
1595
1596 x0_new = Xcell + dirX_sign;
1597 y0_new = Ycell + dirY_sign;
1598 z0_new = Zcell + dirZ_sign;
1599 dxp_seg = x0_new - x0_old;
1600 dyp_seg = y0_new - y0_old;
1601 dzp_seg = z0_new - z0_old;
1602
1603 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
1604 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
1605 Xcell = x0_new;
1606 dyp_seg = dyp/dxp*dxp_seg;
1607 dzp_seg = dzp/dxp*dxp_seg;
1608 y0_new = y0_old + dyp_seg;
1609 z0_new = z0_old + dzp_seg;
1610 }
1611 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
1612 Ycell = y0_new;
1613 dxp_seg = dxp/dyp*dyp_seg;
1614 dzp_seg = dzp/dyp*dyp_seg;
1615 x0_new = x0_old + dxp_seg;
1616 z0_new = z0_old + dzp_seg;
1617 }
1618 else {
1619 Zcell = z0_new;
1620 dxp_seg = dxp/dzp*dzp_seg;
1621 dyp_seg = dyp/dzp*dzp_seg;
1622 x0_new = x0_old + dxp_seg;
1623 y0_new = y0_old + dyp_seg;
1624 }
1625
1626 }
1627
1628 // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, dyp, or dzp)
1629 const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1630 const auto seg_factor_y = static_cast<double>(dyp == 0. ? 1. : dyp_seg/dyp);
1631 const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1632
1633 // compute cell-based weights using the average segment position
1634 double sx_cell[depos_order] = {0.};
1635 double sy_cell[depos_order] = {0.};
1636 double sz_cell[depos_order] = {0.};
1637 double const x0_bar = (x0_new + x0_old)/2.0;
1638 double const y0_bar = (y0_new + y0_old)/2.0;
1639 double const z0_bar = (z0_new + z0_old)/2.0;
1640 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1641 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
1642 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1643
1644 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1645 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1646 double sx_old_cell[depos_order] = {0.};
1647 double sx_new_cell[depos_order] = {0.};
1648 double sy_old_cell[depos_order] = {0.};
1649 double sy_new_cell[depos_order] = {0.};
1650 double sz_old_cell[depos_order] = {0.};
1651 double sz_new_cell[depos_order] = {0.};
1652 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1653 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
1654 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1655 ignore_unused(i0_cell_2, j0_cell_2, k0_cell_2);
1656 for (int m=0; m<depos_order; m++) {
1657 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1658 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
1659 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1660 }
1661 }
1662
1663 // compute node-based weights using the old and new segment positions
1664 double sx_old_node[depos_order+1] = {0.};
1665 double sx_new_node[depos_order+1] = {0.};
1666 double sy_old_node[depos_order+1] = {0.};
1667 double sy_new_node[depos_order+1] = {0.};
1668 double sz_old_node[depos_order+1] = {0.};
1669 double sz_new_node[depos_order+1] = {0.};
1670 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1671 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
1672 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1673
1674 // gather Ex for this segment
1675 amrex::Real weight;
1676 for (int i=0; i<=depos_order-1; i++) {
1677 for (int j=0; j<=depos_order; j++) {
1678 for (int k=0; k<=depos_order; k++) {
1679 weight = sx_cell[i]*( sy_old_node[j]*sz_old_node[k]*one_third
1680 + sy_old_node[j]*sz_new_node[k]*one_sixth
1681 + sy_new_node[j]*sz_old_node[k]*one_sixth
1682 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
1683 Exp += Ex_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k)*weight;
1684 }
1685 }
1686 }
1687
1688 // gather Ey for this segment
1689 for (int i=0; i<=depos_order; i++) {
1690 for (int j=0; j<=depos_order-1; j++) {
1691 for (int k=0; k<=depos_order; k++) {
1692 weight = sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
1693 + sx_old_node[i]*sz_new_node[k]*one_sixth
1694 + sx_new_node[i]*sz_old_node[k]*one_sixth
1695 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1696 Eyp += Ey_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k)*weight;
1697 }
1698 }
1699 }
1700
1701 // gather Ez for this segment
1702 for (int i=0; i<=depos_order; i++) {
1703 for (int j=0; j<=depos_order; j++) {
1704 for (int k=0; k<=depos_order-1; k++) {
1705 weight = sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
1706 + sx_old_node[i]*sy_new_node[j]*one_sixth
1707 + sx_new_node[i]*sy_old_node[j]*one_sixth
1708 + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
1709 Ezp += Ez_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k)*weight;
1710 }
1711 }
1712 }
1713
1714 // update old segment values
1715 if (ns < num_segments-1) {
1716 x0_old = x0_new;
1717 y0_old = y0_new;
1718 z0_old = z0_new;
1719 }
1720
1721 } // end loop over segments
1722
1723#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1724
1725 // compute cell crossings in X-direction
1726 const auto i_old = static_cast<int>(x_old-shift);
1727 const auto i_new = static_cast<int>(x_new-shift);
1728 const int cell_crossings_x = std::abs(i_new-i_old);
1729 num_segments += cell_crossings_x;
1730
1731 // compute cell crossings in Z-direction
1732 const auto k_old = static_cast<int>(z_old-shift);
1733 const auto k_new = static_cast<int>(z_new-shift);
1734 const int cell_crossings_z = std::abs(k_new-k_old);
1735 num_segments += cell_crossings_z;
1736
1737 // need to assert that the number of cell crossings in each direction
1738 // is within the range permitted by the number of guard cells
1739 // e.g., if (num_segments > 5) ...
1740
1741 // compute the initial cell location used to find the cell crossings.
1742 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1743 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1744 double Xcell = 0., Zcell = 0.;
1745 if (num_segments > 1) {
1746 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1747 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1748 }
1749
1750 // loop over the number of segments and deposit
1751 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1752 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1753 double dxp_seg, dzp_seg;
1754 double x0_new, z0_new;
1755 double x0_old = x_old;
1756 double z0_old = z_old;
1757
1758 for (int ns=0; ns<num_segments; ns++) {
1759
1760 if (ns == num_segments-1) { // final segment
1761
1762 x0_new = x_new;
1763 z0_new = z_new;
1764 dxp_seg = x0_new - x0_old;
1765 dzp_seg = z0_new - z0_old;
1766
1767 }
1768 else {
1769
1770 x0_new = Xcell + dirX_sign;
1771 z0_new = Zcell + dirZ_sign;
1772 dxp_seg = x0_new - x0_old;
1773 dzp_seg = z0_new - z0_old;
1774
1775 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1776 Xcell = x0_new;
1777 dzp_seg = dzp/dxp*dxp_seg;
1778 z0_new = z0_old + dzp_seg;
1779 }
1780 else {
1781 Zcell = z0_new;
1782 dxp_seg = dxp/dzp*dzp_seg;
1783 x0_new = x0_old + dxp_seg;
1784 }
1785
1786 }
1787
1788 // compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1789 const auto seg_factor_x = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1790 const auto seg_factor_z = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1791
1792 // compute cell-based weights using the average segment position
1793 double sx_cell[depos_order] = {0.};
1794 double sz_cell[depos_order] = {0.};
1795 double const x0_bar = (x0_new + x0_old)/2.0;
1796 double const z0_bar = (z0_new + z0_old)/2.0;
1797 const int i0_cell = compute_shape_factor_cell(sx_cell, x0_bar-0.5);
1798 const int k0_cell = compute_shape_factor_cell(sz_cell, z0_bar-0.5);
1799
1800 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1801 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1802 double sx_old_cell[depos_order] = {0.};
1803 double sx_new_cell[depos_order] = {0.};
1804 double sz_old_cell[depos_order] = {0.};
1805 double sz_new_cell[depos_order] = {0.};
1806 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1807 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1808 ignore_unused(i0_cell_2, k0_cell_2);
1809 for (int m=0; m<depos_order; m++) {
1810 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1811 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1812 }
1813 }
1814
1815 // compute node-based weights using the old and new segment positions
1816 double sx_old_node[depos_order+1] = {0.};
1817 double sx_new_node[depos_order+1] = {0.};
1818 double sz_old_node[depos_order+1] = {0.};
1819 double sz_new_node[depos_order+1] = {0.};
1820 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1821 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1822
1823 // gather Ex for this segment
1824 amrex::Real weight;
1825 for (int i=0; i<=depos_order-1; i++) {
1826 for (int k=0; k<=depos_order; k++) {
1827 weight = sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1828 Exp += Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 0)*weight;
1829#if defined(WARPX_DIM_RZ)
1830 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{-i m theta}
1831 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1832 const auto dEx = (+ Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode-1)*xy_mid.real()
1833 - Ex_arr(lo.x+i0_cell+i, lo.y+k0_node+k, 0, 2*imode)*xy_mid.imag());
1834 Exp += weight*dEx;
1835 xy_mid = xy_mid*xy_mid0;
1836 }
1837#endif
1838 }
1839 }
1840
1841 // gather out-of-plane Ey for this segment
1842 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1843 for (int i=0; i<=depos_order; i++) {
1844 for (int k=0; k<=depos_order; k++) {
1845 weight = ( sx_old_node[i]*sz_old_node[k]*one_third
1846 + sx_old_node[i]*sz_new_node[k]*one_sixth
1847 + sx_new_node[i]*sz_old_node[k]*one_sixth
1848 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1849 Eyp += Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 0)*weight;
1850#if defined(WARPX_DIM_RZ)
1851 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{-i m theta}
1852 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1853 const auto dEy = (+ Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode-1)*xy_mid.real()
1854 - Ey_arr(lo.x+i0_node+i, lo.y+k0_node+k, 0, 2*imode)*xy_mid.imag());
1855 Eyp += weight*dEy;
1856 xy_mid = xy_mid*xy_mid0;
1857 }
1858#endif
1859 }
1860 }
1861
1862 // gather Ez for this segment
1863 for (int i=0; i<=depos_order; i++) {
1864 for (int k=0; k<=depos_order-1; k++) {
1865 weight = sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1866 Ezp += Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 0)*weight;
1867#if defined(WARPX_DIM_RZ)
1868 Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{-i m theta}
1869 for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1870 const auto dEz = (+ Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode-1)*xy_mid.real()
1871 - Ez_arr(lo.x+i0_node+i, lo.y+k0_cell+k, 0, 2*imode)*xy_mid.imag());
1872 Ezp += weight*dEz;
1873 xy_mid = xy_mid*xy_mid0;
1874 }
1875#endif
1876 }
1877 }
1878
1879 // update old segment values
1880 if (ns < num_segments-1) {
1881 x0_old = x0_new;
1882 z0_old = z0_new;
1883 }
1884
1885 }
1886
1887#ifdef WARPX_DIM_RZ
1888
1889 // Convert Exp and Eyp (which are actually Erp and Ethp) to Exp and Eyp
1890 const amrex::Real Erp = Exp;
1891 const amrex::Real Ethp = Eyp;
1892 Exp = costheta_mid*Erp - sintheta_mid*Ethp;
1893 Eyp = costheta_mid*Ethp + sintheta_mid*Erp;
1894
1895#endif
1896
1897#elif defined(WARPX_DIM_1D_Z)
1898
1899 // compute cell crossings in Z-direction
1900 const auto k_old = static_cast<int>(z_old-shift);
1901 const auto k_new = static_cast<int>(z_new-shift);
1902 const int cell_crossings_z = std::abs(k_new-k_old);
1903 num_segments += cell_crossings_z;
1904
1905 // need to assert that the number of cell crossings in each direction
1906 // is within the range permitted by the number of guard cells
1907 // e.g., if (num_segments > 3) ...
1908
1909 // compute the initial cell location used to find the cell crossings.
1910 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1911 double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1912
1913 // loop over the number of segments and deposit
1914 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1915 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1916 double dzp_seg;
1917 double z0_new;
1918 double z0_old = z_old;
1919
1920 for (int ns=0; ns<num_segments; ns++) {
1921
1922 if (ns == num_segments-1) { // final segment
1923 z0_new = z_new;
1924 dzp_seg = z0_new - z0_old;
1925 }
1926 else {
1927 Zcell = Zcell + dirZ_sign;
1928 z0_new = Zcell;
1929 dzp_seg = z0_new - z0_old;
1930 }
1931
1932 // compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1933 const auto seg_factor = static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1934
1935 // compute cell-based weights using the average segment position
1936 double sz_cell[depos_order] = {0.};
1937 double const z0_bar = (z0_new + z0_old)/2.0;
1938 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1939
1940 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1941 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1942 double sz_old_cell[depos_order] = {0.};
1943 double sz_new_cell[depos_order] = {0.};
1944 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1945 ignore_unused(k0_cell_2);
1946 for (int m=0; m<depos_order; m++) {
1947 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1948 }
1949 }
1950
1951 // compute node-based weights using the old and new segment positions
1952 double sz_old_node[depos_order+1] = {0.};
1953 double sz_new_node[depos_order+1] = {0.};
1954 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1955
1956 // gather out-of-plane Ex and Ey for this segment
1957 for (int k=0; k<=depos_order; k++) {
1958 auto weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1959 Exp += Ex_arr(lo.x+k0_node+k, 0, 0)*weight;
1960 Eyp += Ey_arr(lo.x+k0_node+k, 0, 0)*weight;
1961 }
1962
1963 // gather Ez for this segment
1964 for (int k=0; k<=depos_order-1; k++) {
1965 auto weight = sz_cell[k]*seg_factor;
1966 Ezp += Ez_arr(lo.x+k0_cell+k, 0, 0)*weight;
1967 }
1968
1969 // update old segment values
1970 if (ns < num_segments-1) {
1971 z0_old = z0_new;
1972 }
1973
1974 }
1975
1976#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1977
1978 amrex::ParticleReal Erp = 0.;
1979 amrex::ParticleReal Ethetap = 0.;
1980 // Z for CYLINDER, phi for SPHERE
1981 amrex::ParticleReal E3p = 0.;
1982
1983 // compute cell crossings in X-direction
1984 const auto i_old = static_cast<int>(x_old-shift);
1985 const auto i_new = static_cast<int>(x_new-shift);
1986 const int cell_crossings_x = std::abs(i_new-i_old);
1987 num_segments += cell_crossings_x;
1988
1989 // need to assert that the number of cell crossings in each direction
1990 // is within the range permitted by the number of guard cells
1991 // e.g., if (num_segments > 3) ...
1992
1993 // compute the initial cell location used to find the cell crossings.
1994 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1995 double Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1996
1997 // loop over the number of segments and deposit
1998 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1999 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
2000 double dxp_seg;
2001 double x0_new;
2002 double x0_old = x_old;
2003
2004 for (int ns=0; ns<num_segments; ns++) {
2005
2006 if (ns == num_segments-1) { // final segment
2007 x0_new = x_new;
2008 dxp_seg = x0_new - x0_old;
2009 }
2010 else {
2011 Xcell = Xcell + dirX_sign;
2012 x0_new = Xcell;
2013 dxp_seg = x0_new - x0_old;
2014 }
2015
2016 // compute the segment factor (equal to dt_seg/dt for nonzero dxp)
2017 const auto seg_factor = static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
2018
2019 // compute cell-based weights using the average segment position
2020 double sx_cell[depos_order] = {0.};
2021 double const x0_bar = (x0_new + x0_old)/2.0;
2022 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
2023
2024 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
2025 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
2026 double sx_old_cell[depos_order] = {0.};
2027 double sx_new_cell[depos_order] = {0.};
2028 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
2029 ignore_unused(i0_cell_2);
2030 for (int m=0; m<depos_order; m++) {
2031 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
2032 }
2033 }
2034
2035 // compute node-based weights using the old and new segment positions
2036 double sx_old_node[depos_order+1] = {0.};
2037 double sx_new_node[depos_order+1] = {0.};
2038 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
2039
2040 // gather out-of-plane Ey and Ez for this segment
2041 for (int i=0; i<=depos_order; i++) {
2042 auto weight = 0.5_rt*(sx_old_node[i] + sx_new_node[i])*seg_factor;
2043 Ethetap += Ey_arr(lo.x+i0_node+i, 0, 0)*weight;
2044 E3p += Ez_arr(lo.x+i0_node+i, 0, 0)*weight;
2045 }
2046
2047 // gather Ex for this segment
2048 for (int i=0; i<=depos_order-1; i++) {
2049 auto weight = sx_cell[i]*seg_factor;
2050 Erp += Ex_arr(lo.x+i0_cell+i, 0, 0)*weight;
2051 }
2052
2053 // update old segment values
2054 if (ns < num_segments-1) {
2055 x0_old = x0_new;
2056 }
2057
2058 }
2059
2060#if defined(WARPX_DIM_RCYLINDER)
2061 // Convert Erp and Ethetap to Ex and Ey
2062 Exp += costheta_mid*Erp - sintheta_mid*Ethetap;
2063 Eyp += costheta_mid*Ethetap + sintheta_mid*Erp;
2064 Ezp += E3p;
2065
2066#elif defined(WARPX_DIM_RSPHERE)
2067
2068 // Convert Erp, Ethetap, and Ephi to Ex, Ey, and Ez
2069 Exp += costheta_mid*cosphi_mid*Erp - sintheta_mid*Ethetap - costheta_mid*sinphi_mid*E3p;
2070 Eyp += sintheta_mid*cosphi_mid*Erp + costheta_mid*Ethetap - sintheta_mid*sinphi_mid*E3p;
2071 Ezp += sinphi_mid*Erp + cosphi_mid*E3p;
2072
2073#endif
2074#endif
2075
2076 // Gather magnetic field
2077 // This passes in the uncropped values of the half step position.
2078 // A future fix could pass in the average of the cropped values so
2079 // that the position would be within the domain.
2080 const int depos_order_perp = 1;
2081 const int depos_order_para = 1;
2082 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
2083 xp_nph, yp_nph, zp_nph,
2084 Bxp, Byp, Bzp,
2085 Bx_arr, By_arr, Bz_arr,
2086 Bx_type, By_type, Bz_type,
2087 dinv, xyzmin, lo, n_rz_azimuthal_modes );
2088
2089}
2090
2108AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2109void doGatherShapeN (const amrex::ParticleReal xp,
2110 const amrex::ParticleReal yp,
2111 const amrex::ParticleReal zp,
2112 amrex::ParticleReal& Exp,
2113 amrex::ParticleReal& Eyp,
2114 amrex::ParticleReal& Ezp,
2115 amrex::ParticleReal& Bxp,
2116 amrex::ParticleReal& Byp,
2117 amrex::ParticleReal& Bzp,
2118 amrex::Array4<amrex::Real const> const& ex_arr,
2119 amrex::Array4<amrex::Real const> const& ey_arr,
2120 amrex::Array4<amrex::Real const> const& ez_arr,
2121 amrex::Array4<amrex::Real const> const& bx_arr,
2122 amrex::Array4<amrex::Real const> const& by_arr,
2123 amrex::Array4<amrex::Real const> const& bz_arr,
2124 const amrex::IndexType ex_type,
2125 const amrex::IndexType ey_type,
2126 const amrex::IndexType ez_type,
2127 const amrex::IndexType bx_type,
2128 const amrex::IndexType by_type,
2129 const amrex::IndexType bz_type,
2130 const amrex::XDim3 & dinv,
2131 const amrex::XDim3 & xyzmin,
2132 const amrex::Dim3& lo,
2133 const int n_rz_azimuthal_modes,
2134 const int nox,
2135 const bool galerkin_interpolation)
2136{
2137 if (galerkin_interpolation) {
2138 if (nox == 1) {
2139 doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2140 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2141 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2142 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2143 } else if (nox == 2) {
2144 doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2145 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2146 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2147 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2148 } else if (nox == 3) {
2149 doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2150 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2151 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2152 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2153 } else if (nox == 4) {
2154 doGatherShapeN<4,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2155 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2156 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2157 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2158 }
2159 } else {
2160 if (nox == 1) {
2161 doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2162 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2163 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2164 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2165 } else if (nox == 2) {
2166 doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2167 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2168 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2169 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2170 } else if (nox == 3) {
2171 doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2172 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2173 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2174 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2175 } else if (nox == 4) {
2176 doGatherShapeN<4,0>(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 }
2181 }
2182}
2183
2184
2205AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2206void doGatherShapeNImplicit (
2207 const amrex::ParticleReal xp_n,
2208 const amrex::ParticleReal yp_n,
2209 const amrex::ParticleReal zp_n,
2210 const amrex::ParticleReal xp_nph,
2211 const amrex::ParticleReal yp_nph,
2212 const amrex::ParticleReal zp_nph,
2213 amrex::ParticleReal& Exp,
2214 amrex::ParticleReal& Eyp,
2215 amrex::ParticleReal& Ezp,
2216 amrex::ParticleReal& Bxp,
2217 amrex::ParticleReal& Byp,
2218 amrex::ParticleReal& Bzp,
2219 amrex::Array4<amrex::Real const> const& ex_arr,
2220 amrex::Array4<amrex::Real const> const& ey_arr,
2221 amrex::Array4<amrex::Real const> const& ez_arr,
2222 amrex::Array4<amrex::Real const> const& bx_arr,
2223 amrex::Array4<amrex::Real const> const& by_arr,
2224 amrex::Array4<amrex::Real const> const& bz_arr,
2225 const amrex::IndexType ex_type,
2226 const amrex::IndexType ey_type,
2227 const amrex::IndexType ez_type,
2228 const amrex::IndexType bx_type,
2229 const amrex::IndexType by_type,
2230 const amrex::IndexType bz_type,
2231 const amrex::XDim3 & dinv,
2232 const amrex::XDim3 & xyzmin,
2233 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
2234 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
2235 const amrex::Dim3& lo,
2236 const int n_rz_azimuthal_modes,
2237 const int nox,
2238 const CurrentDepositionAlgo depos_type )
2239{
2240 if (depos_type == CurrentDepositionAlgo::Esirkepov) {
2241 if (nox == 1) {
2242 doGatherShapeNEsirkepovStencilImplicit<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2243 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2244 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2245 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2246 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2247 } else if (nox == 2) {
2248 doGatherShapeNEsirkepovStencilImplicit<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2249 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2250 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2251 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2252 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2253 } else if (nox == 3) {
2254 doGatherShapeNEsirkepovStencilImplicit<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2255 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2256 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2257 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2258 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2259 } else if (nox == 4) {
2260 doGatherShapeNEsirkepovStencilImplicit<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2261 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2262 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2263 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2264 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2265 }
2266 }
2267 else if (depos_type == CurrentDepositionAlgo::Villasenor) {
2268 if (nox == 1) {
2269 doGatherPicnicShapeN<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2270 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2271 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2272 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2273 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2274 } else if (nox == 2) {
2275 doGatherPicnicShapeN<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2276 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2277 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2278 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2279 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2280 } else if (nox == 3) {
2281 doGatherPicnicShapeN<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2282 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2283 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2284 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2285 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2286 } else if (nox == 4) {
2287 doGatherPicnicShapeN<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
2288 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2289 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2290 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2291 dinv, xyzmin, domain_double, do_cropping, lo, n_rz_azimuthal_modes);
2292 }
2293 }
2294 else if (depos_type == CurrentDepositionAlgo::Direct) {
2295 if (nox == 1) {
2296 doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2297 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2298 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2299 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2300 } else if (nox == 2) {
2301 doGatherShapeN<2,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2302 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2303 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2304 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2305 } else if (nox == 3) {
2306 doGatherShapeN<3,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2307 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2308 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2309 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2310 } else if (nox == 4) {
2311 doGatherShapeN<4,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
2312 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
2313 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
2314 dinv, xyzmin, lo, n_rz_azimuthal_modes);
2315 }
2316 }
2317}
2318
2319#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:1388
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:2206
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
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
amrex::Array4
ArrayND< T, 4, true > Array4
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