diff options
Diffstat (limited to 'Source/Particles')
| -rw-r--r-- | Source/Particles/Deposition/CurrentDeposition.H | 83 | ||||
| -rw-r--r-- | Source/Particles/Gather/FieldGather.H | 27 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 129 | ||||
| -rw-r--r-- | Source/Particles/Pusher/GetAndSetPosition.H | 6 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdatePosition.H | 2 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 28 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 2 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 63 |
8 files changed, 136 insertions, 204 deletions
diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 9b9853902..4a392b57e 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -42,7 +42,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real dts2dz = 0.5*dt*dzi; #if (AMREX_SPACEDIM == 2) const amrex::Real invvol = dxi*dzi; -#else // (AMREX_SPACEDIM == 3) +#elif (defined WARPX_DIM_3D) const amrex::Real dyi = 1.0/dx[1]; const amrex::Real dts2dy = 0.5*dt*dyi; const amrex::Real invvol = dxi*dyi*dzi; @@ -66,14 +66,37 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real vy = uyp[ip]*gaminv; const amrex::Real vz = uzp[ip]*gaminv; // wqx, wqy wqz are particle current in each direction +#if (defined WARPX_DIM_RZ) + // In RZ, wqx is actually wqr, and wqy is wqtheta + // Convert to cylinderical at the mid point + const amrex::Real xpmid = xp[ip] - 0.5*dt*vx; + const amrex::Real ypmid = yp[ip] - 0.5*dt*vy; + const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid); + amrex::Real costheta; + amrex::Real sintheta; + if (rpmid > 0.) { + costheta = xpmid/rpmid; + sintheta = ypmid/rpmid; + } else { + costheta = 1.; + sintheta = 0.; + } + const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta); + const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta); +#else const amrex::Real wqx = wq*invvol*vx; const amrex::Real wqy = wq*invvol*vy; +#endif const amrex::Real wqz = wq*invvol*vz; // --- Compute shape factors // x direction // Get particle position after 1/2 push back in position +#if (defined WARPX_DIM_RZ) + const amrex::Real xmid = (rpmid-xmin)*dxi; +#else const amrex::Real xmid = (xp[ip]-xmin)*dxi-dts2dx*vx; +#endif // Compute shape factors for node-centered quantities amrex::Real AMREX_RESTRICT sx [depos_order + 1]; // j: leftmost grid point (node-centered) that the particle touches @@ -83,7 +106,7 @@ void doDepositionShapeN(const amrex::Real * const xp, // j0: leftmost grid point (cell-centered) that the particle touches const int j0 = compute_shape_factor<depos_order>(sx0, xmid-stagger_shift); -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) // y direction const amrex::Real ymid= (yp[ip]-ymin)*dyi-dts2dy*vy; amrex::Real AMREX_RESTRICT sy [depos_order + 1]; @@ -99,7 +122,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const int l0 = compute_shape_factor<depos_order>(sz0, zmid-stagger_shift); // Deposit current into jx_arr, jy_arr and jz_arr -#if (AMREX_SPACEDIM == 2) +#if (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( @@ -113,7 +136,7 @@ void doDepositionShapeN(const amrex::Real * const xp, sx [ix]*sz0[iz]*wqz); } } -#else // (AMREX_SPACEDIM == 3) +#elif (defined WARPX_DIM_3D) for (int iz=0; iz<=depos_order; iz++){ for (int iy=0; iy<=depos_order; iy++){ for (int ix=0; ix<=depos_order; ix++){ @@ -169,7 +192,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dtsdx0 = dt*dxi; const amrex::Real xmin = xyzmin[0]; -#if (AMREX_SPACEDIM == 3) +#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]; @@ -178,11 +201,11 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real dtsdz0 = dt*dzi; const amrex::Real zmin = xyzmin[2]; -#if (AMREX_SPACEDIM == 3) +#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]); -#elif (AMREX_SPACEDIM == 2) +#elif (defined WARPX_DIM_2D) || (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]); @@ -197,34 +220,57 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, // --- Get particle quantities const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq - + uyp[ip]*uyp[ip]*clightsq - + uzp[ip]*uzp[ip]*clightsq); + + uyp[ip]*uyp[ip]*clightsq + + uzp[ip]*uzp[ip]*clightsq); // wqx, wqy wqz are particle current in each direction const amrex::Real wq = q*wp[ip]; const amrex::Real wqx = wq*invdtdx; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) const amrex::Real wqy = wq*invdtdy; #endif const amrex::Real wqz = wq*invdtdz; // computes current and old position in grid units +#if (defined WARPX_DIM_RZ) + const amrex::Real r_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real r_old = std::sqrt((xp[ip] - dt*uxp[ip]*gaminv)*(xp[ip] - dt*uxp[ip]*gaminv) + + (yp[ip] - dt*uyp[ip]*gaminv)*(yp[ip] - dt*uyp[ip]*gaminv)); + const amrex::Real x_new = (r_new - xmin)*dxi; + const amrex::Real x_old = (r_old - xmin)*dxi; +#else const amrex::Real x_new = (xp[ip] - xmin)*dxi; const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv; -#if (AMREX_SPACEDIM == 3) +#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; #endif const amrex::Real z_new = (zp[ip] - zmin)*dzi; const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv; +#if (defined WARPX_DIM_RZ) + amrex::Real costheta; + amrex::Real sintheta; + if (r_new > 0.) { + costheta = xp[ip]/r_new; + sintheta = yp[ip]/r_new; + } else { + costheta = 1.; + sintheta = 0.; + } + const amrex::Real vy = (-uxp[ip]*sintheta + uyp[ip]*costheta)*gaminv; +#elif (defined WARPX_DIM_2D) + const amrex::Real 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 AMREX_RESTRICT sx_new[depos_order + 3] = {0.}; amrex::Real AMREX_RESTRICT sx_old[depos_order + 3] = {0.}; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) amrex::Real AMREX_RESTRICT sy_new[depos_order + 3] = {0.}; amrex::Real AMREX_RESTRICT sy_old[depos_order + 3] = {0.}; #endif @@ -236,7 +282,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, // [ijk]_new: leftmost grid point that the particle touches const int i_new = compute_shape_factor<depos_order>(sx_new+1, x_new); const int i_old = compute_shifted_shape_factor<depos_order>(sx_old, x_old, i_new); -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) const int j_new = compute_shape_factor<depos_order>(sy_new+1, y_new); const int j_old = compute_shifted_shape_factor<depos_order>(sy_old, y_old, j_new); #endif @@ -247,7 +293,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, int dil = 1, diu = 1; if (i_old < i_new) dil = 0; if (i_old > i_new) diu = 0; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) int djl = 1, dju = 1; if (j_old < j_new) djl = 0; if (j_old > j_new) dju = 0; @@ -256,7 +302,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, if (k_old < k_new) dkl = 0; if (k_old > k_new) dku = 0; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) for (int k=dkl; k<=depos_order+2-dku; k++) { for (int j=djl; j<=depos_order+2-dju; j++) { @@ -289,7 +335,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, } } -#elif (AMREX_SPACEDIM == 2) +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) for (int k=dkl; k<=depos_order+2-dku; k++) { amrex::Real sdxi = 0.; @@ -300,8 +346,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * 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*uyp[ip]*gaminv*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])); + 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])); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdyj); } } @@ -313,6 +359,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, } } + #endif } ); diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index be96dd393..8f5e8d4cf 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -42,7 +42,9 @@ void doGatherShapeN(const amrex::Real * const xp, #endif const amrex::Real xmin = xyzmin[0]; +#if (AMREX_SPACEDIM == 3) const amrex::Real ymin = xyzmin[1]; +#endif const amrex::Real zmin = xyzmin[2]; // Loop over particles and gather fields from @@ -53,7 +55,12 @@ void doGatherShapeN(const amrex::Real * const xp, // --- Compute shape factors // x direction // Get particle position +#ifdef WARPX_DIM_RZ + const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real x = (r - xmin)*dxi; +#else const amrex::Real x = (xp[ip]-xmin)*dxi; +#endif // Compute shape factors for node-centered quantities amrex::Real AMREX_RESTRICT sx [depos_order + 1]; // j: leftmost grid point (node-centered) that particle touches @@ -126,6 +133,26 @@ void doGatherShapeN(const amrex::Real * const xp, by_arr(lo.x+j0+ix, lo.y+l0+iz, 0); } } + +#ifdef WARPX_DIM_RZ + // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey + amrex::Real costheta; + amrex::Real sintheta; + if (r > 0.) { + costheta = xp[ip]/r; + sintheta = yp[ip]/r; + } else { + costheta = 1.; + sintheta = 0.; + } + const amrex::Real Exp_save = Exp[ip]; + Exp[ip] = costheta*Exp[ip] - sintheta*Eyp[ip]; + Eyp[ip] = costheta*Eyp[ip] + sintheta*Exp_save; + const amrex::Real Bxp_save = Bxp[ip]; + Bxp[ip] = costheta*Bxp[ip] - sintheta*Byp[ip]; + Byp[ip] = costheta*Byp[ip] + sintheta*Bxp_save; +#endif + #else // (AMREX_SPACEDIM == 3) // Gather field on particle Exp[i] from field on grid ex_arr for (int iz=0; iz<=depos_order; iz++){ diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7d63bd8e7..59d3ce76f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -148,7 +148,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, npart /= 4; } for (long i = 0; i < npart; ++i) { -#if ( AMREX_SPACEDIM == 3 | WARPX_RZ) +#if ( AMREX_SPACEDIM == 3 | WARPX_DIM_RZ) Real weight = q_tot/npart/charge; Real x = distx(mt); Real y = disty(mt); @@ -269,7 +269,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (!part_realbox.ok()) part_realbox = geom.ProbDomain(); int num_ppc = plasma_injector->num_particles_per_cell; -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ Real rmax = std::min(plasma_injector->xmax, part_realbox.hi(0)); #endif @@ -323,7 +323,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) Real density_min = plasma_injector->density_min; Real density_max = plasma_injector->density_max; -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ bool radially_weighted = plasma_injector->radially_weighted; #endif @@ -496,11 +496,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #endif // Save the x and y values to use in the insideBounds checks. - // This is needed with WARPX_RZ since x and y are modified. + // This is needed with WARPX_DIM_RZ since x and y are modified. Real xb = x; Real yb = y; -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ // Replace the x and y, choosing the angle randomly. // These x and y are used to get the momentum and density Real theta = 2.*MathConst::pi*amrex::Random(); @@ -574,7 +574,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) // Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); Real weight = dens * scale_fac; -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ if (radially_weighted) { weight *= 2.*MathConst::pi*xb; } else { @@ -603,7 +603,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = y; p.pos(2) = z; #elif (AMREX_SPACEDIM == 2) -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ pa[PIdx::theta][ip] = theta; #endif p.pos(0) = xb; @@ -893,40 +893,13 @@ PhysicalParticleContainer::FieldGather (int lev, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - const std::array<Real,3>& xyzmin = WarpX::LowerCorner(box, lev); - const int* ixyzmin = box.loVect(); - // // Field Gather // -#ifdef WARPX_RZ - const int ll4symtry = false; - long lvect_fieldgathe = 64; - warpx_geteb_energy_conserving( - &np, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), - Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), - ixyzmin, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); -#endif if (cost) { const Box& tbx = pti.tilebox(); @@ -1195,42 +1168,14 @@ PhysicalParticleContainer::Evolve (int lev, // Field Gather of Aux Data (i.e., the full solution) // BL_PROFILE_VAR_START(blp_pxr_fg); -#ifdef WARPX_RZ - const int ll4symtry = false; - long lvect_fieldgathe = 64; - const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); - const int* ixyzmin_grid = box.loVect(); - warpx_geteb_energy_conserving( - &np_gather, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), - Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), - ixyzmin_grid, - &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(*exfab), - BL_TO_FORTRAN_ANYD(*eyfab), - BL_TO_FORTRAN_ANYD(*ezfab), - BL_TO_FORTRAN_ANYD(*bxfab), - BL_TO_FORTRAN_ANYD(*byfab), - BL_TO_FORTRAN_ANYD(*bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); -#endif if (np_gather < np) { const IntVect& ref_ratio = WarpX::RefRatio(lev-1); const Box& cbox = amrex::coarsen(box,ref_ratio); - const std::array<Real,3>& cxyzmin_grid = WarpX::LowerCorner(cbox, lev-1); - const int* cixyzmin_grid = cbox.loVect(); // Data on the grid FArrayBox const* cexfab = &(*cEx)[pti]; @@ -1293,29 +1238,6 @@ PhysicalParticleContainer::Evolve (int lev, } // Field gather for particles in gather buffers -#ifdef WARPX_RZ - - long ncrse = np - nfine_gather; - warpx_geteb_energy_conserving( - &ncrse, - m_xp[thread_num].dataPtr()+nfine_gather, - m_yp[thread_num].dataPtr()+nfine_gather, - m_zp[thread_num].dataPtr()+nfine_gather, - Exp.dataPtr()+nfine_gather, Eyp.dataPtr()+nfine_gather, Ezp.dataPtr()+nfine_gather, - Bxp.dataPtr()+nfine_gather, Byp.dataPtr()+nfine_gather, Bzp.dataPtr()+nfine_gather, - cixyzmin_grid, - &cxyzmin_grid[0], &cxyzmin_grid[1], &cxyzmin_grid[2], - &cdx[0], &cdx[1], &cdx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(*cexfab), - BL_TO_FORTRAN_ANYD(*ceyfab), - BL_TO_FORTRAN_ANYD(*cezfab), - BL_TO_FORTRAN_ANYD(*cbxfab), - BL_TO_FORTRAN_ANYD(*cbyfab), - BL_TO_FORTRAN_ANYD(*cbzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, @@ -1323,7 +1245,6 @@ PhysicalParticleContainer::Evolve (int lev, cEx->nGrow(), e_is_nodal, nfine_gather, np-nfine_gather, thread_num, lev, lev-1); -#endif } BL_PROFILE_VAR_STOP(blp_pxr_fg); @@ -1657,38 +1578,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); -#ifdef WARPX_RZ - const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); - const int* ixyzmin_grid = box.loVect(); - - const int ll4symtry = false; - long lvect_fieldgathe = 64; - - warpx_geteb_energy_conserving( - &np, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), - Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), - ixyzmin_grid, - &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else - int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, - Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); -#endif + int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 42c61343e..3c74baeb2 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -5,7 +5,7 @@ #include <WarpXParticleContainer.H> #include <AMReX_REAL.H> -#ifndef WARPX_RZ +#ifndef WARPX_DIM_RZ /* \brief Extract the particle's coordinates from the ParticleType struct `p`, * and stores them in the variables `x`, `y`, `z`. */ @@ -42,7 +42,7 @@ void SetPosition( #endif } -# else // if WARPX_RZ is True +# elif defined WARPX_DIM_RZ /* \brief Extract the particle's coordinates from `theta` and the attributes * of the ParticleType struct `p` (which contains the radius), @@ -71,6 +71,6 @@ void SetCylindricalPositionFromCartesian( p.pos(1) = z; } -#endif // WARPX_RZ +#endif // WARPX_DIM_RZ #endif // WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 0a4f579f4..a9df63a30 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -20,7 +20,7 @@ void UpdatePosition( const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2); // Update positions over one time step x += ux * inv_gamma * dt; -#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D +#if (AMREX_SPACEDIM == 3) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * inv_gamma * dt; #endif z += uz * inv_gamma * dt; diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 83556e69f..f049fdb7c 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -421,38 +421,10 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); -#ifdef WARPX_RZ - const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); - const int* ixyzmin_grid = box.loVect(); - - const int ll4symtry = false; - long lvect_fieldgathe = 64; - - warpx_geteb_energy_conserving( - &np, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), - Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), - ixyzmin_grid, - &xyzmin_grid[0], &xyzmin_grid[1], &xyzmin_grid[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); -#else int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); -#endif // Save the position and momenta, making copies auto uxp_save = uxp; diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7cf53260a..a609b4cb3 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -13,7 +13,7 @@ struct PIdx enum { // Particle Attributes stored in amrex::ParticleContainer's struct of array w = 0, // weight ux, uy, uz, Ex, Ey, Ez, Bx, By, Bz, -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ theta, // RZ needs all three position components #endif nattribs diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 89f233b2c..8e4adc528 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -27,7 +27,7 @@ void WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDeviceVector<Real>& y, Cuda::ManagedDeviceVector<Real>& z) const { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ const auto& attribs = GetAttribs(); const auto& theta = attribs[PIdx::theta]; y.resize(x.size()); @@ -44,10 +44,10 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDevi void WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector<Real>& x, const Cuda::ManagedDeviceVector<Real>& y, const Cuda::ManagedDeviceVector<Real>& z) { -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ auto& attribs = GetAttribs(); auto& theta = attribs[PIdx::theta]; - Cuda::DeviceVector<Real> r(x.size()); + Cuda::ManagedDeviceVector<Real> r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { theta[i] = std::atan2(y[i], x[i]); r[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); @@ -80,7 +80,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["Bx"] = PIdx::Bx; particle_comps["By"] = PIdx::By; particle_comps["Bz"] = PIdx::Bz; -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ particle_comps["theta"] = PIdx::theta; #endif @@ -163,7 +163,7 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, p.pos(1) = y; p.pos(2) = z; #elif (AMREX_SPACEDIM == 2) -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ attribs[PIdx::theta] = std::atan2(y, x); x = std::sqrt(x*x + y*y); #endif @@ -209,7 +209,7 @@ WarpXParticleContainer::AddNParticles (int lev, std::size_t np = iend-ibegin; -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ Vector<Real> theta(np); #endif @@ -228,7 +228,7 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = y[i]; p.pos(2) = z[i]; #elif (AMREX_SPACEDIM == 2) -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ theta[i-ibegin] = std::atan2(y[i], x[i]); p.pos(0) = std::sqrt(x[i]*x[i] + y[i]*y[i]); #else @@ -265,7 +265,7 @@ WarpXParticleContainer::AddNParticles (int lev, for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ if (comp == PIdx::theta) { particle_tile.push_back_real(comp, theta.front(), theta.back()); } @@ -394,14 +394,6 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, &WarpX::nox,&WarpX::noy,&WarpX::noz, &j_is_nodal, &lvect,&WarpX::current_deposition_algo); -#ifdef WARPX_RZ - // Rescale current in r-z mode - warpx_current_deposition_rz_volume_scaling( - jx_ptr, &ngJ, jxntot.getVect(), - jy_ptr, &ngJ, jyntot.getVect(), - jz_ptr, &ngJ, jzntot.getVect(), - &xyzmin[0], &dx[0]); -#endif BL_PROFILE_VAR_STOP(blp_pxr_cd); #ifndef AMREX_USE_GPU @@ -503,7 +495,8 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain - const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev);; + // Note that this includes guard cells since it is after tilebox.ngrow + const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); // xyzmin is built on pti.tilebox(), so it does // not include staggering, so the stagger_shift has to be done by hand. // Alternatively, we could define xyzminx from tbx (and the same for 3 @@ -513,35 +506,35 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ - doEsirkepovDepositionShapeN<1>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + doEsirkepovDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doEsirkepovDepositionShapeN<2>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + doEsirkepovDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doEsirkepovDepositionShapeN<3>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + doEsirkepovDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } } else { if (WarpX::nox == 1){ - doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, stagger_shift, q); } else if (WarpX::nox == 2){ - doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, stagger_shift, q); } else if (WarpX::nox == 3){ - doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, stagger_shift, q); } @@ -617,7 +610,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ warpx_charge_deposition_rz_volume_scaling( data_ptr, &ngRho, rholen.getVect(), &xyzmin[0], &dx[0]); @@ -677,7 +670,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ warpx_charge_deposition_rz_volume_scaling( data_ptr, &ngRho, rholen.getVect(), &cxyzmin_tile[0], &cdx[0]); @@ -837,7 +830,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ long ngRho = WarpX::nox; warpx_charge_deposition_rz_volume_scaling( data_ptr, &ngRho, rholen.getVect(), @@ -1022,7 +1015,7 @@ WarpXParticleContainer::PushX (int lev, Real dt) Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ Real* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); #endif // Loop over the particles and update their position @@ -1030,12 +1023,12 @@ WarpXParticleContainer::PushX (int lev, Real dt) [=] AMREX_GPU_DEVICE (long i) { ParticleType& p = pstructs[i]; // Particle object that gets updated Real x, y, z; // Temporary variables -#ifndef WARPX_RZ +#ifndef WARPX_DIM_RZ GetPosition( x, y, z, p ); // Initialize x, y, z UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); SetPosition( p, x, y, z ); // Update the object p #else - // For WARPX_RZ, the particles are still pushed in 3D Cartesian + // For WARPX_DIM_RZ, the particles are still pushed in 3D Cartesian GetCartesianPositionFromCylindrical( x, y, z, p, theta[i] ); UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); SetCylindricalPositionFromCartesian( p, theta[i], x, y, z ); |