8 #ifndef CURRENTDEPOSITION_H_ 9 #define CURRENTDEPOSITION_H_ 16 #include <AMReX_Array4.H> 17 #include <AMReX_REAL.H> 40 template <
int depos_order>
42 const amrex::ParticleReal *
const wp,
43 const amrex::ParticleReal *
const uxp,
44 const amrex::ParticleReal *
const uyp,
45 const amrex::ParticleReal *
const uzp,
46 const int *
const ion_lev,
47 amrex::FArrayBox& jx_fab,
48 amrex::FArrayBox& jy_fab,
49 amrex::FArrayBox& jz_fab,
50 const long np_to_depose,
const amrex::Real dt,
51 const std::array<amrex::Real,3>&
dx,
52 const std::array<amrex::Real,3>& xyzmin,
55 const long n_rz_azimuthal_modes)
57 #if !defined(WARPX_DIM_RZ) 58 amrex::ignore_unused(n_rz_azimuthal_modes);
63 const bool do_ionization = ion_lev;
64 const amrex::Real dxi = 1.0/dx[0];
65 const amrex::Real dzi = 1.0/dx[2];
66 #if !(defined WARPX_DIM_RZ) 67 const amrex::Real dts2dx = 0.5*dt*dxi;
69 const amrex::Real dts2dz = 0.5*dt*dzi;
70 #if (AMREX_SPACEDIM == 2) 71 const amrex::Real invvol = dxi*dzi;
72 #elif (defined WARPX_DIM_3D) 73 const amrex::Real dyi = 1.0/dx[1];
74 const amrex::Real dts2dy = 0.5*dt*dyi;
75 const amrex::Real invvol = dxi*dyi*dzi;
78 const amrex::Real xmin = xyzmin[0];
79 #if (defined WARPX_DIM_3D) 80 const amrex::Real ymin = xyzmin[1];
82 const amrex::Real zmin = xyzmin[2];
84 const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c;
86 amrex::Array4<amrex::Real>
const& jx_arr = jx_fab.array();
87 amrex::Array4<amrex::Real>
const& jy_arr = jy_fab.array();
88 amrex::Array4<amrex::Real>
const& jz_arr = jz_fab.array();
89 amrex::IntVect
const jx_type = jx_fab.box().type();
90 amrex::IntVect
const jy_type = jy_fab.box().type();
91 amrex::IntVect
const jz_type = jz_fab.box().type();
93 constexpr
int zdir = (AMREX_SPACEDIM - 1);
95 constexpr
int CELL = amrex::IndexType::CELL;
100 [=] AMREX_GPU_DEVICE (
long ip) {
102 const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
103 + uyp[ip]*uyp[ip]*clightsq
104 + uzp[ip]*uzp[ip]*clightsq);
105 amrex::Real wq = q*wp[ip];
110 amrex::ParticleReal xp, yp, zp;
111 GetPosition(ip, xp, yp, zp);
113 const amrex::Real vx = uxp[ip]*gaminv;
114 const amrex::Real vy = uyp[ip]*gaminv;
115 const amrex::Real vz = uzp[ip]*gaminv;
117 #if (defined WARPX_DIM_RZ) 120 const amrex::Real xpmid = xp - 0.5*dt*vx;
121 const amrex::Real ypmid = yp - 0.5*dt*vy;
122 const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
123 amrex::Real costheta;
124 amrex::Real sintheta;
126 costheta = xpmid/rpmid;
127 sintheta = ypmid/rpmid;
133 const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
134 const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
136 const amrex::Real wqx = wq*invvol*vx;
137 const amrex::Real wqy = wq*invvol*vy;
139 const amrex::Real wqz = wq*invvol*vz;
144 #if (defined WARPX_DIM_RZ) 146 const double xmid = (rpmid - xmin)*dxi;
148 const double xmid = (xp - xmin)*dxi - dts2dx*vx;
155 double sx_node[depos_order + 1];
156 double sx_cell[depos_order + 1];
160 if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) {
161 j_node = compute_shape_factor(sx_node, xmid);
163 if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
164 j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
167 amrex::Real sx_jx[depos_order + 1] = {0.};
168 amrex::Real sx_jy[depos_order + 1] = {0.};
169 amrex::Real sx_jz[depos_order + 1] = {0.};
170 for (
int ix=0; ix<=depos_order; ix++)
172 sx_jx[ix] = ((jx_type[0] ==
NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
173 sx_jy[ix] = ((jy_type[0] ==
NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
174 sx_jz[ix] = ((jz_type[0] ==
NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
177 int const j_jx = ((jx_type[0] ==
NODE) ? j_node : j_cell);
178 int const j_jy = ((jy_type[0] ==
NODE) ? j_node : j_cell);
179 int const j_jz = ((jz_type[0] ==
NODE) ? j_node : j_cell);
181 #if (defined WARPX_DIM_3D) 184 const double ymid = (yp - ymin)*dyi - dts2dy*vy;
185 double sy_node[depos_order + 1];
186 double sy_cell[depos_order + 1];
189 if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) {
190 k_node = compute_shape_factor(sy_node, ymid);
192 if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
193 k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
195 amrex::Real sy_jx[depos_order + 1] = {0.};
196 amrex::Real sy_jy[depos_order + 1] = {0.};
197 amrex::Real sy_jz[depos_order + 1] = {0.};
198 for (
int iy=0; iy<=depos_order; iy++)
200 sy_jx[iy] = ((jx_type[1] ==
NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
201 sy_jy[iy] = ((jy_type[1] ==
NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
202 sy_jz[iy] = ((jz_type[1] ==
NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
204 int const k_jx = ((jx_type[1] ==
NODE) ? k_node : k_cell);
205 int const k_jy = ((jy_type[1] ==
NODE) ? k_node : k_cell);
206 int const k_jz = ((jz_type[1] ==
NODE) ? k_node : k_cell);
211 const double zmid = (zp - zmin)*dzi - dts2dz*vz;
212 double sz_node[depos_order + 1];
213 double sz_cell[depos_order + 1];
216 if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) {
217 l_node = compute_shape_factor(sz_node, zmid);
219 if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
220 l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
222 amrex::Real sz_jx[depos_order + 1] = {0.};
223 amrex::Real sz_jy[depos_order + 1] = {0.};
224 amrex::Real sz_jz[depos_order + 1] = {0.};
225 for (
int iz=0; iz<=depos_order; iz++)
227 sz_jx[iz] = ((jx_type[zdir] ==
NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
228 sz_jy[iz] = ((jy_type[zdir] ==
NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
229 sz_jz[iz] = ((jz_type[zdir] ==
NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
231 int const l_jx = ((jx_type[zdir] ==
NODE) ? l_node : l_cell);
232 int const l_jy = ((jy_type[zdir] ==
NODE) ? l_node : l_cell);
233 int const l_jz = ((jz_type[zdir] ==
NODE) ? l_node : l_cell);
236 #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) 237 for (
int iz=0; iz<=depos_order; iz++){
238 for (
int ix=0; ix<=depos_order; ix++){
239 amrex::Gpu::Atomic::Add(
240 &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
241 sx_jx[ix]*sz_jx[iz]*wqx);
242 amrex::Gpu::Atomic::Add(
243 &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
244 sx_jy[ix]*sz_jy[iz]*wqy);
245 amrex::Gpu::Atomic::Add(
246 &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
247 sx_jz[ix]*sz_jz[iz]*wqz);
248 #if (defined WARPX_DIM_RZ) 250 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
252 amrex::Gpu::Atomic::Add( &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 2*imode-1), 2.*sx_jx[ix]*sz_jx[iz]*wqx*xy.real());
253 amrex::Gpu::Atomic::Add( &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 2*imode ), 2.*sx_jx[ix]*sz_jx[iz]*wqx*xy.imag());
254 amrex::Gpu::Atomic::Add( &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 2*imode-1), 2.*sx_jy[ix]*sz_jy[iz]*wqy*xy.real());
255 amrex::Gpu::Atomic::Add( &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 2*imode ), 2.*sx_jy[ix]*sz_jy[iz]*wqy*xy.imag());
256 amrex::Gpu::Atomic::Add( &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 2*imode-1), 2.*sx_jz[ix]*sz_jz[iz]*wqz*xy.real());
257 amrex::Gpu::Atomic::Add( &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 2*imode ), 2.*sx_jz[ix]*sz_jz[iz]*wqz*xy.imag());
263 #elif (defined WARPX_DIM_3D) 264 for (
int iz=0; iz<=depos_order; iz++){
265 for (
int iy=0; iy<=depos_order; iy++){
266 for (
int ix=0; ix<=depos_order; ix++){
267 amrex::Gpu::Atomic::Add(
268 &jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz),
269 sx_jx[ix]*sy_jx[iy]*sz_jx[iz]*wqx);
270 amrex::Gpu::Atomic::Add(
271 &jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz),
272 sx_jy[ix]*sy_jy[iy]*sz_jy[iz]*wqy);
273 amrex::Gpu::Atomic::Add(
274 &jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz),
275 sx_jz[ix]*sy_jz[iy]*sz_jz[iz]*wqz);
305 template <
int depos_order>
307 const amrex::ParticleReal *
const wp,
308 const amrex::ParticleReal *
const uxp,
309 const amrex::ParticleReal *
const uyp,
310 const amrex::ParticleReal *
const uzp,
312 const amrex::Array4<amrex::Real>& Jx_arr,
313 const amrex::Array4<amrex::Real>& Jy_arr,
314 const amrex::Array4<amrex::Real>& Jz_arr,
315 const long np_to_depose,
316 const amrex::Real dt,
317 const std::array<amrex::Real,3>&
dx,
318 const std::array<amrex::Real, 3> xyzmin,
319 const amrex::Dim3 lo,
321 const long n_rz_azimuthal_modes)
323 using namespace amrex;
324 #if !defined(WARPX_DIM_RZ) 325 ignore_unused(n_rz_azimuthal_modes);
330 bool const do_ionization = ion_lev;
331 Real
const dxi = 1.0_rt / dx[0];
332 #if !(defined WARPX_DIM_RZ) 333 Real
const dtsdx0 = dt*dxi;
335 Real
const xmin = xyzmin[0];
336 #if (defined WARPX_DIM_3D) 337 Real
const dyi = 1.0_rt / dx[1];
338 Real
const dtsdy0 = dt*dyi;
339 Real
const ymin = xyzmin[1];
341 Real
const dzi = 1.0_rt / dx[2];
342 Real
const dtsdz0 = dt*dzi;
343 Real
const zmin = xyzmin[2];
345 #if (defined WARPX_DIM_3D) 346 Real
const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]);
347 Real
const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]);
348 Real
const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]);
349 #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) 350 Real
const invdtdx = 1.0_rt / (dt*dx[2]);
351 Real
const invdtdz = 1.0_rt / (dt*dx[0]);
352 Real
const invvol = 1.0_rt / (dx[0]*dx[2]);
355 #if (defined WARPX_DIM_RZ) 359 Real
const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c );
364 [=] AMREX_GPU_DEVICE (
long const ip) {
367 Real
const gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
368 + uyp[ip]*uyp[ip]*clightsq
369 + uzp[ip]*uzp[ip]*clightsq);
377 ParticleReal xp, yp, zp;
378 GetPosition(ip, xp, yp, zp);
380 Real
const wqx = wq*invdtdx;
381 #if (defined WARPX_DIM_3D) 382 Real
const wqy = wq*invdtdy;
384 Real
const wqz = wq*invdtdz;
387 #if (defined WARPX_DIM_RZ) 388 Real
const xp_mid = xp - 0.5_rt * dt*uxp[ip]*gaminv;
389 Real
const yp_mid = yp - 0.5_rt * dt*uyp[ip]*gaminv;
390 Real
const xp_old = xp - dt*uxp[ip]*gaminv;
391 Real
const yp_old = yp - dt*uyp[ip]*gaminv;
392 Real
const rp_new = std::sqrt(xp*xp
394 Real
const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
395 Real
const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
396 Real costheta_new, sintheta_new;
397 if (rp_new > 0._rt) {
398 costheta_new = xp/rp_new;
399 sintheta_new = yp/rp_new;
401 costheta_new = 1._rt;
402 sintheta_new = 0._rt;
404 amrex::Real costheta_mid, sintheta_mid;
405 if (rp_mid > 0._rt) {
406 costheta_mid = xp_mid/rp_mid;
407 sintheta_mid = yp_mid/rp_mid;
409 costheta_mid = 1._rt;
410 sintheta_mid = 0._rt;
412 amrex::Real costheta_old, sintheta_old;
413 if (rp_old > 0._rt) {
414 costheta_old = xp_old/rp_old;
415 sintheta_old = yp_old/rp_old;
417 costheta_old = 1._rt;
418 sintheta_old = 0._rt;
424 double const x_new = (rp_new - xmin)*dxi;
425 double const x_old = (rp_old - xmin)*dxi;
428 double const x_new = (xp - xmin)*dxi;
429 double const x_old = x_new - dtsdx0*uxp[ip]*gaminv;
431 #if (defined WARPX_DIM_3D) 433 double const y_new = (yp - ymin)*dyi;
434 double const y_old = y_new - dtsdy0*uyp[ip]*gaminv;
437 double const z_new = (zp - zmin)*dzi;
438 double const z_old = z_new - dtsdz0*uzp[ip]*gaminv;
440 #if (defined WARPX_DIM_RZ) 441 Real
const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
442 #elif (defined WARPX_DIM_XZ) 443 Real
const vy = uyp[ip]*gaminv;
451 double sx_new[depos_order + 3] = {0.};
452 double sx_old[depos_order + 3] = {0.};
453 #if (defined WARPX_DIM_3D) 455 double sy_new[depos_order + 3] = {0.};
456 double sy_old[depos_order + 3] = {0.};
459 double sz_new[depos_order + 3] = {0.};
460 double sz_old[depos_order + 3] = {0.};
468 const int i_new = compute_shape_factor(sx_new+1, x_new);
469 const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
470 #if (defined WARPX_DIM_3D) 471 const int j_new = compute_shape_factor(sy_new+1, y_new);
472 const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
474 const int k_new = compute_shape_factor(sz_new+1, z_new);
475 const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new);
478 int dil = 1, diu = 1;
479 if (i_old < i_new) dil = 0;
480 if (i_old > i_new) diu = 0;
481 #if (defined WARPX_DIM_3D) 482 int djl = 1, dju = 1;
483 if (j_old < j_new) djl = 0;
484 if (j_old > j_new) dju = 0;
486 int dkl = 1, dku = 1;
487 if (k_old < k_new) dkl = 0;
488 if (k_old > k_new) dku = 0;
490 #if (defined WARPX_DIM_3D) 492 for (
int k=dkl; k<=depos_order+2-dku; k++) {
493 for (
int j=djl; j<=depos_order+2-dju; j++) {
494 amrex::Real sdxi = 0._rt;
495 for (
int i=dil;
i<=depos_order+1-diu;
i++) {
496 sdxi += wqx*(sx_old[
i] - sx_new[
i])*((sy_new[j] + 0.5_rt*(sy_old[j] - sy_new[j]))*sz_new[k] +
497 (0.5_rt*sy_new[j] + 1._rt/3._rt*(sy_old[j] - sy_new[j]))*(sz_old[k] - sz_new[k]));
498 amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+
i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdxi);
502 for (
int k=dkl; k<=depos_order+2-dku; k++) {
503 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
504 amrex::Real sdyj = 0._rt;
505 for (
int j=djl; j<=depos_order+1-dju; j++) {
506 sdyj += wqy*(sy_old[j] - sy_new[j])*((sz_new[k] + 0.5_rt*(sz_old[k] - sz_new[k]))*sx_new[
i] +
507 (0.5_rt*sz_new[k] + 1._rt/3._rt*(sz_old[k] - sz_new[k]))*(sx_old[
i] - sx_new[
i]));
508 amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+
i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdyj);
512 for (
int j=djl; j<=depos_order+2-dju; j++) {
513 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
514 amrex::Real sdzk = 0._rt;
515 for (
int k=dkl; k<=depos_order+1-dku; k++) {
516 sdzk += wqz*(sz_old[k] - sz_new[k])*((sx_new[
i] + 0.5_rt*(sx_old[
i] - sx_new[
i]))*sy_new[j] +
517 (0.5_rt*sx_new[
i] + 1._rt/3._rt*(sx_old[
i] - sx_new[
i]))*(sy_old[j] - sy_new[j]));
518 amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdzk);
523 #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) 525 for (
int k=dkl; k<=depos_order+2-dku; k++) {
526 amrex::Real sdxi = 0._rt;
527 for (
int i=dil;
i<=depos_order+1-diu;
i++) {
528 sdxi += wqx*(sx_old[
i] - sx_new[
i])*(sz_new[k] + 0.5_rt*(sz_old[k] - sz_new[k]));
529 amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+
i, lo.y+k_new-1+k, 0, 0), sdxi);
530 #if (defined WARPX_DIM_RZ) 532 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
534 const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
535 amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+
i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real());
536 amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+
i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
537 xy_mid = xy_mid*xy_mid0;
542 for (
int k=dkl; k<=depos_order+2-dku; k++) {
543 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
544 Real
const sdyj = wq*vy*invvol*((sz_new[k] + 0.5_rt * (sz_old[k] - sz_new[k]))*sx_new[
i] +
545 (0.5_rt * sz_new[k] + 1._rt / 3._rt *(sz_old[k] - sz_new[k]))*(sx_old[
i] - sx_new[
i]));
546 amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+
i, lo.y+k_new-1+k, 0, 0), sdyj);
547 #if (defined WARPX_DIM_RZ) 552 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
555 const Complex djt_cmplx = -2._rt * I*(i_new-1 +
i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode*
556 (sx_new[
i]*sz_new[k]*(xy_new - xy_mid) + sx_old[
i]*sz_old[k]*(xy_mid - xy_old));
557 amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+
i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real());
558 amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+
i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
559 xy_new = xy_new*xy_new0;
560 xy_mid = xy_mid*xy_mid0;
561 xy_old = xy_old*xy_old0;
566 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
568 for (
int k=dkl; k<=depos_order+1-dku; k++) {
569 sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[
i] + 0.5_rt * (sx_old[
i] - sx_new[
i]));
570 amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
571 #if (defined WARPX_DIM_RZ) 573 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
575 const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
576 amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real());
577 amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
578 xy_mid = xy_mid*xy_mid0;
614 template <
int depos_order>
616 const amrex::ParticleReal*
const wp,
617 const amrex::ParticleReal*
const uxp,
618 const amrex::ParticleReal*
const uyp,
619 const amrex::ParticleReal*
const uzp,
620 const int*
const ion_lev,
621 amrex::FArrayBox& jx_fab,
622 amrex::FArrayBox& jy_fab,
623 amrex::FArrayBox& jz_fab,
624 const long np_to_depose,
625 const amrex::Real dt,
626 const std::array<amrex::Real,3>&
dx,
627 const std::array<amrex::Real,3>& xyzmin,
628 const amrex::Dim3 lo,
630 const long n_rz_azimuthal_modes)
632 #if (defined WARPX_DIM_RZ) 633 amrex::ignore_unused(GetPosition,
634 wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab,
635 np_to_depose, dt, dx, xyzmin, lo, q, n_rz_azimuthal_modes);
636 amrex::Abort(
"Vay deposition not implemented in RZ geometry");
639 #if !(defined WARPX_DIM_RZ) 640 amrex::ignore_unused(n_rz_azimuthal_modes);
643 const bool do_ionization = ion_lev;
646 const amrex::Real dxi = 1._rt / dx[0];
647 const amrex::Real dzi = 1._rt / dx[2];
648 #if (defined WARPX_DIM_3D) 649 const amrex::Real dyi = 1._rt / dx[1];
653 const amrex::Real invdt = 1._rt / dt;
656 #if (defined WARPX_DIM_XZ) 657 const amrex::Real invvol = dxi * dzi;
658 #elif (defined WARPX_DIM_3D) 659 const amrex::Real invvol = dxi * dyi * dzi;
663 const amrex::Real xmin = xyzmin[0];
664 const amrex::Real zmin = xyzmin[2];
665 #if (defined WARPX_DIM_3D) 666 const amrex::Real ymin = xyzmin[1];
670 #if (defined WARPX_DIM_3D) 671 const amrex::Real onethird = 1._rt / 3._rt;
672 const amrex::Real onesixth = 1._rt / 6._rt;
676 const amrex::Real invcsq = 1._rt / (PhysConst::c * PhysConst::c);
679 amrex::Array4<amrex::Real>
const& jx_arr = jx_fab.array();
680 amrex::Array4<amrex::Real>
const& jy_arr = jy_fab.array();
681 amrex::Array4<amrex::Real>
const& jz_arr = jz_fab.array();
684 amrex::ParallelFor(np_to_depose, [=] AMREX_GPU_DEVICE (
long ip)
687 const amrex::Real invgam = 1._rt / std::sqrt(1._rt + uxp[ip] * uxp[ip] * invcsq
688 + uyp[ip] * uyp[ip] * invcsq
689 + uzp[ip] * uzp[ip] * invcsq);
691 amrex::Real wq = q * wp[ip];
692 if (do_ionization) wq *= ion_lev[ip];
695 amrex::ParticleReal xp, yp, zp;
696 GetPosition(ip, xp, yp, zp);
699 const amrex::Real vx = uxp[ip] * invgam;
700 const amrex::Real vy = uyp[ip] * invgam;
701 const amrex::Real vz = uzp[ip] * invgam;
704 #if (defined WARPX_DIM_XZ) 705 const amrex::Real wqy = wq * vy * invvol;
710 double const x_new = (xp - xmin) * dxi;
711 double const x_old = x_new - vx * dt * dxi;
712 #if (defined WARPX_DIM_3D) 714 double const y_new = (yp - ymin) * dyi;
715 double const y_old = y_new - vy * dt * dyi;
718 double const z_new = (zp - zmin) * dzi;
719 double const z_old = z_new - vz * dt * dzi;
723 double sx_new[depos_order+1] = {0.};
724 double sx_old[depos_order+1] = {0.};
725 #if (defined WARPX_DIM_3D) 727 double sy_new[depos_order+1] = {0.};
728 double sy_old[depos_order+1] = {0.};
731 double sz_new[depos_order+1] = {0.};
732 double sz_old[depos_order+1] = {0.};
739 const int i_new = compute_shape_factor(sx_new, x_new);
740 #if (defined WARPX_DIM_3D) 743 const int j_new = compute_shape_factor(sy_new, y_new);
747 const int k_new = compute_shape_factor(sz_new, z_new);
753 const int i_old = compute_shape_factor(sx_old, x_old);
754 #if (defined WARPX_DIM_3D) 757 const int j_old = compute_shape_factor(sy_old, y_old);
761 const int k_old = compute_shape_factor(sz_old, z_old);
764 #if (defined WARPX_DIM_XZ) 766 for (
int k=0; k<=depos_order; k++) {
767 for (
int i=0;
i<=depos_order;
i++) {
771 amrex::Real
const sxn_szn = sx_new[
i] * sz_new[k];
772 amrex::Real
const sxo_szn = sx_old[
i] * sz_new[k];
773 amrex::Real
const sxn_szo = sx_new[
i] * sz_old[k];
774 amrex::Real
const sxo_szo = sx_old[
i] * sz_old[k];
777 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new +
i, lo.y + k_new + k, 0, 0),
778 wq * invvol * invdt * 0.5_rt * sxn_szn);
780 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old +
i, lo.y + k_new + k, 0, 0),
781 - wq * invvol * invdt * 0.5_rt * sxo_szn);
783 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new +
i, lo.y + k_old + k, 0, 0),
784 wq * invvol * invdt * 0.5_rt * sxn_szo);
786 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old +
i, lo.y + k_old + k, 0, 0),
787 - wq * invvol * invdt * 0.5_rt * sxo_szo);
790 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new +
i, lo.y + k_new + k, 0, 0),
791 wqy * 0.25_rt * sxn_szn);
793 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new +
i, lo.y + k_old + k, 0, 0),
794 wqy * 0.25_rt * sxn_szo);
796 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old +
i, lo.y + k_new + k, 0, 0),
797 wqy * 0.25_rt * sxo_szn);
799 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old +
i, lo.y + k_old + k, 0, 0),
800 wqy * 0.25_rt * sxo_szo);
803 amrex::Gpu::Atomic::Add(&jz_arr(lo.x + i_new +
i, lo.y + k_new + k, 0, 0),
804 wq * invvol * invdt * 0.5_rt * sxn_szn);
806 amrex::Gpu::Atomic::Add(&jz_arr(lo.x+i_new+
i,lo.y+k_old+k,0,0),
807 - wq * invvol * invdt * 0.5_rt * sxn_szo);
809 amrex::Gpu::Atomic::Add(&jz_arr(lo.x+i_old+
i,lo.y+k_new+k,0,0),
810 wq * invvol * invdt * 0.5_rt * sxo_szn);
812 amrex::Gpu::Atomic::Add(&jz_arr(lo.x + i_old +
i, lo.y + k_old + k, 0, 0),
813 - wq * invvol * invdt * 0.5_rt * sxo_szo);
817 #elif (defined WARPX_DIM_3D) 819 for (
int k=0; k<=depos_order; k++) {
820 for (
int j=0; j<=depos_order; j++) {
822 amrex::Real
const syn_szn = sy_new[j] * sz_new[k];
823 amrex::Real
const syo_szn = sy_old[j] * sz_new[k];
824 amrex::Real
const syn_szo = sy_new[j] * sz_old[k];
825 amrex::Real
const syo_szo = sy_old[j] * sz_old[k];
827 for (
int i=0;
i<=depos_order;
i++) {
829 amrex::Real
const sxn_syn_szn = sx_new[
i] * syn_szn;
830 amrex::Real
const sxo_syn_szn = sx_old[
i] * syn_szn;
831 amrex::Real
const sxn_syo_szn = sx_new[
i] * syo_szn;
832 amrex::Real
const sxo_syo_szn = sx_old[
i] * syo_szn;
833 amrex::Real
const sxn_syn_szo = sx_new[
i] * syn_szo;
834 amrex::Real
const sxo_syn_szo = sx_old[
i] * syn_szo;
835 amrex::Real
const sxn_syo_szo = sx_new[
i] * syo_szo;
836 amrex::Real
const sxo_syo_szo = sx_old[
i] * syo_szo;
839 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new +
i, lo.y + j_new + j, lo.z + k_new + k),
840 wq * invvol * invdt * onethird * sxn_syn_szn);
842 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old +
i, lo.y + j_new + j, lo.z + k_new + k),
843 - wq * invvol * invdt * onethird * sxo_syn_szn);
845 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new +
i, lo.y + j_old + j, lo.z + k_new + k),
846 wq * invvol * invdt * onesixth * sxn_syo_szn);
848 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old +
i, lo.y + j_old + j,lo.z + k_new + k),
849 - wq * invvol * invdt * onesixth * sxo_syo_szn);
851 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new +
i, lo.y + j_new + j, lo.z + k_old + k),
852 wq * invvol * invdt * onesixth * sxn_syn_szo);
854 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old +
i, lo.y + j_new + j, lo.z + k_old + k),
855 - wq * invvol * invdt * onesixth * sxo_syn_szo);
857 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new +
i, lo.y + j_old + j, lo.z + k_old + k),
858 wq * invvol * invdt * onethird * sxn_syo_szo);
860 amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old +
i, lo.y + j_old + j, lo.z + k_old + k),
861 - wq * invvol * invdt * onethird * sxo_syo_szo);
864 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new +
i, lo.y + j_new + j, lo.z + k_new + k),
865 wq * invvol * invdt * onethird * sxn_syn_szn);
867 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new +
i, lo.y + j_old + j, lo.z + k_new + k),
868 - wq * invvol * invdt * onethird * sxn_syo_szn);
870 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old +
i, lo.y + j_new + j, lo.z + k_new + k),
871 wq * invvol * invdt * onesixth * sxo_syn_szn);
873 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old +
i, lo.y + j_old + j, lo.z + k_new + k),
874 - wq * invvol * invdt * onesixth * sxo_syo_szn);
876 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new +
i, lo.y + j_new + j, lo.z + k_old + k),
877 wq * invvol * invdt * onesixth * sxn_syn_szo);
879 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new +
i, lo.y + j_old + j, lo.z + k_old + k),
880 - wq * invvol * invdt * onesixth * sxn_syo_szo);
882 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old +
i, lo.y + j_new + j, lo.z + k_old + k),
883 wq * invvol * invdt * onethird * sxo_syn_szo);
885 amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old +
i, lo.y + j_old + j, lo.z + k_old + k),
886 - wq * invvol * invdt * onethird * sxo_syo_szo);
889 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new +
i, lo.y + j_new + j, lo.z + k_new + k),
890 wq * invvol * invdt * onethird * sxn_syn_szn);
892 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new +
i, lo.y + j_new + j, lo.z + k_old + k),
893 - wq * invvol * invdt * onethird * sxn_syn_szo);
895 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old +
i, lo.y + j_new + j, lo.z + k_new + k),
896 wq * invvol * invdt * onesixth * sxo_syn_szn);
898 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old +
i, lo.y + j_new + j, lo.z + k_old + k),
899 - wq * invvol * invdt * onesixth * sxo_syn_szo);
901 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new +
i, lo.y + j_old + j, lo.z + k_new + k),
902 wq * invvol * invdt * onesixth * sxn_syo_szn);
904 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new +
i, lo.y + j_old + j, lo.z + k_old + k),
905 - wq * invvol * invdt * onesixth * sxn_syo_szo);
907 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old +
i, lo.y + j_old + j, lo.z + k_new + k),
908 wq * invvol * invdt * onethird * sxo_syo_szn);
910 amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old +
i, lo.y + j_old + j, lo.z + k_old + k),
911 - wq * invvol * invdt * onethird * sxo_syo_szo);
917 #endif // #if !(defined WARPX_DIM_RZ) 919 #endif // CURRENTDEPOSITION_H_ amrex::GpuComplex< amrex::Real > Complex
Definition: WarpX_Complex.H:22
Definition: wp_parser.tab.cpp:115
void doVayDepositionShapeN(const GetParticlePosition &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp, const amrex::ParticleReal *const uyp, const amrex::ParticleReal *const uzp, const int *const ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, const long np_to_depose, const amrex::Real dt, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const long n_rz_azimuthal_modes)
Vay current deposition (Vay et al, 2013) for thread thread_num: deposit D in real space and store the...
Definition: CurrentDeposition.H:615
int dx
Definition: compute_domain.py:35
Definition: ShapeFactors.H:118
void doDepositionShapeN(const GetParticlePosition &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp, const amrex::ParticleReal *const uyp, const amrex::ParticleReal *const uzp, const int *const ion_lev, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, const long np_to_depose, const amrex::Real dt, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const long n_rz_azimuthal_modes)
Current Deposition for thread thread_num /param GetPosition : A functor for returning the particle po...
Definition: CurrentDeposition.H:41
Definition: ShapeFactors.H:22
i
Definition: check_interp_points_and_weights.py:171
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel...
Definition: GetAndSetPosition.H:25
void doEsirkepovDepositionShapeN(const GetParticlePosition &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp, const amrex::ParticleReal *const uyp, const amrex::ParticleReal *const uzp, const int *ion_lev, const amrex::Array4< amrex::Real > &Jx_arr, const amrex::Array4< amrex::Real > &Jy_arr, const amrex::Array4< amrex::Real > &Jz_arr, const long np_to_depose, const amrex::Real dt, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > xyzmin, const amrex::Dim3 lo, const amrex::Real q, const long n_rz_azimuthal_modes)
Esirkepov Current Deposition for thread thread_num.
Definition: CurrentDeposition.H:306