diff options
Diffstat (limited to 'Source/Particles/WarpXParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 109 |
1 files changed, 75 insertions, 34 deletions
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index cd0c37bd8..99596f8c3 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -366,9 +366,15 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, BL_PROFILE_VAR_START(blp_pxr_cd); Real dxi = 1.0/dx[0]; Real dzi = 1.0/dx[2]; - Real invvol = dxi*dzi; Real dts2dx = 0.5*dt*dxi; Real dts2dz = 0.5*dt*dzi; +#if (AMREX_SPACEDIM == 2) + Real invvol = dxi*dzi; +#else // (AMREX_SPACEDIM == 3) + Real dyi = 1.0/dx[1]; + Real dts2dy = 0.5*dt*dyi; + Real invvol = dxi*dyi*dzi; +#endif Real clightsq = 1.0/PhysConst::c/PhysConst::c; Real q = this->charge; @@ -392,59 +398,94 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; +#if (AMREX_SPACEDIM == 3) + Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const Real ymin = xyzmin[1]; +#endif + ParallelFor( np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { + + // macro-particle current in all 3D (even for 2D runs) + // wqx, wqy and wqz are the particle current in each + // direction Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + uyp[ip]*uyp[ip]*clightsq + uzp[ip]*uzp[ip]*clightsq); - + Real wq = q*wp[ip]; Real vx = uxp[ip]*gaminv; + Real wqx = wq*invvol*vx; Real vy = uyp[ip]*gaminv; + Real wqy = wq*invvol*vy; Real vz = uzp[ip]*gaminv; + Real wqz = wq*invvol*vz; + + // --- x direction + // particle position in cell unit (1/2 timestep back) + // Real x = (Real) (xp[ip]-xmin)*dxi; + // particle position in cell unit + // Real xmid=x-dts2dx*vx; + Real xmid= (Real) (xp[ip]-xmin)*dxi-dts2dx*vx; + // x index of cell containing particle (cell-centered) + int j = (int) xmid; + // x index of cell containing particle (node-centered) + int j0= (int) (xmid-stagger_shift); + // for the cell-centered quantities, compute current on + // each point + Real xint = xmid-j; + Real sx[2] = {1.0 - xint, xint}; + // for the node-centered quantities, compute current on + // each point + xint = xmid-stagger_shift-j0; + Real sx0[2] = {1.0 - xint, xint}; - Real wq=q*wp[ip]; - Real wqx=wq*invvol*vx; - Real wqy=wq*invvol*vy; - Real wqz=wq*invvol*vz; - - Real x = (Real) (xp[ip]-xmin)*dxi; +#if (AMREX_SPACEDIM == 3) + // --- y direction + Real y = (Real) (yp[ip]-ymin)*dyi; + Real ymid=y-dts2dy*vy; + int k = (int) ymid; + int k0= (int) (ymid-stagger_shift); + Real yint = ymid-k; + Real sy[2] = {1.0 - yint, yint}; + yint = ymid-stagger_shift-k0; + Real sy0[2] = {1.0 - yint, yint}; +#endif + // --- z direction Real z = (Real) (zp[ip]-zmin)*dzi; - Real xmid=x-dts2dx*vx; Real zmid=z-dts2dz*vz; - - int j = (int) xmid; int l = (int) zmid; - int j0= (int) (xmid-stagger_shift); int l0= (int) (zmid-stagger_shift); - - Real xint = xmid-j; Real zint = zmid-l; - - Real sx[2] = {1.0 - xint, xint}; Real sz[2] = {1.0 - zint, zint}; - - xint = xmid-stagger_shift-j0; zint = zmid-stagger_shift-l0; - - Real sx0[2] = {1.0 - xint, xint}; Real sz0[2] = {1.0 - zint, zint}; - jx_arr(lo.x+j0 ,lo.y+l ,0) += sx0[0]*sz[0]*wqx; - jx_arr(lo.x+j0+1,lo.y+l ,0) += sx0[1]*sz[0]*wqx; - jx_arr(lo.x+j0 ,lo.y+l+1,0) += sx0[0]*sz[1]*wqx; - jx_arr(lo.x+j0+1,lo.y+l+1,0) += sx0[1]*sz[1]*wqx; - - jy_arr(lo.x+j ,lo.y+l ,0) += sx[0]*sz[0]*wqy; - jy_arr(lo.x+j+1,lo.y+l ,0) += sx[1]*sz[0]*wqy; - jy_arr(lo.x+j ,lo.y+l+1,0) += sx[0]*sz[1]*wqy; - jy_arr(lo.x+j+1,lo.y+l+1,0) += sx[1]*sz[1]*wqy; - - jz_arr(lo.x+j ,lo.y+l0 ,0) += sx[0]*sz0[0]*wqz; - jz_arr(lo.x+j+1,lo.y+l0 ,0) += sx[1]*sz0[0]*wqz; - jz_arr(lo.x+j ,lo.y+l0+1,0) += sx[0]*sz0[1]*wqz; - jz_arr(lo.x+j+1,lo.y+l0+1,0) += sx[1]*sz0[1]*wqz; + // Real xbounds[2], ybounds[2], zbounds[2]; + // compute_shape_factor<WarpX::nox>(xbounds, sx, xint); + // deposit_current() + +#if (AMREX_SPACEDIM == 2) + for (int iz=0; iz<2; iz++){ + for (int ix=0; ix<2; ix++){ + jx_arr(lo.x+j0+ix ,lo.y+l +iz ,0) += sx0[ix]*sz [iz]*wqx; + jy_arr(lo.x+j +ix ,lo.y+l +iz ,0) += sx [ix]*sz [iz]*wqy; + jz_arr(lo.x+j +ix ,lo.y+l0+iz ,0) += sx [ix]*sz0[iz]*wqz; + } + } +#else // (AMREX_SPACEDIM == 3) + for (int iz=0; iz<2; iz++){ + for (int iy=0; iy<2; iy++){ + for (int ix=0; ix<2; ix++){ + jx_arr(lo.x+j0+ix ,lo.y+k +iy ,lo.z+l +iz ) += sx0[ix]*sy [iy]*sz [iz]*wqx; + jy_arr(lo.x+j +ix ,lo.y+k0+iy ,lo.z+l +iz ) += sx [ix]*sy0[iy]*sz [iz]*wqy; + jz_arr(lo.x+j +ix ,lo.y+k +iy ,lo.z+l0+iz ) += sx [ix]*sy [iy]*sz0[iz]*wqz; + } + } + } +#endif } ); + #ifdef WARPX_RZ // Rescale current in r-z mode warpx_current_deposition_rz_volume_scaling( |