8 #ifndef WARPX_FIELDGATHER_H_
9 #define WARPX_FIELDGATHER_H_
36 template <
int depos_order,
int galerkin_
interpolation>
39 const amrex::ParticleReal yp,
40 const amrex::ParticleReal zp,
41 amrex::ParticleReal& Exp,
42 amrex::ParticleReal& Eyp,
43 amrex::ParticleReal& Ezp,
44 amrex::ParticleReal& Bxp,
45 amrex::ParticleReal& Byp,
46 amrex::ParticleReal& Bzp,
62 const int n_rz_azimuthal_modes)
64 using namespace amrex;
66 #if defined(WARPX_DIM_XZ)
70 #if defined(WARPX_DIM_1D_Z)
78 #if (AMREX_SPACEDIM >= 2)
79 const amrex::Real dxi = 1.0_rt/
dx[0];
81 const amrex::Real dzi = 1.0_rt/
dx[2];
82 #if defined(WARPX_DIM_3D)
83 const amrex::Real dyi = 1.0_rt/
dx[1];
86 #if (AMREX_SPACEDIM >= 2)
87 const amrex::Real
xmin = xyzmin[0];
89 #if defined(WARPX_DIM_3D)
90 const amrex::Real ymin = xyzmin[1];
92 const amrex::Real zmin = xyzmin[2];
94 constexpr
int zdir = WARPX_ZINDEX;
101 Compute_shape_factor<depos_order - galerkin_interpolation >
const compute_shape_factor_galerkin;
103 #if (AMREX_SPACEDIM >= 2)
107 const amrex::Real rp = std::sqrt(xp*xp + yp*yp);
108 const amrex::Real
x = (rp -
xmin)*dxi;
110 const amrex::Real
x = (xp-
xmin)*dxi;
117 amrex::Real sx_node[depos_order + 1];
118 amrex::Real sx_cell[depos_order + 1];
119 amrex::Real sx_node_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt};
120 amrex::Real sx_cell_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt};
126 if ((ey_type[0] ==
NODE) || (ez_type[0] ==
NODE) || (bx_type[0] ==
NODE)) {
127 j_node = compute_shape_factor(sx_node,
x);
129 if ((ey_type[0] == CELL) || (ez_type[0] == CELL) || (bx_type[0] == CELL)) {
130 j_cell = compute_shape_factor(sx_cell,
x - 0.5_rt);
132 if ((ex_type[0] ==
NODE) || (by_type[0] ==
NODE) || (bz_type[0] ==
NODE)) {
133 j_node_v = compute_shape_factor_galerkin(sx_node_galerkin,
x);
135 if ((ex_type[0] == CELL) || (by_type[0] == CELL) || (bz_type[0] == CELL)) {
136 j_cell_v = compute_shape_factor_galerkin(sx_cell_galerkin,
x - 0.5_rt);
138 const amrex::Real (&sx_ex)[depos_order + 1 - galerkin_interpolation] = ((ex_type[0] ==
NODE) ? sx_node_galerkin : sx_cell_galerkin);
139 const amrex::Real (&sx_ey)[depos_order + 1 ] = ((ey_type[0] ==
NODE) ? sx_node : sx_cell );
140 const amrex::Real (&sx_ez)[depos_order + 1 ] = ((ez_type[0] ==
NODE) ? sx_node : sx_cell );
141 const amrex::Real (&sx_bx)[depos_order + 1 ] = ((bx_type[0] ==
NODE) ? sx_node : sx_cell );
142 const amrex::Real (&sx_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[0] ==
NODE) ? sx_node_galerkin : sx_cell_galerkin);
143 const amrex::Real (&sx_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[0] ==
NODE) ? sx_node_galerkin : sx_cell_galerkin);
144 int const j_ex = ((ex_type[0] ==
NODE) ? j_node_v : j_cell_v);
145 int const j_ey = ((ey_type[0] ==
NODE) ? j_node : j_cell );
146 int const j_ez = ((ez_type[0] ==
NODE) ? j_node : j_cell );
147 int const j_bx = ((bx_type[0] ==
NODE) ? j_node : j_cell );
148 int const j_by = ((by_type[0] ==
NODE) ? j_node_v : j_cell_v);
149 int const j_bz = ((bz_type[0] ==
NODE) ? j_node_v : j_cell_v);
152 #if defined(WARPX_DIM_3D)
154 const amrex::Real y = (yp-ymin)*dyi;
155 amrex::Real sy_node[depos_order + 1];
156 amrex::Real sy_cell[depos_order + 1];
157 amrex::Real sy_node_v[depos_order + 1 - galerkin_interpolation];
158 amrex::Real sy_cell_v[depos_order + 1 - galerkin_interpolation];
163 if ((ex_type[1] ==
NODE) || (ez_type[1] ==
NODE) || (by_type[1] ==
NODE)) {
164 k_node = compute_shape_factor(sy_node, y);
166 if ((ex_type[1] == CELL) || (ez_type[1] == CELL) || (by_type[1] == CELL)) {
167 k_cell = compute_shape_factor(sy_cell, y - 0.5_rt);
169 if ((ey_type[1] ==
NODE) || (bx_type[1] ==
NODE) || (bz_type[1] ==
NODE)) {
170 k_node_v = compute_shape_factor_galerkin(sy_node_v, y);
172 if ((ey_type[1] == CELL) || (bx_type[1] == CELL) || (bz_type[1] == CELL)) {
173 k_cell_v = compute_shape_factor_galerkin(sy_cell_v, y - 0.5_rt);
175 const amrex::Real (&sy_ex)[depos_order + 1 ] = ((ex_type[1] ==
NODE) ? sy_node : sy_cell );
176 const amrex::Real (&sy_ey)[depos_order + 1 - galerkin_interpolation] = ((ey_type[1] ==
NODE) ? sy_node_v : sy_cell_v);
177 const amrex::Real (&sy_ez)[depos_order + 1 ] = ((ez_type[1] ==
NODE) ? sy_node : sy_cell );
178 const amrex::Real (&sy_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[1] ==
NODE) ? sy_node_v : sy_cell_v);
179 const amrex::Real (&sy_by)[depos_order + 1 ] = ((by_type[1] ==
NODE) ? sy_node : sy_cell );
180 const amrex::Real (&sy_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[1] ==
NODE) ? sy_node_v : sy_cell_v);
181 int const k_ex = ((ex_type[1] ==
NODE) ? k_node : k_cell );
182 int const k_ey = ((ey_type[1] ==
NODE) ? k_node_v : k_cell_v);
183 int const k_ez = ((ez_type[1] ==
NODE) ? k_node : k_cell );
184 int const k_bx = ((bx_type[1] ==
NODE) ? k_node_v : k_cell_v);
185 int const k_by = ((by_type[1] ==
NODE) ? k_node : k_cell );
186 int const k_bz = ((bz_type[1] ==
NODE) ? k_node_v : k_cell_v);
190 const amrex::Real z = (zp-zmin)*dzi;
191 amrex::Real sz_node[depos_order + 1];
192 amrex::Real sz_cell[depos_order + 1];
193 amrex::Real sz_node_v[depos_order + 1 - galerkin_interpolation];
194 amrex::Real sz_cell_v[depos_order + 1 - galerkin_interpolation];
199 if ((ex_type[zdir] ==
NODE) || (ey_type[zdir] ==
NODE) || (bz_type[zdir] ==
NODE)) {
200 l_node = compute_shape_factor(sz_node, z);
202 if ((ex_type[zdir] == CELL) || (ey_type[zdir] == CELL) || (bz_type[zdir] == CELL)) {
203 l_cell = compute_shape_factor(sz_cell, z - 0.5_rt);
205 if ((ez_type[zdir] ==
NODE) || (bx_type[zdir] ==
NODE) || (by_type[zdir] ==
NODE)) {
206 l_node_v = compute_shape_factor_galerkin(sz_node_v, z);
208 if ((ez_type[zdir] == CELL) || (bx_type[zdir] == CELL) || (by_type[zdir] == CELL)) {
209 l_cell_v = compute_shape_factor_galerkin(sz_cell_v, z - 0.5_rt);
211 const amrex::Real (&sz_ex)[depos_order + 1 ] = ((ex_type[zdir] ==
NODE) ? sz_node : sz_cell );
212 const amrex::Real (&sz_ey)[depos_order + 1 ] = ((ey_type[zdir] ==
NODE) ? sz_node : sz_cell );
213 const amrex::Real (&sz_ez)[depos_order + 1 - galerkin_interpolation] = ((ez_type[zdir] ==
NODE) ? sz_node_v : sz_cell_v);
214 const amrex::Real (&sz_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[zdir] ==
NODE) ? sz_node_v : sz_cell_v);
215 const amrex::Real (&sz_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[zdir] ==
NODE) ? sz_node_v : sz_cell_v);
216 const amrex::Real (&sz_bz)[depos_order + 1 ] = ((bz_type[zdir] ==
NODE) ? sz_node : sz_cell );
217 int const l_ex = ((ex_type[zdir] ==
NODE) ? l_node : l_cell );
218 int const l_ey = ((ey_type[zdir] ==
NODE) ? l_node : l_cell );
219 int const l_ez = ((ez_type[zdir] ==
NODE) ? l_node_v : l_cell_v);
220 int const l_bx = ((bx_type[zdir] ==
NODE) ? l_node_v : l_cell_v);
221 int const l_by = ((by_type[zdir] ==
NODE) ? l_node_v : l_cell_v);
222 int const l_bz = ((bz_type[zdir] ==
NODE) ? l_node : l_cell );
230 #if defined(WARPX_DIM_1D_Z)
234 for (
int iz=0;
iz<=depos_order;
iz++){
236 ey_arr(lo.
x+l_ey+
iz, 0, 0, 0);
238 ex_arr(lo.
x+l_ex+
iz, 0, 0, 0);
240 bz_arr(lo.
x+l_bz+
iz, 0, 0, 0);
246 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
248 ez_arr(lo.
x+l_ez+
iz, 0, 0, 0);
250 bx_arr(lo.
x+l_bx+
iz, 0, 0, 0);
252 by_arr(lo.
x+l_by+
iz, 0, 0, 0);
255 #elif defined(WARPX_DIM_XZ)
257 for (
int iz=0;
iz<=depos_order;
iz++){
258 for (
int ix=0;
ix<=depos_order;
ix++){
259 Eyp += sx_ey[
ix]*sz_ey[
iz]*
260 ey_arr(lo.
x+j_ey+
ix, lo.
y+l_ey+
iz, 0, 0);
265 for (
int iz=0;
iz<=depos_order;
iz++){
266 for (
int ix=0;
ix<=depos_order-galerkin_interpolation;
ix++){
267 Exp += sx_ex[
ix]*sz_ex[
iz]*
268 ex_arr(lo.
x+j_ex+
ix, lo.
y+l_ex+
iz, 0, 0);
269 Bzp += sx_bz[
ix]*sz_bz[
iz]*
270 bz_arr(lo.
x+j_bz+
ix, lo.
y+l_bz+
iz, 0, 0);
275 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
276 for (
int ix=0;
ix<=depos_order;
ix++){
277 Ezp += sx_ez[
ix]*sz_ez[
iz]*
278 ez_arr(lo.
x+j_ez+
ix, lo.
y+l_ez+
iz, 0, 0);
279 Bxp += sx_bx[
ix]*sz_bx[
iz]*
280 bx_arr(lo.
x+j_bx+
ix, lo.
y+l_bx+
iz, 0, 0);
284 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
285 for (
int ix=0;
ix<=depos_order-galerkin_interpolation;
ix++){
286 Byp += sx_by[
ix]*sz_by[
iz]*
287 by_arr(lo.
x+j_by+
ix, lo.
y+l_by+
iz, 0, 0);
291 #elif defined(WARPX_DIM_RZ)
293 amrex::ParticleReal Erp = 0.;
294 amrex::ParticleReal Ethetap = 0.;
295 amrex::ParticleReal Brp = 0.;
296 amrex::ParticleReal Bthetap = 0.;
299 for (
int iz=0;
iz<=depos_order;
iz++){
300 for (
int ix=0;
ix<=depos_order;
ix++){
301 Ethetap += sx_ey[
ix]*sz_ey[
iz]*
302 ey_arr(lo.
x+j_ey+
ix, lo.
y+l_ey+
iz, 0, 0);
307 for (
int iz=0;
iz<=depos_order;
iz++){
308 for (
int ix=0;
ix<=depos_order-galerkin_interpolation;
ix++){
309 Erp += sx_ex[
ix]*sz_ex[
iz]*
310 ex_arr(lo.
x+j_ex+
ix, lo.
y+l_ex+
iz, 0, 0);
311 Bzp += sx_bz[
ix]*sz_bz[
iz]*
312 bz_arr(lo.
x+j_bz+
ix, lo.
y+l_bz+
iz, 0, 0);
317 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
318 for (
int ix=0;
ix<=depos_order;
ix++){
319 Ezp += sx_ez[
ix]*sz_ez[
iz]*
320 ez_arr(lo.
x+j_ez+
ix, lo.
y+l_ez+
iz, 0, 0);
321 Brp += sx_bx[
ix]*sz_bx[
iz]*
322 bx_arr(lo.
x+j_bx+
ix, lo.
y+l_bx+
iz, 0, 0);
326 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
327 for (
int ix=0;
ix<=depos_order-galerkin_interpolation;
ix++){
328 Bthetap += sx_by[
ix]*sz_by[
iz]*
329 by_arr(lo.
x+j_by+
ix, lo.
y+l_by+
iz, 0, 0);
333 amrex::Real costheta;
334 amrex::Real sintheta;
345 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
348 for (
int iz=0;
iz<=depos_order;
iz++){
349 for (
int ix=0;
ix<=depos_order;
ix++){
350 const amrex::Real dEy = (+ ey_arr(lo.
x+j_ey+
ix, lo.
y+l_ey+
iz, 0, 2*imode-1)*xy.
real()
351 - ey_arr(lo.
x+j_ey+
ix, lo.
y+l_ey+
iz, 0, 2*imode)*xy.
imag());
352 Ethetap += sx_ey[
ix]*sz_ey[
iz]*dEy;
357 for (
int iz=0;
iz<=depos_order;
iz++){
358 for (
int ix=0;
ix<=depos_order-galerkin_interpolation;
ix++){
359 const amrex::Real dEx = (+ ex_arr(lo.
x+j_ex+
ix, lo.
y+l_ex+
iz, 0, 2*imode-1)*xy.
real()
360 - ex_arr(lo.
x+j_ex+
ix, lo.
y+l_ex+
iz, 0, 2*imode)*xy.
imag());
361 Erp += sx_ex[
ix]*sz_ex[
iz]*dEx;
362 const amrex::Real dBz = (+ bz_arr(lo.
x+j_bz+
ix, lo.
y+l_bz+
iz, 0, 2*imode-1)*xy.
real()
363 - bz_arr(lo.
x+j_bz+
ix, lo.
y+l_bz+
iz, 0, 2*imode)*xy.
imag());
364 Bzp += sx_bz[
ix]*sz_bz[
iz]*dBz;
369 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
370 for (
int ix=0;
ix<=depos_order;
ix++){
371 const amrex::Real dEz = (+ ez_arr(lo.
x+j_ez+
ix, lo.
y+l_ez+
iz, 0, 2*imode-1)*xy.
real()
372 - ez_arr(lo.
x+j_ez+
ix, lo.
y+l_ez+
iz, 0, 2*imode)*xy.
imag());
373 Ezp += sx_ez[
ix]*sz_ez[
iz]*dEz;
374 const amrex::Real dBx = (+ bx_arr(lo.
x+j_bx+
ix, lo.
y+l_bx+
iz, 0, 2*imode-1)*xy.
real()
375 - bx_arr(lo.
x+j_bx+
ix, lo.
y+l_bx+
iz, 0, 2*imode)*xy.
imag());
376 Brp += sx_bx[
ix]*sz_bx[
iz]*dBx;
380 for (
int iz=0;
iz<=depos_order-galerkin_interpolation;
iz++){
381 for (
int ix=0;
ix<=depos_order-galerkin_interpolation;
ix++){
382 const amrex::Real dBy = (+ by_arr(lo.
x+j_by+
ix, lo.
y+l_by+
iz, 0, 2*imode-1)*xy.
real()
383 - by_arr(lo.
x+j_by+
ix, lo.
y+l_by+
iz, 0, 2*imode)*xy.
imag());
384 Bthetap += sx_by[
ix]*sz_by[
iz]*dBy;
391 Exp += costheta*Erp - sintheta*Ethetap;
392 Eyp += costheta*Ethetap + sintheta*Erp;
393 Bxp += costheta*Brp - sintheta*Bthetap;
394 Byp += costheta*Bthetap + sintheta*Brp;
398 for (
int iz=0;
iz<=depos_order;
iz++){
399 for (
int iy=0;
iy<=depos_order;
iy++){
400 for (
int ix=0;
ix<= depos_order - galerkin_interpolation;
ix++){
401 Exp += sx_ex[
ix]*sy_ex[
iy]*sz_ex[
iz]*
402 ex_arr(lo.
x+j_ex+
ix, lo.
y+k_ex+
iy, lo.
z+l_ex+
iz);
407 for (
int iz=0;
iz<=depos_order;
iz++){
408 for (
int iy=0;
iy<= depos_order - galerkin_interpolation;
iy++){
409 for (
int ix=0;
ix<=depos_order;
ix++){
410 Eyp += sx_ey[
ix]*sy_ey[
iy]*sz_ey[
iz]*
411 ey_arr(lo.
x+j_ey+
ix, lo.
y+k_ey+
iy, lo.
z+l_ey+
iz);
416 for (
int iz=0;
iz<= depos_order - galerkin_interpolation;
iz++){
417 for (
int iy=0;
iy<=depos_order;
iy++){
418 for (
int ix=0;
ix<=depos_order;
ix++){
419 Ezp += sx_ez[
ix]*sy_ez[
iy]*sz_ez[
iz]*
420 ez_arr(lo.
x+j_ez+
ix, lo.
y+k_ez+
iy, lo.
z+l_ez+
iz);
425 for (
int iz=0;
iz<=depos_order;
iz++){
426 for (
int iy=0;
iy<= depos_order - galerkin_interpolation;
iy++){
427 for (
int ix=0;
ix<= depos_order - galerkin_interpolation;
ix++){
428 Bzp += sx_bz[
ix]*sy_bz[
iy]*sz_bz[
iz]*
429 bz_arr(lo.
x+j_bz+
ix, lo.
y+k_bz+
iy, lo.
z+l_bz+
iz);
434 for (
int iz=0;
iz<= depos_order - galerkin_interpolation;
iz++){
435 for (
int iy=0;
iy<=depos_order;
iy++){
436 for (
int ix=0;
ix<= depos_order - galerkin_interpolation;
ix++){
437 Byp += sx_by[
ix]*sy_by[
iy]*sz_by[
iz]*
438 by_arr(lo.
x+j_by+
ix, lo.
y+k_by+
iy, lo.
z+l_by+
iz);
443 for (
int iz=0;
iz<= depos_order - galerkin_interpolation;
iz++){
444 for (
int iy=0;
iy<= depos_order - galerkin_interpolation;
iy++){
445 for (
int ix=0;
ix<=depos_order;
ix++){
446 Bxp += sx_bx[
ix]*sy_bx[
iy]*sz_bx[
iz]*
447 bx_arr(lo.
x+j_bx+
ix, lo.
y+k_bx+
iy, lo.
z+l_bx+
iz);
472 template <
int depos_order>
475 [[maybe_unused]]
const amrex::ParticleReal xp_n,
476 [[maybe_unused]]
const amrex::ParticleReal yp_n,
477 const amrex::ParticleReal zp_n,
478 [[maybe_unused]]
const amrex::ParticleReal xp_nph,
479 [[maybe_unused]]
const amrex::ParticleReal yp_nph,
480 const amrex::ParticleReal zp_nph,
481 amrex::ParticleReal& Exp,
482 amrex::ParticleReal& Eyp,
483 amrex::ParticleReal& Ezp,
484 amrex::ParticleReal& Bxp,
485 amrex::ParticleReal& Byp,
486 amrex::ParticleReal& Bzp,
502 const int n_rz_azimuthal_modes)
504 using namespace amrex;
505 #if !defined(WARPX_DIM_RZ)
509 #if !defined(WARPX_DIM_1D_Z)
510 Real
const dxi = 1.0_rt /
dx[0];
512 #if !defined(WARPX_DIM_1D_Z)
513 Real
const xmin = xyzmin[0];
515 #if defined(WARPX_DIM_3D)
516 Real
const dyi = 1.0_rt /
dx[1];
517 Real
const ymin = xyzmin[1];
519 Real
const dzi = 1.0_rt /
dx[2];
520 Real
const zmin = xyzmin[2];
522 #if !defined(WARPX_DIM_1D_Z)
523 Real constexpr one_third = 1.0_rt / 3.0_rt;
524 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
527 #if !defined(WARPX_DIM_1D_Z)
528 const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
530 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
531 const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
533 const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
536 #if defined(WARPX_DIM_RZ)
537 amrex::Real
const xp_new = xp_np1;
538 amrex::Real
const yp_new = yp_np1;
539 amrex::Real
const xp_mid = xp_nph;
540 amrex::Real
const yp_mid = yp_nph;
541 amrex::Real
const xp_old = xp_n;
542 amrex::Real
const yp_old = yp_n;
543 amrex::Real
const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
544 amrex::Real
const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
545 amrex::Real
const rp_mid = (rp_new + rp_old)/2._rt;
546 amrex::Real costheta_mid, sintheta_mid;
547 if (rp_mid > 0._rt) {
548 costheta_mid = xp_mid/rp_mid;
549 sintheta_mid = yp_mid/rp_mid;
551 costheta_mid = 1._rt;
552 sintheta_mid = 0._rt;
556 double const x_new = (rp_new -
xmin)*dxi;
557 double const x_old = (rp_old -
xmin)*dxi;
559 #if !defined(WARPX_DIM_1D_Z)
561 double const x_new = (xp_np1 -
xmin)*dxi;
562 double const x_old = (xp_n -
xmin)*dxi;
565 #if defined(WARPX_DIM_3D)
567 double const y_new = (yp_np1 - ymin)*dyi;
568 double const y_old = (yp_n - ymin)*dyi;
571 double const z_new = (zp_np1 - zmin)*dzi;
572 double const z_old = (zp_n - zmin)*dzi;
579 #if !defined(WARPX_DIM_1D_Z)
580 double sx_E_new[depos_order + 3] = {0.};
581 double sx_E_old[depos_order + 3] = {0.};
583 #if defined(WARPX_DIM_3D)
585 double sy_E_new[depos_order + 3] = {0.};
586 double sy_E_old[depos_order + 3] = {0.};
589 double sz_E_new[depos_order + 3] = {0.};
590 double sz_E_old[depos_order + 3] = {0.};
592 #if defined(WARPX_DIM_3D)
593 double sx_B_new[depos_order + 3] = {0.};
594 double sx_B_old[depos_order + 3] = {0.};
595 double sy_B_new[depos_order + 3] = {0.};
596 double sy_B_old[depos_order + 3] = {0.};
597 double sz_B_new[depos_order + 3] = {0.};
598 double sz_B_old[depos_order + 3] = {0.};
601 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
604 double sx_By_new[depos_order + 2] = {0.};
605 double sx_By_old[depos_order + 2] = {0.};
606 double sz_By_new[depos_order + 2] = {0.};
607 double sz_By_old[depos_order + 2] = {0.};
616 #if !defined(WARPX_DIM_1D_Z)
617 const int i_E_new = compute_shape_factor(sx_E_new+1, x_new);
618 const int i_E_old = compute_shifted_shape_factor(sx_E_old, x_old, i_E_new);
620 #if defined(WARPX_DIM_3D)
621 const int j_E_new = compute_shape_factor(sy_E_new+1, y_new);
622 const int j_E_old = compute_shifted_shape_factor(sy_E_old, y_old, j_E_new);
624 const int k_E_new = compute_shape_factor(sz_E_new+1, z_new);
625 const int k_E_old = compute_shifted_shape_factor(sz_E_old, z_old, k_E_new);
627 #if defined(WARPX_DIM_3D)
628 const int i_B_new = compute_shape_factor(sx_B_new+1, x_new - 0.5_rt);
629 const int i_B_old = compute_shifted_shape_factor(sx_B_old, x_old - 0.5_rt, i_B_new);
630 const int j_B_new = compute_shape_factor(sy_B_new+1, y_new - 0.5_rt);
631 const int j_B_old = compute_shifted_shape_factor(sy_B_old, y_old - 0.5_rt, j_B_new);
632 const int k_B_new = compute_shape_factor(sz_B_new+1, z_new - 0.5_rt);
633 const int k_B_old = compute_shifted_shape_factor(sz_B_old, z_old - 0.5_rt, k_B_new);
637 #if !defined(WARPX_DIM_1D_Z)
638 int dil_E = 1, diu_E = 1;
639 if (i_E_old < i_E_new) { dil_E = 0; }
640 if (i_E_old > i_E_new) { diu_E = 0; }
642 #if defined(WARPX_DIM_3D)
643 int djl_E = 1, dju_E = 1;
644 if (j_E_old < j_E_new) { djl_E = 0; }
645 if (j_E_old > j_E_new) { dju_E = 0; }
647 int dkl_E = 1, dku_E = 1;
648 if (k_E_old < k_E_new) { dkl_E = 0; }
649 if (k_E_old > k_E_new) { dku_E = 0; }
651 #if defined(WARPX_DIM_3D)
652 int dil_B = 1, diu_B = 1;
653 if (i_B_old < i_B_new) { dil_B = 0; }
654 if (i_B_old > i_B_new) { diu_B = 0; }
655 int djl_B = 1, dju_B = 1;
656 if (j_B_old < j_B_new) { djl_B = 0; }
657 if (j_B_old > j_B_new) { dju_B = 0; }
658 int dkl_B = 1, dku_B = 1;
659 if (k_B_old < k_B_new) { dkl_B = 0; }
660 if (k_B_old > k_B_new) { dku_B = 0; }
663 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
666 const int i_By_new = compute_shape_factor_By(sx_By_new+1, x_new - 0.5_rt);
667 const int i_By_old = compute_shifted_shape_factor_By(sx_By_old, x_old - 0.5_rt, i_By_new);
668 const int k_By_new = compute_shape_factor_By(sz_By_new+1, z_new - 0.5_rt);
669 const int k_By_old = compute_shifted_shape_factor_By(sz_By_old, z_old - 0.5_rt, k_By_new);
670 int dil_By = 1, diu_By = 1;
671 if (i_By_old < i_By_new) { dil_By = 0; }
672 if (i_By_old > i_By_new) { diu_By = 0; }
673 int dkl_By = 1, dku_By = 1;
674 if (k_By_old < k_By_new) { dkl_By = 0; }
675 if (k_By_old > k_By_new) { dku_By = 0; }
678 #if defined(WARPX_DIM_3D)
680 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
681 for (
int j=djl_E; j<=depos_order+2-dju_E; j++) {
682 const amrex::Real sdzjk = one_third*(sy_E_new[j]*sz_E_new[k] + sy_E_old[j]*sz_E_old[k])
683 +one_sixth*(sy_E_new[j]*sz_E_old[k] + sy_E_old[j]*sz_E_new[k]);
684 amrex::Real sdxi = 0._rt;
685 for (
int i=dil_E;
i<=depos_order+1-diu_E;
i++) {
686 sdxi += (sx_E_old[
i] - sx_E_new[
i]);
687 auto sdxiov =
static_cast<amrex::Real
>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
688 Exp += Ex_arr(lo.
x+i_E_new-1+
i, lo.
y+j_E_new-1+j, lo.
z+k_E_new-1+k)*sdxiov*sdzjk;
692 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
693 for (
int i=dil_E;
i<=depos_order+2-diu_E;
i++) {
694 const amrex::Real sdyik = one_third*(sx_E_new[
i]*sz_E_new[k] + sx_E_old[
i]*sz_E_old[k])
695 +one_sixth*(sx_E_new[
i]*sz_E_old[k] + sx_E_old[
i]*sz_E_new[k]);
696 amrex::Real sdyj = 0._rt;
697 for (
int j=djl_E; j<=depos_order+1-dju_E; j++) {
698 sdyj += (sy_E_old[j] - sy_E_new[j]);
699 auto sdyjov =
static_cast<amrex::Real
>((y_new - y_old) == 0. ? 1. : sdyj/(y_new - y_old));
700 Eyp += Ey_arr(lo.
x+i_E_new-1+
i, lo.
y+j_E_new-1+j, lo.
z+k_E_new-1+k)*sdyjov*sdyik;
704 for (
int j=djl_E; j<=depos_order+2-dju_E; j++) {
705 for (
int i=dil_E;
i<=depos_order+2-diu_E;
i++) {
706 const amrex::Real sdzij = one_third*(sx_E_new[
i]*sy_E_new[j] + sx_E_old[
i]*sy_E_old[j])
707 +one_sixth*(sx_E_new[
i]*sy_E_old[j] + sx_E_old[
i]*sy_E_new[j]);
708 amrex::Real sdzk = 0._rt;
709 for (
int k=dkl_E; k<=depos_order+1-dku_E; k++) {
710 sdzk += (sz_E_old[k] - sz_E_new[k]);
711 auto sdzkov =
static_cast<amrex::Real
>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
712 Ezp += Ez_arr(lo.
x+i_E_new-1+
i, lo.
y+j_E_new-1+j, lo.
z+k_E_new-1+k)*sdzkov*sdzij;
716 for (
int k=dkl_B; k<=depos_order+2-dku_B; k++) {
717 for (
int j=djl_B; j<=depos_order+2-dju_B; j++) {
718 const amrex::Real sdzjk = one_third*(sy_B_new[j]*sz_B_new[k] + sy_B_old[j]*sz_B_old[k])
719 +one_sixth*(sy_B_new[j]*sz_B_old[k] + sy_B_old[j]*sz_B_new[k]);
720 amrex::Real sdxi = 0._rt;
721 for (
int i=dil_B;
i<=depos_order+1-diu_B;
i++) {
722 sdxi += (sx_B_old[
i] - sx_B_new[
i]);
723 auto sdxiov =
static_cast<amrex::Real
>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
724 Bxp += Bx_arr(lo.
x+i_B_new-1+
i, lo.
y+j_B_new-1+j, lo.
z+k_B_new-1+k)*sdxiov*sdzjk;
728 for (
int k=dkl_B; k<=depos_order+2-dku_B; k++) {
729 for (
int i=dil_B;
i<=depos_order+2-diu_B;
i++) {
730 const amrex::Real sdyik = one_third*(sx_B_new[
i]*sz_B_new[k] + sx_B_old[
i]*sz_B_old[k])
731 +one_sixth*(sx_B_new[
i]*sz_B_old[k] + sx_B_old[
i]*sz_B_new[k]);
732 amrex::Real sdyj = 0._rt;
733 for (
int j=djl_B; j<=depos_order+1-dju_B; j++) {
734 sdyj += (sy_B_old[j] - sy_B_new[j]);
735 auto sdyjov =
static_cast<amrex::Real
>((y_new - y_old) == 0. ? 1. : sdyj/(y_new - y_old));
736 Byp += By_arr(lo.
x+i_B_new-1+
i, lo.
y+j_B_new-1+j, lo.
z+k_B_new-1+k)*sdyjov*sdyik;
740 for (
int j=djl_B; j<=depos_order+2-dju_B; j++) {
741 for (
int i=dil_B;
i<=depos_order+2-diu_B;
i++) {
742 const amrex::Real sdzij = one_third*(sx_B_new[
i]*sy_B_new[j] + sx_B_old[
i]*sy_B_old[j])
743 +one_sixth*(sx_B_new[
i]*sy_B_old[j] + sx_B_old[
i]*sy_B_new[j]);
744 amrex::Real sdzk = 0._rt;
745 for (
int k=dkl_B; k<=depos_order+1-dku_B; k++) {
746 sdzk += (sz_B_old[k] - sz_B_new[k]);
747 auto sdzkov =
static_cast<amrex::Real
>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
748 Bzp += Bz_arr(lo.
x+i_B_new-1+
i, lo.
y+j_B_new-1+j, lo.
z+k_E_new-1+k)*sdzkov*sdzij;
753 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
755 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
756 const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
757 amrex::Real sdxi = 0._rt;
758 for (
int i=dil_E;
i<=depos_order+1-diu_E;
i++) {
759 sdxi += (sx_E_old[
i] - sx_E_new[
i]);
760 auto sdxiov =
static_cast<amrex::Real
>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
761 Exp += Ex_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
762 Bzp += Bz_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 0)*sdxiov*sdzk;
765 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
766 for (
int i=dil_E;
i<=depos_order+2-diu_E;
i++) {
768 one_third*(sx_E_new[
i]*sz_E_new[k] + sx_E_old[
i]*sz_E_old[k])
769 +one_sixth*(sx_E_new[
i]*sz_E_old[k] + sx_E_old[
i]*sz_E_new[k]));
770 Eyp += Ey_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 0)*sdyj;
773 for (
int i=dil_E;
i<=depos_order+2-diu_E;
i++) {
774 const amrex::Real sdxi = 0.5_rt*(sx_E_new[
i] + sx_E_old[
i]);
775 amrex::Real sdzk = 0._rt;
776 for (
int k=dkl_E; k<=depos_order+1-dku_E; k++) {
777 sdzk += (sz_E_old[k] - sz_E_new[k]);
778 auto sdzkov =
static_cast<amrex::Real
>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
779 Ezp += Ez_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
780 Bxp += Bx_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 0)*sdzkov*sdxi;
783 for (
int k=dkl_By; k<=depos_order+1-dku_By; k++) {
784 for (
int i=dil_By;
i<=depos_order+1-diu_By;
i++) {
786 one_third*(sx_By_new[
i]*sz_By_new[k] + sx_By_old[
i]*sz_By_old[k])
787 +one_sixth*(sx_By_new[
i]*sz_By_old[k] + sx_By_old[
i]*sz_By_new[k]));
788 Byp += By_arr(lo.
x+i_By_new-1+
i, lo.
y+k_By_new-1+k, 0, 0)*sdyj;
795 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
799 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
800 const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]);
801 amrex::Real sdxi = 0._rt;
802 for (
int i=dil_E;
i<=depos_order+1-diu_E;
i++) {
803 sdxi += (sx_E_old[
i] - sx_E_new[
i]);
804 auto sdxiov =
static_cast<amrex::Real
>((x_new - x_old) == 0. ? 1. : sdxi/(x_new - x_old));
805 const amrex::Real dEx = (+ Ex_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.
real()
806 - Ex_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode)*xy_mid.
imag());
807 const amrex::Real dBz = (+ Bz_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.
real()
808 - Bz_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode)*xy_mid.
imag());
809 Exp += dEx*sdxiov*sdzk;
810 Bzp += dBz*sdxiov*sdzk;
814 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
815 for (
int i=dil_E;
i<=depos_order+2-diu_E;
i++) {
817 one_third*(sx_E_new[
i]*sz_E_new[k] + sx_E_old[
i]*sz_E_old[k])
818 +one_sixth*(sx_E_new[
i]*sz_E_old[k] + sx_E_old[
i]*sz_E_new[k]));
819 const amrex::Real dEy = (+ Ey_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.
real()
820 - Ey_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode)*xy_mid.
imag());
826 for (
int i=dil_E;
i<=depos_order+2-diu_E;
i++) {
827 const amrex::Real sdxi = 0.5_rt*(sx_E_new[
i] + sx_E_old[
i]);
828 amrex::Real sdzk = 0._rt;
829 for (
int k=dkl_E; k<=depos_order+1-dku_E; k++) {
830 sdzk += (sz_E_old[k] - sz_E_new[k]);
831 auto sdzkov =
static_cast<amrex::Real
>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
832 const amrex::Real dEz = (+ Ez_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.
real()
833 - Ez_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode)*xy_mid.
imag());
834 const amrex::Real dBx = (+ Bx_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode-1)*xy_mid.
real()
835 - Bx_arr(lo.
x+i_E_new-1+
i, lo.
y+k_E_new-1+k, 0, 2*imode)*xy_mid.
imag());
836 Ezp += dEz*sdzkov*sdxi;
837 Bxp += dBx*sdzkov*sdxi;
841 for (
int k=dkl_By; k<=depos_order+1-dku_By; k++) {
842 for (
int i=dil_By;
i<=depos_order+1-diu_By;
i++) {
844 one_third*(sx_By_new[
i]*sz_By_new[k] + sx_By_old[
i]*sz_By_old[k])
845 +one_sixth*(sx_By_new[
i]*sz_By_old[k] + sx_By_old[
i]*sz_By_new[k]));
846 const amrex::Real dBy = (+ By_arr(lo.
x+i_By_new-1+
i, lo.
y+k_By_new-1+k, 0, 2*imode-1)*xy_mid.
real()
847 - By_arr(lo.
x+i_By_new-1+
i, lo.
y+k_By_new-1+k, 0, 2*imode)*xy_mid.
imag());
851 xy_mid = xy_mid*xy_mid0;
855 const amrex::Real Exp_save = Exp;
856 Exp = costheta_mid*Exp - sintheta_mid*Eyp;
857 Eyp = costheta_mid*Eyp + sintheta_mid*Exp_save;
858 const amrex::Real Bxp_save = Bxp;
859 Bxp = costheta_mid*Bxp - sintheta_mid*Byp;
860 Byp = costheta_mid*Byp + sintheta_mid*Bxp_save;
864 #elif defined(WARPX_DIM_1D_Z)
866 for (
int k=dkl_E; k<=depos_order+2-dku_E; k++) {
867 amrex::Real
const sdzk = 0.5_rt*(sz_E_old[k] + sz_E_new[k]);
868 Exp += Ex_arr(lo.
x+k_E_new-1+k, 0, 0, 0)*sdzk;
869 Eyp += Ey_arr(lo.
x+k_E_new-1+k, 0, 0, 0)*sdzk;
870 Bzp += Bz_arr(lo.
x+k_E_new-1+k, 0, 0, 0)*sdzk;
872 amrex::Real sdzk = 0._rt;
873 for (
int k=dkl_E; k<=depos_order+1-dku_E; k++) {
874 sdzk += (sz_E_old[k] - sz_E_new[k]);
875 auto sdzkov =
static_cast<amrex::Real
>((z_new - z_old) == 0. ? 1. : sdzk/(z_new - z_old));
876 Bxp += Bx_arr(lo.
x+k_E_new-1+k, 0, 0, 0)*sdzkov;
877 Byp += By_arr(lo.
x+k_E_new-1+k, 0, 0, 0)*sdzkov;
878 Ezp += Ez_arr(lo.
x+k_E_new-1+k, 0, 0, 0)*sdzkov;
902 template <
int depos_order>
905 [[maybe_unused]]
const amrex::ParticleReal xp_n,
906 [[maybe_unused]]
const amrex::ParticleReal yp_n,
907 const amrex::ParticleReal zp_n,
908 [[maybe_unused]]
const amrex::ParticleReal xp_nph,
909 [[maybe_unused]]
const amrex::ParticleReal yp_nph,
910 const amrex::ParticleReal zp_nph,
911 amrex::ParticleReal& Exp,
912 amrex::ParticleReal& Eyp,
913 amrex::ParticleReal& Ezp,
914 amrex::ParticleReal& Bxp,
915 amrex::ParticleReal& Byp,
916 amrex::ParticleReal& Bzp,
932 const int n_rz_azimuthal_modes)
934 using namespace amrex;
935 #if !defined(WARPX_DIM_RZ)
939 #if !defined(WARPX_DIM_1D_Z)
940 Real
const dxi = 1.0_rt /
dx[0];
942 #if !defined(WARPX_DIM_1D_Z)
943 Real
const xmin = xyzmin[0];
945 #if defined(WARPX_DIM_3D)
946 Real
const dyi = 1.0_rt /
dx[1];
947 Real
const ymin = xyzmin[1];
949 Real
const dzi = 1.0_rt /
dx[2];
950 Real
const zmin = xyzmin[2];
952 #if !defined(WARPX_DIM_1D_Z)
953 const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n;
955 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
956 const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n;
958 const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n;
960 #if !defined(WARPX_DIM_1D_Z)
961 Real constexpr one_third = 1.0_rt / 3.0_rt;
962 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
966 #if defined(WARPX_DIM_RZ)
967 amrex::Real
const xp_new = xp_np1;
968 amrex::Real
const yp_new = yp_np1;
969 amrex::Real
const xp_mid = xp_nph;
970 amrex::Real
const yp_mid = yp_nph;
971 amrex::Real
const xp_old = xp_n;
972 amrex::Real
const yp_old = yp_n;
973 amrex::Real
const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
974 amrex::Real
const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
975 amrex::Real
const rp_mid = (rp_new + rp_old)/2._rt;
976 amrex::Real costheta_mid, sintheta_mid;
977 if (rp_mid > 0._rt) {
978 costheta_mid = xp_mid/rp_mid;
979 sintheta_mid = yp_mid/rp_mid;
981 costheta_mid = 1._rt;
982 sintheta_mid = 0._rt;
986 double const x_new = (rp_new -
xmin)*dxi;
987 double const x_old = (rp_old -
xmin)*dxi;
988 double const x_bar = (rp_mid -
xmin)*dxi;
989 #elif !defined(WARPX_DIM_1D_Z)
991 double const x_new = (xp_np1 -
xmin)*dxi;
992 double const x_old = (xp_n -
xmin)*dxi;
993 double const x_bar = (xp_nph -
xmin)*dxi;
995 #if defined(WARPX_DIM_3D)
997 double const y_new = (yp_np1 - ymin)*dyi;
998 double const y_old = (yp_n - ymin)*dyi;
999 double const y_bar = (yp_nph - ymin)*dyi;
1002 double const z_new = (zp_np1 - zmin)*dzi;
1003 double const z_old = (zp_n - zmin)*dzi;
1004 double const z_bar = (zp_nph - zmin)*dzi;
1013 int num_segments = 1;
1015 if ( (depos_order % 2) == 0 ) {
shift = 0.5; }
1017 #if defined(WARPX_DIM_3D)
1020 const auto i_old =
static_cast<int>(x_old-
shift);
1021 const auto i_new =
static_cast<int>(x_new-
shift);
1022 const int cell_crossings_x = std::abs(i_new-i_old);
1023 num_segments += cell_crossings_x;
1026 const auto j_old =
static_cast<int>(y_old-
shift);
1027 const auto j_new =
static_cast<int>(y_new-
shift);
1028 const int cell_crossings_y = std::abs(j_new-j_old);
1029 num_segments += cell_crossings_y;
1032 const auto k_old =
static_cast<int>(z_old-
shift);
1033 const auto k_new =
static_cast<int>(z_new-
shift);
1034 const int cell_crossings_z = std::abs(k_new-k_old);
1035 num_segments += cell_crossings_z;
1043 const double dxp = x_new - x_old;
1044 const double dyp = y_new - y_old;
1045 const double dzp = z_new - z_old;
1046 const auto dirX_sign =
static_cast<double>(dxp < 0. ? -1. : 1.);
1047 const auto dirY_sign =
static_cast<double>(dyp < 0. ? -1. : 1.);
1048 const auto dirZ_sign =
static_cast<double>(dzp < 0. ? -1. : 1.);
1049 double Xcell = 0., Ycell = 0., Zcell = 0.;
1050 if (num_segments > 1) {
1051 Xcell =
static_cast<double>(i_old) +
shift + 0.5*(1.-dirX_sign);
1052 Ycell =
static_cast<double>(j_old) +
shift + 0.5*(1.-dirY_sign);
1053 Zcell =
static_cast<double>(k_old) +
shift + 0.5*(1.-dirZ_sign);
1059 double dxp_seg, dyp_seg, dzp_seg;
1060 double x0_new, y0_new, z0_new;
1061 double x0_old = x_old;
1062 double y0_old = y_old;
1063 double z0_old = z_old;
1065 for (
int ns=0; ns<num_segments; ns++) {
1067 if (ns == num_segments-1) {
1072 dxp_seg = x0_new - x0_old;
1073 dyp_seg = y0_new - y0_old;
1074 dzp_seg = z0_new - z0_old;
1079 x0_new = Xcell + dirX_sign;
1080 y0_new = Ycell + dirY_sign;
1081 z0_new = Zcell + dirZ_sign;
1082 dxp_seg = x0_new - x0_old;
1083 dyp_seg = y0_new - y0_old;
1084 dzp_seg = z0_new - z0_old;
1086 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
1087 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
1089 dyp_seg = dyp/dxp*dxp_seg;
1090 dzp_seg = dzp/dxp*dxp_seg;
1091 y0_new = y0_old + dyp_seg;
1092 z0_new = z0_old + dzp_seg;
1094 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
1096 dxp_seg = dxp/dyp*dyp_seg;
1097 dzp_seg = dzp/dyp*dyp_seg;
1098 x0_new = x0_old + dxp_seg;
1099 z0_new = z0_old + dzp_seg;
1103 dxp_seg = dxp/dzp*dzp_seg;
1104 dyp_seg = dyp/dzp*dzp_seg;
1105 x0_new = x0_old + dxp_seg;
1106 y0_new = y0_old + dyp_seg;
1112 const auto seg_factor_x =
static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1113 const auto seg_factor_y =
static_cast<double>(dyp == 0. ? 1. : dyp_seg/dyp);
1114 const auto seg_factor_z =
static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1117 double sx_cell[depos_order] = {0.};
1118 double sy_cell[depos_order] = {0.};
1119 double sz_cell[depos_order] = {0.};
1120 double const x0_bar = (x0_new + x0_old)/2.0;
1121 double const y0_bar = (y0_new + y0_old)/2.0;
1122 double const z0_bar = (z0_new + z0_old)/2.0;
1123 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1124 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
1125 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1127 if constexpr (depos_order >= 3) {
1129 double sx_old_cell[depos_order] = {0.};
1130 double sx_new_cell[depos_order] = {0.};
1131 double sy_old_cell[depos_order] = {0.};
1132 double sy_new_cell[depos_order] = {0.};
1133 double sz_old_cell[depos_order] = {0.};
1134 double sz_new_cell[depos_order] = {0.};
1135 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1136 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
1137 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1139 for (
int m=0; m<depos_order; m++) {
1140 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1141 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
1142 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1147 double sx_old_node[depos_order+1] = {0.};
1148 double sx_new_node[depos_order+1] = {0.};
1149 double sy_old_node[depos_order+1] = {0.};
1150 double sy_new_node[depos_order+1] = {0.};
1151 double sz_old_node[depos_order+1] = {0.};
1152 double sz_new_node[depos_order+1] = {0.};
1153 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1154 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
1155 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1159 for (
int i=0;
i<=depos_order-1;
i++) {
1160 for (
int j=0; j<=depos_order; j++) {
1161 for (
int k=0; k<=depos_order; k++) {
1162 weight = sx_cell[
i]*( sy_old_node[j]*sz_old_node[k]*one_third
1163 + sy_old_node[j]*sz_new_node[k]*one_sixth
1164 + sy_new_node[j]*sz_old_node[k]*one_sixth
1165 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
1166 Exp += Ex_arr(lo.
x+i0_cell+
i, lo.
y+j0_node+j, lo.
z+k0_node+k)*weight;
1172 for (
int i=0;
i<=depos_order;
i++) {
1173 for (
int j=0; j<=depos_order-1; j++) {
1174 for (
int k=0; k<=depos_order; k++) {
1175 weight = sy_cell[j]*( sx_old_node[
i]*sz_old_node[k]*one_third
1176 + sx_old_node[
i]*sz_new_node[k]*one_sixth
1177 + sx_new_node[
i]*sz_old_node[k]*one_sixth
1178 + sx_new_node[
i]*sz_new_node[k]*one_third )*seg_factor_y;
1179 Eyp += Ey_arr(lo.
x+i0_node+
i, lo.
y+j0_cell+j, lo.
z+k0_node+k)*weight;
1185 for (
int i=0;
i<=depos_order;
i++) {
1186 for (
int j=0; j<=depos_order; j++) {
1187 for (
int k=0; k<=depos_order-1; k++) {
1188 weight = sz_cell[k]*( sx_old_node[
i]*sy_old_node[j]*one_third
1189 + sx_old_node[
i]*sy_new_node[j]*one_sixth
1190 + sx_new_node[
i]*sy_old_node[j]*one_sixth
1191 + sx_new_node[
i]*sy_new_node[j]*one_third )*seg_factor_z;
1192 Ezp += Ez_arr(lo.
x+i0_node+
i, lo.
y+j0_node+j, lo.
z+k0_cell+k)*weight;
1198 if (ns < num_segments-1) {
1207 const int depos_order_B = 1;
1209 double sz_bar_node[depos_order_B+1] = {0.};
1210 double sz_bar_cell[depos_order_B+1] = {0.};
1211 const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar);
1212 const int k_bar_cell = compute_shape_factor_B(sz_bar_cell, z_bar-0.5);
1213 double sy_bar_node[depos_order_B+1] = {0.};
1214 double sy_bar_cell[depos_order_B+1] = {0.};
1215 const int j_bar_node = compute_shape_factor_B(sy_bar_node, y_bar);
1216 const int j_bar_cell = compute_shape_factor_B(sy_bar_cell, y_bar-0.5);
1217 double sx_bar_node[depos_order_B+1] = {0.};
1218 double sx_bar_cell[depos_order_B+1] = {0.};
1219 const int i_bar_node = compute_shape_factor_B(sx_bar_node, x_bar);
1220 const int i_bar_cell = compute_shape_factor_B(sx_bar_cell, x_bar-0.5);
1223 for (
int i=0;
i<=depos_order_B;
i++) {
1224 for (
int j=0; j<=depos_order_B; j++) {
1225 for (
int k=0; k<=depos_order_B; k++) {
1226 weight =
static_cast<amrex::Real
>(sx_bar_node[
i]*sy_bar_cell[j]*sz_bar_cell[k]);
1227 Bxp += Bx_arr(lo.
x+i_bar_node+
i, lo.
y+j_bar_cell+j, lo.
z+k_bar_cell+k)*weight;
1229 weight =
static_cast<amrex::Real
>(sx_bar_cell[
i]*sy_bar_node[j]*sz_bar_cell[k]);
1230 Byp += By_arr(lo.
x+i_bar_cell+
i, lo.
y+j_bar_node+j, lo.
z+k_bar_cell+k)*weight;
1232 weight =
static_cast<amrex::Real
>(sx_bar_cell[
i]*sy_bar_cell[j]*sz_bar_node[k]);
1233 Bzp += Bz_arr(lo.
x+i_bar_cell+
i, lo.
y+j_bar_cell+j, lo.
z+k_bar_node+k)*weight;
1238 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1241 const auto i_old =
static_cast<int>(x_old-
shift);
1242 const auto i_new =
static_cast<int>(x_new-
shift);
1243 const int cell_crossings_x = std::abs(i_new-i_old);
1244 num_segments += cell_crossings_x;
1247 const auto k_old =
static_cast<int>(z_old-
shift);
1248 const auto k_new =
static_cast<int>(z_new-
shift);
1249 const int cell_crossings_z = std::abs(k_new-k_old);
1250 num_segments += cell_crossings_z;
1258 const double dxp = x_new - x_old;
1259 const double dzp = z_new - z_old;
1260 const auto dirX_sign =
static_cast<double>(dxp < 0. ? -1. : 1.);
1261 const auto dirZ_sign =
static_cast<double>(dzp < 0. ? -1. : 1.);
1262 double Xcell = 0., Zcell = 0.;
1263 if (num_segments > 1) {
1264 Xcell =
static_cast<double>(i_old) +
shift + 0.5*(1.-dirX_sign);
1265 Zcell =
static_cast<double>(k_old) +
shift + 0.5*(1.-dirZ_sign);
1271 double dxp_seg, dzp_seg;
1272 double x0_new, z0_new;
1273 double x0_old = x_old;
1274 double z0_old = z_old;
1276 for (
int ns=0; ns<num_segments; ns++) {
1278 if (ns == num_segments-1) {
1282 dxp_seg = x0_new - x0_old;
1283 dzp_seg = z0_new - z0_old;
1288 x0_new = Xcell + dirX_sign;
1289 z0_new = Zcell + dirZ_sign;
1290 dxp_seg = x0_new - x0_old;
1291 dzp_seg = z0_new - z0_old;
1293 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1295 dzp_seg = dzp/dxp*dxp_seg;
1296 z0_new = z0_old + dzp_seg;
1300 dxp_seg = dxp/dzp*dzp_seg;
1301 x0_new = x0_old + dxp_seg;
1307 const auto seg_factor_x =
static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1308 const auto seg_factor_z =
static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1311 double sx_cell[depos_order] = {0.};
1312 double sz_cell[depos_order] = {0.};
1313 double const x0_bar = (x0_new + x0_old)/2.0;
1314 double const z0_bar = (z0_new + z0_old)/2.0;
1315 const int i0_cell = compute_shape_factor_cell(sx_cell, x0_bar-0.5);
1316 const int k0_cell = compute_shape_factor_cell(sz_cell, z0_bar-0.5);
1318 if constexpr (depos_order >= 3) {
1320 double sx_old_cell[depos_order] = {0.};
1321 double sx_new_cell[depos_order] = {0.};
1322 double sz_old_cell[depos_order] = {0.};
1323 double sz_new_cell[depos_order] = {0.};
1324 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1325 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1327 for (
int m=0; m<depos_order; m++) {
1328 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1329 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1334 double sx_old_node[depos_order+1] = {0.};
1335 double sx_new_node[depos_order+1] = {0.};
1336 double sz_old_node[depos_order+1] = {0.};
1337 double sz_new_node[depos_order+1] = {0.};
1338 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1339 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1343 for (
int i=0;
i<=depos_order-1;
i++) {
1344 for (
int k=0; k<=depos_order; k++) {
1345 weight = sx_cell[
i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1346 Exp += Ex_arr(lo.
x+i0_cell+
i, lo.
y+k0_node+k, 0, 0)*weight;
1347 #if defined(WARPX_DIM_RZ)
1349 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1350 const auto dEx = (+ Ex_arr(lo.
x+i0_cell+
i, lo.
y+k0_node+k, 0, 2*imode-1)*xy_mid.
real()
1351 - Ex_arr(lo.
x+i0_cell+
i, lo.
y+k0_node+k, 0, 2*imode)*xy_mid.
imag());
1353 xy_mid = xy_mid*xy_mid0;
1360 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1361 for (
int i=0;
i<=depos_order;
i++) {
1362 for (
int k=0; k<=depos_order; k++) {
1363 weight = ( sx_old_node[
i]*sz_old_node[k]*one_third
1364 + sx_old_node[
i]*sz_new_node[k]*one_sixth
1365 + sx_new_node[
i]*sz_old_node[k]*one_sixth
1366 + sx_new_node[
i]*sz_new_node[k]*one_third )*seg_factor_y;
1367 Eyp += Ey_arr(lo.
x+i0_node+
i, lo.
y+k0_node+k, 0, 0)*weight;
1368 #if defined(WARPX_DIM_RZ)
1370 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1371 const auto dEy = (+ Ey_arr(lo.
x+i0_node+
i, lo.
y+k0_node+k, 0, 2*imode-1)*xy_mid.
real()
1372 - Ey_arr(lo.
x+i0_node+
i, lo.
y+k0_node+k, 0, 2*imode)*xy_mid.
imag());
1374 xy_mid = xy_mid*xy_mid0;
1381 for (
int i=0;
i<=depos_order;
i++) {
1382 for (
int k=0; k<=depos_order-1; k++) {
1383 weight = sz_cell[k]*(sx_old_node[
i] + sx_new_node[
i])/2.0_rt*seg_factor_z;
1384 Ezp += Ez_arr(lo.
x+i0_node+
i, lo.
y+k0_cell+k, 0, 0)*weight;
1385 #if defined(WARPX_DIM_RZ)
1387 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1388 const auto dEz = (+ Ez_arr(lo.
x+i0_node+
i, lo.
y+k0_cell+k, 0, 2*imode-1)*xy_mid.
real()
1389 - Ez_arr(lo.
x+i0_node+
i, lo.
y+k0_cell+k, 0, 2*imode)*xy_mid.
imag());
1391 xy_mid = xy_mid*xy_mid0;
1398 if (ns < num_segments-1) {
1406 const int depos_order_B = 1;
1408 double sz_bar_node[depos_order_B+1] = {0.};
1409 double sz_bar_cell[depos_order_B+1] = {0.};
1410 const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar);
1411 const int k_bar_cell = compute_shape_factor_B(sz_bar_cell, z_bar-0.5);
1412 double sx_bar_node[depos_order_B+1] = {0.};
1413 double sx_bar_cell[depos_order_B+1] = {0.};
1414 const int i_bar_node = compute_shape_factor_B(sx_bar_node, x_bar);
1415 const int i_bar_cell = compute_shape_factor_B(sx_bar_cell, x_bar-0.5);
1417 for (
int i=0;
i<=depos_order_B;
i++) {
1418 for (
int k=0; k<=depos_order_B; k++) {
1419 const auto weight_Bz =
static_cast<amrex::Real
>(sx_bar_cell[
i]*sz_bar_node[k]);
1420 Bzp += Bz_arr(lo.
x+i_bar_cell+
i, lo.
y+k_bar_node+k, 0, 0)*weight_Bz;
1422 const auto weight_Bx =
static_cast<amrex::Real
>(sx_bar_node[
i]*sz_bar_cell[k]);
1423 Bxp += Bx_arr(lo.
x+i_bar_node+
i, lo.
y+k_bar_cell+k, 0, 0)*weight_Bx;
1425 const auto weight_By =
static_cast<amrex::Real
>(sx_bar_cell[
i]*sz_bar_cell[k]);
1426 Byp += By_arr(lo.
x+i_bar_cell+
i, lo.
y+k_bar_cell+k, 0, 0)*weight_By;
1427 #if defined(WARPX_DIM_RZ)
1429 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1430 const auto dBx = (+ Bx_arr(lo.
x+i_bar_node+
i, lo.
y+k_bar_cell+k, 0, 2*imode-1)*xy_mid.
real()
1431 - Bx_arr(lo.
x+i_bar_node+
i, lo.
y+k_bar_cell+k, 0, 2*imode)*xy_mid.
imag());
1432 const auto dBy = (+ By_arr(lo.
x+i_bar_cell+
i, lo.
y+k_bar_cell+k, 0, 2*imode-1)*xy_mid.
real()
1433 - By_arr(lo.
x+i_bar_cell+
i, lo.
y+k_bar_cell+k, 0, 2*imode)*xy_mid.
imag());
1434 const auto dBz = (+ Bz_arr(lo.
x+i_bar_cell+
i, lo.
y+k_bar_node+k, 0, 2*imode-1)*xy_mid.
real()
1435 - Bz_arr(lo.
x+i_bar_cell+
i, lo.
y+k_bar_node+k, 0, 2*imode)*xy_mid.
imag());
1436 Bxp += weight_Bx*dBx;
1437 Byp += weight_By*dBy;
1438 Bzp += weight_Bz*dBz;
1439 xy_mid = xy_mid*xy_mid0;
1448 const amrex::Real Exp_save = Exp;
1449 Exp = costheta_mid*Exp - sintheta_mid*Eyp;
1450 Eyp = costheta_mid*Eyp + sintheta_mid*Exp_save;
1451 const amrex::Real Bxp_save = Bxp;
1452 Bxp = costheta_mid*Bxp - sintheta_mid*Byp;
1453 Byp = costheta_mid*Byp + sintheta_mid*Bxp_save;
1457 #elif defined(WARPX_DIM_1D_Z)
1460 const auto k_old =
static_cast<int>(z_old-
shift);
1461 const auto k_new =
static_cast<int>(z_new-
shift);
1462 const int cell_crossings_z = std::abs(k_new-k_old);
1463 num_segments += cell_crossings_z;
1470 double const dzp = z_new - z_old;
1471 const auto dirZ_sign =
static_cast<double>(dzp < 0. ? -1. : 1.);
1472 double Zcell =
static_cast<double>(k_old) +
shift + 0.5*(1.-dirZ_sign);
1479 double z0_old = z_old;
1481 for (
int ns=0; ns<num_segments; ns++) {
1483 if (ns == num_segments-1) {
1485 dzp_seg = z0_new - z0_old;
1488 Zcell = Zcell + dirZ_sign;
1490 dzp_seg = z0_new - z0_old;
1494 const auto seg_factor =
static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1497 double sz_cell[depos_order] = {0.};
1498 double const z0_bar = (z0_new + z0_old)/2.0;
1499 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1501 if constexpr (depos_order >= 3) {
1503 double sz_old_cell[depos_order] = {0.};
1504 double sz_new_cell[depos_order] = {0.};
1505 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1507 for (
int m=0; m<depos_order; m++) {
1508 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1513 double sz_old_node[depos_order+1] = {0.};
1514 double sz_new_node[depos_order+1] = {0.};
1515 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1518 for (
int k=0; k<=depos_order; k++) {
1519 auto weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1520 Exp += Ex_arr(lo.
x+k0_node+k, 0, 0)*weight;
1521 Eyp += Ey_arr(lo.
x+k0_node+k, 0, 0)*weight;
1525 for (
int k=0; k<=depos_order-1; k++) {
1526 auto weight = sz_cell[k]*seg_factor;
1527 Ezp += Ez_arr(lo.
x+k0_cell+k, 0, 0)*weight;
1531 if (ns < num_segments-1) {
1538 const int depos_order_B = 1;
1540 double sz_bar_node[depos_order_B+1] = {0.};
1541 double sz_bar_cell[depos_order_B+1] = {0.};
1542 const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar);
1543 const int k_bar_cell = compute_shape_factor_B(sz_bar_cell, z_bar-0.5_rt);
1546 for (
int k=0; k<=depos_order_B; k++) {
1547 weight =
static_cast<amrex::Real
>(sz_bar_node[k]);
1548 Bzp += Bz_arr(lo.
x+k_bar_node+k, 0, 0)*weight;
1550 weight =
static_cast<amrex::Real
>(sz_bar_cell[k]);
1551 Bxp += Bx_arr(lo.
x+k_bar_cell+k, 0, 0)*weight;
1552 Byp += By_arr(lo.
x+k_bar_cell+k, 0, 0)*weight;
1575 template <
int depos_order,
int lower_in_v>
1578 amrex::ParticleReal *
const Exp, amrex::ParticleReal *
const Eyp,
1579 amrex::ParticleReal *
const Ezp, amrex::ParticleReal *
const Bxp,
1580 amrex::ParticleReal *
const Byp, amrex::ParticleReal *
const Bzp,
1587 const long np_to_gather,
1588 const std::array<amrex::Real, 3>&
dx,
1589 const std::array<amrex::Real, 3> xyzmin,
1591 const int n_rz_azimuthal_modes)
1617 amrex::ParticleReal xp, yp, zp;
1618 getPosition(ip, xp, yp, zp);
1619 getExternalEB(ip, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip]);
1621 doGatherShapeN<depos_order, lower_in_v>(
1622 xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip],
1623 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1624 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1625 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1649 const amrex::ParticleReal yp,
1650 const amrex::ParticleReal zp,
1651 amrex::ParticleReal& Exp,
1652 amrex::ParticleReal& Eyp,
1653 amrex::ParticleReal& Ezp,
1654 amrex::ParticleReal& Bxp,
1655 amrex::ParticleReal& Byp,
1656 amrex::ParticleReal& Bzp,
1672 const int n_rz_azimuthal_modes,
1674 const bool galerkin_interpolation)
1676 if (galerkin_interpolation) {
1678 doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1679 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1680 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1681 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1682 }
else if (nox == 2) {
1683 doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1684 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1685 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1686 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1687 }
else if (nox == 3) {
1688 doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1689 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1690 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1691 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1692 }
else if (nox == 4) {
1693 doGatherShapeN<4,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1694 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1695 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1696 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1700 doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1701 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1702 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1703 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1704 }
else if (nox == 2) {
1705 doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1706 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1707 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1708 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1709 }
else if (nox == 3) {
1710 doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1711 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1712 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1713 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1714 }
else if (nox == 4) {
1715 doGatherShapeN<4,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1716 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1717 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1718 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1745 const amrex::ParticleReal xp_n,
1746 const amrex::ParticleReal yp_n,
1747 const amrex::ParticleReal zp_n,
1748 const amrex::ParticleReal xp_nph,
1749 const amrex::ParticleReal yp_nph,
1750 const amrex::ParticleReal zp_nph,
1751 amrex::ParticleReal& Exp,
1752 amrex::ParticleReal& Eyp,
1753 amrex::ParticleReal& Ezp,
1754 amrex::ParticleReal& Bxp,
1755 amrex::ParticleReal& Byp,
1756 amrex::ParticleReal& Bzp,
1772 const int n_rz_azimuthal_modes,
1774 const int depos_type )
1776 if (depos_type==0) {
1778 doGatherShapeNEsirkepovStencilImplicit<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1779 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1780 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1781 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1782 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1783 }
else if (nox == 2) {
1784 doGatherShapeNEsirkepovStencilImplicit<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1785 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1786 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1787 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1788 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1789 }
else if (nox == 3) {
1790 doGatherShapeNEsirkepovStencilImplicit<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1791 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1792 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1793 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1794 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1795 }
else if (nox == 4) {
1796 doGatherShapeNEsirkepovStencilImplicit<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1797 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1798 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1799 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1800 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1803 else if (depos_type==3) {
1805 doGatherPicnicShapeN<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1806 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1807 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1808 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1809 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1810 }
else if (nox == 2) {
1811 doGatherPicnicShapeN<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1812 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1813 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1814 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1815 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1816 }
else if (nox == 3) {
1817 doGatherPicnicShapeN<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1818 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1819 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1820 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1821 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1822 }
else if (nox == 4) {
1823 doGatherPicnicShapeN<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
1824 Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1825 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1826 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1827 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1830 else if (depos_type==1) {
1832 doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1833 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1834 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1835 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1836 }
else if (nox == 2) {
1837 doGatherShapeN<2,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1838 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1839 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1840 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1841 }
else if (nox == 3) {
1842 doGatherShapeN<3,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1843 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1844 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1845 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
1846 }
else if (nox == 4) {
1847 doGatherShapeN<4,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
1848 ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
1849 ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
1850 dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeN(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &ex_arr, amrex::Array4< amrex::Real const > const &ey_arr, amrex::Array4< amrex::Real const > const &ez_arr, amrex::Array4< amrex::Real const > const &bx_arr, amrex::Array4< amrex::Real const > const &by_arr, amrex::Array4< amrex::Real const > const &bz_arr, const amrex::IndexType ex_type, const amrex::IndexType ey_type, const amrex::IndexType ez_type, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx, const amrex::GpuArray< amrex::Real, 3 > &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Field gather for a single particle.
Definition: FieldGather.H:38
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeNImplicit(const amrex::ParticleReal xp_n, const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, const amrex::ParticleReal xp_nph, const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &ex_arr, amrex::Array4< amrex::Real const > const &ey_arr, amrex::Array4< amrex::Real const > const &ez_arr, amrex::Array4< amrex::Real const > const &bx_arr, amrex::Array4< amrex::Real const > const &by_arr, amrex::Array4< amrex::Real const > const &bz_arr, const amrex::IndexType ex_type, const amrex::IndexType ey_type, const amrex::IndexType ez_type, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx_arr, const amrex::GpuArray< amrex::Real, 3 > &xyzmin_arr, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes, const int nox, const int depos_type)
Field gather for a single particle.
Definition: FieldGather.H:1744
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherPicnicShapeN([[maybe_unused]] const amrex::ParticleReal xp_n, [[maybe_unused]] const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, [[maybe_unused]] const amrex::ParticleReal xp_nph, [[maybe_unused]] const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &Ex_arr, amrex::Array4< amrex::Real const > const &Ey_arr, amrex::Array4< amrex::Real const > const &Ez_arr, amrex::Array4< amrex::Real const > const &Bx_arr, amrex::Array4< amrex::Real const > const &By_arr, amrex::Array4< amrex::Real const > const &Bz_arr, [[maybe_unused]] const amrex::IndexType Ex_type, [[maybe_unused]] const amrex::IndexType Ey_type, [[maybe_unused]] const amrex::IndexType Ez_type, [[maybe_unused]] const amrex::IndexType Bx_type, [[maybe_unused]] const amrex::IndexType By_type, [[maybe_unused]] const amrex::IndexType Bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx, const amrex::GpuArray< amrex::Real, 3 > &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Energy conserving field gather for thread thread_num for the implicit scheme This uses the same stenc...
Definition: FieldGather.H:904
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherShapeNEsirkepovStencilImplicit([[maybe_unused]] const amrex::ParticleReal xp_n, [[maybe_unused]] const amrex::ParticleReal yp_n, const amrex::ParticleReal zp_n, [[maybe_unused]] const amrex::ParticleReal xp_nph, [[maybe_unused]] const amrex::ParticleReal yp_nph, const amrex::ParticleReal zp_nph, amrex::ParticleReal &Exp, amrex::ParticleReal &Eyp, amrex::ParticleReal &Ezp, amrex::ParticleReal &Bxp, amrex::ParticleReal &Byp, amrex::ParticleReal &Bzp, amrex::Array4< amrex::Real const > const &Ex_arr, amrex::Array4< amrex::Real const > const &Ey_arr, amrex::Array4< amrex::Real const > const &Ez_arr, amrex::Array4< amrex::Real const > const &Bx_arr, amrex::Array4< amrex::Real const > const &By_arr, amrex::Array4< amrex::Real const > const &Bz_arr, [[maybe_unused]] const amrex::IndexType Ex_type, [[maybe_unused]] const amrex::IndexType Ey_type, [[maybe_unused]] const amrex::IndexType Ez_type, [[maybe_unused]] const amrex::IndexType Bx_type, [[maybe_unused]] const amrex::IndexType By_type, [[maybe_unused]] const amrex::IndexType Bz_type, const amrex::GpuArray< amrex::Real, 3 > &dx, const amrex::GpuArray< amrex::Real, 3 > &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Energy conserving field gather for thread thread_num for the implicit scheme This uses the same stenc...
Definition: FieldGather.H:474
AMREX_FORCE_INLINE Array4< Real const > array() const noexcept
const Box & box() const noexcept
AMREX_GPU_HOST_DEVICE IndexType ixType() const noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box shift(const Box &b, int dir, int nzones) noexcept
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
i
Definition: check_interp_points_and_weights.py:174
tuple dx
lab frame
Definition: stencil.py:429
xmin
Definition: stencil.py:424
Definition: ShapeFactors.H:168
Definition: ShapeFactors.H:29
Definition: ShapeFactors.H:95
Functor class that assigns external field values (E and B) to particles.
Definition: GetExternalFields.H:25
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T real() const noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T imag() const noexcept