aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles/Deposition/CurrentDeposition.H
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles/Deposition/CurrentDeposition.H')
-rw-r--r--Source/Particles/Deposition/CurrentDeposition.H121
1 files changed, 63 insertions, 58 deletions
diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H
index 6da0f1155..2737eb008 100644
--- a/Source/Particles/Deposition/CurrentDeposition.H
+++ b/Source/Particles/Deposition/CurrentDeposition.H
@@ -4,6 +4,9 @@
#include "ShapeFactors.H"
#include <WarpX_Complex.H>
+#include <AMReX_Array4.H>
+#include <AMReX_REAL.H>
+
/* \brief Current Deposition for thread thread_num
* /param xp, yp, zp : Pointer to arrays of particle positions.
* \param wp : Pointer to array of particle weights.
@@ -208,69 +211,71 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
const amrex::Real q,
const long n_rz_azimuthal_modes)
{
+ using namespace amrex;
+
// Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
// (do_ionization=1)
- const bool do_ionization = ion_lev;
- const amrex::Real dxi = 1.0/dx[0];
- const amrex::Real dtsdx0 = dt*dxi;
- const amrex::Real xmin = xyzmin[0];
+ bool const do_ionization = ion_lev;
+ Real const dxi = 1.0_rt / dx[0];
+ Real const dtsdx0 = dt*dxi;
+ Real const xmin = xyzmin[0];
#if (defined WARPX_DIM_3D)
- const amrex::Real dyi = 1.0/dx[1];
- const amrex::Real dtsdy0 = dt*dyi;
- const amrex::Real ymin = xyzmin[1];
+ Real const dyi = 1.0_rt / dx[1];
+ Real const dtsdy0 = dt*dyi;
+ Real const ymin = xyzmin[1];
#endif
- const amrex::Real dzi = 1.0/dx[2];
- const amrex::Real dtsdz0 = dt*dzi;
- const amrex::Real zmin = xyzmin[2];
+ Real const dzi = 1.0_rt / dx[2];
+ Real const dtsdz0 = dt*dzi;
+ Real const zmin = xyzmin[2];
#if (defined WARPX_DIM_3D)
- const amrex::Real invdtdx = 1.0/(dt*dx[1]*dx[2]);
- const amrex::Real invdtdy = 1.0/(dt*dx[0]*dx[2]);
- const amrex::Real invdtdz = 1.0/(dt*dx[0]*dx[1]);
+ Real const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]);
+ Real const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]);
+ Real const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]);
#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
- const amrex::Real invdtdx = 1.0/(dt*dx[2]);
- const amrex::Real invdtdz = 1.0/(dt*dx[0]);
- const amrex::Real invvol = 1.0/(dx[0]*dx[2]);
+ Real const invdtdx = 1.0_rt / (dt*dx[2]);
+ Real const invdtdz = 1.0_rt / (dt*dx[0]);
+ Real const invvol = 1.0_rt / (dx[0]*dx[2]);
#endif
#if (defined WARPX_DIM_RZ)
- const Complex I = Complex{0., 1.};
+ Complex const I = Complex{0., 1.};
#endif
- const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c;
+ Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c );
// Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
amrex::ParallelFor(
np_to_depose,
- [=] AMREX_GPU_DEVICE (long ip) {
+ [=] AMREX_GPU_DEVICE (long const ip) {
// --- Get particle quantities
- const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
+ Real const gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
+ uyp[ip]*uyp[ip]*clightsq
+ uzp[ip]*uzp[ip]*clightsq);
// wqx, wqy wqz are particle current in each direction
- amrex::Real wq = q*wp[ip];
+ Real wq = q*wp[ip];
if (do_ionization){
wq *= ion_lev[ip];
}
- const amrex::Real wqx = wq*invdtdx;
+ Real const wqx = wq*invdtdx;
#if (defined WARPX_DIM_3D)
- const amrex::Real wqy = wq*invdtdy;
+ Real const wqy = wq*invdtdy;
#endif
- const amrex::Real wqz = wq*invdtdz;
+ Real const wqz = wq*invdtdz;
// computes current and old position in grid units
#if (defined WARPX_DIM_RZ)
- const amrex::Real xp_mid = xp[ip] - 0.5*dt*uxp[ip]*gaminv;
- const amrex::Real yp_mid = yp[ip] - 0.5*dt*uyp[ip]*gaminv;
- const amrex::Real xp_old = xp[ip] - dt*uxp[ip]*gaminv;
- const amrex::Real yp_old = yp[ip] - dt*uyp[ip]*gaminv;
- const amrex::Real rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
- const amrex::Real rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
- const amrex::Real rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
- amrex::Real costheta_new, sintheta_new;
- if (rp_new > 0.) {
+ Real const xp_mid = xp[ip] - 0.5_rt * dt*uxp[ip]*gaminv;
+ Real const yp_mid = yp[ip] - 0.5_rt * dt*uyp[ip]*gaminv;
+ Real const xp_old = xp[ip] - dt*uxp[ip]*gaminv;
+ Real const yp_old = yp[ip] - dt*uyp[ip]*gaminv;
+ Real const rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
+ Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
+ Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
+ Real costheta_new, sintheta_new;
+ if (rp_new > 0._rt) {
costheta_new = xp[ip]/rp_new;
sintheta_new = yp[ip]/rp_new;
} else {
@@ -278,7 +283,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
sintheta_new = 0.;
}
amrex::Real costheta_mid, sintheta_mid;
- if (rp_mid > 0.) {
+ if (rp_mid > 0._rt) {
costheta_mid = xp_mid/rp_mid;
sintheta_mid = yp_mid/rp_mid;
} else {
@@ -286,7 +291,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
sintheta_mid = 0.;
}
amrex::Real costheta_old, sintheta_old;
- if (rp_old > 0.) {
+ if (rp_old > 0._rt) {
costheta_old = xp_old/rp_old;
sintheta_old = yp_old/rp_old;
} else {
@@ -296,37 +301,37 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
const Complex xy_new0 = Complex{costheta_new, sintheta_new};
const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
const Complex xy_old0 = Complex{costheta_old, sintheta_old};
- const amrex::Real x_new = (rp_new - xmin)*dxi;
- const amrex::Real x_old = (rp_old - xmin)*dxi;
+ Real const x_new = (rp_new - xmin)*dxi;
+ Real const x_old = (rp_old - xmin)*dxi;
#else
- const amrex::Real x_new = (xp[ip] - xmin)*dxi;
- const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv;
+ Real const x_new = (xp[ip] - xmin)*dxi;
+ Real const x_old = x_new - dtsdx0*uxp[ip]*gaminv;
#endif
#if (defined WARPX_DIM_3D)
- const amrex::Real y_new = (yp[ip] - ymin)*dyi;
- const amrex::Real y_old = y_new - dtsdy0*uyp[ip]*gaminv;
+ Real const y_new = (yp[ip] - ymin)*dyi;
+ Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv;
#endif
- const amrex::Real z_new = (zp[ip] - zmin)*dzi;
- const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv;
+ Real const z_new = (zp[ip] - zmin)*dzi;
+ Real const z_old = z_new - dtsdz0*uzp[ip]*gaminv;
#if (defined WARPX_DIM_RZ)
- const amrex::Real vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
+ Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
#elif (defined WARPX_DIM_XZ)
- const amrex::Real vy = uyp[ip]*gaminv;
+ Real const vy = uyp[ip]*gaminv;
#endif
// Shape factor arrays
// Note that there are extra values above and below
// to possibly hold the factor for the old particle
// which can be at a different grid location.
- amrex::Real sx_new[depos_order + 3] = {0.};
- amrex::Real sx_old[depos_order + 3] = {0.};
+ Real sx_new[depos_order + 3] = {0.};
+ Real sx_old[depos_order + 3] = {0.};
#if (defined WARPX_DIM_3D)
- amrex::Real sy_new[depos_order + 3] = {0.};
- amrex::Real sy_old[depos_order + 3] = {0.};
+ Real sy_new[depos_order + 3] = {0.};
+ Real sy_old[depos_order + 3] = {0.};
#endif
- amrex::Real sz_new[depos_order + 3] = {0.};
- amrex::Real sz_old[depos_order + 3] = {0.};
+ Real sz_new[depos_order + 3] = {0.};
+ Real sz_old[depos_order + 3] = {0.};
// --- Compute shape factors
// Compute shape factors for position as they are now and at old positions
@@ -397,7 +402,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 comes from the normalization of the modes
- const Complex djr_cmplx = amrex::Real(2.)*sdxi*xy_mid;
+ const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
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());
amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
xy_mid = xy_mid*xy_mid0;
@@ -407,8 +412,8 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
}
for (int k=dkl; k<=depos_order+2-dku; k++) {
for (int i=dil; i<=depos_order+2-diu; i++) {
- const amrex::Real sdyj = wq*vy*invvol*((sz_new[k] + 0.5*(sz_old[k] - sz_new[k]))*sx_new[i] +
- (0.5*sz_new[k] + 1./3.*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
+ Real const sdyj = wq*vy*invvol*((sz_new[k] + 0.5_rt * (sz_old[k] - sz_new[k]))*sx_new[i] +
+ (0.5_rt * sz_new[k] + 1._rt / 3._rt *(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
#if (defined WARPX_DIM_RZ)
Complex xy_new = xy_new0;
@@ -418,7 +423,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 comes from the normalization of the modes
// The minus sign comes from the different convention with respect to Davidson et al.
- const Complex djt_cmplx = -amrex::Real(2.)*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode*
+ const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode*
(sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old));
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());
amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
@@ -430,15 +435,15 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
}
}
for (int i=dil; i<=depos_order+2-diu; i++) {
- amrex::Real sdzk = 0.;
+ Real sdzk = 0.;
for (int k=dkl; k<=depos_order+1-dku; k++) {
- sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5*(sx_old[i] - sx_new[i]));
+ sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5_rt * (sx_old[i] - sx_new[i]));
amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
#if (defined WARPX_DIM_RZ)
Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 comes from the normalization of the modes
- const Complex djz_cmplx = amrex::Real(2.)*sdzk*xy_mid;
+ const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
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());
amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
xy_mid = xy_mid*xy_mid0;