diff options
| -rwxr-xr-x | Examples/Tests/PML/analysis_pml_ckc.py | 5 | ||||
| -rw-r--r-- | Source/Diagnostics/WarpXIO.cpp | 1 | ||||
| -rw-r--r-- | Source/Particles/Deposition/CurrentDeposition.H | 179 | ||||
| -rw-r--r-- | Source/Particles/Deposition/Make.package | 3 | ||||
| -rw-r--r-- | Source/Particles/Make.package | 1 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 33 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 19 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 144 | ||||
| -rw-r--r-- | Source/WarpX.H | 1 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 11 |
10 files changed, 374 insertions, 23 deletions
diff --git a/Examples/Tests/PML/analysis_pml_ckc.py b/Examples/Tests/PML/analysis_pml_ckc.py index 90d13cecc..e095864c0 100755 --- a/Examples/Tests/PML/analysis_pml_ckc.py +++ b/Examples/Tests/PML/analysis_pml_ckc.py @@ -30,5 +30,8 @@ energy_end = energyE + energyB Reflectivity = energy_end/energy_start Reflectivity_theory = 1.8015e-06 -assert( abs(Reflectivity-Reflectivity_theory) < 5./100 * Reflectivity_theory ) +print("Reflectivity", Reflectivity) +print("Reflectivity_theory", Reflectivity_theory) + +assert( Reflectivity < 105./100 * Reflectivity_theory ) diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 24272c588..38399bf9e 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -755,7 +755,6 @@ WarpX::WriteSlicePlotFile () const // writing plotfile // const std::string& slice_plotfilename = amrex::Concatenate(slice_plot_file,istep[0]); amrex::Print() << " Writing slice plotfile " << slice_plotfilename << "\n"; - const int ngrow = 0; Vector<std::string> rfs; VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H new file mode 100644 index 000000000..efda97892 --- /dev/null +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -0,0 +1,179 @@ +#ifndef CURRENTDEPOSITION_H_ +#define CURRENTDEPOSITION_H_ + +using namespace amrex; + +// Compute shape factor and return index of leftmost cell where +// particle writes. +// Specialized templates are defined below for orders 1, 2 and 3. +template <int depos_order> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int compute_shape_factor(Real* const sx, Real xint) {return 0;}; + +// Compute shape factor for order 1. +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int compute_shape_factor <1> (Real* const sx, Real xmid){ + int j = (int) xmid; + Real xint = xmid-j; + sx[0] = 1.0 - xint; + sx[1] = xint; + return j; +} + +// Compute shape factor for order 2. +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int compute_shape_factor <2> (Real* const sx, Real xmid){ + int j = (int) (xmid+0.5); + Real xint = xmid-j; + sx[0] = 0.5*(0.5-xint)*(0.5-xint); + sx[1] = 0.75-xint*xint; + sx[2] = 0.5*(0.5+xint)*(0.5+xint); + // index of the leftmost cell where particle deposits + return j-1; +} + +// Compute shape factor for order 3. +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int compute_shape_factor <3> (Real* const sx, Real xmid){ + int j = (int) xmid; + Real xint = xmid-j; + sx[0] = 1.0/6.0*(1.0-xint)*(1.0-xint)*(1.0-xint); + sx[1] = 2.0/3.0-xint*xint*(1-xint/2.0); + sx[2] = 2.0/3.0-(1-xint)*(1-xint)*(1.0-0.5*(1-xint)); + sx[3] = 1.0/6.0*xint*xint*xint; + // index of the leftmost cell where particle deposits + return j-1; +} + +/* \brief Current Deposition for thread thread_num + * /param xp, yp, zp : Pointer to arrays of particle positions. + * \param wp : Pointer to array of particle weights. + * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param jx_arr : Array4 of current density, either full array or tile. + * \param jy_arr : Array4 of current density, either full array or tile. + * \param jz_arr : Array4 of current density, either full array or tile. + * \param offset : Index of first particle for which current is deposited + * \param np_to_depose : Number of particles for which current is deposited. + Particles [offset,offset+np_tp_depose] deposit current. + * \param dt : Time step for particle level + * \param dx : 3D cell size + * \param xyzmin : Physical lower bounds of domain. + * \param lo : Index lower bounds of domain. + * \param stagger_shift: 0 if nodal, 0.5 if staggered. + * /param q : species charge. + */ +template <int depos_order> +void doDepositionShapeN(const Real * const xp, const Real * const yp, const Real * const zp, + const Real * const wp, const Real * const uxp, + const Real * const uyp, const Real * const uzp, + const amrex::Array4<amrex::Real>& jx_arr, + const amrex::Array4<amrex::Real>& jy_arr, + const amrex::Array4<amrex::Real>& jz_arr, + const long offset, const long np_to_depose, + const amrex::Real dt, const std::array<amrex::Real,3>& dx, + const std::array<Real, 3> xyzmin, + const Dim3 lo, + const amrex::Real stagger_shift, + const amrex::Real q) +{ + const Real dxi = 1.0/dx[0]; + const Real dzi = 1.0/dx[2]; + const Real dts2dx = 0.5*dt*dxi; + const Real dts2dz = 0.5*dt*dzi; +#if (AMREX_SPACEDIM == 2) + const Real invvol = dxi*dzi; +#else // (AMREX_SPACEDIM == 3) + const Real dyi = 1.0/dx[1]; + const Real dts2dy = 0.5*dt*dyi; + const Real invvol = dxi*dyi*dzi; +#endif + + const Real xmin = xyzmin[0]; + const Real ymin = xyzmin[1]; + const Real zmin = xyzmin[2]; + const Real clightsq = 1.0/PhysConst::c/PhysConst::c; + + // Loop over particles and deposit into jx_arr, jy_arr and jz_arr + ParallelFor( np_to_depose, + [=] AMREX_GPU_DEVICE (long ip) { + // --- Get particle quantities + const Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + + uyp[ip]*uyp[ip]*clightsq + + uzp[ip]*uzp[ip]*clightsq); + const Real wq = q*wp[ip]; + const Real vx = uxp[ip]*gaminv; + const Real vy = uyp[ip]*gaminv; + const Real vz = uzp[ip]*gaminv; + // wqx, wqy wqz are particle current in each direction + const Real wqx = wq*invvol*vx; + const Real wqy = wq*invvol*vy; + const Real wqz = wq*invvol*vz; + + // --- Compute shape factors + // x direction + // Get particle position after 1/2 push back in position + const Real xmid = (xp[ip]-xmin)*dxi-dts2dx*vx; + // Compute shape factors for node-centered quantities + Real AMREX_RESTRICT sx [depos_order + 1]; + // j: leftmost grid point (node-centered) that the particle touches + const int j = compute_shape_factor<depos_order>(sx, xmid); + // Compute shape factors for cell-centered quantities + Real AMREX_RESTRICT sx0[depos_order + 1]; + // j0: leftmost grid point (cell-centered) that the particle touches + const int j0 = compute_shape_factor<depos_order>(sx0, xmid-stagger_shift); + +#if (AMREX_SPACEDIM == 3) + // y direction + const Real ymid= (yp[ip]-ymin)*dyi-dts2dy*vy; + Real AMREX_RESTRICT sy [depos_order + 1]; + const int k = compute_shape_factor<depos_order>(sy, ymid); + Real AMREX_RESTRICT sy0[depos_order + 1]; + const int k0 = compute_shape_factor<depos_order>(sy0, ymid-stagger_shift); +#endif + // z direction + const Real zmid= (zp[ip]-zmin)*dzi-dts2dz*vz; + Real AMREX_RESTRICT sz [depos_order + 1]; + const int l = compute_shape_factor<depos_order>(sz, zmid); + Real AMREX_RESTRICT sz0[depos_order + 1]; + const int l0 = compute_shape_factor<depos_order>(sz0, zmid-stagger_shift); + + // Deposit current into jx_arr, jy_arr and jz_arr +#if (AMREX_SPACEDIM == 2) + for (int iz=0; iz<=depos_order; iz++){ + for (int ix=0; ix<=depos_order; ix++){ + amrex::Gpu::Atomic::Add( + &jx_arr(lo.x+j0+ix, lo.y+l +iz, 0), + sx0[ix]*sz [iz]*wqx); + amrex::Gpu::Atomic::Add( + &jy_arr(lo.x+j +ix, lo.y+l +iz, 0), + sx [ix]*sz [iz]*wqy); + amrex::Gpu::Atomic::Add( + &jz_arr(lo.x+j +ix, lo.y+l0+iz, 0), + sx [ix]*sz0[iz]*wqz); + } + } +#else // (AMREX_SPACEDIM == 3) + for (int iz=0; iz<=depos_order; iz++){ + for (int iy=0; iy<=depos_order; iy++){ + for (int ix=0; ix<=depos_order; ix++){ + amrex::Gpu::Atomic::Add( + &jx_arr(lo.x+j0+ix, lo.y+k +iy, lo.z+l +iz), + sx0[ix]*sy [iy]*sz [iz]*wqx); + amrex::Gpu::Atomic::Add( + &jy_arr(lo.x+j +ix, lo.y+k0+iy, lo.z+l +iz), + sx [ix]*sy0[iy]*sz [iz]*wqy); + amrex::Gpu::Atomic::Add( + &jz_arr(lo.x+j +ix, lo.y+k +iy, lo.z+l0+iz), + sx [ix]*sy [iy]*sz0[iz]*wqz); + } + } + } +#endif + } + ); +} + +#endif // CURRENTDEPOSITION_H_ diff --git a/Source/Particles/Deposition/Make.package b/Source/Particles/Deposition/Make.package new file mode 100644 index 000000000..0d5ebe2a7 --- /dev/null +++ b/Source/Particles/Deposition/Make.package @@ -0,0 +1,3 @@ +CEXE_headers += CurrentDeposition.H +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index f33095738..2038472a1 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -11,6 +11,7 @@ CEXE_headers += RigidInjectedParticleContainer.H CEXE_headers += PhysicalParticleContainer.H include $(WARPX_HOME)/Source/Particles/Pusher/Make.package +include $(WARPX_HOME)/Source/Particles/Deposition/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d6c4973c3..28d611020 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1529,17 +1529,30 @@ PhysicalParticleContainer::Evolve (int lev, // // Current Deposition // - // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, - 0, np_current, thread_num, - lev, lev, dt); - if (has_buffer){ - // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, - np_current, np-np_current, thread_num, - lev, lev-1, dt); + if (WarpX::use_picsar_deposition) { + // Deposit inside domains + DepositCurrentFortran(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt); + if (has_buffer){ + // Deposit in buffers + DepositCurrentFortran(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt); + } + } else { + // Deposit inside domains + DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt); + if (has_buffer){ + // Deposit in buffers + DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt); + } } - + // // copy particle data back // diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 0e800bf1d..662b2e1b8 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -189,6 +189,21 @@ public: int depos_lev, amrex::Real dt); + virtual void DepositCurrentFortran(WarpXParIter& pti, + RealVector& wp, + RealVector& uxp, + RealVector& uyp, + RealVector& uzp, + amrex::MultiFab* jx, + amrex::MultiFab* jy, + amrex::MultiFab* jz, + const long offset, + const long np_to_depose, + int thread_num, + int lev, + int depos_lev, + amrex::Real dt); + // If particles start outside of the domain, ContinuousInjection // makes sure that they are initialized when they enter the domain, and // NOT before. Virtual function, overriden by derived classes. @@ -253,8 +268,8 @@ public: AddIntComp(comm); } - int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } - + int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + protected: std::map<std::string, int> particle_comps; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 317f46fd4..7316dcc95 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -10,6 +10,7 @@ // Import low-level single-particle kernels #include <GetAndSetPosition.H> #include <UpdatePosition.H> +#include <CurrentDeposition.H> using namespace amrex; @@ -279,7 +280,7 @@ WarpXParticleContainer::AddNParticles (int lev, Redistribute(); } -/* \brief Current Deposition for thread thread_num +/* \brief Current Deposition for thread thread_num using PICSAR * \param pti : Particle iterator * \param wp : Array of particle weights * \param uxp uyp uzp : Array of particle @@ -292,16 +293,15 @@ WarpXParticleContainer::AddNParticles (int lev, * \param lev : Level of box that contains particles * \param depos_lev : Level on which particles deposit (if buffers are used) * \param dt : Time step for particle level - * \param ngJ : Number of ghosts cells (of lev) */ void -WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, - RealVector& wp, RealVector& uxp, - RealVector& uyp, RealVector& uzp, - MultiFab* jx, MultiFab* jy, MultiFab* jz, - const long offset, const long np_to_depose, - int thread_num, int lev, int depos_lev, - Real dt) +WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, + RealVector& wp, RealVector& uxp, + RealVector& uyp, RealVector& uzp, + MultiFab* jx, MultiFab* jy, MultiFab* jz, + const long offset, const long np_to_depose, + int thread_num, int lev, int depos_lev, + Real dt) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || (depos_lev==(lev )), @@ -414,6 +414,130 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #endif } +/* \brief Current Deposition for thread thread_num + * \param pti : Particle iterator + * \param wp : Array of particle weights + * \param uxp uyp uzp : Array of particle + * \param jx jy jz : Full array of current density + * \param offset : Index of first particle for which current is deposited + * \param np_to_depose: Number of particles for which current is deposited. + Particles [offset,offset+np_tp_depose] deposit current + * \param thread_num : Thread number (if tiling) + * \param lev : Level of box that contains particles + * \param depos_lev : Level on which particles deposit (if buffers are used) + * \param dt : Time step for particle level + */ +void +WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, + RealVector& wp, RealVector& uxp, + RealVector& uyp, RealVector& uzp, + MultiFab* jx, MultiFab* jy, MultiFab* jz, + const long offset, const long np_to_depose, + int thread_num, int lev, int depos_lev, + Real dt) +{ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || + (depos_lev==(lev )), + "Deposition buffers only work for lev-1"); + // If no particles, do not do anything + if (np_to_depose == 0) return; + + const long ngJ = jx->nGrow(); + const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0)); + int j_is_nodal = jx->is_nodal() and jy->is_nodal() and jz->is_nodal(); + Real q = this->charge; + const Real stagger_shift = j_is_nodal ? 0.0 : 0.5; + + BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); + + // Get tile box where current is deposited. + // The tile box is different when depositing in the buffers (depos_lev<lev) + // or when depositing inside the level (depos_lev=lev) + Box tilebox; + if (lev == depos_lev) { + tilebox = pti.tilebox(); + } else { + const IntVect& ref_ratio = WarpX::RefRatio(depos_lev); + tilebox = amrex::coarsen(pti.tilebox(),ref_ratio); + } + + // Staggered tile boxes (different in each direction) + Box tbx = convert(tilebox, WarpX::jx_nodal_flag); + Box tby = convert(tilebox, WarpX::jy_nodal_flag); + Box tbz = convert(tilebox, WarpX::jz_nodal_flag); + tilebox.grow(ngJ); + +#ifdef AMREX_USE_GPU + // No tiling on GPU: jx_ptr points to the full + // jx array (same for jy_ptr and jz_ptr). + Array4<Real> const& jx_arr = jx->array(pti); + Array4<Real> const& jy_arr = jy->array(pti); + Array4<Real> const& jz_arr = jz->array(pti); +#else + // Tiling is on: jx_ptr points to local_jx[thread_num] + // (same for jy_ptr and jz_ptr) + tbx.grow(ngJ); + tby.grow(ngJ); + tbz.grow(ngJ); + + local_jx[thread_num].resize(tbx); + local_jy[thread_num].resize(tby); + local_jz[thread_num].resize(tbz); + + // local_jx[thread_num] is set to zero + local_jx[thread_num].setVal(0.0); + local_jy[thread_num].setVal(0.0); + local_jz[thread_num].setVal(0.0); + + Array4<Real> const& jx_arr = local_jx[thread_num].array(); + Array4<Real> const& jy_arr = local_jy[thread_num].array(); + Array4<Real> const& jz_arr = local_jz[thread_num].array(); +#endif + // GPU, no tiling: deposit directly in jx + // CPU, tiling: deposit into local_jx + // (same for jx and jz) + + Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + + // Lower corner of tile box physical domain + const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev);; + // xyzmin is built on pti.tilebox(), so it does + // not include staggering, so the stagger_shift has to be done by hand. + // Alternatively, we could define xyzminx from tbx (and the same for 3 + // directions and for jx, jy, jz). This way, sx0 would not be needed. + // Better for memory? worth trying? + const Dim3 lo = lbound(tilebox); + + if (WarpX::nox == 1){ + doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), + uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + jz_arr, offset, np_to_depose, dt, dx, + xyzmin, lo, stagger_shift, q); + } else if (WarpX::nox == 2){ + doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), + uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + jz_arr, offset, np_to_depose, dt, dx, + xyzmin, lo, stagger_shift, q); + } else if (WarpX::nox == 3){ + doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), + uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + jz_arr, offset, np_to_depose, dt, dx, + xyzmin, lo, stagger_shift, q); + } + +#ifndef AMREX_USE_GPU + BL_PROFILE_VAR_START(blp_accumulate); + // CPU, tiling: atomicAdd local_jx into jx + // (same for jx and jz) + (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); + (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1); + (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1); + BL_PROFILE_VAR_STOP(blp_accumulate); +#endif +} + void WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, MultiFab* rhomf, MultiFab* crhomf, int icomp, @@ -932,3 +1056,5 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, if (pld.m_lev == lev-1){ } } + + diff --git a/Source/WarpX.H b/Source/WarpX.H index 4ad3d119f..7d34178f5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -81,6 +81,7 @@ public: static amrex::Vector<amrex::Real> B_external; // Algorithms + static long use_picsar_deposition; static long current_deposition_algo; static long charge_deposition_algo; static long field_gathering_algo; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 877882037..b20c6211c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -37,6 +37,7 @@ Vector<int> WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; +long WarpX::use_picsar_deposition = 1; long WarpX::current_deposition_algo; long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; @@ -484,7 +485,17 @@ WarpX::ReadParameters () { ParmParse pp("algo"); + // If not in RZ mode, read use_picsar_deposition + // In RZ mode, use_picsar_deposition is on, as the C++ version + // of the deposition does not support RZ +#ifndef WARPX_RZ + pp.query("use_picsar_deposition", use_picsar_deposition); +#endif current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); + if (!use_picsar_deposition){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo >= 2, + "if not use_picsar_deposition, cannot use Esirkepov deposition."); + } charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); particle_pusher_algo = GetAlgorithmInteger(pp, "particle_pusher"); |