diff options
Diffstat (limited to 'Source')
22 files changed, 978 insertions, 530 deletions
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index bb1300562..c0ec3cd90 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -12,6 +12,8 @@ #include <limits> #include <WarpX.H> +#include <WarpX_QED_K.H> +#include <WarpX_QED_Field_Pushers.cpp> #include <WarpXConst.H> #include <WarpX_f.H> #include <WarpXUtil.H> @@ -358,10 +360,24 @@ WarpX::OneStep_nosub (Real cur_time) // Push E and B from {n} to {n+1} // (And update guard cells immediately afterwards) #ifdef WARPX_USE_PSATD + if (use_hybrid_QED) + { + WarpX::Hybrid_QED_Push(dt); + FillBoundaryE(guard_cells.ng_alloc_EB, guard_cells.ng_Extra); + } PushPSATD(dt[0]); - if (do_pml) DampPML(); FillBoundaryE(guard_cells.ng_alloc_EB, guard_cells.ng_Extra); FillBoundaryB(guard_cells.ng_alloc_EB, guard_cells.ng_Extra); + + if (use_hybrid_QED) + { + WarpX::Hybrid_QED_Push(dt); + FillBoundaryE(guard_cells.ng_alloc_EB, guard_cells.ng_Extra); + + } + if (do_pml) DampPML(); + +; #else EvolveF(0.5*dt[0], DtType::FirstHalf); FillBoundaryF(guard_cells.ng_FieldSolverF); diff --git a/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp b/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp new file mode 100644 index 000000000..afc205aa2 --- /dev/null +++ b/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp @@ -0,0 +1,175 @@ +#include <cmath> +#include <limits> + +#include <WarpX.H> +#include <WarpXConst.H> +#include <WarpX_f.H> +#include <WarpX_K.H> +#include <WarpX_PML_kernels.H> +#include <WarpX_FDTD.H> +#ifdef WARPX_USE_PY +#include <WarpX_py.H> +#endif + +#include <WarpX_QED_K.H> + +#include <PML_current.H> + +#ifdef BL_USE_SENSEI_INSITU +#include <AMReX_AmrMeshInSituBridge.H> +#endif + +using namespace amrex; + + +void +WarpX::Hybrid_QED_Push (amrex::Vector<amrex::Real> dt) +{ + if (WarpX::do_nodal == 0) { + Print()<<"The do_nodal flag is tripped.\n"; + try{ + throw "Error: The Hybrid QED method is currently only compatible with the nodal scheme.\n"; + } + catch (const char* msg) { + std::cerr << msg << std::endl; + exit(0); + } + } + for (int lev = 0; lev <= finest_level; ++lev) { + Hybrid_QED_Push(lev, dt[lev]); + } +} + +void +WarpX::Hybrid_QED_Push (int lev, Real a_dt) +{ + BL_PROFILE("WarpX::Hybrid_QED_Push()"); + Hybrid_QED_Push(lev, PatchType::fine, a_dt); + if (lev > 0) + { + Hybrid_QED_Push(lev, PatchType::coarse, a_dt); + } +} + +void +WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt) +{ + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array<Real,3>& dx_vec= WarpX::CellSize(patch_level); + const Real dx = dx_vec[0]; + const Real dy = dx_vec[1]; + const Real dz = dx_vec[2]; + + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; + if (patch_type == PatchType::fine) + { + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + } + + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + + // xmin is only used by the kernel for cylindrical geometry, + // in which case it is actually rmin. + const Real xmin = Geom(0).ProbLo(0); + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Real wt = amrex::second(); + + // Get boxes for E and B + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + + // Get field arrays + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); + + // Define grown box with 1 ghost cell for finite difference stencil + const Box& gex = amrex::grow(tex,1); + const Box& gey = amrex::grow(tey,1); + const Box& gez = amrex::grow(tez,1); + + // Temporary arrays for electric field, protected by Elixir on GPU + FArrayBox tmpEx_fab(gex,1); + Elixir tmpEx_eli = tmpEx_fab.elixir(); + auto const& tmpEx = tmpEx_fab.array(); + + FArrayBox tmpEy_fab(gey,1); + Elixir tmpEy_eli = tmpEy_fab.elixir(); + auto const& tmpEy = tmpEy_fab.array(); + + FArrayBox tmpEz_fab(gez,1); + Elixir tmpEz_eli = tmpEz_fab.elixir(); + auto const& tmpEz = tmpEz_fab.array(); + + // Copy electric field to temporary arrays + AMREX_PARALLEL_FOR_4D( + gex, 1, i, j, k, n, + { tmpEx(i,j,k,n) = Exfab(i,j,k,n); } + ); + + AMREX_PARALLEL_FOR_4D( + gey, 1, i, j, k, n, + { tmpEy(i,j,k,n) = Eyfab(i,j,k,n); } + ); + + AMREX_PARALLEL_FOR_4D( + gez, 1, i, j, k, n, + { tmpEz(i,j,k,n) = Ezfab(i,j,k,n); } + ); + + // Apply QED correction to electric field, using temporary arrays. + amrex::ParallelFor( + tbx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_hybrid_QED_push(j,k,l, Exfab, Eyfab, Ezfab, Bxfab, Byfab, + Bzfab, tmpEx, tmpEy, tmpEz, dx, dy, dz, + a_dt); + } + ); + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (amrex::second() - wt) / cbx.d_numPts(); + auto costfab = cost->array(mfi); + + amrex::ParallelFor( + cbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + costfab(i,j,k) += wt; + } + ); + } + } +} diff --git a/Source/FieldSolver/WarpX_QED_K.H b/Source/FieldSolver/WarpX_QED_K.H new file mode 100644 index 000000000..7c55ce4df --- /dev/null +++ b/Source/FieldSolver/WarpX_QED_K.H @@ -0,0 +1,322 @@ +#ifndef WarpX_QED_K_h +#define WarpX_QED_K_h + +#include <AMReX_FArrayBox.H> +#include <WarpXConst.H> +#include <cmath> + +using namespace amrex; + +/** + * calc_M calculates the Magnetization field of the vacuum at a specific point and returns it as a three component vector + * \param[in] arr This is teh empty array that will be filled with the components of the M-field + * \param[in] ex The x-component of the E-field at the point at whicht the M-field is to be calculated + * \param[in] ey The y-component of the E-field at the point at whicht the M-field is to be calculated + * \param[in] ez The z-component of the E-field at the point at whicht the M-field is to be calculated + * \param[in] bx The x-component of the B-field at the point at whicht the M-field is to be calculated + * \param[in] by The y-component of the B-field at the point at whicht the M-field is to be calculated + * \param[in] bz The z-component of the B-field at the point at whicht the M-field is to be calculated + * \param[in] xi The quantum parameter being used for the simulation + * \param[in] c2 The speed of light squared + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void calc_M(Real arr [], Real ex, Real ey, Real ez, Real bx, Real by, Real bz, Real xi, Real c2) +{ + const Real ee = ex*ex+ey*ey+ez*ez; + const Real bb = bx*bx+by*by+bz*bz; + const Real eb = ex*bx+ey*by+ez*bz; + arr[0] = -2*xi*c2*( 2*bx*(ee-c2*bb) - 7*ex*eb ); + arr[1] = -2*xi*c2*( 2*by*(ee-c2*bb) - 7*ey*eb ); + arr[2] = -2*xi*c2*( 2*bz*(ee-c2*bb) - 7*ez*eb ); +}; + + +/** + * warpx_hybrid_QED_push uses an FDTD scheme to calculate QED corrections to + * Maxwell's equations and preforms a half timestep correction to the E-fields + * + * \param[in] Ex This function modifies the Ex field at the end + * \param[in] Ey This function modifies the Ey field at the end + * \param[in] Ez This function modifies the Ez field at the end + * \param[in] Bx The QED corrections are non-linear functions of B. However, + * they do not modify B itself + * \param[in] By The QED corrections are non-linear functions of B. However, + * they do not modify B itself + * \param[in] Bz The QED corrections are non-linear functions of B. However, + * they do not modify B itself + * \param[in] tempEx Since the corrections to E at a given node are non-linear functions + * of the values of E on the surronding nodes, temp arrays are used so that + * modifications to one node do not influence the corrections to the surronding nodes + * \param[in] tempEy Since the corrections to E at a given node are non-linear functions + * of the values of E on the surronding nodes, temp arrays are used so that modifications to + * one node do not influence the corrections to the surronding nodes + * \param[in] tempEz Since the corrections to E at a given node are non-linear functions + * of the values of E on the surronding nodes, temp arrays are used so that modifications to + * one node do not influence the corrections to the surronding nodes + * \param[in] dx The x spatial step, used for calculating curls + * \param[in] dy The y spatial step, used for calculating curls + * \param[in] dz The z spatial step, used for calulating curls + * \param[in] dt The temporal step, used for the half push/correction to the E-fields at the end of the function + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_hybrid_QED_push (int j, int k, int l, Array4<Real> const& Ex, Array4<Real> + const& Ey, Array4<Real> const& Ez, Array4<Real> const& Bx, + Array4<Real> const& By, Array4<Real const> const& Bz, + Array4<Real> const& tmpEx, Array4<Real> const& tmpEy, + Array4<Real> const& tmpEz, Real dx, Real dy, Real dz, Real dt) +{ + + +// Defining constants to be used in the calculations + +constexpr amrex::Real c2 = PhysConst::c * PhysConst::c; +const amrex::Real xi = WarpX::quantum_xi; +const amrex::Real dxi = 1./dx; +const amrex::Real dzi = 1./dz; + +#if (AMREX_SPACEDIM == 3) +const amrex::Real dyi = 1./dy; + + // Picking out points for stencil to be used in curl function of M + + amrex::Real Mpx [3] = {0.,0.,0.}; + amrex::Real Mnx [3] = {0.,0.,0.}; + amrex::Real Mpy [3] = {0.,0.,0.}; + amrex::Real Mny [3] = {0.,0.,0.}; + amrex::Real Mpz [3] = {0.,0.,0.}; + amrex::Real Mnz [3] = {0.,0.,0.}; + + // Calcualting the M-field at the chosen stencil points + + calc_M(Mpx, tmpEx(j+1,k,l), tmpEy(j+1,k,l), tmpEz(j+1,k,l), + Bx(j+1,k,l), By(j+1,k,l), Bz(j+1,k,l), xi, c2); + calc_M(Mnx, tmpEx(j-1,k,l), tmpEy(j-1,k,l), tmpEz(j-1,k,l), + Bx(j-1,k,l), By(j-1,k,l), Bz(j-1,k,l), xi, c2); + calc_M(Mpy, tmpEx(j,k+1,l), tmpEy(j,k+1,l), tmpEz(j,k+1,l), + Bx(j,k+1,l), By(j,k+1,l), Bz(j,k+1,l), xi, c2); + calc_M(Mny, tmpEx(j,k-1,l), tmpEy(j,k-1,l), tmpEz(j,k-1,l), + Bx(j,k-1,l), By(j,k-1,l), Bz(j,k-1,l), xi, c2); + calc_M(Mpz, tmpEx(j,k,l+1), tmpEy(j,k,l+1), tmpEz(j,k,l+1), + Bx(j,k,l+1), By(j,k,l+1), Bz(j,k,l+1), xi, c2); + calc_M(Mnz, tmpEx(j,k,l-1), tmpEy(j,k,l-1), tmpEz(j,k,l-1), + Bx(j,k,l-1), By(j,k,l-1), Bz(j,k,l-1), xi, c2); + + // Calculating necessary curls + + const amrex::Real VxM[3] = { + 0.5*( (Mpy[2]-Mny[2])*dyi - (Mpz[1]-Mnz[1])*dzi ), + 0.5*( (Mpz[0]-Mnz[0])*dzi - (Mpx[2]-Mnx[2])*dxi ), + 0.5*( (Mpx[1]-Mnx[1])*dxi - (Mpy[0]-Mny[0])*dyi ), + }; + + const amrex::Real VxE[3] = { + 0.5*( (tmpEz(j,k+1,l)-tmpEz(j,k-1,l) )*dyi - (tmpEy(j,k,l+1)-tmpEy(j,k,l-1) )*dzi ), + 0.5*( (tmpEx(j,k,l+1)-tmpEx(j,k,l-1) )*dzi - (tmpEz(j+1,k,l)-tmpEz(j-1,k,l) )*dxi ), + 0.5*( (tmpEy(j+1,k,l)-tmpEy(j-1,k,l) )*dxi - (tmpEx(j,k+1,l)-tmpEx(j,k-1,l) )*dyi ), + }; + + const amrex::Real VxB[3] = { + 0.5*( (Bz(j,k+1,l)-Bz(j,k-1,l) )*dyi - (By(j,k,l+1)-By(j,k,l-1) )*dzi ), + 0.5*( (Bx(j,k,l+1)-Bx(j,k,l-1) )*dzi - (Bz(j+1,k,l)-Bz(j-1,k,l) )*dxi ), + 0.5*( (By(j+1,k,l)-By(j-1,k,l) )*dxi - (Bx(j,k+1,l)-Bx(j,k-1,l) )*dyi ), + }; + + // Defining comapct values for QED corrections + + const amrex::Real ex = tmpEx(j,k,l); + const amrex::Real ey = tmpEy(j,k,l); + const amrex::Real ez = tmpEz(j,k,l); + const amrex::Real bx = Bx(j,k,l); + const amrex::Real by = By(j,k,l); + const amrex::Real bz = Bz(j,k,l); + const amrex::Real ee = ex*ex + ey*ey + ez*ez; + const amrex::Real bb = bx*bx + by*by + bz*bz; + const amrex::Real eb = ex*bx + ey*by + ez*bz; + const amrex::Real EVxE = ex*VxE[0] + ey*VxE[1] + ez*VxE[2]; + const amrex::Real BVxE = bx*VxE[0] + by*VxE[1] + bz*VxE[2]; + const amrex::Real EVxB = ex*VxB[0] + ey*VxB[1] + ez*VxB[2]; + const amrex::Real BVxB = bx*VxB[0] + by*VxB[1] + bz*VxB[2]; + + const amrex::Real beta = 4*xi*( ee - c2*bb ) + PhysConst::ep0; + + const amrex::Real Alpha[3] = { + 2*xi*c2*( -7*bx*EVxE - 7*VxE[0]*eb + 4*ex*BVxE ) + VxM[0], + 2*xi*c2*( -7*by*EVxE - 7*VxE[1]*eb + 4*ey*BVxE ) + VxM[1], + 2*xi*c2*( -7*bz*EVxE - 7*VxE[2]*eb + 4*ez*BVxE ) + VxM[2] + }; + + const amrex::Real Omega[3] = { + Alpha[0] + 2*xi*c2*( 4*ex*EVxB + 2*VxB[0]*( ee - c2*bb ) + 7*c2*bx*BVxB ), + Alpha[1] + 2*xi*c2*( 4*ey*EVxB + 2*VxB[1]*( ee - c2*bb ) + 7*c2*by*BVxB ), + Alpha[2] + 2*xi*c2*( 4*ez*EVxB + 2*VxB[2]*( ee - c2*bb ) + 7*c2*bz*BVxB ) + }; + + // Calcualting matrix values for the QED correction algorithm + + const amrex::Real a00 = beta + xi*( 8*ex*ex + 14*c2*bx*bx ); + + const amrex::Real a11 = beta + xi*( 8*ey*ey + 14*c2*by*by ); + + const amrex::Real a22 = beta + xi*( 8*ez*ez + 14*c2*bz*bz ); + + const amrex::Real a01 = xi*( 2*ex*ey + 14*c2*bx*by ); + + const amrex::Real a02 = xi*( 2*ex*ez + 14*c2*bx*bz ); + + const amrex::Real a12 = xi*( 2*ez*ey + 14*c2*bz*by ); + + const amrex::Real detA = a00*( a11*a22 - a12*a12 ) - a01*( a01*a22 - a02*a12 )+a02*( a01*a12 - a02*a11 ); + + // Calcualting the rows of the inverse matrix using the general 3x3 inverse form + + const amrex::Real invAx[3] = {a22*a11 - a12*a12, a12*a02 - a22*a01, a12*a01 - a11*a02}; + + const amrex::Real invAy[3] = {a02*a12 - a22*a01, a00*a22 - a02*a02, a01*a02 - a12*a00}; + + const amrex::Real invAz[3] = {a12*a01 - a02*a11, a02*a01 - a12*a00, a11*a00 - a01*a01}; + + // Calcualting the final QED corrections by mutliplying the Omega vector with the inverse matrix + + const amrex::Real dEx = (-1/detA)*(invAx[0]*Omega[0] + + invAx[1]*Omega[1] + + invAx[2]*Omega[2]); + + const amrex::Real dEy = (-1/detA)*(invAy[0]*Omega[0] + + invAy[1]*Omega[1] + + invAy[2]*Omega[2]); + + const amrex::Real dEz = (-1/detA)*(invAz[0]*Omega[0] + + invAz[1]*Omega[1] + + invAz[2]*Omega[2]); + + // Adding the QED corrections to the origional fields + + Ex(j,k,l) = Ex(j,k,l) + 0.5*dt*dEx; + + Ey(j,k,l) = Ey(j,k,l) + 0.5*dt*dEy; + + Ez(j,k,l) = Ez(j,k,l) + 0.5*dt*dEz; + + +// 2D case - follows naturally from 3D case +#else + + // Picking out points for stencil to be used in curl function of M + + amrex::Real Mpx [3] = {0.,0.,0.}; + amrex::Real Mnx [3] = {0.,0.,0.}; + amrex::Real Mpz [3] = {0.,0.,0.}; + amrex::Real Mnz [3] = {0.,0.,0.}; + + // Calcualting the M-field at the chosen stencil points + + calc_M(Mpx, tmpEx(j+1,k,0), tmpEy(j+1,k,0), tmpEz(j+1,k,0), + Bx(j+1,k,0), By(j+1,k,0), Bz(j+1,k,0), xi, c2); + calc_M(Mnx, tmpEx(j-1,k,0), tmpEy(j-1,k,0), tmpEz(j-1,k,0), + Bx(j-1,k,0), By(j-1,k,0), Bz(j-1,k,0), xi, c2); + calc_M(Mpz, tmpEx(j,k+1,0), tmpEy(j,k+1,0), tmpEz(j,k+1,0), + Bx(j,k+1,0), By(j,k+1,0), Bz(j,k+1,0), xi, c2); + calc_M(Mnz, tmpEx(j,k-1,0), tmpEy(j,k-1,0), tmpEz(j,k-1,0), + Bx(j,k-1,0), By(j,k-1,0), Bz(j,k-1,0), xi, c2); + + // Calculating necessary curls + + const amrex::Real VxM[3] = { + -0.5*(Mpz[1]-Mnz[1])*dzi, + 0.5*( (Mpz[0]-Mnz[0])*dzi - (Mpx[2]-Mnx[2])*dxi ), + 0.5*(Mpx[1]-Mnx[1])*dxi, + }; + + const amrex::Real VxE[3] = { + -0.5*(tmpEy(j,k+1,0)-tmpEy(j,k-1,0) )*dzi, + 0.5*( (tmpEx(j,k+1,0)-tmpEx(j,k-1,0) )*dzi - (tmpEz(j+1,k,0)-tmpEz(j-1,k,0) )*dxi ), + 0.5*(tmpEy(j+1,k,0)-tmpEy(j-1,k,0) )*dxi, + }; + + const amrex::Real VxB[3] = { + -0.5*(By(j,k+1,0)-By(j,k-1,0) )*dzi, + 0.5*( (Bx(j,k+1,0)-Bx(j,k-1,0) )*dzi - (Bz(j+1,k,0)-Bz(j-1,k,0) )*dxi ), + 0.5*(By(j+1,k,0)-By(j-1,k,0) )*dxi, + }; + + // Defining comapct values for QED corrections + + const amrex::Real ex = tmpEx(j,k,0); + const amrex::Real ey = tmpEy(j,k,0); + const amrex::Real ez = tmpEz(j,k,0); + const amrex::Real bx = Bx(j,k,0); + const amrex::Real by = By(j,k,0); + const amrex::Real bz = Bz(j,k,0); + const amrex::Real ee = ex*ex + ey*ey + ez*ez; + const amrex::Real bb = bx*bx + by*by + bz*bz; + const amrex::Real eb = ex*bx + ey*by + ez*bz; + const amrex::Real EVxE = ex*VxE[0] + ey*VxE[1] + ez*VxE[2]; + const amrex::Real BVxE = bx*VxE[0] + by*VxE[1] + bz*VxE[2]; + const amrex::Real EVxB = ex*VxB[0] + ey*VxB[1] + ez*VxB[2]; + const amrex::Real BVxB = bx*VxB[0] + by*VxB[1] + bz*VxB[2]; + + const amrex::Real beta = 4*xi*( ee - c2*bb ) + PhysConst::ep0; + + const amrex::Real Alpha[3] = { + 2*xi*c2*( -7*bx*EVxE - 7*VxE[0]*eb + 4*ex*BVxE ) + VxM[0], + 2*xi*c2*( -7*by*EVxE - 7*VxE[1]*eb + 4*ey*BVxE ) + VxM[1], + 2*xi*c2*( -7*bz*EVxE - 7*VxE[2]*eb + 4*ez*BVxE ) + VxM[2] + }; + + const amrex::Real Omega[3] = { + Alpha[0] + 2*xi*c2*( 4*ex*EVxB + 2*VxB[0]*( ee - c2*bb ) + 7*c2*bx*BVxB ), + Alpha[1] + 2*xi*c2*( 4*ey*EVxB + 2*VxB[1]*( ee - c2*bb ) + 7*c2*by*BVxB ), + Alpha[2] + 2*xi*c2*( 4*ez*EVxB + 2*VxB[2]*( ee - c2*bb ) + 7*c2*bz*BVxB ) + }; + + // Calcualting matrix values for the QED correction algorithm + + const amrex::Real a00 = beta + xi*( 8*ex*ex + 14*c2*bx*bx ); + + const amrex::Real a11 = beta + xi*( 8*ey*ey + 14*c2*by*by ); + + const amrex::Real a22 = beta + xi*( 8*ez*ez + 14*c2*bz*bz ); + + const amrex::Real a01 = xi*( 2*ex*ey + 14*c2*bx*by ); + + const amrex::Real a02 = xi*( 2*ex*ez + 14*c2*bx*bz ); + + const amrex::Real a12 = xi*( 2*ez*ey + 14*c2*bz*by ); + + const amrex::Real detA = a00*( a11*a22 - a12*a12 ) - a01*( a01*a22 - a02*a12 ) + a02*( a01*a12 - a02*a11 ); + + // Calcualting inverse matrix values using the general 3x3 form + + const amrex::Real invAx[3] = {a22*a11 - a12*a12, a12*a02 - a22*a01, a12*a01 - a11*a02}; + + const amrex::Real invAy[3] = {a02*a12 - a22*a01, a00*a22 - a02*a02, a01*a02 - a12*a00}; + + const amrex::Real invAz[3] = {a12*a01 - a02*a11, a02*a01 - a12*a00, a11*a00 - a01*a01}; + + // Calcualting the final QED corrections by mutliplying the Omega vector with the inverse matrix + + const amrex::Real dEx = (-1/detA)*(invAx[0]*Omega[0] + + invAx[1]*Omega[1] + + invAx[2]*Omega[2]); + + const amrex::Real dEy = (-1/detA)*(invAy[0]*Omega[0] + + invAy[1]*Omega[1] + + invAy[2]*Omega[2]); + + const amrex::Real dEz = (-1/detA)*(invAz[0]*Omega[0] + + invAz[1]*Omega[1] + + invAz[2]*Omega[2]); + + // Adding the QED corrections to the origional fields + + Ex(j,k,0) = Ex(j,k,0) + 0.5*dt*dEx; + + Ey(j,k,0) = Ey(j,k,0) + 0.5*dt*dEy; + + Ez(j,k,0) = Ez(j,k,0) + 0.5*dt*dEz; + +#endif + +} + +#endif diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 0d56783ed..80d5e7f16 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -64,16 +64,16 @@ public: virtual void PostRestart () final; - void calculate_laser_plane_coordinates (const int np, const int thread_num, + void calculate_laser_plane_coordinates (const WarpXParIter& pti, const int np, amrex::Real * AMREX_RESTRICT const pplane_Xp, amrex::Real * AMREX_RESTRICT const pplane_Yp); - void update_laser_particle (const int np, amrex::ParticleReal * AMREX_RESTRICT const puxp, + void update_laser_particle (WarpXParIter& pti, const int np, amrex::ParticleReal * AMREX_RESTRICT const puxp, amrex::ParticleReal * AMREX_RESTRICT const puyp, amrex::ParticleReal * AMREX_RESTRICT const puzp, amrex::ParticleReal const * AMREX_RESTRICT const pwp, amrex::Real const * AMREX_RESTRICT const amplitude, - const amrex::Real dt, const int thread_num); + const amrex::Real dt); protected: diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index de192e65e..2bc522a8e 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -16,6 +16,7 @@ #include <WarpX_Complex.H> #include <WarpX_f.H> #include <MultiParticleContainer.H> +#include <GetAndSetPosition.H> using namespace amrex; using namespace WarpXLaserProfiles; @@ -432,13 +433,6 @@ LaserParticleContainer::Evolve (int lev, plane_Yp.resize(np); amplitude_E.resize(np); - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); - if (rho) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 0, 0, @@ -454,7 +448,7 @@ LaserParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_pp); // Find the coordinates of the particles in the emission plane - calculate_laser_plane_coordinates(np, thread_num, + calculate_laser_plane_coordinates(pti, np, plane_Xp.dataPtr(), plane_Yp.dataPtr()); @@ -465,9 +459,9 @@ LaserParticleContainer::Evolve (int lev, t_lab, amplitude_E.dataPtr()); // Calculate the corresponding momentum and position for the particles - update_laser_particle(np, uxp.dataPtr(), uyp.dataPtr(), + update_laser_particle(pti, np, uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), wp.dataPtr(), - amplitude_E.dataPtr(), dt, thread_num); + amplitude_E.dataPtr(), dt); BL_PROFILE_VAR_STOP(blp_pp); // @@ -489,13 +483,6 @@ LaserParticleContainer::Evolve (int lev, } } - // - // copy particle data back - // - BL_PROFILE_VAR_START(blp_copy); - pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); - if (rho) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 1, 0, @@ -585,14 +572,12 @@ LaserParticleContainer::PushP (int lev, Real dt, * in laser plane coordinate. */ void -LaserParticleContainer::calculate_laser_plane_coordinates ( - const int np, const int thread_num, - Real * AMREX_RESTRICT const pplane_Xp, - Real * AMREX_RESTRICT const pplane_Yp) +LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& pti, const int np, + Real * AMREX_RESTRICT const pplane_Xp, + Real * AMREX_RESTRICT const pplane_Yp) { - ParticleReal const * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr(); - ParticleReal const * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr(); - ParticleReal const * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); + Real tmp_u_X_0 = u_X[0]; Real tmp_u_X_2 = u_X[2]; Real tmp_position_0 = position[0]; @@ -608,19 +593,21 @@ LaserParticleContainer::calculate_laser_plane_coordinates ( amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) { + ParticleReal x, y, z; + GetPosition(i, x, y, z); #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) pplane_Xp[i] = - tmp_u_X_0 * (xp[i] - tmp_position_0) + - tmp_u_X_1 * (yp[i] - tmp_position_1) + - tmp_u_X_2 * (zp[i] - tmp_position_2); + tmp_u_X_0 * (x - tmp_position_0) + + tmp_u_X_1 * (y - tmp_position_1) + + tmp_u_X_2 * (z - tmp_position_2); pplane_Yp[i] = - tmp_u_Y_0 * (xp[i] - tmp_position_0) + - tmp_u_Y_1 * (yp[i] - tmp_position_1) + - tmp_u_Y_2 * (zp[i] - tmp_position_2); + tmp_u_Y_0 * (x - tmp_position_0) + + tmp_u_Y_1 * (y - tmp_position_1) + + tmp_u_Y_2 * (z - tmp_position_2); #elif (AMREX_SPACEDIM == 2) pplane_Xp[i] = - tmp_u_X_0 * (xp[i] - tmp_position_0) + - tmp_u_X_2 * (zp[i] - tmp_position_2); + tmp_u_X_0 * (x - tmp_position_0) + + tmp_u_X_2 * (z - tmp_position_2); pplane_Yp[i] = 0.; #endif } @@ -629,23 +616,26 @@ LaserParticleContainer::calculate_laser_plane_coordinates ( /* \brief push laser particles, in simulation coordinates. * + * \param pti: Particle iterator * \param np: number of laser particles * \param puxp, puyp, puzp: pointers to arrays of particle momenta. * \param pwp: pointer to array of particle weights. * \param amplitude: Electric field amplitude at the position of each particle. * \param dt: time step. - * \param thread_num: thread number */ void -LaserParticleContainer::update_laser_particle( - const int np, ParticleReal * AMREX_RESTRICT const puxp, ParticleReal * AMREX_RESTRICT const puyp, - ParticleReal * AMREX_RESTRICT const puzp, ParticleReal const * AMREX_RESTRICT const pwp, - Real const * AMREX_RESTRICT const amplitude, const Real dt, - const int thread_num) +LaserParticleContainer::update_laser_particle(WarpXParIter& pti, + const int np, + ParticleReal * AMREX_RESTRICT const puxp, + ParticleReal * AMREX_RESTRICT const puyp, + ParticleReal * AMREX_RESTRICT const puzp, + ParticleReal const * AMREX_RESTRICT const pwp, + Real const * AMREX_RESTRICT const amplitude, + const Real dt) { - ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr(); - ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr(); - ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + Real tmp_p_X_0 = p_X[0]; Real tmp_p_X_1 = p_X[1]; Real tmp_p_X_2 = p_X[2]; @@ -678,12 +668,16 @@ LaserParticleContainer::update_laser_particle( puxp[i] = gamma * vx; puyp[i] = gamma * vy; puzp[i] = gamma * vz; + // Push the the particle positions - xp[i] += vx * dt; + ParticleReal x, y, z; + GetPosition(i, x, y, z); + x += vx * dt; #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) - yp[i] += vy * dt; + y += vy * dt; #endif - zp[i] += vz * dt; + z += vz * dt; + SetPosition(i, x, y, z); } ); } diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 669fc4c57..b03e4224f 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -8,10 +8,11 @@ #ifndef CHARGEDEPOSITION_H_ #define CHARGEDEPOSITION_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" /* \brief Charge Deposition for thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param wp : Pointer to array of particle weights. * \param ion_lev : Pointer to array of particle ionization level. This is required to have the charge of each macroparticle @@ -25,9 +26,7 @@ * /param q : species charge. */ template <int depos_order> -void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doChargeDepositionShapeN(const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, const int * const ion_lev, const amrex::Array4<amrex::Real>& rho_arr, @@ -63,14 +62,17 @@ void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, wq *= ion_lev[ip]; } + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + // --- Compute shape factors // x direction // Get particle position in grid coordinates #if (defined WARPX_DIM_RZ) - const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real r = std::sqrt(xp*xp + yp*yp); const amrex::Real x = (r - xmin)*dxi; #else - const amrex::Real x = (xp[ip] - xmin)*dxi; + const amrex::Real x = (xp - xmin)*dxi; #endif // Compute shape factors for node-centered quantities amrex::Real sx[depos_order + 1]; @@ -79,12 +81,12 @@ void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_3D) // y direction - const amrex::Real y = (yp[ip] - ymin)*dyi; + const amrex::Real y = (yp - ymin)*dyi; amrex::Real sy[depos_order + 1]; const int j = compute_shape_factor<depos_order>(sy, y); #endif // z direction - const amrex::Real z = (zp[ip] - zmin)*dzi; + const amrex::Real z = (zp - zmin)*dzi; amrex::Real sz[depos_order + 1]; const int k = compute_shape_factor<depos_order>(sz, z); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 90039dea2..af3b0006b 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -8,6 +8,7 @@ #ifndef CURRENTDEPOSITION_H_ #define CURRENTDEPOSITION_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" #include <WarpX_Complex.H> @@ -16,7 +17,7 @@ /** * \brief Current Deposition for thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. * \param ion_lev : Pointer to array of particle ionization level. This is @@ -34,9 +35,7 @@ * /param q : species charge. */ template <int depos_order> -void doDepositionShapeN(const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doDepositionShapeN(const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, const amrex::ParticleReal * const uxp, const amrex::ParticleReal * const uyp, @@ -95,6 +94,10 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, if (do_ionization){ wq *= ion_lev[ip]; } + + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + const amrex::Real vx = uxp[ip]*gaminv; const amrex::Real vy = uyp[ip]*gaminv; const amrex::Real vz = uzp[ip]*gaminv; @@ -102,8 +105,8 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #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 xpmid = xp - 0.5*dt*vx; + const amrex::Real ypmid = yp - 0.5*dt*vy; const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid); amrex::Real costheta; amrex::Real sintheta; @@ -128,7 +131,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_RZ) const amrex::Real xmid = (rpmid - xmin)*dxi; #else - const amrex::Real xmid = (xp[ip] - xmin)*dxi - dts2dx*vx; + const amrex::Real xmid = (xp - xmin)*dxi - dts2dx*vx; #endif // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current // sx_j[xyz] shape factor along x for the centering of each current @@ -153,7 +156,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #if (defined WARPX_DIM_3D) // y direction - const amrex::Real ymid = (yp[ip] - ymin)*dyi - dts2dy*vy; + const amrex::Real ymid = (yp - ymin)*dyi - dts2dy*vy; amrex::Real sy_node[depos_order + 1]; amrex::Real sy_cell[depos_order + 1]; int k_node; @@ -173,7 +176,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, #endif // z direction - const amrex::Real zmid = (zp[ip] - zmin)*dzi - dts2dz*vz; + const amrex::Real zmid = (zp - zmin)*dzi - dts2dz*vz; amrex::Real sz_node[depos_order + 1]; amrex::Real sz_cell[depos_order + 1]; int l_node; @@ -230,7 +233,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, /** * \brief Esirkepov Current Deposition for thread thread_num * - * \param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. * \param ion_lev : Pointer to array of particle ionization level. This is @@ -249,9 +252,7 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order> -void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, const amrex::ParticleReal * const uxp, const amrex::ParticleReal * const uyp, @@ -308,14 +309,18 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, // --- Get particle quantities Real const 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 Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; } + + ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + Real const wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) Real const wqy = wq*invdtdy; @@ -324,17 +329,18 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, // computes current and old position in grid units #if (defined WARPX_DIM_RZ) - 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 xp_mid = xp - 0.5_rt * dt*uxp[ip]*gaminv; + Real const yp_mid = yp - 0.5_rt * dt*uyp[ip]*gaminv; + Real const xp_old = xp - dt*uxp[ip]*gaminv; + Real const yp_old = yp - dt*uyp[ip]*gaminv; + Real const rp_new = std::sqrt(xp*xp + + yp*yp); 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; + costheta_new = xp/rp_new; + sintheta_new = yp/rp_new; } else { costheta_new = 1.; sintheta_new = 0.; @@ -361,14 +367,14 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, Real const x_new = (rp_new - xmin)*dxi; Real const x_old = (rp_old - xmin)*dxi; #else - Real const x_new = (xp[ip] - xmin)*dxi; + Real const x_new = (xp - xmin)*dxi; Real const x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif #if (defined WARPX_DIM_3D) - Real const y_new = (yp[ip] - ymin)*dyi; + Real const y_new = (yp - ymin)*dyi; Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv; #endif - Real const z_new = (zp[ip] - zmin)*dzi; + Real const z_new = (zp - zmin)*dzi; Real const z_old = z_new - dtsdz0*uzp[ip]*gaminv; #if (defined WARPX_DIM_RZ) diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 0a58f8425..12d9b6291 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -8,12 +8,13 @@ #ifndef FIELDGATHER_H_ #define FIELDGATHER_H_ +#include "GetAndSetPosition.H" #include "ShapeFactors.H" #include <WarpX_Complex.H> /** * \brief Field gather for particles handled by thread thread_num - * \param xp, yp, zp : Pointer to arrays of particle positions. + * /param GetPosition : A functor for returning the particle position. * \param Exp, Eyp, Ezp: Pointer to array of electric field on particles. * \param Bxp, Byp, Bzp: Pointer to array of magnetic field on particles. * \param ex_arr ey_arr: Array4 of current density, either full array or tile. @@ -27,9 +28,7 @@ * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order, int lower_in_v> -void doGatherShapeN(const amrex::ParticleReal * const xp, - const amrex::ParticleReal * const yp, - const amrex::ParticleReal * const zp, +void doGatherShapeN(const GetParticlePosition& GetPosition, amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, @@ -80,14 +79,18 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, amrex::ParallelFor( np_to_gather, [=] AMREX_GPU_DEVICE (long ip) { + + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + // --- Compute shape factors // x direction // Get particle position #ifdef WARPX_DIM_RZ - const amrex::Real rp = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real rp = std::sqrt(xp*xp + yp*yp); const amrex::Real x = (rp - xmin)*dxi; #else - const amrex::Real x = (xp[ip]-xmin)*dxi; + const amrex::Real x = (xp-xmin)*dxi; #endif // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current @@ -129,7 +132,7 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, #if (AMREX_SPACEDIM == 3) // y direction - const amrex::Real y = (yp[ip]-ymin)*dyi; + const amrex::Real y = (yp-ymin)*dyi; amrex::Real sy_node[depos_order + 1]; amrex::Real sy_cell[depos_order + 1]; amrex::Real sy_node_v[depos_order + 1 - lower_in_v]; @@ -165,7 +168,7 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, #endif // z direction - const amrex::Real z = (zp[ip]-zmin)*dzi; + const amrex::Real z = (zp-zmin)*dzi; amrex::Real sz_node[depos_order + 1]; amrex::Real sz_cell[depos_order + 1]; amrex::Real sz_node_v[depos_order + 1 - lower_in_v]; @@ -245,8 +248,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, amrex::Real costheta; amrex::Real sintheta; if (rp > 0.) { - costheta = xp[ip]/rp; - sintheta = yp[ip]/rp; + costheta = xp/rp; + sintheta = yp/rp; } else { costheta = 1.; sintheta = 0.; diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 7750d5ce8..ef32ebcb0 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -55,9 +55,6 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; // Do nothing diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index ab85170ac..d007d3c1c 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -21,7 +21,7 @@ // Import low-level single-particle kernels #include <UpdatePositionPhoton.H> - +#include <GetAndSetPosition.H> using namespace amrex; @@ -57,19 +57,13 @@ void PhotonParticleContainer::InitData() } void -PhotonParticleContainer::PushPX(WarpXParIter& pti, - Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp, - Real dt, DtType a_dt_type) +PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); + // Extract pointers to the different particle quantities - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); @@ -82,15 +76,19 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { - copy_attribs(pti, x, y, z); + copy_attribs(pti); } + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - - UpdatePositionPhoton( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePositionPhoton( x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 0500d091a..421dcd842 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -72,10 +72,7 @@ public: void AssignExternalFieldOnParticles ( WarpXParIter& pti, RealVector& Exp, RealVector& Eyp, RealVector& Ezp, RealVector& Bxp, - RealVector& Byp, RealVector& Bzp, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, int lev); + RealVector& Byp, RealVector& Bzp, int lev); virtual void FieldGather (int lev, const amrex::MultiFab& Ex, @@ -101,7 +98,6 @@ public: const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, - int thread_num, int lev, int depos_lev); @@ -163,11 +159,7 @@ public: amrex::Real dt, DtType a_dt_type=DtType::Full) override; - virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, - amrex::Real dt, DtType a_dt_type=DtType::Full); + virtual void PushPX (WarpXParIter& pti, amrex::Real dt, DtType a_dt_type=DtType::Full); virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, @@ -190,8 +182,7 @@ public: RealVector& uzp, RealVector& wp ); - void copy_attribs(WarpXParIter& pti,const amrex::ParticleReal* xp, - const amrex::ParticleReal* yp, const amrex::ParticleReal* zp); + void copy_attribs (WarpXParIter& pti); virtual void PostRestart () final {} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4c98f436e..6572657ff 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -20,6 +20,7 @@ #include <WarpXWrappers.h> #include <IonizationEnergiesTable.H> #include <FieldGather.H> +#include <GetAndSetPosition.H> #include <WarpXAlgorithmSelection.H> @@ -967,10 +968,7 @@ PhysicalParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<Mul void PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, RealVector& Exp, RealVector& Eyp, RealVector& Ezp, - RealVector& Bxp, RealVector& Byp, RealVector& Bzp, - const Gpu::ManagedDeviceVector<ParticleReal>& xp, - const Gpu::ManagedDeviceVector<ParticleReal>& yp, - const Gpu::ManagedDeviceVector<ParticleReal>& zp, int lev) + RealVector& Bxp, RealVector& Byp, RealVector& Bzp, int lev) { const long np = pti.numParticles(); /// get WarpX class object @@ -990,9 +988,7 @@ PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, Bzp.assign(np,mypc.m_B_external_particle[2]); } if (mypc.m_E_ext_particle_s=="parse_e_ext_particle_function") { - const Real* const AMREX_RESTRICT xp_data = xp.dataPtr(); - const Real* const AMREX_RESTRICT yp_data = yp.dataPtr(); - const Real* const AMREX_RESTRICT zp_data = zp.dataPtr(); + const auto GetPosition = GetParticlePosition(pti); Real* const AMREX_RESTRICT Exp_data = Exp.dataPtr(); Real* const AMREX_RESTRICT Eyp_data = Eyp.dataPtr(); Real* const AMREX_RESTRICT Ezp_data = Ezp.dataPtr(); @@ -1001,20 +997,20 @@ PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, ParserWrapper *zfield_partparser = mypc.m_Ez_particle_parser.get(); Real time = warpx.gett_new(lev); amrex::ParallelFor(pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Exp_data[i] = xfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Eyp_data[i] = yfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Ezp_data[i] = zfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - }, + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal x, y, z; + GetPosition(i, x, y, z); + Exp_data[i] = xfield_partparser->getField(x, y, z, time); + Eyp_data[i] = yfield_partparser->getField(x, y, z, time); + Ezp_data[i] = zfield_partparser->getField(x, y, z, time); + }, /* To allocate shared memory for the GPU threads. */ /* But, for now only 4 doubles (x,y,z,t) are allocated. */ amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 ); } if (mypc.m_B_ext_particle_s=="parse_b_ext_particle_function") { - const Real* const AMREX_RESTRICT xp_data = xp.dataPtr(); - const Real* const AMREX_RESTRICT yp_data = yp.dataPtr(); - const Real* const AMREX_RESTRICT zp_data = zp.dataPtr(); + const auto GetPosition = GetParticlePosition(pti); Real* const AMREX_RESTRICT Bxp_data = Bxp.dataPtr(); Real* const AMREX_RESTRICT Byp_data = Byp.dataPtr(); Real* const AMREX_RESTRICT Bzp_data = Bzp.dataPtr(); @@ -1024,10 +1020,12 @@ PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, Real time = warpx.gett_new(lev); amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - Bxp_data[i] = xfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Byp_data[i] = yfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - Bzp_data[i] = zfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); - }, + ParticleReal x, y, z; + GetPosition(i, x, y, z); + Bxp_data[i] = xfield_partparser->getField(x, y, z, time); + Byp_data[i] = yfield_partparser->getField(x, y, z, time); + Bzp_data[i] = zfield_partparser->getField(x, y, z, time); + }, /* To allocate shared memory for the GPU threads. */ /* But, for now only 4 doubles (x,y,z,t) are allocated. */ amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 @@ -1056,11 +1054,6 @@ PhysicalParticleContainer::FieldGather (int lev, #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { Real wt = amrex::second(); @@ -1087,18 +1080,13 @@ PhysicalParticleContainer::FieldGather (int lev, const FArrayBox& bzfab = Bz[pti]; // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - - // // Field Gather // 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); + 0, np, lev, lev); if (cost) { const Box& tbx = pti.tilebox(); @@ -1110,9 +1098,6 @@ PhysicalParticleContainer::FieldGather (int lev, costarr(i,j,k) += wt; }); } - // synchronize avoids cudaStreams from over-writing the temporary arrays used to - // store positions - Gpu::synchronize(); } } } @@ -1238,13 +1223,6 @@ PhysicalParticleContainer::Evolve (int lev, const long np_current = (cjx) ? nfine_current : np; - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); - if (rho) { // Deposit charge before particle push, in component 0 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1274,7 +1252,7 @@ PhysicalParticleContainer::Evolve (int lev, 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); + 0, np_gather, lev, lev); if (np_gather < np) { @@ -1310,7 +1288,7 @@ PhysicalParticleContainer::Evolve (int lev, cbxfab, cbyfab, cbzfab, cEx->nGrow(), e_is_nodal, nfine_gather, np-nfine_gather, - thread_num, lev, lev-1); + lev, lev-1); } BL_PROFILE_VAR_STOP(blp_fg); @@ -1328,7 +1306,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], dt, a_dt_type); + PushPX(pti, dt, a_dt_type); BL_PROFILE_VAR_STOP(blp_ppc_pp); // @@ -1352,14 +1330,6 @@ PhysicalParticleContainer::Evolve (int lev, np_current, np-np_current, thread_num, lev, lev-1, dt); } - - - // - // copy particle data back - // - BL_PROFILE_VAR_START(blp_copy); - pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); } if (rho) { @@ -1478,7 +1448,6 @@ PhysicalParticleContainer::SplitParticles(int lev) { auto& mypc = WarpX::GetInstance().GetPartContainer(); auto& pctmp_split = mypc.GetPCtmp(); - Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; RealVector psplit_x, psplit_y, psplit_z, psplit_w; RealVector psplit_ux, psplit_uy, psplit_uz; long np_split_to_add = 0; @@ -1493,7 +1462,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // Loop over particle interator for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - pti.GetPosition(xp, yp, zp); + const auto GetPosition = GetParticlePosition(pti); const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); @@ -1518,6 +1487,8 @@ PhysicalParticleContainer::SplitParticles(int lev) auto& uzp = attribs[PIdx::uz]; const long np = pti.numParticles(); for(int i=0; i<np; i++){ + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); auto& p = particles[i]; if (p.id() == DoSplitParticleID){ // If particle is tagged, split it and put the @@ -1530,9 +1501,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int ishift = -1; ishift < 2; ishift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x and z - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + kshift*split_offset[2] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1544,17 +1515,17 @@ PhysicalParticleContainer::SplitParticles(int lev) // 4 particles in 2d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in z - psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*split_offset[2] ); + psplit_x.push_back( xp ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1569,9 +1540,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int jshift = -1; jshift < 2; jshift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x, y and z - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] + jshift*split_offset[1] ); - psplit_z.push_back( zp[i] + kshift*split_offset[2] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp + jshift*split_offset[1] ); + psplit_z.push_back( zp + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1584,25 +1555,25 @@ PhysicalParticleContainer::SplitParticles(int lev) // 6 particles in 3d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*split_offset[0] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] ); + psplit_x.push_back( xp + ishift*split_offset[0] ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in y - psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] + ishift*split_offset[1] ); - psplit_z.push_back( zp[i] ); + psplit_x.push_back( xp ); + psplit_y.push_back( yp + ishift*split_offset[1] ); + psplit_z.push_back( zp ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in z - psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*split_offset[2] ); + psplit_x.push_back( xp ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1639,19 +1610,16 @@ PhysicalParticleContainer::SplitParticles(int lev) } void -PhysicalParticleContainer::PushPX(WarpXParIter& pti, - Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp, - Real dt, DtType a_dt_type) +PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); @@ -1664,7 +1632,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) { - copy_attribs(pti, x, y, z); + copy_attribs(pti); } int* AMREX_RESTRICT ion_lev = nullptr; @@ -1697,8 +1665,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); } - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); }else{ @@ -1708,8 +1678,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); } @@ -1723,8 +1695,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); #endif @@ -1737,8 +1711,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumBoris( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { @@ -1750,8 +1726,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumVay( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { @@ -1763,8 +1741,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt ); + SetPosition(i, x, y, z); } ); } else { @@ -1830,11 +1810,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt, #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1858,16 +1833,11 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - 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); + 0, np, lev, lev); // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities @@ -1944,8 +1914,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, - const ParticleReal* yp, const ParticleReal* zp) +void PhysicalParticleContainer::copy_attribs (WarpXParIter& pti) { auto& attribs = pti.GetAttribs(); ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); @@ -1962,11 +1931,15 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleRea ParticleReal* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); ParticleReal* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { - xpold[i]=xp[i]; - ypold[i]=yp[i]; - zpold[i]=zp[i]; + ParticleReal x, y, z; + GetPosition(i, x, y, z); + xpold[i]=x; + ypold[i]=y; + zpold[i]=z; uxpold[i]=uxp[i]; uypold[i]=uyp[i]; @@ -2024,11 +1997,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -2037,10 +2005,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real if ( !slice_box.intersects(tile_real_box) ) continue; - pti.GetPosition(m_xp[thread_num],m_yp[thread_num],m_zp[thread_num]); - Real *const AMREX_RESTRICT xpnew = m_xp[thread_num].dataPtr(); - Real *const AMREX_RESTRICT ypnew = m_yp[thread_num].dataPtr(); - Real *const AMREX_RESTRICT zpnew = m_zp[thread_num].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); auto& attribs = pti.GetAttribs(); Real* const AMREX_RESTRICT wpnew = attribs[PIdx::w].dataPtr(); @@ -2078,12 +2043,14 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) { - Flag[i] = 0; - if ( (((zpnew[i] >= z_new) && (zpold[i] <= z_old)) || - ((zpnew[i] <= z_new) && (zpold[i] >= z_old))) ) - { - Flag[i] = 1; - } + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + Flag[i] = 0; + if ( (((zp >= z_new) && (zpold[i] <= z_old)) || + ((zp <= z_new) && (zpold[i] >= z_old))) ) + { + Flag[i] = 1; + } }); // exclusive scan to obtain location indices using flag values @@ -2120,6 +2087,8 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) { + ParticleReal xp_new, yp_new, zp_new; + GetPosition(i, xp_new, yp_new, zp_new); if (Flag[i] == 1) { // Lorentz Transform particles to lab-frame @@ -2127,12 +2096,10 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real (uxpnew[i]*uxpnew[i] + uypnew[i]*uypnew[i] + uzpnew[i]*uzpnew[i])); - const Real t_new_p = gammaboost*t_boost - - uzfrm*zpnew[i]*inv_c2; - const Real z_new_p = gammaboost*(zpnew[i] - + betaboost*Phys_c*t_boost); - const Real uz_new_p = gammaboost*uzpnew[i] - - gamma_new_p*uzfrm; + const Real t_new_p = gammaboost*t_boost - uzfrm*zp_new*inv_c2; + const Real z_new_p = gammaboost*(zp_new + betaboost*Phys_c*t_boost); + const Real uz_new_p = gammaboost*uzpnew[i] - gamma_new_p*uzfrm; + const Real gamma_old_p = std::sqrt(1.0 + inv_c2* (uxpold[i]*uxpold[i] + uypold[i]*uypold[i] @@ -2150,12 +2117,10 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real const Real weight_new = (t_lab - t_old_p) / (t_new_p - t_old_p); - const Real xp = xpold[i]*weight_old - + xpnew[i]*weight_new; - const Real yp = ypold[i]*weight_old - + ypnew[i]*weight_new; - const Real zp = z_old_p*weight_old - + z_new_p*weight_new; + const Real xp = xpold[i]*weight_old + xp_new*weight_new; + const Real yp = ypold[i]*weight_old + yp_new*weight_new; + const Real zp = z_old_p*weight_old + z_new_p*weight_new; + const Real uxp = uxpold[i]*weight_old + uxpnew[i]*weight_new; const Real uyp = uypold[i]*weight_old @@ -2197,7 +2162,6 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param e_is_nodal: 0 if E is staggered, 1 if E is nodal * \param offset: index of first particle for which fields are gathered * \param np_to_gather: number of particles onto which fields are gathered - * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. @@ -2219,7 +2183,6 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, - int thread_num, int lev, int gather_lev) { @@ -2231,9 +2194,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // initializing the field value to the externally applied field before // gathering fields from the grid to the particles. - AssignExternalFieldOnParticles(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - m_xp[thread_num], m_yp[thread_num], - m_zp[thread_num], lev); + AssignExternalFieldOnParticles(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, lev); // Get cell size on gather_lev @@ -2252,9 +2213,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // Add guard cells to the box. box.grow(ngE); - const ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - const ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - const ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const auto GetPosition = GetParticlePosition(pti, offset); // Lower corner of tile box physical domain const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev); @@ -2265,7 +2224,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ if (WarpX::nox == 1){ - doGatherShapeN<1,1>(xp, yp, zp, + doGatherShapeN<1,1>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2273,7 +2232,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doGatherShapeN<2,1>(xp, yp, zp, + doGatherShapeN<2,1>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2281,7 +2240,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doGatherShapeN<3,1>(xp, yp, zp, + doGatherShapeN<3,1>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2291,7 +2250,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } else { if (WarpX::nox == 1){ - doGatherShapeN<1,0>(xp, yp, zp, + doGatherShapeN<1,0>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2299,7 +2258,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doGatherShapeN<2,0>(xp, yp, zp, + doGatherShapeN<2,0>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, @@ -2307,7 +2266,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, np_to_gather, dx, xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doGatherShapeN<3,0>(xp, yp, zp, + doGatherShapeN<3,0>(GetPosition, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 7180f3c55..594260703 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -12,72 +12,99 @@ #include <WarpXParticleContainer.H> #include <AMReX_REAL.H> -#ifndef WARPX_DIM_RZ - -/** \brief Extract the particle's coordinates from the ParticleType struct `p`, - * and stores them in the variables `x`, `y`, `z`. */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void GetPosition( - amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const WarpXParticleContainer::ParticleType& p) +/** \brief Functor that can be used to extract the positions of the macroparticles + * inside a ParallelFor kernel + * + * \param a_pti iterator to the tile containing the macroparticles + * \param a_offset offset to apply to the particle indices +*/ +struct GetParticlePosition { -#if (AMREX_SPACEDIM==3) - x = p.pos(0); - y = p.pos(1); - z = p.pos(2); -#else - x = p.pos(0); - y = std::numeric_limits<amrex::ParticleReal>::quiet_NaN(); - z = p.pos(1); + using PType = WarpXParticleContainer::ParticleType; + using RType = amrex::ParticleReal; + + const PType* AMREX_RESTRICT m_structs; +#if (defined WARPX_DIM_RZ) + const RType* m_theta; +#elif (AMREX_SPACEDIM == 2) + static constexpr RType m_snan = std::numeric_limits<RType>::quiet_NaN(); #endif -} + GetParticlePosition (const WarpXParIter& a_pti, int a_offset = 0) noexcept + { + const auto& aos = a_pti.GetArrayOfStructs(); + m_structs = aos().dataPtr() + a_offset; +#if (defined WARPX_DIM_RZ) + const auto& soa = a_pti.GetStructOfArrays(); + m_theta = soa.GetRealData(PIdx::theta).dataPtr() + a_offset; +#endif + } -/** \brief Set the particle's coordinates in the ParticleType struct `p`, - * from their values in the variables `x`, `y`, `z`. */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void SetPosition( - WarpXParticleContainer::ParticleType& p, - const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z) -{ -#if (AMREX_SPACEDIM==3) - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; + /** \brief Extract the cartesian position coordinates of the particle + * located at index `i + a_offset` and store them in the variables + * `x`, `y`, `z` */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (const int i, RType& x, RType& y, RType& z) const noexcept + { +#ifdef WARPX_DIM_RZ + RType r = m_structs[i].pos(0); + x = r*std::cos(m_theta[i]); + y = r*std::sin(m_theta[i]); + z = m_structs[i].pos(1); +#elif WARPX_DIM_3D + x = m_structs[i].pos(0); + y = m_structs[i].pos(1); + z = m_structs[i].pos(2); #else - p.pos(0) = x; - p.pos(1) = z; + x = m_structs[i].pos(0); + y = m_snan; + z = m_structs[i].pos(1); #endif -} - -# 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), - * and store them in the variables `x`, `y`, `z` */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void GetCartesianPositionFromCylindrical( - amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const WarpXParticleContainer::ParticleType& p, const amrex::ParticleReal theta) +/** \brief Functor that can be used to modify the positions of the macroparticles, + * inside a ParallelFor kernel. + * + * \param a_pti iterator to the tile being modified + * \param a_offset offset to apply to the particle indices +*/ +struct SetParticlePosition { - const amrex::ParticleReal r = p.pos(0); - x = r*std::cos(theta); - y = r*std::sin(theta); - z = p.pos(1); -} + using PType = WarpXParticleContainer::ParticleType; + using RType = amrex::ParticleReal; -/** \brief Set the particle's cylindrical coordinates by setting `theta` - * and the attributes of the ParticleType struct `p` (which stores the radius), - * from the values of `x`, `y`, `z` */ -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void SetCylindricalPositionFromCartesian( - WarpXParticleContainer::ParticleType& p, amrex::ParticleReal& theta, - const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z ) -{ - theta = std::atan2(y, x); - p.pos(0) = std::sqrt(x*x + y*y); - p.pos(1) = z; -} + PType* AMREX_RESTRICT m_structs; +#if (defined WARPX_DIM_RZ) + RType* AMREX_RESTRICT m_theta; +#endif + SetParticlePosition (WarpXParIter& a_pti, int a_offset = 0) noexcept + { + auto& aos = a_pti.GetArrayOfStructs(); + m_structs = aos().dataPtr() + a_offset; +#if (defined WARPX_DIM_RZ) + auto& soa = a_pti.GetStructOfArrays(); + m_theta = soa.GetRealData(PIdx::theta).dataPtr() + a_offset; +#endif + } -#endif // WARPX_DIM_RZ + /** \brief Set the position of the particle at index `i + a_offset` + * to `x`, `y`, `z` */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (const int i, RType x, RType y, RType z) const noexcept + { +#ifdef WARPX_DIM_RZ + m_theta[i] = std::atan2(y, x); + m_structs[i].pos(0) = std::sqrt(x*x + y*y); + m_structs[i].pos(1) = z; +#elif WARPX_DIM_3D + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = y; + m_structs[i].pos(2) = z; +#else + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = z; +#endif + } +}; #endif // WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 7f86a106d..1968a1439 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -15,10 +15,9 @@ /** \brief Push the particle's positions over one timestep, * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE -void UpdatePosition( - amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, - const amrex::Real dt ) +void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, + const amrex::Real dt ) { constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); @@ -31,7 +30,6 @@ void UpdatePosition( y += uy * inv_gamma * dt; #endif z += uz * inv_gamma * dt; - } #endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H index 3a28e87a1..44c0afda9 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -32,7 +32,6 @@ void UpdatePositionPhoton( y += uy * c_over_umod * dt; #endif z += uz * c_over_umod * dt; - } #endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index 5e5749093..7ad7fd82e 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -68,11 +68,7 @@ public: amrex::Real dt, DtType a_dt_type=DtType::Full) override; - virtual void PushPX(WarpXParIter& pti, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, - amrex::Real dt, DtType a_dt_type=DtType::Full) override; + virtual void PushPX (WarpXParIter& pti, amrex::Real dt, DtType a_dt_type=DtType::Full) override; virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index f9db682d7..953874c0f 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -24,6 +24,7 @@ #include <UpdateMomentumVay.H> #include <UpdateMomentumBorisWithRadiationReaction.H> #include <UpdateMomentumHigueraCary.H> +#include <GetAndSetPosition.H> using namespace amrex; @@ -87,8 +88,6 @@ RigidInjectedParticleContainer::RemapParticles() // Note that the particles are already in the boosted frame. // This value is saved to advance the particles not injected yet - Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { auto& attribs = pti.GetAttribs(); @@ -96,31 +95,32 @@ RigidInjectedParticleContainer::RemapParticles() auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; - // Copy data from particle container to temp arrays - pti.GetPosition(xp, yp, zp); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); // Loop over particles const long np = pti.numParticles(); for (int i=0 ; i < np ; i++) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + const Real gammapr = std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/csq); const Real vzpr = uzp[i]/gammapr; // Back out the value of z_lab - const Real z_lab = (zp[i] + uz_boost*t_lab + WarpX::gamma_boost*t_lab*vzpr)/(WarpX::gamma_boost + uz_boost*vzpr/csq); + const Real z_lab = (zp + uz_boost*t_lab + WarpX::gamma_boost*t_lab*vzpr)/(WarpX::gamma_boost + uz_boost*vzpr/csq); // Time of the particle in the boosted frame given its position in the lab frame at t=0. const Real tpr = WarpX::gamma_boost*t_lab - uz_boost*z_lab/csq; // Adjust the position, taking away its motion from its own velocity and adding // the motion from the average velocity - zp[i] = zp[i] + tpr*vzpr - tpr*vzbeam_ave_boosted; - - } + zp += tpr*vzpr - tpr*vzbeam_ave_boosted; - // Copy the data back to the particle container - pti.SetPosition(xp, yp, zp); + SetPosition(i, xp, yp, zp); + } } } } @@ -147,8 +147,6 @@ RigidInjectedParticleContainer::BoostandRemapParticles() #pragma omp parallel #endif { - Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; - for (WarpXParIter pti(*this, 0); pti.isValid(); ++pti) { @@ -157,13 +155,16 @@ RigidInjectedParticleContainer::BoostandRemapParticles() auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; - // Copy data from particle container to temp arrays - pti.GetPosition(xp, yp, zp); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); // Loop over particles const long np = pti.numParticles(); for (int i=0 ; i < np ; i++) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + const Real gamma_lab = std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/(PhysConst::c*PhysConst::c)); const Real vx_lab = uxp[i]/gamma_lab; @@ -173,24 +174,24 @@ RigidInjectedParticleContainer::BoostandRemapParticles() // t0_lab is the time in the lab frame that the particles reaches z=0 // The location and time (z=0, t=0) is a synchronization point between the // lab and boosted frames. - const Real t0_lab = -zp[i]/vz_lab; + const Real t0_lab = -zp/vz_lab; if (!projected) { - xp[i] += t0_lab*vx_lab; - yp[i] += t0_lab*vy_lab; + xp += t0_lab*vx_lab; + yp += t0_lab*vy_lab; } if (focused) { // Correct for focusing effect from shift from z=0 to zinject const Real tfocus = -zinject_plane*WarpX::gamma_boost/vz_lab; - xp[i] -= tfocus*vx_lab; - yp[i] -= tfocus*vy_lab; + xp -= tfocus*vx_lab; + yp -= tfocus*vy_lab; } // Time of the particle in the boosted frame given its position in the lab frame at t=0. - const Real tpr = -WarpX::gamma_boost*WarpX::beta_boost*zp[i]/PhysConst::c; + const Real tpr = -WarpX::gamma_boost*WarpX::beta_boost*zp/PhysConst::c; // Position of the particle in the boosted frame given its position in the lab frame at t=0. - const Real zpr = WarpX::gamma_boost*zp[i]; + const Real zpr = WarpX::gamma_boost*zp; // Momentum of the particle in the boosted frame (assuming that it is fixed). uzp[i] = WarpX::gamma_boost*(uzp[i] - WarpX::beta_boost*PhysConst::c*gamma_lab); @@ -198,30 +199,23 @@ RigidInjectedParticleContainer::BoostandRemapParticles() // Put the particle at the location in the boosted frame at boost frame t=0, if (rigid_advance) { // with the particle moving at the average velocity - zp[i] = zpr - vzbeam_ave_boosted*tpr; + zp = zpr - vzbeam_ave_boosted*tpr; } else { // with the particle moving with its own velocity const Real gammapr = std::sqrt(1. + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])/(PhysConst::c*PhysConst::c)); const Real vzpr = uzp[i]/gammapr; - zp[i] = zpr - vzpr*tpr; + zp = zpr - vzpr*tpr; } + SetPosition(i, xp, yp, zp); } - - // Copy the data back to the particle container - pti.SetPosition(xp, yp, zp); - } } } void -RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, - Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp, - Real dt, DtType a_dt_type) +RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. @@ -234,9 +228,9 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Gpu::ManagedDeviceVector<ParticleReal> xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + ParticleReal* const AMREX_RESTRICT ux = uxp.dataPtr(); ParticleReal* const AMREX_RESTRICT uy = uyp.dataPtr(); ParticleReal* const AMREX_RESTRICT uz = uzp.dataPtr(); @@ -247,14 +241,38 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, ParticleReal* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr(); ParticleReal* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); - if (!done_injecting_lev) { + if (!done_injecting_lev) + { // If the old values are not already saved, create copies here. - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; + const auto np = pti.numParticles(); + + xp_save.resize(np); + yp_save.resize(np); + zp_save.resize(np); + + uxp_save.resize(np); + uyp_save.resize(np); + uzp_save.resize(np); + + amrex::Real* const AMREX_RESTRICT xp_save_ptr = xp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT yp_save_ptr = yp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT zp_save_ptr = zp_save.dataPtr(); + + amrex::Real* const AMREX_RESTRICT uxp_save_ptr = uxp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT uyp_save_ptr = uyp_save.dataPtr(); + amrex::Real* const AMREX_RESTRICT uzp_save_ptr = uzp_save.dataPtr(); + + amrex::ParallelFor( np, + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + xp_save_ptr[i] = xp; + yp_save_ptr[i] = yp; + zp_save_ptr[i] = zp; + uxp_save_ptr[i] = ux[i]; + uyp_save_ptr[i] = uy[i]; + uzp_save_ptr[i] = uz[i]; + }); // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -265,20 +283,22 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, const Real vz_ave_boosted = vzbeam_ave_boosted; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const Real dtscale = dt - (z_plane_previous - z[i])/(vz_ave_boosted + v_boost); - if (0. < dtscale && dtscale < dt) { - Exp[i] *= dtscale; - Eyp[i] *= dtscale; - Ezp[i] *= dtscale; - Bxp[i] *= dtscale; - Byp[i] *= dtscale; - Bzp[i] *= dtscale; - } - } - ); + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + const Real dtscale = dt - (z_plane_previous - zp)/(vz_ave_boosted + v_boost); + if (0. < dtscale && dtscale < dt) { + Exp[i] *= dtscale; + Eyp[i] *= dtscale; + Ezp[i] *= dtscale; + Bxp[i] *= dtscale; + Byp[i] *= dtscale; + Bzp[i] *= dtscale; + } + } + ); } - PhysicalParticleContainer::PushPX(pti, xp, yp, zp, dt, a_dt_type); + PhysicalParticleContainer::PushPX(pti, dt, a_dt_type); if (!done_injecting_lev) { @@ -296,23 +316,25 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, const bool rigid = rigid_advance; const Real inv_csq = 1./(PhysConst::c*PhysConst::c); amrex::ParallelFor( pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - if (z[i] <= z_plane_lev) { - ux[i] = ux_save[i]; - uy[i] = uy_save[i]; - uz[i] = uz_save[i]; - x[i] = x_save[i]; - y[i] = y_save[i]; - if (rigid) { - z[i] = z_save[i] + dt*vz_ave_boosted; - } - else { - const Real gi = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_csq); - z[i] = z_save[i] + dt*uz[i]*gi; - } - } - } - ); + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + if (zp <= z_plane_lev) { + ux[i] = ux_save[i]; + uy[i] = uy_save[i]; + uz[i] = uz_save[i]; + xp = x_save[i]; + yp = y_save[i]; + if (rigid) { + zp = z_save[i] + dt*vz_ave_boosted; + } + else { + const Real gi = 1./std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_csq); + zp = z_save[i] + dt*uz[i]*gi; + } + SetPosition(i, xp, yp, zp); + } + }); } } @@ -370,11 +392,6 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, #pragma omp parallel #endif { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -401,16 +418,11 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - 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); + 0, np, lev, lev); // Save the position and momenta, making copies auto uxp_save = uxp; @@ -419,7 +431,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - const ParticleReal* const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + const auto GetPosition = GetParticlePosition(pti); ParticleReal* const AMREX_RESTRICT uxpp = uxp.dataPtr(); ParticleReal* const AMREX_RESTRICT uypp = uyp.dataPtr(); ParticleReal* const AMREX_RESTRICT uzpp = uzp.dataPtr(); @@ -485,15 +497,16 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const ParticleReal* const AMREX_RESTRICT uz_save = uzp_save.dataPtr(); const ParticleReal zz = zinject_plane_levels[lev]; amrex::ParallelFor( pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - if (zp[i] <= zz) { - uxpp[i] = ux_save[i]; - uypp[i] = uy_save[i]; - uzpp[i] = uz_save[i]; - } - } - ); - + [=] AMREX_GPU_DEVICE (long i) { + ParticleReal xp, yp, zp; + GetPosition(i, xp, yp, zp); + if (zp <= zz) { + uxpp[i] = ux_save[i]; + uypp[i] = uy_save[i]; + uzpp[i] = uz_save[i]; + } + } + ); } } } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index cac1d9883..b1def39d3 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -81,14 +81,6 @@ public: WarpXParIter (ContainerType& pc, int level); -#if (AMREX_SPACEDIM == 2) - void GetPosition (amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& x, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& y, - amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& z) const; - void SetPosition (const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& x, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& y, - const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& z); -#endif const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { return GetStructOfArrays().GetRealData(); } @@ -379,8 +371,6 @@ protected: using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>; using PairIndex = std::pair<int, int>; - amrex::Vector<DataContainer> m_xp, m_yp, m_zp; - // Whether to dump particle quantities. // If true, particle position is always dumped. int plot_species = 1; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 85cac0dbe..cd5bcaaf9 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -30,56 +30,6 @@ WarpXParIter::WarpXParIter (ContainerType& pc, int level) { } -#if (AMREX_SPACEDIM == 2) -void -WarpXParIter::GetPosition (Gpu::ManagedDeviceVector<ParticleReal>& xp, - Gpu::ManagedDeviceVector<ParticleReal>& yp, - Gpu::ManagedDeviceVector<ParticleReal>& zp) const -{ - amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(xp, zp); -#ifdef WARPX_DIM_RZ - const auto& attribs = GetAttribs(); - const auto& thetap = attribs[PIdx::theta]; - yp.resize(xp.size()); - ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT theta = thetap.dataPtr(); - amrex::ParallelFor( xp.size(), - [=] AMREX_GPU_DEVICE (long i) { - // The x stored in the particles is actually the radius - y[i] = x[i]*std::sin(theta[i]); - x[i] = x[i]*std::cos(theta[i]); - }); -#else - yp.resize(xp.size(), std::numeric_limits<ParticleReal>::quiet_NaN()); -#endif -} - -void -WarpXParIter::SetPosition (const Gpu::ManagedDeviceVector<ParticleReal>& xp, - const Gpu::ManagedDeviceVector<ParticleReal>& yp, - const Gpu::ManagedDeviceVector<ParticleReal>& zp) -{ -#ifdef WARPX_DIM_RZ - auto& attribs = GetAttribs(); - auto& thetap = attribs[PIdx::theta]; - Gpu::ManagedDeviceVector<ParticleReal> rp(xp.size()); - const ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); - ParticleReal* const AMREX_RESTRICT r = rp.dataPtr(); - ParticleReal* const AMREX_RESTRICT theta = thetap.dataPtr(); - amrex::ParallelFor( xp.size(), - [=] AMREX_GPU_DEVICE (long i) { - theta[i] = std::atan2(y[i], x[i]); - r[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); - }); - amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(rp, zp); -#else - amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(xp, zp); -#endif -} -#endif - WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) : ParticleContainer<0,0,PIdx::nattribs>(amr_core->GetParGDB()) , species_id(ispecies) @@ -116,10 +66,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) local_jx.resize(num_threads); local_jy.resize(num_threads); local_jz.resize(num_threads); - m_xp.resize(num_threads); - m_yp.resize(num_threads); - m_zp.resize(num_threads); - } void @@ -349,9 +295,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // CPU, tiling: deposit into local_jx // (same for jx and jz) - ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const auto GetPosition = GetParticlePosition(pti, offset); // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -362,19 +306,19 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doEsirkepovDepositionShapeN<2>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doEsirkepovDepositionShapeN<3>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); @@ -382,19 +326,19 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } else { if (WarpX::nox == 1){ doDepositionShapeN<1>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ doDepositionShapeN<2>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ doDepositionShapeN<3>( - xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q); @@ -488,9 +432,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // GPU, no tiling: deposit directly in rho // CPU, tiling: deposit into local_rho - ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const auto GetPosition = GetParticlePosition(pti, offset); // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -500,13 +442,13 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<1>(GetPosition, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<2>(GetPosition, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + doChargeDepositionShapeN<3>(GetPosition, wp.dataPtr()+offset, ion_lev, rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -545,8 +487,6 @@ WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::Mult const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -615,8 +555,6 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -772,32 +710,27 @@ WarpXParticleContainer::PushX (int lev, amrex::Real dt) // // Particle Push // - // Extract pointers to particle position and momenta, for this particle tile - // - positions are stored as an array of struct, in `ParticleType` - ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); + + const auto GetPosition = GetParticlePosition(pti); + auto SetPosition = SetParticlePosition(pti); + // - momenta are stored as a struct of array, in `attribs` auto& attribs = pti.GetAttribs(); ParticleReal* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); ParticleReal* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); #ifdef WARPX_DIM_RZ + auto& aos = pti.GetArrayOfStructs(); + ParticleType* AMREX_RESTRICT const pstruct = aos().dataPtr(); ParticleReal* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); #endif // Loop over the particles and update their position amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - ParticleType& p = pstructs[i]; // Particle object that gets updated - ParticleReal x, y, z; // Temporary variables -#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_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 ); -#endif + ParticleReal x, y, z; + GetPosition(i, x, y, z); + UpdatePosition(x, y, z, ux[i], uy[i], uz[i], dt); + SetPosition(i, x, y, z); } ); diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index 34e08118d..18a77f6ca 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -9,6 +9,7 @@ #define WARPX_CONST_H_ #include <AMReX_REAL.H> +#include <cmath> // Math constants namespace MathConst @@ -29,6 +30,8 @@ namespace PhysConst static constexpr amrex::Real hbar = 1.054571817e-34; static constexpr amrex::Real alpha = mu0/(4*MathConst::pi)*q_e*q_e*c/hbar; static constexpr amrex::Real r_e = 1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c); + static constexpr amrex::Real xi = (2.*alpha*alpha*ep0*ep0*hbar*hbar*hbar)/ + (45.*m_e*m_e*m_e*m_e*c*c*c*c*c); // SI value is 1.3050122.e-52 } #endif diff --git a/Source/WarpX.H b/Source/WarpX.H index ae25a8168..81893175f 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -236,6 +236,26 @@ public: void EvolveE (int lev, PatchType patch_type, amrex::Real dt); void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type); + /** \brief apply QED correction on electric field + * \param dt vector of time steps (for all levels) + */ + void Hybrid_QED_Push ( amrex::Vector<amrex::Real> dt); + + /** \brief apply QED correction on electric field for level lev + * \param lev mesh refinement level + * \param dt time step + */ + void Hybrid_QED_Push (int lev, amrex::Real dt); + + /** \brief apply QED correction on electric field for level lev and patch type patch_type + * \param lev mesh refinement level + * \param dt patch_type which MR patch: PatchType::fine or PatchType::coarse + * \param dt time step + */ + void Hybrid_QED_Push (int lev, PatchType patch_type, amrex::Real dt); + + static amrex::Real quantum_xi; + #ifdef WARPX_DIM_RZ void ApplyInverseVolumeScalingToCurrentDensity(amrex::MultiFab* Jx, amrex::MultiFab* Jy, @@ -634,6 +654,8 @@ private: // Other runtime parameters int verbose = 1; + bool use_hybrid_QED = 0; + int max_step = std::numeric_limits<int>::max(); amrex::Real stop_time = std::numeric_limits<amrex::Real>::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d93fab7df..3fe61274a 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -55,6 +55,7 @@ int WarpX::do_moving_window = 0; int WarpX::moving_window_dir = -1; Real WarpX::moving_window_v = std::numeric_limits<amrex::Real>::max(); +Real WarpX::quantum_xi = PhysConst::xi; Real WarpX::gamma_boost = 1.; Real WarpX::beta_boost = 0.; Vector<int> WarpX::boost_direction = {0,0,0}; @@ -341,6 +342,7 @@ WarpX::ReadParameters () pp.query("verbose", verbose); pp.query("regrid_int", regrid_int); pp.query("do_subcycling", do_subcycling); + pp.query("use_hybrid_QED", use_hybrid_QED); pp.query("exchange_all_guard_cells", exchange_all_guard_cells); pp.query("override_sync_int", override_sync_int); @@ -456,6 +458,8 @@ WarpX::ReadParameters () pp.query("n_current_deposition_buffer", n_current_deposition_buffer); pp.query("sort_int", sort_int); + pp.query("quantum_xi", quantum_xi); + pp.query("do_pml", do_pml); pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); |