8 #ifndef WARPX_CURRENTDEPOSITION_H_
9 #define WARPX_CURRENTDEPOSITION_H_
47 template <
int depos_order>
50 const amrex::ParticleReal yp,
51 const amrex::ParticleReal zp,
52 const amrex::ParticleReal wq,
53 const amrex::ParticleReal vx,
54 const amrex::ParticleReal vy,
55 const amrex::ParticleReal vz,
62 const amrex::Real relative_time,
64 const amrex::Real dxi,
65 const amrex::Real dyi),
67 const amrex::Real
xmin,
68 const amrex::Real ymin),
69 const amrex::Real invvol,
71 const int n_rz_azimuthal_modes)
73 using namespace amrex::literals;
74 #if !defined(WARPX_DIM_RZ)
77 #if defined(WARPX_DIM_1D_Z)
80 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
84 constexpr
int zdir = WARPX_ZINDEX;
89 #if defined(WARPX_DIM_RZ)
92 const amrex::Real xpmid = xp + relative_time*vx;
93 const amrex::Real ypmid = yp + relative_time*vy;
94 const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid);
98 costheta = xpmid/rpmid;
99 sintheta = ypmid/rpmid;
105 const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta);
106 const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta);
108 const amrex::Real wqx = wq*invvol*vx;
109 const amrex::Real wqy = wq*invvol*vy;
111 const amrex::Real wqz = wq*invvol*vz;
115 #if (AMREX_SPACEDIM >= 2)
118 #if defined(WARPX_DIM_RZ)
120 const double xmid = (rpmid -
xmin)*dxi;
122 const double xmid = ((xp -
xmin) + relative_time*vx)*dxi;
129 double sx_node[depos_order + 1] = {0.};
130 double sx_cell[depos_order + 1] = {0.};
133 if (jx_type[0] ==
NODE || jy_type[0] ==
NODE || jz_type[0] ==
NODE) {
134 j_node = compute_shape_factor(sx_node, xmid);
136 if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
137 j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
140 amrex::Real sx_jx[depos_order + 1] = {0._rt};
141 amrex::Real sx_jy[depos_order + 1] = {0._rt};
142 amrex::Real sx_jz[depos_order + 1] = {0._rt};
143 for (
int ix=0;
ix<=depos_order;
ix++)
145 sx_jx[
ix] = ((jx_type[0] ==
NODE) ? amrex::Real(sx_node[
ix]) : amrex::Real(sx_cell[
ix]));
146 sx_jy[
ix] = ((jy_type[0] ==
NODE) ? amrex::Real(sx_node[
ix]) : amrex::Real(sx_cell[
ix]));
147 sx_jz[
ix] = ((jz_type[0] ==
NODE) ? amrex::Real(sx_node[
ix]) : amrex::Real(sx_cell[
ix]));
150 int const j_jx = ((jx_type[0] ==
NODE) ? j_node : j_cell);
151 int const j_jy = ((jy_type[0] ==
NODE) ? j_node : j_cell);
152 int const j_jz = ((jz_type[0] ==
NODE) ? j_node : j_cell);
155 #if defined(WARPX_DIM_3D)
158 const double ymid = ((yp - ymin) + relative_time*vy)*dyi;
159 double sy_node[depos_order + 1] = {0.};
160 double sy_cell[depos_order + 1] = {0.};
163 if (jx_type[1] ==
NODE || jy_type[1] ==
NODE || jz_type[1] ==
NODE) {
164 k_node = compute_shape_factor(sy_node, ymid);
166 if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
167 k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
169 amrex::Real sy_jx[depos_order + 1] = {0._rt};
170 amrex::Real sy_jy[depos_order + 1] = {0._rt};
171 amrex::Real sy_jz[depos_order + 1] = {0._rt};
172 for (
int iy=0;
iy<=depos_order;
iy++)
174 sy_jx[
iy] = ((jx_type[1] ==
NODE) ? amrex::Real(sy_node[
iy]) : amrex::Real(sy_cell[
iy]));
175 sy_jy[
iy] = ((jy_type[1] ==
NODE) ? amrex::Real(sy_node[
iy]) : amrex::Real(sy_cell[
iy]));
176 sy_jz[
iy] = ((jz_type[1] ==
NODE) ? amrex::Real(sy_node[
iy]) : amrex::Real(sy_cell[
iy]));
178 int const k_jx = ((jx_type[1] ==
NODE) ? k_node : k_cell);
179 int const k_jy = ((jy_type[1] ==
NODE) ? k_node : k_cell);
180 int const k_jz = ((jz_type[1] ==
NODE) ? k_node : k_cell);
185 const double zmid = ((zp - zmin) + relative_time*vz)*dzi;
186 double sz_node[depos_order + 1] = {0.};
187 double sz_cell[depos_order + 1] = {0.};
190 if (jx_type[zdir] ==
NODE || jy_type[zdir] ==
NODE || jz_type[zdir] ==
NODE) {
191 l_node = compute_shape_factor(sz_node, zmid);
193 if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
194 l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
196 amrex::Real sz_jx[depos_order + 1] = {0._rt};
197 amrex::Real sz_jy[depos_order + 1] = {0._rt};
198 amrex::Real sz_jz[depos_order + 1] = {0._rt};
199 for (
int iz=0;
iz<=depos_order;
iz++)
201 sz_jx[
iz] = ((jx_type[zdir] ==
NODE) ? amrex::Real(sz_node[
iz]) : amrex::Real(sz_cell[
iz]));
202 sz_jy[
iz] = ((jy_type[zdir] ==
NODE) ? amrex::Real(sz_node[
iz]) : amrex::Real(sz_cell[
iz]));
203 sz_jz[
iz] = ((jz_type[zdir] ==
NODE) ? amrex::Real(sz_node[
iz]) : amrex::Real(sz_cell[
iz]));
205 int const l_jx = ((jx_type[zdir] ==
NODE) ? l_node : l_cell);
206 int const l_jy = ((jy_type[zdir] ==
NODE) ? l_node : l_cell);
207 int const l_jz = ((jz_type[zdir] ==
NODE) ? l_node : l_cell);
210 #if defined(WARPX_DIM_1D_Z)
211 for (
int iz=0;
iz<=depos_order;
iz++){
213 &jx_arr(lo.
x+l_jx+
iz, 0, 0, 0),
216 &jy_arr(lo.
x+l_jy+
iz, 0, 0, 0),
219 &jz_arr(lo.
x+l_jz+
iz, 0, 0, 0),
223 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
224 for (
int iz=0;
iz<=depos_order;
iz++){
225 for (
int ix=0;
ix<=depos_order;
ix++){
227 &jx_arr(lo.
x+j_jx+
ix, lo.
y+l_jx+
iz, 0, 0),
228 sx_jx[
ix]*sz_jx[
iz]*wqx);
230 &jy_arr(lo.
x+j_jy+
ix, lo.
y+l_jy+
iz, 0, 0),
231 sx_jy[
ix]*sz_jy[
iz]*wqy);
233 &jz_arr(lo.
x+j_jz+
ix, lo.
y+l_jz+
iz, 0, 0),
234 sx_jz[
ix]*sz_jz[
iz]*wqz);
235 #if defined(WARPX_DIM_RZ)
237 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
250 #elif defined(WARPX_DIM_3D)
251 for (
int iz=0;
iz<=depos_order;
iz++){
252 for (
int iy=0;
iy<=depos_order;
iy++){
253 for (
int ix=0;
ix<=depos_order;
ix++){
255 &jx_arr(lo.
x+j_jx+
ix, lo.
y+k_jx+
iy, lo.
z+l_jx+
iz),
256 sx_jx[
ix]*sy_jx[
iy]*sz_jx[
iz]*wqx);
258 &jy_arr(lo.
x+j_jy+
ix, lo.
y+k_jy+
iy, lo.
z+l_jy+
iz),
259 sx_jy[
ix]*sy_jy[
iy]*sz_jy[
iz]*wqy);
261 &jz_arr(lo.
x+j_jz+
ix, lo.
y+k_jz+
iy, lo.
z+l_jz+
iz),
262 sx_jz[
ix]*sy_jz[
iy]*sz_jz[
iz]*wqz);
291 template <
int depos_order>
293 const amrex::ParticleReal *
const wp,
294 const amrex::ParticleReal *
const uxp,
295 const amrex::ParticleReal *
const uyp,
296 const amrex::ParticleReal *
const uzp,
302 amrex::Real relative_time,
303 const std::array<amrex::Real,3>&
dx,
304 const std::array<amrex::Real,3>& xyzmin,
307 int n_rz_azimuthal_modes)
309 using namespace amrex::literals;
311 #if !defined(WARPX_DIM_RZ)
317 const bool do_ionization = ion_lev;
318 const amrex::Real dzi = 1.0_rt/
dx[2];
319 #if defined(WARPX_DIM_1D_Z)
320 const amrex::Real invvol = dzi;
322 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
323 const amrex::Real dxi = 1.0_rt/
dx[0];
324 const amrex::Real invvol = dxi*dzi;
325 #elif defined(WARPX_DIM_3D)
326 const amrex::Real dxi = 1.0_rt/
dx[0];
327 const amrex::Real dyi = 1.0_rt/
dx[1];
328 const amrex::Real invvol = dxi*dyi*dzi;
331 #if (AMREX_SPACEDIM >= 2)
332 const amrex::Real
xmin = xyzmin[0];
334 #if defined(WARPX_DIM_3D)
335 const amrex::Real ymin = xyzmin[1];
337 const amrex::Real zmin = xyzmin[2];
352 amrex::ParticleReal xp, yp, zp;
353 GetPosition(ip, xp, yp, zp);
356 const amrex::Real gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
357 + uyp[ip]*uyp[ip]*clightsq
358 + uzp[ip]*uzp[ip]*clightsq);
359 const amrex::Real vx = uxp[ip]*gaminv;
360 const amrex::Real vy = uyp[ip]*gaminv;
361 const amrex::Real vz = uzp[ip]*gaminv;
363 amrex::Real wq = q*wp[ip];
368 doDepositionShapeNKernel<depos_order>(xp, yp, zp, wq, vx, vy, vz, jx_arr, jy_arr, jz_arr,
369 jx_type, jy_type, jz_type,
373 invvol, lo, n_rz_azimuthal_modes);
400 template <
int depos_order>
402 const amrex::ParticleReal *
const wp,
403 const amrex::ParticleReal *
const uxp_n,
404 const amrex::ParticleReal *
const uyp_n,
405 const amrex::ParticleReal *
const uzp_n,
406 const amrex::ParticleReal *
const uxp,
407 const amrex::ParticleReal *
const uyp,
408 const amrex::ParticleReal *
const uzp,
409 const int *
const ion_lev,
413 const long np_to_deposit,
414 const std::array<amrex::Real,3>&
dx,
415 const std::array<amrex::Real,3>& xyzmin,
418 const int n_rz_azimuthal_modes)
420 using namespace amrex::literals;
421 #if !defined(WARPX_DIM_RZ)
427 const bool do_ionization = ion_lev;
428 const amrex::Real dzi = 1.0_rt/
dx[2];
429 #if defined(WARPX_DIM_1D_Z)
430 const amrex::Real invvol = dzi;
432 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
433 const amrex::Real dxi = 1.0_rt/
dx[0];
434 const amrex::Real invvol = dxi*dzi;
435 #elif defined(WARPX_DIM_3D)
436 const amrex::Real dxi = 1.0_rt/
dx[0];
437 const amrex::Real dyi = 1.0_rt/
dx[1];
438 const amrex::Real invvol = dxi*dyi*dzi;
441 #if (AMREX_SPACEDIM >= 2)
442 const amrex::Real
xmin = xyzmin[0];
444 #if defined(WARPX_DIM_3D)
445 const amrex::Real ymin = xyzmin[1];
447 const amrex::Real zmin = xyzmin[2];
460 amrex::ParticleReal xp, yp, zp;
461 GetPosition(ip, xp, yp, zp);
467 const amrex::ParticleReal uxp_np1 = 2._prt*uxp[ip] - uxp_n[ip];
468 const amrex::ParticleReal uyp_np1 = 2._prt*uyp[ip] - uyp_n[ip];
469 const amrex::ParticleReal uzp_np1 = 2._prt*uzp[ip] - uzp_n[ip];
470 const amrex::ParticleReal gamma_n = std::sqrt(1._prt + (uxp_n[ip]*uxp_n[ip] + uyp_n[ip]*uyp_n[ip] + uzp_n[ip]*uzp_n[ip])*inv_c2);
471 const amrex::ParticleReal gamma_np1 = std::sqrt(1._prt + (uxp_np1*uxp_np1 + uyp_np1*uyp_np1 + uzp_np1*uzp_np1)*inv_c2);
472 const amrex::ParticleReal gaminv = 2.0_prt/(gamma_n + gamma_np1);
474 const amrex::Real vx = uxp[ip]*gaminv;
475 const amrex::Real vy = uyp[ip]*gaminv;
476 const amrex::Real vz = uzp[ip]*gaminv;
478 amrex::Real wq = q*wp[ip];
483 const amrex::Real relative_time = 0._rt;
484 doDepositionShapeNKernel<depos_order>(xp, yp, zp, wq, vx, vy, vz, jx_arr, jy_arr, jz_arr,
485 jx_type, jy_type, jz_type,
489 invvol, lo, n_rz_azimuthal_modes);
518 template <
int depos_order>
520 const amrex::ParticleReal *
const wp,
521 const amrex::ParticleReal *
const uxp,
522 const amrex::ParticleReal *
const uyp,
523 const amrex::ParticleReal *
const uzp,
529 const amrex::Real relative_time,
530 const std::array<amrex::Real,3>&
dx,
531 const std::array<amrex::Real,3>& xyzmin,
534 int n_rz_azimuthal_modes,
540 using namespace amrex::literals;
542 #if defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)
543 using namespace amrex;
554 constexpr
int zdir = WARPX_ZINDEX;
561 const auto domain = geom.
Domain();
564 sample_tbox.
grow(depos_order);
572 const int nblocks = a_bins.
numBins();
576 const std::size_t shared_mem_bytes = npts*
sizeof(amrex::Real);
580 "Tile size too big for GPU shared memory current deposition");
587 const int bin_id = blockIdx.x;
588 const unsigned int bin_start = offsets_ptr[bin_id];
589 const unsigned int bin_stop = offsets_ptr[bin_id+1];
591 if (bin_start == bin_stop) {
return; }
596 ParticleReal xp, yp, zp;
597 GetPosition(permutation[bin_start], xp, yp, zp);
598 #if defined(WARPX_DIM_3D)
599 IntVect iv =
IntVect(
int( amrex::Math::floor((xp-plo[0]) * dxiarr[0]) ),
600 int( amrex::Math::floor((yp-plo[1]) * dxiarr[1]) ),
601 int( amrex::Math::floor((zp-plo[2]) * dxiarr[2]) ));
602 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
603 IntVect iv =
IntVect(
int( amrex::Math::floor((xp-plo[0]) * dxiarr[0]) ),
604 int( amrex::Math::floor((zp-plo[1]) * dxiarr[1]) ));
605 #elif defined(WARPX_DIM_1D_Z)
606 IntVect iv =
IntVect(
int( amrex::Math::floor((zp-plo[0]) * dxiarr[0]) ));
608 iv += domain.smallEnd();
612 buffer_box.
grow(depos_order);
618 amrex::Real*
const shared = gsm.
dataPtr();
628 volatile amrex::Real* vs = shared;
629 for (
int i = threadIdx.x;
i < npts;
i += blockDim.x){
633 for (
unsigned int ip_orig = bin_start+threadIdx.x; ip_orig<bin_stop; ip_orig += blockDim.x)
635 const unsigned int ip = permutation[ip_orig];
636 depositComponent<depos_order>(GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_buff, jx_type,
637 relative_time,
dx, xyzmin, lo, q, n_rz_azimuthal_modes,
638 ip, zdir,
NODE, CELL, 0);
642 addLocalToGlobal(tbox_x, jx_arr, jx_buff);
643 for (
int i = threadIdx.x;
i < npts;
i += blockDim.x){
648 for (
unsigned int ip_orig = bin_start+threadIdx.x; ip_orig<bin_stop; ip_orig += blockDim.x)
650 const unsigned int ip = permutation[ip_orig];
651 depositComponent<depos_order>(GetPosition, wp, uxp, uyp, uzp, ion_lev, jy_buff, jy_type,
652 relative_time,
dx, xyzmin, lo, q, n_rz_azimuthal_modes,
653 ip, zdir,
NODE, CELL, 1);
657 addLocalToGlobal(tbox_y, jy_arr, jy_buff);
658 for (
int i = threadIdx.x;
i < npts;
i += blockDim.x){
663 for (
unsigned int ip_orig = bin_start+threadIdx.x; ip_orig<bin_stop; ip_orig += blockDim.x)
665 const unsigned int ip = permutation[ip_orig];
666 depositComponent<depos_order>(GetPosition, wp, uxp, uyp, uzp, ion_lev, jz_buff, jz_type,
667 relative_time,
dx, xyzmin, lo, q, n_rz_azimuthal_modes,
668 ip, zdir,
NODE, CELL, 2);
672 addLocalToGlobal(tbox_z, jz_arr, jz_buff);
678 ignore_unused(GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time,
dx, xyzmin, lo, q, n_rz_azimuthal_modes, a_bins, box, geom, a_tbox_max_size);
707 template <
int depos_order>
709 const amrex::ParticleReal *
const wp,
710 const amrex::ParticleReal *
const uxp,
711 const amrex::ParticleReal *
const uyp,
712 const amrex::ParticleReal *
const uzp,
719 amrex::Real relative_time,
720 const std::array<amrex::Real,3>&
dx,
721 std::array<amrex::Real, 3> xyzmin,
724 int n_rz_azimuthal_modes)
726 using namespace amrex;
727 using namespace amrex::literals;
729 #if !defined(WARPX_DIM_RZ)
735 bool const do_ionization = ion_lev;
736 #if !defined(WARPX_DIM_1D_Z)
737 Real
const dxi = 1.0_rt /
dx[0];
739 #if !defined(WARPX_DIM_1D_Z)
740 Real
const xmin = xyzmin[0];
742 #if defined(WARPX_DIM_3D)
743 Real
const dyi = 1.0_rt /
dx[1];
744 Real
const ymin = xyzmin[1];
746 Real
const dzi = 1.0_rt /
dx[2];
747 Real
const zmin = xyzmin[2];
749 #if defined(WARPX_DIM_3D)
750 Real
const invdtdx = 1.0_rt / (
dt*
dx[1]*
dx[2]);
751 Real
const invdtdy = 1.0_rt / (
dt*
dx[0]*
dx[2]);
752 Real
const invdtdz = 1.0_rt / (
dt*
dx[0]*
dx[1]);
753 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
754 Real
const invdtdx = 1.0_rt / (
dt*
dx[2]);
755 Real
const invdtdz = 1.0_rt / (
dt*
dx[0]);
756 Real
const invvol = 1.0_rt / (
dx[0]*
dx[2]);
757 #elif defined(WARPX_DIM_1D_Z)
758 Real
const invdtdz = 1.0_rt / (
dt*
dx[0]);
759 Real
const invvol = 1.0_rt / (
dx[2]);
762 #if defined(WARPX_DIM_RZ)
767 #if !defined(WARPX_DIM_1D_Z)
768 Real constexpr one_third = 1.0_rt / 3.0_rt;
769 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
777 Real
const gaminv = 1.0_rt/std::sqrt(1.0_rt + uxp[ip]*uxp[ip]*clightsq
778 + uyp[ip]*uyp[ip]*clightsq
779 + uzp[ip]*uzp[ip]*clightsq);
787 ParticleReal xp, yp, zp;
788 GetPosition(ip, xp, yp, zp);
790 #if !defined(WARPX_DIM_1D_Z)
791 Real
const wqx = wq*invdtdx;
793 #if defined(WARPX_DIM_3D)
794 Real
const wqy = wq*invdtdy;
796 Real
const wqz = wq*invdtdz;
799 #if defined(WARPX_DIM_RZ)
800 Real
const xp_new = xp + (relative_time + 0.5_rt*
dt)*uxp[ip]*gaminv;
801 Real
const yp_new = yp + (relative_time + 0.5_rt*
dt)*uyp[ip]*gaminv;
802 Real
const xp_mid = xp_new - 0.5_rt*
dt*uxp[ip]*gaminv;
803 Real
const yp_mid = yp_new - 0.5_rt*
dt*uyp[ip]*gaminv;
804 Real
const xp_old = xp_new -
dt*uxp[ip]*gaminv;
805 Real
const yp_old = yp_new -
dt*uyp[ip]*gaminv;
806 Real
const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
807 Real
const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
808 Real
const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
809 Real costheta_new, sintheta_new;
810 if (rp_new > 0._rt) {
811 costheta_new = xp_new/rp_new;
812 sintheta_new = yp_new/rp_new;
814 costheta_new = 1._rt;
815 sintheta_new = 0._rt;
817 amrex::Real costheta_mid, sintheta_mid;
818 if (rp_mid > 0._rt) {
819 costheta_mid = xp_mid/rp_mid;
820 sintheta_mid = yp_mid/rp_mid;
822 costheta_mid = 1._rt;
823 sintheta_mid = 0._rt;
825 amrex::Real costheta_old, sintheta_old;
826 if (rp_old > 0._rt) {
827 costheta_old = xp_old/rp_old;
828 sintheta_old = yp_old/rp_old;
830 costheta_old = 1._rt;
831 sintheta_old = 0._rt;
837 double const x_new = (rp_new -
xmin)*dxi;
838 double const x_old = (rp_old -
xmin)*dxi;
840 #if !defined(WARPX_DIM_1D_Z)
842 double const x_new = (xp -
xmin + (relative_time + 0.5_rt*
dt)*uxp[ip]*gaminv)*dxi;
843 double const x_old = x_new -
dt*dxi*uxp[ip]*gaminv;
846 #if defined(WARPX_DIM_3D)
848 double const y_new = (yp - ymin + (relative_time + 0.5_rt*
dt)*uyp[ip]*gaminv)*dyi;
849 double const y_old = y_new -
dt*dyi*uyp[ip]*gaminv;
852 double const z_new = (zp - zmin + (relative_time + 0.5_rt*
dt)*uzp[ip]*gaminv)*dzi;
853 double const z_old = z_new -
dt*dzi*uzp[ip]*gaminv;
855 #if defined(WARPX_DIM_RZ)
856 Real
const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
857 #elif defined(WARPX_DIM_XZ)
858 Real
const vy = uyp[ip]*gaminv;
859 #elif defined(WARPX_DIM_1D_Z)
860 Real
const vx = uxp[ip]*gaminv;
861 Real
const vy = uyp[ip]*gaminv;
869 #if !defined(WARPX_DIM_1D_Z)
870 double sx_new[depos_order + 3] = {0.};
871 double sx_old[depos_order + 3] = {0.};
873 #if defined(WARPX_DIM_3D)
875 double sy_new[depos_order + 3] = {0.};
876 double sy_old[depos_order + 3] = {0.};
879 double sz_new[depos_order + 3] = {0.};
880 double sz_old[depos_order + 3] = {0.};
888 #if !defined(WARPX_DIM_1D_Z)
889 const int i_new = compute_shape_factor(sx_new+1, x_new);
890 const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
892 #if defined(WARPX_DIM_3D)
893 const int j_new = compute_shape_factor(sy_new+1, y_new);
894 const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
896 const int k_new = compute_shape_factor(sz_new+1, z_new);
897 const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new);
900 #if !defined(WARPX_DIM_1D_Z)
901 int dil = 1, diu = 1;
902 if (i_old < i_new) { dil = 0; }
903 if (i_old > i_new) { diu = 0; }
905 #if defined(WARPX_DIM_3D)
906 int djl = 1, dju = 1;
907 if (j_old < j_new) { djl = 0; }
908 if (j_old > j_new) { dju = 0; }
910 int dkl = 1, dku = 1;
911 if (k_old < k_new) { dkl = 0; }
912 if (k_old > k_new) { dku = 0; }
914 #if defined(WARPX_DIM_3D)
916 for (
int k=dkl; k<=depos_order+2-dku; k++) {
917 for (
int j=djl; j<=depos_order+2-dju; j++) {
918 amrex::Real sdxi = 0._rt;
919 for (
int i=dil;
i<=depos_order+1-diu;
i++) {
920 sdxi += wqx*(sx_old[
i] - sx_new[
i])*(
921 one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k])
922 +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k]));
927 for (
int k=dkl; k<=depos_order+2-dku; k++) {
928 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
929 amrex::Real sdyj = 0._rt;
930 for (
int j=djl; j<=depos_order+1-dju; j++) {
931 sdyj += wqy*(sy_old[j] - sy_new[j])*(
932 one_third*(sx_new[
i]*sz_new[k] + sx_old[
i]*sz_old[k])
933 +one_sixth*(sx_new[
i]*sz_old[k] + sx_old[
i]*sz_new[k]));
938 for (
int j=djl; j<=depos_order+2-dju; j++) {
939 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
940 amrex::Real sdzk = 0._rt;
941 for (
int k=dkl; k<=depos_order+1-dku; k++) {
942 sdzk += wqz*(sz_old[k] - sz_new[k])*(
943 one_third*(sx_new[
i]*sy_new[j] + sx_old[
i]*sy_old[j])
944 +one_sixth*(sx_new[
i]*sy_old[j] + sx_old[
i]*sy_new[j]));
950 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
952 for (
int k=dkl; k<=depos_order+2-dku; k++) {
953 amrex::Real sdxi = 0._rt;
954 for (
int i=dil;
i<=depos_order+1-diu;
i++) {
955 sdxi += wqx*(sx_old[
i] - sx_new[
i])*0.5_rt*(sz_new[k] + sz_old[k]);
957 #if defined(WARPX_DIM_RZ)
959 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
961 const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
964 xy_mid = xy_mid*xy_mid0;
969 for (
int k=dkl; k<=depos_order+2-dku; k++) {
970 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
971 Real
const sdyj = wq*vy*invvol*(
972 one_third*(sx_new[
i]*sz_new[k] + sx_old[
i]*sz_old[k])
973 +one_sixth*(sx_new[
i]*sz_old[k] + sx_old[
i]*sz_new[k]));
975 #if defined(WARPX_DIM_RZ)
980 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
983 const Complex djt_cmplx = -2._rt * I*(i_new-1 +
i +
xmin*dxi)*wq*invdtdx/(amrex::Real)imode
984 *(
Complex(sx_new[
i]*sz_new[k], 0._rt)*(xy_new - xy_mid)
985 +
Complex(sx_old[
i]*sz_old[k], 0._rt)*(xy_mid - xy_old));
988 xy_new = xy_new*xy_new0;
989 xy_mid = xy_mid*xy_mid0;
990 xy_old = xy_old*xy_old0;
995 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
997 for (
int k=dkl; k<=depos_order+1-dku; k++) {
998 sdzk += wqz*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[
i] + sx_old[
i]);
1000 #if defined(WARPX_DIM_RZ)
1002 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1004 const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
1007 xy_mid = xy_mid*xy_mid0;
1012 #elif defined(WARPX_DIM_1D_Z)
1014 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1015 amrex::Real
const sdxi = wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1018 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1019 amrex::Real
const sdyj = wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1022 amrex::Real sdzk = 0._rt;
1023 for (
int k=dkl; k<=depos_order+1-dku; k++) {
1024 sdzk += wqz*(sz_old[k] - sz_new[k]);
1056 template <
int depos_order>
1058 const amrex::ParticleReal *
const yp_n,
1059 const amrex::ParticleReal *
const zp_n,
1061 const amrex::ParticleReal *
const wp,
1062 [[maybe_unused]]
const amrex::ParticleReal *
const uxp_n,
1063 [[maybe_unused]]
const amrex::ParticleReal *
const uyp_n,
1064 [[maybe_unused]]
const amrex::ParticleReal *
const uzp_n,
1065 [[maybe_unused]]
const amrex::ParticleReal *
const uxp_nph,
1066 [[maybe_unused]]
const amrex::ParticleReal *
const uyp_nph,
1067 [[maybe_unused]]
const amrex::ParticleReal *
const uzp_nph,
1068 const int *
const ion_lev,
1072 const long np_to_deposit,
1073 const amrex::Real
dt,
1074 const std::array<amrex::Real, 3>&
dx,
1075 const std::array<amrex::Real, 3> xyzmin,
1077 const amrex::Real q,
1078 const int n_rz_azimuthal_modes)
1080 using namespace amrex;
1081 #if !defined(WARPX_DIM_RZ)
1087 bool const do_ionization = ion_lev;
1088 #if !defined(WARPX_DIM_1D_Z)
1089 Real
const dxi = 1.0_rt /
dx[0];
1091 #if !defined(WARPX_DIM_1D_Z)
1092 Real
const xmin = xyzmin[0];
1094 #if defined(WARPX_DIM_3D)
1095 Real
const dyi = 1.0_rt /
dx[1];
1096 Real
const ymin = xyzmin[1];
1098 Real
const dzi = 1.0_rt /
dx[2];
1099 Real
const zmin = xyzmin[2];
1101 #if defined(WARPX_DIM_3D)
1102 Real
const invdtdx = 1.0_rt / (
dt*
dx[1]*
dx[2]);
1103 Real
const invdtdy = 1.0_rt / (
dt*
dx[0]*
dx[2]);
1104 Real
const invdtdz = 1.0_rt / (
dt*
dx[0]*
dx[1]);
1105 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1106 Real
const invdtdx = 1.0_rt / (
dt*
dx[2]);
1107 Real
const invdtdz = 1.0_rt / (
dt*
dx[0]);
1108 Real
const invvol = 1.0_rt / (
dx[0]*
dx[2]);
1109 #elif defined(WARPX_DIM_1D_Z)
1110 Real
const invdtdz = 1.0_rt / (
dt*
dx[0]);
1111 Real
const invvol = 1.0_rt / (
dx[2]);
1114 #if defined(WARPX_DIM_RZ)
1118 #if !defined(WARPX_DIM_1D_Z)
1119 Real constexpr one_third = 1.0_rt / 3.0_rt;
1120 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
1127 #if !defined(WARPX_DIM_3D)
1132 const amrex::ParticleReal uxp_np1 = 2._prt*uxp_nph[ip] - uxp_n[ip];
1133 const amrex::ParticleReal uyp_np1 = 2._prt*uyp_nph[ip] - uyp_n[ip];
1134 const amrex::ParticleReal uzp_np1 = 2._prt*uzp_nph[ip] - uzp_n[ip];
1135 const amrex::ParticleReal gamma_n = std::sqrt(1._prt + (uxp_n[ip]*uxp_n[ip] + uyp_n[ip]*uyp_n[ip] + uzp_n[ip]*uzp_n[ip])*inv_c2);
1136 const amrex::ParticleReal gamma_np1 = std::sqrt(1._prt + (uxp_np1*uxp_np1 + uyp_np1*uyp_np1 + uzp_np1*uzp_np1)*inv_c2);
1137 const amrex::ParticleReal gaminv = 2.0_prt/(gamma_n + gamma_np1);
1146 ParticleReal xp_nph, yp_nph, zp_nph;
1147 GetPosition(ip, xp_nph, yp_nph, zp_nph);
1149 #if !defined(WARPX_DIM_1D_Z)
1150 ParticleReal
const xp_np1 = 2._prt*xp_nph - xp_n[ip];
1154 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
1155 ParticleReal
const yp_np1 = 2._prt*yp_nph - yp_n[ip];
1159 ParticleReal
const zp_np1 = 2._prt*zp_nph - zp_n[ip];
1161 #if !defined(WARPX_DIM_1D_Z)
1162 amrex::Real
const wqx = wq*invdtdx;
1164 #if defined(WARPX_DIM_3D)
1165 amrex::Real
const wqy = wq*invdtdy;
1167 amrex::Real
const wqz = wq*invdtdz;
1170 #if defined(WARPX_DIM_RZ)
1171 amrex::Real
const xp_new = xp_np1;
1172 amrex::Real
const yp_new = yp_np1;
1173 amrex::Real
const xp_mid = xp_nph;
1174 amrex::Real
const yp_mid = yp_nph;
1175 amrex::Real
const xp_old = xp_n[ip];
1176 amrex::Real
const yp_old = yp_n[ip];
1177 amrex::Real
const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
1178 amrex::Real
const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
1179 amrex::Real
const rp_mid = (rp_new + rp_old)/2._rt;
1180 amrex::Real costheta_new, sintheta_new;
1181 if (rp_new > 0._rt) {
1182 costheta_new = xp_new/rp_new;
1183 sintheta_new = yp_new/rp_new;
1185 costheta_new = 1._rt;
1186 sintheta_new = 0._rt;
1188 amrex::Real costheta_mid, sintheta_mid;
1189 if (rp_mid > 0._rt) {
1190 costheta_mid = xp_mid/rp_mid;
1191 sintheta_mid = yp_mid/rp_mid;
1193 costheta_mid = 1._rt;
1194 sintheta_mid = 0._rt;
1196 amrex::Real costheta_old, sintheta_old;
1197 if (rp_old > 0._rt) {
1198 costheta_old = xp_old/rp_old;
1199 sintheta_old = yp_old/rp_old;
1201 costheta_old = 1._rt;
1202 sintheta_old = 0._rt;
1208 double const x_new = (rp_new -
xmin)*dxi;
1209 double const x_old = (rp_old -
xmin)*dxi;
1211 #if !defined(WARPX_DIM_1D_Z)
1213 double const x_new = (xp_np1 -
xmin)*dxi;
1214 double const x_old = (xp_n[ip] -
xmin)*dxi;
1217 #if defined(WARPX_DIM_3D)
1219 double const y_new = (yp_np1 - ymin)*dyi;
1220 double const y_old = (yp_n[ip] - ymin)*dyi;
1223 double const z_new = (zp_np1 - zmin)*dzi;
1224 double const z_old = (zp_n[ip] - zmin)*dzi;
1226 #if defined(WARPX_DIM_RZ)
1227 amrex::Real
const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
1228 #elif defined(WARPX_DIM_XZ)
1229 amrex::Real
const vy = uyp_nph[ip]*gaminv;
1230 #elif defined(WARPX_DIM_1D_Z)
1231 amrex::Real
const vx = uxp_nph[ip]*gaminv;
1232 amrex::Real
const vy = uyp_nph[ip]*gaminv;
1240 #if !defined(WARPX_DIM_1D_Z)
1241 double sx_new[depos_order + 3] = {0.};
1242 double sx_old[depos_order + 3] = {0.};
1244 #if defined(WARPX_DIM_3D)
1246 double sy_new[depos_order + 3] = {0.};
1247 double sy_old[depos_order + 3] = {0.};
1250 double sz_new[depos_order + 3] = {0.};
1251 double sz_old[depos_order + 3] = {0.};
1259 #if !defined(WARPX_DIM_1D_Z)
1260 const int i_new = compute_shape_factor(sx_new+1, x_new);
1261 const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new);
1263 #if defined(WARPX_DIM_3D)
1264 const int j_new = compute_shape_factor(sy_new+1, y_new);
1265 const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new);
1267 const int k_new = compute_shape_factor(sz_new+1, z_new);
1268 const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new);
1271 #if !defined(WARPX_DIM_1D_Z)
1272 int dil = 1, diu = 1;
1273 if (i_old < i_new) { dil = 0; }
1274 if (i_old > i_new) { diu = 0; }
1276 #if defined(WARPX_DIM_3D)
1277 int djl = 1, dju = 1;
1278 if (j_old < j_new) { djl = 0; }
1279 if (j_old > j_new) { dju = 0; }
1281 int dkl = 1, dku = 1;
1282 if (k_old < k_new) { dkl = 0; }
1283 if (k_old > k_new) { dku = 0; }
1285 #if defined(WARPX_DIM_3D)
1287 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1288 for (
int j=djl; j<=depos_order+2-dju; j++) {
1289 amrex::Real sdxi = 0._rt;
1290 for (
int i=dil;
i<=depos_order+1-diu;
i++) {
1291 sdxi += wqx*(sx_old[
i] - sx_new[
i])*(
1292 one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k])
1293 +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k]));
1298 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1299 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
1300 amrex::Real sdyj = 0._rt;
1301 for (
int j=djl; j<=depos_order+1-dju; j++) {
1302 sdyj += wqy*(sy_old[j] - sy_new[j])*(
1303 one_third*(sx_new[
i]*sz_new[k] + sx_old[
i]*sz_old[k])
1304 +one_sixth*(sx_new[
i]*sz_old[k] + sx_old[
i]*sz_new[k]));
1309 for (
int j=djl; j<=depos_order+2-dju; j++) {
1310 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
1311 amrex::Real sdzk = 0._rt;
1312 for (
int k=dkl; k<=depos_order+1-dku; k++) {
1313 sdzk += wqz*(sz_old[k] - sz_new[k])*(
1314 one_third*(sx_new[
i]*sy_new[j] + sx_old[
i]*sy_old[j])
1315 +one_sixth*(sx_new[
i]*sy_old[j] + sx_old[
i]*sy_new[j]));
1321 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1323 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1324 amrex::Real sdxi = 0._rt;
1325 for (
int i=dil;
i<=depos_order+1-diu;
i++) {
1326 sdxi += wqx*(sx_old[
i] - sx_new[
i])*0.5_rt*(sz_new[k] + sz_old[k]);
1328 #if defined(WARPX_DIM_RZ)
1330 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1332 const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
1335 xy_mid = xy_mid*xy_mid0;
1340 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1341 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
1342 Real
const sdyj = wq*vy*invvol*(
1343 one_third*(sx_new[
i]*sz_new[k] + sx_old[
i]*sz_old[k])
1344 +one_sixth*(sx_new[
i]*sz_old[k] + sx_old[
i]*sz_new[k]));
1346 #if defined(WARPX_DIM_RZ)
1351 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1354 const Complex djt_cmplx = -2._rt * I*(i_new-1 +
i +
xmin*dxi)*wq*invdtdx/(amrex::Real)imode
1355 *(
Complex(sx_new[
i]*sz_new[k], 0._rt)*(xy_new - xy_mid)
1356 +
Complex(sx_old[
i]*sz_old[k], 0._rt)*(xy_mid - xy_old));
1359 xy_new = xy_new*xy_new0;
1360 xy_mid = xy_mid*xy_mid0;
1361 xy_old = xy_old*xy_old0;
1366 for (
int i=dil;
i<=depos_order+2-diu;
i++) {
1368 for (
int k=dkl; k<=depos_order+1-dku; k++) {
1369 sdzk += wqz*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[
i] + sx_old[
i]);
1371 #if defined(WARPX_DIM_RZ)
1373 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1375 const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
1378 xy_mid = xy_mid*xy_mid0;
1383 #elif defined(WARPX_DIM_1D_Z)
1385 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1386 amrex::Real
const sdxi = wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1389 for (
int k=dkl; k<=depos_order+2-dku; k++) {
1390 amrex::Real
const sdyj = wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]);
1393 amrex::Real sdzk = 0._rt;
1394 for (
int k=dkl; k<=depos_order+1-dku; k++) {
1395 sdzk += wqz*(sz_old[k] - sz_new[k]);
1429 template <
int depos_order>
1431 const amrex::ParticleReal *
const yp_n,
1432 const amrex::ParticleReal *
const zp_n,
1434 const amrex::ParticleReal *
const wp,
1435 [[maybe_unused]]
const amrex::ParticleReal *
const uxp_n,
1436 [[maybe_unused]]
const amrex::ParticleReal *
const uyp_n,
1437 [[maybe_unused]]
const amrex::ParticleReal *
const uzp_n,
1438 [[maybe_unused]]
const amrex::ParticleReal *
const uxp_nph,
1439 [[maybe_unused]]
const amrex::ParticleReal *
const uyp_nph,
1440 [[maybe_unused]]
const amrex::ParticleReal *
const uzp_nph,
1441 const int *
const ion_lev,
1445 const long np_to_deposit,
1446 const amrex::Real
dt,
1447 const std::array<amrex::Real, 3>&
dx,
1448 const std::array<amrex::Real, 3> xyzmin,
1450 const amrex::Real q,
1451 const int n_rz_azimuthal_modes)
1453 using namespace amrex;
1454 #if !defined(WARPX_DIM_RZ)
1460 bool const do_ionization = ion_lev;
1461 #if !defined(WARPX_DIM_1D_Z)
1462 Real
const dxi = 1.0_rt /
dx[0];
1464 #if !defined(WARPX_DIM_1D_Z)
1465 Real
const xmin = xyzmin[0];
1467 #if defined(WARPX_DIM_3D)
1468 Real
const dyi = 1.0_rt /
dx[1];
1469 Real
const ymin = xyzmin[1];
1471 Real
const dzi = 1.0_rt /
dx[2];
1472 Real
const zmin = xyzmin[2];
1474 #if defined(WARPX_DIM_3D)
1475 Real
const invvol = 1.0_rt / (
dx[0]*
dx[1]*
dx[2]);
1476 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1477 Real
const invvol = 1.0_rt / (
dx[0]*
dx[2]);
1478 #elif defined(WARPX_DIM_1D_Z)
1479 Real
const invvol = 1.0_rt / (
dx[2]);
1482 #if !defined(WARPX_DIM_1D_Z)
1483 Real constexpr one_third = 1.0_rt / 3.0_rt;
1484 Real constexpr one_sixth = 1.0_rt / 6.0_rt;
1492 #if !defined(WARPX_DIM_3D)
1497 const amrex::ParticleReal uxp_np1 = 2._prt*uxp_nph[ip] - uxp_n[ip];
1498 const amrex::ParticleReal uyp_np1 = 2._prt*uyp_nph[ip] - uyp_n[ip];
1499 const amrex::ParticleReal uzp_np1 = 2._prt*uzp_nph[ip] - uzp_n[ip];
1500 const amrex::ParticleReal gamma_n = std::sqrt(1._prt + (uxp_n[ip]*uxp_n[ip] + uyp_n[ip]*uyp_n[ip] + uzp_n[ip]*uzp_n[ip])*inv_c2);
1501 const amrex::ParticleReal gamma_np1 = std::sqrt(1._prt + (uxp_np1*uxp_np1 + uyp_np1*uyp_np1 + uzp_np1*uzp_np1)*inv_c2);
1502 const amrex::ParticleReal gaminv = 2.0_prt/(gamma_n + gamma_np1);
1511 ParticleReal xp_nph, yp_nph, zp_nph;
1512 GetPosition(ip, xp_nph, yp_nph, zp_nph);
1514 #if !defined(WARPX_DIM_1D_Z)
1515 ParticleReal
const xp_np1 = 2._prt*xp_nph - xp_n[ip];
1519 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ)
1520 ParticleReal
const yp_np1 = 2._prt*yp_nph - yp_n[ip];
1524 ParticleReal
const zp_np1 = 2._prt*zp_nph - zp_n[ip];
1527 #if defined(WARPX_DIM_RZ)
1528 amrex::Real
const xp_new = xp_np1;
1529 amrex::Real
const yp_new = yp_np1;
1530 amrex::Real
const xp_mid = xp_nph;
1531 amrex::Real
const yp_mid = yp_nph;
1532 amrex::Real
const xp_old = xp_n[ip];
1533 amrex::Real
const yp_old = yp_n[ip];
1534 amrex::Real
const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
1535 amrex::Real
const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
1536 amrex::Real
const rp_mid = (rp_new + rp_old)/2._rt;
1537 amrex::Real costheta_mid, sintheta_mid;
1538 if (rp_mid > 0._rt) {
1539 costheta_mid = xp_mid/rp_mid;
1540 sintheta_mid = yp_mid/rp_mid;
1542 costheta_mid = 1._rt;
1543 sintheta_mid = 0._rt;
1548 double const x_new = (rp_new -
xmin)*dxi;
1549 double const x_old = (rp_old -
xmin)*dxi;
1550 amrex::Real
const vx = (rp_new - rp_old)/
dt;
1551 amrex::Real
const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
1552 #elif defined(WARPX_DIM_XZ)
1554 double const x_new = (xp_np1 -
xmin)*dxi;
1555 double const x_old = (xp_n[ip] -
xmin)*dxi;
1556 amrex::Real
const vx = (xp_np1 - xp_n[ip])/
dt;
1557 amrex::Real
const vy = uyp_nph[ip]*gaminv;
1558 #elif defined(WARPX_DIM_1D_Z)
1559 amrex::Real
const vx = uxp_nph[ip]*gaminv;
1560 amrex::Real
const vy = uyp_nph[ip]*gaminv;
1561 #elif defined(WARPX_DIM_3D)
1563 double const x_new = (xp_np1 -
xmin)*dxi;
1564 double const x_old = (xp_n[ip] -
xmin)*dxi;
1565 double const y_new = (yp_np1 - ymin)*dyi;
1566 double const y_old = (yp_n[ip] - ymin)*dyi;
1567 amrex::Real
const vx = (xp_np1 - xp_n[ip])/
dt;
1568 amrex::Real
const vy = (yp_np1 - yp_n[ip])/
dt;
1572 double const z_new = (zp_np1 - zmin)*dzi;
1573 double const z_old = (zp_n[ip] - zmin)*dzi;
1574 amrex::Real
const vz = (zp_np1 - zp_n[ip])/
dt;
1577 amrex::Real
const wqx = wq*vx*invvol;
1578 amrex::Real
const wqy = wq*vy*invvol;
1579 amrex::Real
const wqz = wq*vz*invvol;
1587 int num_segments = 1;
1589 if ( (depos_order % 2) == 0 ) {
shift = 0.5; }
1591 #if defined(WARPX_DIM_3D)
1594 const auto i_old =
static_cast<int>(x_old-
shift);
1595 const auto i_new =
static_cast<int>(x_new-
shift);
1596 const int cell_crossings_x = std::abs(i_new-i_old);
1597 num_segments += cell_crossings_x;
1600 const auto j_old =
static_cast<int>(y_old-
shift);
1601 const auto j_new =
static_cast<int>(y_new-
shift);
1602 const int cell_crossings_y = std::abs(j_new-j_old);
1603 num_segments += cell_crossings_y;
1606 const auto k_old =
static_cast<int>(z_old-
shift);
1607 const auto k_new =
static_cast<int>(z_new-
shift);
1608 const int cell_crossings_z = std::abs(k_new-k_old);
1609 num_segments += cell_crossings_z;
1617 const double dxp = x_new - x_old;
1618 const double dyp = y_new - y_old;
1619 const double dzp = z_new - z_old;
1620 const auto dirX_sign =
static_cast<double>(dxp < 0. ? -1. : 1.);
1621 const auto dirY_sign =
static_cast<double>(dyp < 0. ? -1. : 1.);
1622 const auto dirZ_sign =
static_cast<double>(dzp < 0. ? -1. : 1.);
1623 double Xcell = 0., Ycell = 0., Zcell = 0.;
1624 if (num_segments > 1) {
1625 Xcell =
static_cast<double>(i_old) +
shift + 0.5*(1.-dirX_sign);
1626 Ycell =
static_cast<double>(j_old) +
shift + 0.5*(1.-dirY_sign);
1627 Zcell =
static_cast<double>(k_old) +
shift + 0.5*(1.-dirZ_sign);
1633 double dxp_seg, dyp_seg, dzp_seg;
1634 double x0_new, y0_new, z0_new;
1635 double x0_old = x_old;
1636 double y0_old = y_old;
1637 double z0_old = z_old;
1639 for (
int ns=0; ns<num_segments; ns++) {
1641 if (ns == num_segments-1) {
1646 dxp_seg = x0_new - x0_old;
1647 dyp_seg = y0_new - y0_old;
1648 dzp_seg = z0_new - z0_old;
1653 x0_new = Xcell + dirX_sign;
1654 y0_new = Ycell + dirY_sign;
1655 z0_new = Zcell + dirZ_sign;
1656 dxp_seg = x0_new - x0_old;
1657 dyp_seg = y0_new - y0_old;
1658 dzp_seg = z0_new - z0_old;
1660 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
1661 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
1663 dyp_seg = dyp/dxp*dxp_seg;
1664 dzp_seg = dzp/dxp*dxp_seg;
1665 y0_new = y0_old + dyp_seg;
1666 z0_new = z0_old + dzp_seg;
1668 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
1670 dxp_seg = dxp/dyp*dyp_seg;
1671 dzp_seg = dzp/dyp*dyp_seg;
1672 x0_new = x0_old + dxp_seg;
1673 z0_new = z0_old + dzp_seg;
1677 dxp_seg = dxp/dzp*dzp_seg;
1678 dyp_seg = dyp/dzp*dzp_seg;
1679 x0_new = x0_old + dxp_seg;
1680 y0_new = y0_old + dyp_seg;
1686 const auto seg_factor_x =
static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1687 const auto seg_factor_y =
static_cast<double>(dyp == 0. ? 1. : dyp_seg/dyp);
1688 const auto seg_factor_z =
static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1691 double sx_cell[depos_order] = {0.};
1692 double sy_cell[depos_order] = {0.};
1693 double sz_cell[depos_order] = {0.};
1694 double const x0_bar = (x0_new + x0_old)/2.0;
1695 double const y0_bar = (y0_new + y0_old)/2.0;
1696 double const z0_bar = (z0_new + z0_old)/2.0;
1697 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1698 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
1699 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1701 if constexpr (depos_order >= 3) {
1703 double sx_old_cell[depos_order] = {0.};
1704 double sx_new_cell[depos_order] = {0.};
1705 double sy_old_cell[depos_order] = {0.};
1706 double sy_new_cell[depos_order] = {0.};
1707 double sz_old_cell[depos_order] = {0.};
1708 double sz_new_cell[depos_order] = {0.};
1709 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1710 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
1711 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1713 for (
int m=0; m<depos_order; m++) {
1714 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1715 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
1716 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1721 double sx_old_node[depos_order+1] = {0.};
1722 double sx_new_node[depos_order+1] = {0.};
1723 double sy_old_node[depos_order+1] = {0.};
1724 double sy_new_node[depos_order+1] = {0.};
1725 double sz_old_node[depos_order+1] = {0.};
1726 double sz_new_node[depos_order+1] = {0.};
1727 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1728 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
1729 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1732 amrex::Real this_Jx;
1733 for (
int i=0;
i<=depos_order-1;
i++) {
1734 for (
int j=0; j<=depos_order; j++) {
1735 for (
int k=0; k<=depos_order; k++) {
1736 this_Jx = wqx*sx_cell[
i]*( sy_old_node[j]*sz_old_node[k]*one_third
1737 + sy_old_node[j]*sz_new_node[k]*one_sixth
1738 + sy_new_node[j]*sz_old_node[k]*one_sixth
1739 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
1746 amrex::Real this_Jy;
1747 for (
int i=0;
i<=depos_order;
i++) {
1748 for (
int j=0; j<=depos_order-1; j++) {
1749 for (
int k=0; k<=depos_order; k++) {
1750 this_Jy = wqy*sy_cell[j]*( sx_old_node[
i]*sz_old_node[k]*one_third
1751 + sx_old_node[
i]*sz_new_node[k]*one_sixth
1752 + sx_new_node[
i]*sz_old_node[k]*one_sixth
1753 + sx_new_node[
i]*sz_new_node[k]*one_third )*seg_factor_y;
1760 amrex::Real this_Jz;
1761 for (
int i=0;
i<=depos_order;
i++) {
1762 for (
int j=0; j<=depos_order; j++) {
1763 for (
int k=0; k<=depos_order-1; k++) {
1764 this_Jz = wqz*sz_cell[k]*( sx_old_node[
i]*sy_old_node[j]*one_third
1765 + sx_old_node[
i]*sy_new_node[j]*one_sixth
1766 + sx_new_node[
i]*sy_old_node[j]*one_sixth
1767 + sx_new_node[
i]*sy_new_node[j]*one_third )*seg_factor_z;
1774 if (ns < num_segments-1) {
1782 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1785 const auto i_old =
static_cast<int>(x_old-
shift);
1786 const auto i_new =
static_cast<int>(x_new-
shift);
1787 const int cell_crossings_x = std::abs(i_new-i_old);
1788 num_segments += cell_crossings_x;
1791 const auto k_old =
static_cast<int>(z_old-
shift);
1792 const auto k_new =
static_cast<int>(z_new-
shift);
1793 const int cell_crossings_z = std::abs(k_new-k_old);
1794 num_segments += cell_crossings_z;
1802 const double dxp = x_new - x_old;
1803 const double dzp = z_new - z_old;
1804 const auto dirX_sign =
static_cast<double>(dxp < 0. ? -1. : 1.);
1805 const auto dirZ_sign =
static_cast<double>(dzp < 0. ? -1. : 1.);
1806 double Xcell = 0., Zcell = 0.;
1807 if (num_segments > 1) {
1808 Xcell =
static_cast<double>(i_old) +
shift + 0.5*(1.-dirX_sign);
1809 Zcell =
static_cast<double>(k_old) +
shift + 0.5*(1.-dirZ_sign);
1815 double dxp_seg, dzp_seg;
1816 double x0_new, z0_new;
1817 double x0_old = x_old;
1818 double z0_old = z_old;
1820 for (
int ns=0; ns<num_segments; ns++) {
1822 if (ns == num_segments-1) {
1826 dxp_seg = x0_new - x0_old;
1827 dzp_seg = z0_new - z0_old;
1832 x0_new = Xcell + dirX_sign;
1833 z0_new = Zcell + dirZ_sign;
1834 dxp_seg = x0_new - x0_old;
1835 dzp_seg = z0_new - z0_old;
1837 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1839 dzp_seg = dzp/dxp*dxp_seg;
1840 z0_new = z0_old + dzp_seg;
1844 dxp_seg = dxp/dzp*dzp_seg;
1845 x0_new = x0_old + dxp_seg;
1851 const auto seg_factor_x =
static_cast<double>(dxp == 0. ? 1. : dxp_seg/dxp);
1852 const auto seg_factor_z =
static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1855 double sx_cell[depos_order] = {0.};
1856 double sz_cell[depos_order] = {0.};
1857 double const x0_bar = (x0_new + x0_old)/2.0;
1858 double const z0_bar = (z0_new + z0_old)/2.0;
1859 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1860 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1862 if constexpr (depos_order >= 3) {
1864 double sx_old_cell[depos_order] = {0.};
1865 double sx_new_cell[depos_order] = {0.};
1866 double sz_old_cell[depos_order] = {0.};
1867 double sz_new_cell[depos_order] = {0.};
1868 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1869 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1871 for (
int m=0; m<depos_order; m++) {
1872 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1873 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1878 double sx_old_node[depos_order+1] = {0.};
1879 double sx_new_node[depos_order+1] = {0.};
1880 double sz_old_node[depos_order+1] = {0.};
1881 double sz_new_node[depos_order+1] = {0.};
1882 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1883 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1886 amrex::Real this_Jx;
1887 for (
int i=0;
i<=depos_order-1;
i++) {
1888 for (
int k=0; k<=depos_order; k++) {
1889 this_Jx = wqx*sx_cell[
i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1891 #if defined(WARPX_DIM_RZ)
1893 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1895 const Complex djr_cmplx = 2._rt*this_Jx*xy_mid;
1898 xy_mid = xy_mid*xy_mid0;
1905 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1906 amrex::Real this_Jy;
1907 for (
int i=0;
i<=depos_order;
i++) {
1908 for (
int k=0; k<=depos_order; k++) {
1909 this_Jy = wqy*( sx_old_node[
i]*sz_old_node[k]*one_third
1910 + sx_old_node[
i]*sz_new_node[k]*one_sixth
1911 + sx_new_node[
i]*sz_old_node[k]*one_sixth
1912 + sx_new_node[
i]*sz_new_node[k]*one_third )*seg_factor_y;
1914 #if defined(WARPX_DIM_RZ)
1917 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1919 const Complex djy_cmplx = 2._rt*this_Jy*xy_mid;
1922 xy_mid = xy_mid*xy_mid0;
1929 amrex::Real this_Jz;
1930 for (
int i=0;
i<=depos_order;
i++) {
1931 for (
int k=0; k<=depos_order-1; k++) {
1932 this_Jz = wqz*sz_cell[k]*(sx_old_node[
i] + sx_new_node[
i])/2.0_rt*seg_factor_z;
1934 #if defined(WARPX_DIM_RZ)
1936 for (
int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
1938 const Complex djz_cmplx = 2._rt*this_Jz*xy_mid;
1941 xy_mid = xy_mid*xy_mid0;
1948 if (ns < num_segments-1) {
1955 #elif defined(WARPX_DIM_1D_Z)
1958 const auto k_old =
static_cast<int>(z_old-
shift);
1959 const auto k_new =
static_cast<int>(z_new-
shift);
1960 const int cell_crossings_z = std::abs(k_new-k_old);
1961 num_segments += cell_crossings_z;
1968 double const dzp = z_new - z_old;
1969 const auto dirZ_sign =
static_cast<double>(dzp < 0. ? -1. : 1.);
1970 double Zcell =
static_cast<double>(k_old) +
shift + 0.5*(1.-dirZ_sign);
1977 double z0_old = z_old;
1979 for (
int ns=0; ns<num_segments; ns++) {
1981 if (ns == num_segments-1) {
1983 dzp_seg = z0_new - z0_old;
1986 Zcell = Zcell + dirZ_sign;
1988 dzp_seg = z0_new - z0_old;
1992 const auto seg_factor =
static_cast<double>(dzp == 0. ? 1. : dzp_seg/dzp);
1995 double sz_cell[depos_order] = {0.};
1996 double const z0_bar = (z0_new + z0_old)/2.0;
1997 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1999 if constexpr (depos_order >= 3) {
2001 double sz_old_cell[depos_order] = {0.};
2002 double sz_new_cell[depos_order] = {0.};
2003 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
2005 for (
int m=0; m<depos_order; m++) {
2006 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
2011 double sz_old_node[depos_order+1] = {0.};
2012 double sz_new_node[depos_order+1] = {0.};
2013 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
2016 for (
int k=0; k<=depos_order; k++) {
2017 const amrex::Real weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
2023 for (
int k=0; k<=depos_order-1; k++) {
2024 const amrex::Real this_Jz = wqz*sz_cell[k]*seg_factor;
2029 if (ns < num_segments-1) {
2066 template <
int depos_order>
2068 const amrex::ParticleReal*
const wp,
2069 const amrex::ParticleReal*
const uxp,
2070 const amrex::ParticleReal*
const uyp,
2071 const amrex::ParticleReal*
const uzp,
2072 const int*
const ion_lev,
2078 amrex::Real relative_time,
2079 const std::array<amrex::Real,3>&
dx,
2080 const std::array<amrex::Real,3>& xyzmin,
2083 int n_rz_azimuthal_modes)
2085 using namespace amrex::literals;
2087 #if defined(WARPX_DIM_RZ)
2089 wp, uxp, uyp, uzp, ion_lev, Dx_fab, Dy_fab, Dz_fab,
2090 np_to_deposit,
dt, relative_time,
dx, xyzmin, lo, q, n_rz_azimuthal_modes);
2094 #if defined(WARPX_DIM_1D_Z)
2096 wp, uxp, uyp, uzp, ion_lev, Dx_fab, Dy_fab, Dz_fab,
2097 np_to_deposit,
dt, relative_time,
dx, xyzmin, lo, q, n_rz_azimuthal_modes);
2101 #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z)
2105 const bool do_ionization = ion_lev;
2108 const amrex::Real dxi = 1._rt /
dx[0];
2109 const amrex::Real dzi = 1._rt /
dx[2];
2110 #if defined(WARPX_DIM_3D)
2111 const amrex::Real dyi = 1._rt /
dx[1];
2115 const amrex::Real invdt = 1._rt /
dt;
2118 #if defined(WARPX_DIM_XZ)
2119 const amrex::Real invvol = dxi * dzi;
2120 #elif defined(WARPX_DIM_3D)
2121 const amrex::Real invvol = dxi * dyi * dzi;
2125 const amrex::Real
xmin = xyzmin[0];
2126 const amrex::Real zmin = xyzmin[2];
2127 #if defined(WARPX_DIM_3D)
2128 const amrex::Real ymin = xyzmin[1];
2132 #if defined(WARPX_DIM_3D)
2135 #elif defined(WARPX_DIM_XZ)
2154 const amrex::Real invgam = 1._rt / std::sqrt(1._rt + uxp[ip] * uxp[ip] * invcsq
2155 + uyp[ip] * uyp[ip] * invcsq
2156 + uzp[ip] * uzp[ip] * invcsq);
2158 amrex::Real wq = q * wp[ip];
2159 if (do_ionization) { wq *= ion_lev[ip]; }
2162 amrex::ParticleReal xp, yp, zp;
2163 GetPosition(ip, xp, yp, zp);
2166 const amrex::Real vx = uxp[ip] * invgam;
2167 const amrex::Real vy = uyp[ip] * invgam;
2168 const amrex::Real vz = uzp[ip] * invgam;
2171 xp += relative_time * vx;
2172 yp += relative_time * vy;
2173 zp += relative_time * vz;
2176 #if defined(WARPX_DIM_XZ)
2177 const amrex::Real wqy = wq * vy * invvol;
2182 double const x_new = (xp -
xmin + 0.5_rt*
dt*vx) * dxi;
2183 double const x_old = (xp -
xmin - 0.5_rt*
dt*vx) * dxi;
2184 #if defined(WARPX_DIM_3D)
2186 double const y_new = (yp - ymin + 0.5_rt*
dt*vy) * dyi;
2187 double const y_old = (yp - ymin - 0.5_rt*
dt*vy) * dyi;
2190 double const z_new = (zp - zmin + 0.5_rt*
dt*vz) * dzi;
2191 double const z_old = (zp - zmin - 0.5_rt*
dt*vz) * dzi;
2195 double sx_new[depos_order+1] = {0.};
2196 double sx_old[depos_order+1] = {0.};
2197 #if defined(WARPX_DIM_3D)
2199 double sy_new[depos_order+1] = {0.};
2200 double sy_old[depos_order+1] = {0.};
2203 double sz_new[depos_order+1] = {0.};
2204 double sz_old[depos_order+1] = {0.};
2211 const int i_new = compute_shape_factor(sx_new, x_new);
2212 #if defined(WARPX_DIM_3D)
2215 const int j_new = compute_shape_factor(sy_new, y_new);
2219 const int k_new = compute_shape_factor(sz_new, z_new);
2225 const int i_old = compute_shape_factor(sx_old, x_old);
2226 #if defined(WARPX_DIM_3D)
2229 const int j_old = compute_shape_factor(sy_old, y_old);
2233 const int k_old = compute_shape_factor(sz_old, z_old);
2236 #if defined(WARPX_DIM_XZ)
2238 for (
int k=0; k<=depos_order; k++) {
2239 for (
int i=0;
i<=depos_order;
i++) {
2243 auto const sxn_szn =
static_cast<amrex::Real
>(sx_new[
i] * sz_new[k]);
2244 auto const sxo_szn =
static_cast<amrex::Real
>(sx_old[
i] * sz_new[k]);
2245 auto const sxn_szo =
static_cast<amrex::Real
>(sx_new[
i] * sz_old[k]);
2246 auto const sxo_szo =
static_cast<amrex::Real
>(sx_old[
i] * sz_old[k]);
2248 if (i_new == i_old && k_new == k_old) {
2251 wq * invvol * invdt * (sxn_szn - sxo_szo));
2254 wq * invvol * invdt * (sxn_szo - sxo_szn));
2258 wqy * 0.25_rt * (sxn_szn + sxn_szo + sxo_szn + sxo_szo));
2262 wq * invvol * invdt * sxn_szn);
2265 - wq * invvol * invdt * sxo_szo);
2268 wq * invvol * invdt * sxn_szo);
2271 - wq * invvol * invdt * sxo_szn);
2275 wqy * 0.25_rt * sxn_szn);
2278 wqy * 0.25_rt * sxn_szo);
2281 wqy * 0.25_rt * sxo_szn);
2284 wqy * 0.25_rt * sxo_szo);
2290 #elif defined(WARPX_DIM_3D)
2292 for (
int k=0; k<=depos_order; k++) {
2293 for (
int j=0; j<=depos_order; j++) {
2295 auto const syn_szn =
static_cast<amrex::Real
>(sy_new[j] * sz_new[k]);
2296 auto const syo_szn =
static_cast<amrex::Real
>(sy_old[j] * sz_new[k]);
2297 auto const syn_szo =
static_cast<amrex::Real
>(sy_new[j] * sz_old[k]);
2298 auto const syo_szo =
static_cast<amrex::Real
>(sy_old[j] * sz_old[k]);
2300 for (
int i=0;
i<=depos_order;
i++) {
2302 auto const sxn_syn_szn =
static_cast<amrex::Real
>(sx_new[
i]) * syn_szn;
2303 auto const sxo_syn_szn =
static_cast<amrex::Real
>(sx_old[
i]) * syn_szn;
2304 auto const sxn_syo_szn =
static_cast<amrex::Real
>(sx_new[
i]) * syo_szn;
2305 auto const sxo_syo_szn =
static_cast<amrex::Real
>(sx_old[
i]) * syo_szn;
2306 auto const sxn_syn_szo =
static_cast<amrex::Real
>(sx_new[
i]) * syn_szo;
2307 auto const sxo_syn_szo =
static_cast<amrex::Real
>(sx_old[
i]) * syn_szo;
2308 auto const sxn_syo_szo =
static_cast<amrex::Real
>(sx_new[
i]) * syo_szo;
2309 auto const sxo_syo_szo =
static_cast<amrex::Real
>(sx_old[
i]) * syo_szo;
2311 if (i_new == i_old && j_new == j_old && k_new == k_old) {
2314 wq * invvol * invdt * (sxn_syn_szn - sxo_syo_szo));
2317 wq * invvol * invdt * (sxn_syn_szo - sxo_syo_szn));
2320 wq * invvol * invdt * (sxn_syo_szn - sxo_syn_szo));
2323 wq * invvol * invdt * (sxo_syn_szn - sxn_syo_szo));
2327 wq * invvol * invdt * sxn_syn_szn);
2330 - wq * invvol * invdt * sxo_syo_szo);
2333 wq * invvol * invdt * sxn_syn_szo);
2336 - wq * invvol * invdt * sxo_syo_szn);
2339 wq * invvol * invdt * sxn_syo_szn);
2342 - wq * invvol * invdt * sxo_syn_szo);
2345 wq * invvol * invdt * sxo_syn_szn);
2348 - wq * invvol * invdt * sxn_syo_szo);
2356 #if defined(WARPX_DIM_3D)
2359 const amrex::Real t_a = temp_arr(
i,j,k,0);
2360 const amrex::Real t_b = temp_arr(
i,j,k,1);
2361 const amrex::Real t_c = temp_arr(
i,j,k,2);
2362 const amrex::Real t_d = temp_arr(
i,j,k,3);
2363 Dx_arr(
i,j,k) += (1._rt/6._rt)*(2_rt*t_a + t_b + t_c - 2._rt*t_d);
2364 Dy_arr(
i,j,k) += (1._rt/6._rt)*(2_rt*t_a + t_b - 2._rt*t_c + t_d);
2365 Dz_arr(
i,j,k) += (1._rt/6._rt)*(2_rt*t_a - 2._rt*t_b + t_c + t_d);
2367 #elif defined(WARPX_DIM_XZ)
2370 const amrex::Real t_a = temp_arr(
i,j,0,0);
2371 const amrex::Real t_b = temp_arr(
i,j,0,1);
2372 Dx_arr(
i,j,0) += (0.5_rt)*(t_a + t_b);
2373 Dz_arr(
i,j,0) += (0.5_rt)*(t_a - t_b);
#define AMREX_ALWAYS_ASSERT(EX)
#define AMREX_GPU_HOST_DEVICE
#define AMREX_D_DECL(a, b, c)
AMREX_GPU_HOST_DEVICE AMREX_INLINE void doDepositionShapeNKernel(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, const amrex::ParticleReal wq, const amrex::ParticleReal vx, const amrex::ParticleReal vy, const amrex::ParticleReal vz, amrex::Array4< amrex::Real > const &jx_arr, amrex::Array4< amrex::Real > const &jy_arr, amrex::Array4< amrex::Real > const &jz_arr, amrex::IntVect const &jx_type, amrex::IntVect const &jy_type, amrex::IntVect const &jz_type, const amrex::Real relative_time, AMREX_D_DECL(const amrex::Real dzi, const amrex::Real dxi, const amrex::Real dyi), AMREX_D_DECL(const amrex::Real zmin, const amrex::Real xmin, const amrex::Real ymin), const amrex::Real invvol, const amrex::Dim3 lo, const int n_rz_azimuthal_modes)
Kernel for the direct current deposition for thread thread_num.
Definition: CurrentDeposition.H:49
void doDepositionShapeNImplicit(const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, const amrex::ParticleReal *const uxp_n, const amrex::ParticleReal *const uyp_n, const amrex::ParticleReal *const uzp_n, 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_deposit, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, const amrex::Dim3 lo, const amrex::Real q, const int n_rz_azimuthal_modes)
Direct current deposition for thread thread_num for the implicit scheme The only difference from doDe...
Definition: CurrentDeposition.H:401
void doEsirkepovDepositionShapeN(const GetParticlePosition< PIdx > &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, long np_to_deposit, amrex::Real dt, amrex::Real relative_time, const std::array< amrex::Real, 3 > &dx, std::array< amrex::Real, 3 > xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes)
Esirkepov Current Deposition for thread thread_num.
Definition: CurrentDeposition.H:708
void doVayDepositionShapeN(const GetParticlePosition< PIdx > &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 &Dx_fab, amrex::FArrayBox &Dy_fab, amrex::FArrayBox &Dz_fab, long np_to_deposit, amrex::Real dt, amrex::Real relative_time, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, amrex::Dim3 lo, amrex::Real q, int 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:2067
void doDepositionSharedShapeN(const GetParticlePosition< PIdx > &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, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, long np_to_deposit, const amrex::Real relative_time, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes, const amrex::DenseBins< WarpXParticleContainer::ParticleTileType::ParticleTileDataType > &a_bins, const amrex::Box &box, const amrex::Geometry &geom, const amrex::IntVect &a_tbox_max_size)
Current Deposition for thread thread_num using shared memory.
Definition: CurrentDeposition.H:519
void doDepositionShapeN(const GetParticlePosition< PIdx > &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, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, long np_to_deposit, amrex::Real relative_time, const std::array< amrex::Real, 3 > &dx, const std::array< amrex::Real, 3 > &xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes)
Current Deposition for thread thread_num.
Definition: CurrentDeposition.H:292
void doVillasenorDepositionShapeNImplicit(const amrex::ParticleReal *const xp_n, const amrex::ParticleReal *const yp_n, const amrex::ParticleReal *const zp_n, const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, [[maybe_unused]]const amrex::ParticleReal *const uxp_n, [[maybe_unused]]const amrex::ParticleReal *const uyp_n, [[maybe_unused]]const amrex::ParticleReal *const uzp_n, [[maybe_unused]]const amrex::ParticleReal *const uxp_nph, [[maybe_unused]]const amrex::ParticleReal *const uyp_nph, [[maybe_unused]]const amrex::ParticleReal *const uzp_nph, const int *const 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_deposit, 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 int n_rz_azimuthal_modes)
Villasenor and Buneman Current Deposition for thread thread_num for implicit scheme....
Definition: CurrentDeposition.H:1430
void doChargeConservingDepositionShapeNImplicit(const amrex::ParticleReal *const xp_n, const amrex::ParticleReal *const yp_n, const amrex::ParticleReal *const zp_n, const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *const wp, [[maybe_unused]]const amrex::ParticleReal *const uxp_n, [[maybe_unused]]const amrex::ParticleReal *const uyp_n, [[maybe_unused]]const amrex::ParticleReal *const uzp_n, [[maybe_unused]]const amrex::ParticleReal *const uxp_nph, [[maybe_unused]]const amrex::ParticleReal *const uyp_nph, [[maybe_unused]]const amrex::ParticleReal *const uzp_nph, const int *const 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_deposit, 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 int n_rz_azimuthal_modes)
Esirkepov Current Deposition for thread thread_num for implicit scheme The difference from doEsirkepo...
Definition: CurrentDeposition.H:1057
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition: TextMsg.H:15
#define WARPX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition: TextMsg.H:13
amrex::GpuComplex< amrex::Real > Complex
Definition: WarpX_Complex.H:22
static amrex::IntVect shared_tilesize
tileSize to use for shared current deposition operations
Definition: WarpX.H:237
static int shared_mem_current_tpb
number of threads to use per block in shared deposition
Definition: WarpX.H:234
AMREX_FORCE_INLINE Array4< Real const > array() const noexcept
const Box & box() const noexcept
AMREX_GPU_HOST_DEVICE Long numPts() const noexcept
AMREX_GPU_HOST_DEVICE Box & grow(int i) noexcept
GpuArray< Real, AMREX_SPACEDIM > InvCellSizeArray() const noexcept
index_type * offsetsPtr() noexcept
index_type * permutationPtr() noexcept
Long numBins() const noexcept
GpuArray< Real, AMREX_SPACEDIM > ProbLoArray() const noexcept
const Box & Domain() const noexcept
static std::size_t sharedMemPerBlock() noexcept
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
void streamSynchronize() noexcept
gpuStream_t gpuStream() noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 end(Box const &box) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box convert(const Box &b, const IntVect &typ) noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 begin(Box const &box) noexcept
void launch(T const &n, L &&f) 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)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int getTileIndex(const IntVect &iv, const Box &box, const bool a_do_tiling, const IntVect &a_tile_size, Box &tbx)
i
Definition: check_interp_points_and_weights.py:174
float dt
Definition: stencil.py:442
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
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T real() const noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T imag() const noexcept
AMREX_GPU_DEVICE T * dataPtr() noexcept