From 4ee2de44648757a767f48d76393b816d9b7f4715 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Tue, 3 Sep 2019 16:25:25 -0700 Subject: In current deposition, allowed more flexible centering --- Source/Particles/WarpXParticleContainer.cpp | 30 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 16 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 80d1caa0f..29b86be49 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -298,9 +298,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, const long ngJ = jx->nGrow(); const std::array& 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); BL_PROFILE_VAR_NS("PPC::CurrentDeposition", blp_deposit); @@ -357,18 +355,13 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, 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 - // Note that this includes guard cells since it is after tilebox.ngrow - const std::array& 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); BL_PROFILE_VAR_START(blp_deposit); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { + // Lower corner of tile box physical domain + // Note that this includes guard cells since it is after tilebox.ngrow + const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>( xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, @@ -386,24 +379,29 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } } else { + // Lower corner of tile box physical domain + // Note that this includes guard cells since it is after tilebox.ngrow + std::array xyzminx = WarpX::LowerCornerWithCentering(tbx, depos_lev); + std::array xyzminy = WarpX::LowerCornerWithCentering(tby, depos_lev); + std::array xyzminz = WarpX::LowerCornerWithCentering(tbz, depos_lev); if (WarpX::nox == 1){ doDepositionShapeN<1>( xp, yp, zp, 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, - stagger_shift, q); + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, + xyzminx, xyzminy, xyzminz, lo, q); } else if (WarpX::nox == 2){ doDepositionShapeN<2>( xp, yp, zp, 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, - stagger_shift, q); + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, + xyzminx, xyzminy, xyzminz, lo, q); } else if (WarpX::nox == 3){ doDepositionShapeN<3>( xp, yp, zp, 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, - stagger_shift, q); + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, + xyzminx, xyzminy, xyzminz, lo, q); } } BL_PROFILE_VAR_STOP(blp_deposit); -- cgit v1.2.3 From 9d567430aeb46181f81e3e1f1dbc48b6291b1e95 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 30 Sep 2019 11:41:23 -0700 Subject: Generalized direct current deposition allowing different centerings for Jx, y, and z --- Source/Particles/Deposition/CurrentDeposition.H | 129 +++++++++++++++--------- Source/Particles/WarpXParticleContainer.cpp | 11 +- 2 files changed, 82 insertions(+), 58 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index a047a2ed2..548cadeb6 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -18,9 +18,10 @@ * \param np_to_depose : Number of particles for which current is deposited. * \param dt : Time step for particle level * \param dx : 3D cell size - * \param xyzminx : Physical lower bounds of domain for jx. - * \param xyzminy : Physical lower bounds of domain for jy. - * \param xyzminz : Physical lower bounds of domain for jz. + * \param tbx : Box of domain for jx. + * \param tby : Box of domain for jy. + * \param tbz : Box of domain for jz. + * \param depos_lev : Refinement level being used for the deposition. * \param lo : Index lower bounds of domain. * /param q : species charge. */ @@ -38,9 +39,10 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, const amrex::Array4& jz_arr, const long np_to_depose, const amrex::Real dt, const std::array& dx, - const std::array& xyzminx, - const std::array& xyzminy, - const std::array& xyzminz, + const amrex::Box& tbx, + const amrex::Box& tby, + const amrex::Box& tbz, + const int depos_lev, const amrex::Dim3 lo, const amrex::Real q) { @@ -59,16 +61,20 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, const amrex::Real invvol = dxi*dyi*dzi; #endif - // This assumes that the centering in the directions transverse - // to the direction of the current are the same. - const amrex::Real xmin_x = xyzminx[0]; - const amrex::Real xmin_yz = xyzminz[0]; + const std::array xyzmin_jx = WarpX::LowerCornerWithCentering(tbx, depos_lev); + const std::array xyzmin_jy = WarpX::LowerCornerWithCentering(tby, depos_lev); + const std::array xyzmin_jz = WarpX::LowerCornerWithCentering(tbz, depos_lev); + const amrex::Real xmin_jx = xyzmin_jx[0]; + const amrex::Real xmin_jy = xyzmin_jy[0]; + const amrex::Real xmin_jz = xyzmin_jz[0]; #if (defined WARPX_DIM_3D) - const amrex::Real ymin_y = xyzminy[1]; - const amrex::Real ymin_xz = xyzminx[1]; + const amrex::Real ymin_jx = xyzmin_jx[1]; + const amrex::Real ymin_jy = xyzmin_jy[1]; + const amrex::Real ymin_jz = xyzmin_jz[1]; #endif - const amrex::Real zmin_z = xyzminz[2]; - const amrex::Real zmin_xy = xyzminx[2]; + const amrex::Real zmin_jx = xyzmin_jx[2]; + const amrex::Real zmin_jy = xyzmin_jy[2]; + const amrex::Real zmin_jz = xyzmin_jz[2]; const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c; @@ -115,51 +121,74 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, // x direction // Get particle position after 1/2 push back in position #if (defined WARPX_DIM_RZ) - const amrex::Real xmid_yz = (rpmid - xmin_yz)*dxi; - const amrex::Real xmid_x = (rpmid - xmin_x )*dxi; + const amrex::Real xmid_jx = (rpmid - xmin_jx)*dxi; + const amrex::Real xmid_jz = (rpmid - xmin_jz)*dxi; #else - const amrex::Real xmid_yz = (xp[ip] - xmin_yz)*dxi - dts2dx*vx; - const amrex::Real xmid_x = (xp[ip] - xmin_x )*dxi - dts2dx*vx; + const amrex::Real xmid_jx = (xp[ip] - xmin_jx)*dxi - dts2dx*vx; + const amrex::Real xmid_jy = (xp[ip] - xmin_jy)*dxi - dts2dx*vx; + const amrex::Real xmid_jz = (xp[ip] - xmin_jz)*dxi - dts2dx*vx; #endif - // Compute shape factors for jy, jz, perpendicular to the x-axis - amrex::Real sx_yz[depos_order + 1]; - // j_yz: leftmost grid point (node-centered) that the particle touches - const int j_yz = compute_shape_factor(sx_yz, xmid_yz); - // Compute shape factors for jx, parallel to the x-axis - amrex::Real sx_x[depos_order + 1]; - // j_x: leftmost grid point (cell-centered) that the particle touches - const int j_x = compute_shape_factor(sx_x, xmid_x); + // j_[xyz]: leftmost grid point in x that the particle touches for the centering of each current + // sx_[xyz] shape factor along x for the centering of each current + // There are only two possible centerings, node or cell centered, so at most only two shape factor + // arrays will be needed. + // sx_jx is always calculated. sx_2 only if needed (when Jy or Jz have different centerings than Jx). + amrex::Real sx_jx[depos_order + 1]; + amrex::Real sx_2[depos_order + 1]; + const int j_jx = compute_shape_factor(sx_jx, xmid_jx); + // If Jy has the same centering as Jx, use the same index and shape factor, otherwise + // calculate and use the different shape factor. + const int j_jy = (tby.type(0) == tbx.type(0) ? j_jx : compute_shape_factor(sx_2, xmid_jy)); + const amrex::Real (&sx_jy)[depos_order + 1] = (tby.type(0) == tbx.type(0) ? sx_jx : sx_2); + // If Jz has the same centering as Jx, use the same index and shape factor, otherwise + // calculate and use the different shape factor, unless it has been already calculated for Jy. + const int j_jz = (tbz.type(0) == tbx.type(0) ? j_jx : + (tbz.type(0) == tby.type(0) ? j_jy : + compute_shape_factor(sx_2, xmid_jz))); + const amrex::Real (&sx_jz)[depos_order + 1] = (tbz.type(0) == tbx.type(0) ? sx_jx : sx_2); #if (defined WARPX_DIM_3D) // y direction - const amrex::Real ymid_xz = (yp[ip] - ymin_xz)*dyi - dts2dy*vy; - const amrex::Real ymid_y = (yp[ip] - ymin_y )*dyi - dts2dy*vy; - amrex::Real sy_xz[depos_order + 1]; - const int k_xz = compute_shape_factor(sy_xz, ymid_xz); - amrex::Real sy_y[depos_order + 1]; - const int k_y = compute_shape_factor(sy_y, ymid_y); + const amrex::Real ymid_jx = (yp[ip] - ymin_jx)*dxi - dts2dx*vx; + const amrex::Real ymid_jy = (yp[ip] - ymin_jy)*dxi - dts2dx*vx; + const amrex::Real ymid_jz = (yp[ip] - ymin_jz)*dxi - dts2dx*vx; + amrex::Real sy_jx[depos_order + 1]; + amrex::Real sy_2[depos_order + 1]; + const int k_jx = compute_shape_factor(sy_jx, ymid_jx); + const int k_jy = (tby.type(1) == tbx.type(1) ? k_jx : compute_shape_factor(sy_2, ymid_jy)); + const amrex::Real (&sy_jy)[depos_order + 1] = (tby.type(1) == tbx.type(1) ? sy_jx : sy_2); + const int k_jz = (tbz.type(1) == tbx.type(1) ? k_jx : + (tbz.type(1) == tby.type(1) ? k_jy : + compute_shape_factor(sy_2, ymid_jz))); + const amrex::Real (&sy_jz)[depos_order + 1] = (tbz.type(1) == tbx.type(1) ? sy_jx : sy_2); #endif // z direction - const amrex::Real zmid_xy = (zp[ip] - zmin_xy)*dzi - dts2dz*vz; - const amrex::Real zmid_z = (zp[ip] - zmin_z )*dzi - dts2dz*vz; - amrex::Real sz_xy[depos_order + 1]; - const int l_xy = compute_shape_factor(sz_xy, zmid_xy); - amrex::Real sz_z[depos_order + 1]; - const int l_z = compute_shape_factor(sz_z, zmid_z); + const amrex::Real zmid_jx = (zp[ip] - zmin_jx)*dxi - dts2dx*vx; + const amrex::Real zmid_jy = (zp[ip] - zmin_jy)*dxi - dts2dx*vx; + const amrex::Real zmid_jz = (zp[ip] - zmin_jz)*dxi - dts2dx*vx; + amrex::Real sz_jx[depos_order + 1]; + amrex::Real sz_2[depos_order + 1]; + const int l_jx = compute_shape_factor(sz_jx, zmid_jx); + const int l_jy = (tby.type(2) == tbx.type(2) ? l_jx : compute_shape_factor(sz_2, zmid_jy)); + const amrex::Real (&sz_jy)[depos_order + 1] = (tby.type(2) == tbx.type(2) ? sz_jx : sz_2); + const int l_jz = (tbz.type(2) == tbx.type(2) ? l_jx : + (tbz.type(2) == tby.type(2) ? l_jy : + compute_shape_factor(sz_2, zmid_jz))); + const amrex::Real (&sz_jz)[depos_order + 1] = (tbz.type(2) == tbx.type(2) ? sz_jx : sz_2); // Deposit current into jx_arr, jy_arr and jz_arr #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( - &jx_arr(lo.x+j_x+ix, lo.y+l_xy+iz, 0, 0), - sx_x[ix]*sz_xy[iz]*wqx); + &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0), + sx_jx[ix]*sz_jx[iz]*wqx); amrex::Gpu::Atomic::Add( - &jy_arr(lo.x+j_yz+ix, lo.y+l_xy+iz, 0, 0), - sx_yz[ix]*sz_xy[iz]*wqy); + &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0), + sx_jy[ix]*sz_jy[iz]*wqy); amrex::Gpu::Atomic::Add( - &jz_arr(lo.x+j_yz+ix, lo.y+l_z+iz, 0, 0), - sx_yz[ix]*sz_z[iz]*wqz); + &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0), + sx_jz[ix]*sz_jz[iz]*wqz); } } #elif (defined WARPX_DIM_3D) @@ -167,14 +196,14 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, for (int iy=0; iy<=depos_order; iy++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( - &jx_arr(lo.x+j_x+ix, lo.y+k_xz+iy, lo.z+l_xy+iz), - sx_x[ix]*sy_xz[iy]*sz_xy[iz]*wqx); + &jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz), + sx_jx[ix]*sy_jx[iy]*sz_jx[iz]*wqx); amrex::Gpu::Atomic::Add( - &jy_arr(lo.x+j_yz+ix, lo.y+k_y+iy, lo.z+l_xy+iz), - sx_yz[ix]*sy_y[iy]*sz_xy[iz]*wqy); + &jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz), + sx_jy[ix]*sy_jy[iy]*sz_jy[iz]*wqy); amrex::Gpu::Atomic::Add( - &jz_arr(lo.x+j_yz+ix, lo.y+k_xz+iy, lo.z+l_z+iz), - sx_yz[ix]*sy_xz[iy]*sz_z[iz]*wqz); + &jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz), + sx_jz[ix]*sy_jz[iy]*sz_jz[iz]*wqz); } } } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9fd92ceb4..c7dbd09ba 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -395,29 +395,24 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, WarpX::n_rz_azimuthal_modes); } } else { - // Lower corner of tile box physical domain - // Note that this includes guard cells since it is after tilebox.ngrow - std::array xyzminx = WarpX::LowerCornerWithCentering(tbx, depos_lev); - std::array xyzminy = WarpX::LowerCornerWithCentering(tby, depos_lev); - std::array xyzminz = WarpX::LowerCornerWithCentering(tbz, depos_lev); if (WarpX::nox == 1){ doDepositionShapeN<1>( xp, yp, zp, 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, - xyzminx, xyzminy, xyzminz, lo, q); + tbx, tby, tbz, depos_lev, lo, q); } else if (WarpX::nox == 2){ doDepositionShapeN<2>( xp, yp, zp, 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, - xyzminx, xyzminy, xyzminz, lo, q); + tbx, tby, tbz, depos_lev, lo, q); } else if (WarpX::nox == 3){ doDepositionShapeN<3>( xp, yp, zp, 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, - xyzminx, xyzminy, xyzminz, lo, q); + tbx, tby, tbz, depos_lev, lo, q); } } BL_PROFILE_VAR_STOP(blp_deposit); -- cgit v1.2.3 From 6d69f086406831458de3eae58a941003d66993d6 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 4 Nov 2019 13:35:29 -0800 Subject: Reworked generalized current deposition centering --- Source/Particles/Deposition/CurrentDeposition.H | 157 ++++++++++++------------ Source/Particles/WarpXParticleContainer.cpp | 24 ++-- 2 files changed, 95 insertions(+), 86 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index f8e0c285d..c1502e311 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -15,16 +15,13 @@ required to have the charge of each macroparticle since q is a scalar. For non-ionizable species, ion_lev is a null pointer. - * \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 jx_fab : FArrayBox of current density, either full array or tile. + * \param jy_fab : FArrayBox of current density, either full array or tile. + * \param jz_fab : FArrayBox of current density, either full array or tile. * \param np_to_depose : Number of particles for which current is deposited. * \param dt : Time step for particle level * \param dx : 3D cell size - * \param tbx : Box of domain for jx. - * \param tby : Box of domain for jy. - * \param tbz : Box of domain for jz. - * \param depos_lev : Refinement level being used for the deposition. + * \param xyzmin : Physical lower bounds of domain. * \param lo : Index lower bounds of domain. * /param q : species charge. */ @@ -37,15 +34,12 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, const amrex::ParticleReal * const uyp, const amrex::ParticleReal * const uzp, const int * const ion_lev, - const amrex::Array4& jx_arr, - const amrex::Array4& jy_arr, - const amrex::Array4& jz_arr, + amrex::FArrayBox& jx_fab, + amrex::FArrayBox& jy_fab, + amrex::FArrayBox& jz_fab, const long np_to_depose, const amrex::Real dt, const std::array& dx, - const amrex::Box& tbx, - const amrex::Box& tby, - const amrex::Box& tbz, - const int depos_lev, + const std::array& xyzmin, const amrex::Dim3 lo, const amrex::Real q) { @@ -64,31 +58,31 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, const amrex::Real invvol = dxi*dyi*dzi; #endif - const std::array xyzmin_jx = WarpX::LowerCornerWithCentering(tbx, depos_lev); - const std::array xyzmin_jy = WarpX::LowerCornerWithCentering(tby, depos_lev); - const std::array xyzmin_jz = WarpX::LowerCornerWithCentering(tbz, depos_lev); - const amrex::Real xmin_jx = xyzmin_jx[0]; - const amrex::Real xmin_jy = xyzmin_jy[0]; - const amrex::Real xmin_jz = xyzmin_jz[0]; -#if (defined WARPX_DIM_3D) - const amrex::Real ymin_jx = xyzmin_jx[1]; - const amrex::Real ymin_jy = xyzmin_jy[1]; - const amrex::Real ymin_jz = xyzmin_jz[1]; -#endif - const amrex::Real zmin_jx = xyzmin_jx[2]; - const amrex::Real zmin_jy = xyzmin_jy[2]; - const amrex::Real zmin_jz = xyzmin_jz[2]; + const amrex::Real xmin = xyzmin[0]; + const amrex::Real ymin = xyzmin[1]; + const amrex::Real zmin = xyzmin[2]; const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c; - // Loop over particles and deposit into jx_arr, jy_arr and jz_arr + amrex::Array4 const& jx_arr = jx_fab.array(); + amrex::Array4 const& jy_arr = jy_fab.array(); + amrex::Array4 const& jz_arr = jz_fab.array(); + amrex::IntVect const jx_type = jx_fab.box().type(); + amrex::IntVect const jy_type = jy_fab.box().type(); + amrex::IntVect const jz_type = jz_fab.box().type(); + + constexpr int zdir = (AMREX_SPACEDIM - 1); + constexpr int NODE = amrex::IndexType::NODE; + constexpr int CELL = amrex::IndexType::CELL; + + // Loop over particles and deposit into jx_fab, jy_fab and jz_fab amrex::ParallelFor( np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { // --- Get particle quantities const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq - + uyp[ip]*uyp[ip]*clightsq - + uzp[ip]*uzp[ip]*clightsq); + + uyp[ip]*uyp[ip]*clightsq + + uzp[ip]*uzp[ip]*clightsq); amrex::Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; @@ -124,61 +118,70 @@ void doDepositionShapeN(const amrex::ParticleReal * const xp, // x direction // Get particle position after 1/2 push back in position #if (defined WARPX_DIM_RZ) - const amrex::Real xmid_jx = (rpmid - xmin_jx)*dxi; - const amrex::Real xmid_jy = (rpmid - xmin_jy)*dxi; - const amrex::Real xmid_jz = (rpmid - xmin_jz)*dxi; + const amrex::Real xmid = (rpmid - xmin)*dxi; #else - const amrex::Real xmid_jx = (xp[ip] - xmin_jx)*dxi - dts2dx*vx; - const amrex::Real xmid_jy = (xp[ip] - xmin_jy)*dxi - dts2dx*vx; - const amrex::Real xmid_jz = (xp[ip] - xmin_jz)*dxi - dts2dx*vx; + const amrex::Real xmid = (xp[ip] - xmin)*dxi - dts2dx*vx; #endif - // j_[xyz]: leftmost grid point in x that the particle touches for the centering of each current - // sx_[xyz] shape factor along x for the centering of each current + // 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 // There are only two possible centerings, node or cell centered, so at most only two shape factor // arrays will be needed. - // sx_jx is always calculated. sx_2 only if needed (when Jy or Jz have different centerings than Jx). - amrex::Real sx_jx[depos_order + 1]; - amrex::Real sx_2[depos_order + 1]; - const int j_jx = compute_shape_factor(sx_jx, xmid_jx); - // If Jy has the same centering as Jx, use the same index and shape factor, otherwise - // calculate and use the different shape factor. - const int j_jy = (tby.type(0) == tbx.type(0) ? j_jx : compute_shape_factor(sx_2, xmid_jy)); - const amrex::Real (&sx_jy)[depos_order + 1] = (tby.type(0) == tbx.type(0) ? sx_jx : sx_2); - // If Jz has the same centering as Jx, use the same index and shape factor, otherwise - // calculate and use the different shape factor, unless it has been already calculated for Jy. - const int j_jz = (tbz.type(0) == tbx.type(0) ? j_jx : - (tbz.type(0) == tby.type(0) ? j_jy : - compute_shape_factor(sx_2, xmid_jz))); - const amrex::Real (&sx_jz)[depos_order + 1] = (tbz.type(0) == tbx.type(0) ? sx_jx : sx_2); + amrex::Real sx_node[depos_order + 1]; + amrex::Real sx_cell[depos_order + 1]; + int j_node; + int j_cell; + if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) { + j_node = compute_shape_factor(sx_node, xmid); + } + if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) { + j_cell = compute_shape_factor(sx_cell, xmid - 0.5); + } + const amrex::Real (&sx_jx)[depos_order + 1] = ((jx_type[0] == NODE) ? sx_node : sx_cell); + const amrex::Real (&sx_jy)[depos_order + 1] = ((jy_type[0] == NODE) ? sx_node : sx_cell); + const amrex::Real (&sx_jz)[depos_order + 1] = ((jz_type[0] == NODE) ? sx_node : sx_cell); + int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell); + int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell); + int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell); #if (defined WARPX_DIM_3D) // y direction - const amrex::Real ymid_jx = (yp[ip] - ymin_jx)*dxi - dts2dx*vx; - const amrex::Real ymid_jy = (yp[ip] - ymin_jy)*dxi - dts2dx*vx; - const amrex::Real ymid_jz = (yp[ip] - ymin_jz)*dxi - dts2dx*vx; - amrex::Real sy_jx[depos_order + 1]; - amrex::Real sy_2[depos_order + 1]; - const int k_jx = compute_shape_factor(sy_jx, ymid_jx); - const int k_jy = (tby.type(1) == tbx.type(1) ? k_jx : compute_shape_factor(sy_2, ymid_jy)); - const amrex::Real (&sy_jy)[depos_order + 1] = (tby.type(1) == tbx.type(1) ? sy_jx : sy_2); - const int k_jz = (tbz.type(1) == tbx.type(1) ? k_jx : - (tbz.type(1) == tby.type(1) ? k_jy : - compute_shape_factor(sy_2, ymid_jz))); - const amrex::Real (&sy_jz)[depos_order + 1] = (tbz.type(1) == tbx.type(1) ? sy_jx : sy_2); + const amrex::Real ymid = (yp[ip] - ymin)*dyi - dts2dy*vy; + amrex::Real sy_node[depos_order + 1]; + amrex::Real sy_cell[depos_order + 1]; + int k_node; + int k_cell; + if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) { + k_node = compute_shape_factor(sy_node, ymid); + } + if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) { + k_cell = compute_shape_factor(sy_cell, ymid - 0.5); + } + const amrex::Real (&sy_jx)[depos_order + 1] = ((jx_type[1] == NODE) ? sy_node : sy_cell); + const amrex::Real (&sy_jy)[depos_order + 1] = ((jy_type[1] == NODE) ? sy_node : sy_cell); + const amrex::Real (&sy_jz)[depos_order + 1] = ((jz_type[1] == NODE) ? sy_node : sy_cell); + int const k_jx = ((jx_type[1] == NODE) ? k_node : k_cell); + int const k_jy = ((jy_type[1] == NODE) ? k_node : k_cell); + int const k_jz = ((jz_type[1] == NODE) ? k_node : k_cell); #endif + // z direction - const amrex::Real zmid_jx = (zp[ip] - zmin_jx)*dxi - dts2dx*vx; - const amrex::Real zmid_jy = (zp[ip] - zmin_jy)*dxi - dts2dx*vx; - const amrex::Real zmid_jz = (zp[ip] - zmin_jz)*dxi - dts2dx*vx; - amrex::Real sz_jx[depos_order + 1]; - amrex::Real sz_2[depos_order + 1]; - const int l_jx = compute_shape_factor(sz_jx, zmid_jx); - const int l_jy = (tby.type(2) == tbx.type(2) ? l_jx : compute_shape_factor(sz_2, zmid_jy)); - const amrex::Real (&sz_jy)[depos_order + 1] = (tby.type(2) == tbx.type(2) ? sz_jx : sz_2); - const int l_jz = (tbz.type(2) == tbx.type(2) ? l_jx : - (tbz.type(2) == tby.type(2) ? l_jy : - compute_shape_factor(sz_2, zmid_jz))); - const amrex::Real (&sz_jz)[depos_order + 1] = (tbz.type(2) == tbx.type(2) ? sz_jx : sz_2); + const amrex::Real zmid = (zp[ip] - zmin)*dzi - dts2dz*vz; + amrex::Real sz_node[depos_order + 1]; + amrex::Real sz_cell[depos_order + 1]; + int l_node; + int l_cell; + if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) { + l_node = compute_shape_factor(sz_node, zmid); + } + if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) { + l_cell = compute_shape_factor(sz_cell, zmid - 0.5); + } + const amrex::Real (&sz_jx)[depos_order + 1] = ((jx_type[zdir] == NODE) ? sz_node : sz_cell); + const amrex::Real (&sz_jy)[depos_order + 1] = ((jy_type[zdir] == NODE) ? sz_node : sz_cell); + const amrex::Real (&sz_jz)[depos_order + 1] = ((jz_type[zdir] == NODE) ? sz_node : sz_cell); + int const l_jx = ((jx_type[zdir] == NODE) ? l_node : l_cell); + int const l_jy = ((jy_type[zdir] == NODE) ? l_node : l_cell); + int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell); // Deposit current into jx_arr, jy_arr and jz_arr #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c55642065..e31a9030c 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -298,6 +298,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #ifdef AMREX_USE_GPU // No tiling on GPU: jx_ptr points to the full // jx array (same for jy_ptr and jz_ptr). + auto & jx_fab = jx[pti]; + auto & jy_fab = jy[pti]; + auto & jz_fab = jz[pti]; Array4 const& jx_arr = jx->array(pti); Array4 const& jy_arr = jy->array(pti); Array4 const& jz_arr = jz->array(pti); @@ -317,6 +320,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, local_jy[thread_num].setVal(0.0); local_jz[thread_num].setVal(0.0); + auto & jx_fab = local_jx[thread_num]; + auto & jy_fab = local_jy[thread_num]; + auto & jz_fab = local_jz[thread_num]; Array4 const& jx_arr = local_jx[thread_num].array(); Array4 const& jy_arr = local_jy[thread_num].array(); Array4 const& jz_arr = local_jz[thread_num].array(); @@ -329,13 +335,13 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + // Lower corner of tile box physical domain + // Note that this includes guard cells since it is after tilebox.ngrow const Dim3 lo = lbound(tilebox); + const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); BL_PROFILE_VAR_START(blp_deposit); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { - // Lower corner of tile box physical domain - // Note that this includes guard cells since it is after tilebox.ngrow - const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>( xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, @@ -360,20 +366,20 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, doDepositionShapeN<1>( xp, yp, zp, 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, - tbx, tby, tbz, depos_lev, lo, q); + 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, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, - tbx, tby, tbz, depos_lev, lo, q); + 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, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, - tbx, tby, tbz, depos_lev, lo, q); + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, + xyzmin, lo, q); } } BL_PROFILE_VAR_STOP(blp_deposit); -- cgit v1.2.3 From 647ab2871103d86c4d29888e1bbe19de8b127680 Mon Sep 17 00:00:00 2001 From: grote Date: Tue, 5 Nov 2019 09:25:04 -0800 Subject: Fix for deposition generalization for GPU --- Source/Particles/WarpXParticleContainer.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index e31a9030c..d5520ba06 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -298,9 +298,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #ifdef AMREX_USE_GPU // No tiling on GPU: jx_ptr points to the full // jx array (same for jy_ptr and jz_ptr). - auto & jx_fab = jx[pti]; - auto & jy_fab = jy[pti]; - auto & jz_fab = jz[pti]; + auto & jx_fab = jx->get(pti); + auto & jy_fab = jy->get(pti); + auto & jz_fab = jz->get(pti); Array4 const& jx_arr = jx->array(pti); Array4 const& jy_arr = jy->array(pti); Array4 const& jz_arr = jz->array(pti); -- cgit v1.2.3 dd5e1856555e7d10e3e54684caa52c6&follow=1'>[ci] yarn formatGravatar matthewp 1-8/+8 2022-01-07[ci] release (#2339)astro@0.22.9Gravatar github-actions[bot] 28-39/+40 2022-01-07[ci] yarn formatGravatar matthewp 3-7/+8 2022-01-07Handle loading the Code package in the static build (#2337)Gravatar Matthew Phillips 8-4/+87 2022-01-07[ci] update lockfile (#2334)Gravatar Fred K. Schott 1-154/+154 2022-01-07[ci] yarn formatGravatar matthewp 1-8/+8 2022-01-07[ci] release (#2333)astro@0.22.8Gravatar github-actions[bot] 28-39/+40 2022-01-07[ci] collect statsGravatar FredKSchott 1-0/+1 2022-01-06[ci] yarn formatGravatar matthewp 4-54/+54 2022-01-06[ci] update lockfile (#2327)Gravatar Fred K. Schott 1-58/+64 2022-01-06Fix subpath support regressions (#2330)Gravatar Matthew Phillips 12-22/+566 2022-01-06[ci] yarn formatGravatar natemoo-re 1-2/+2 2022-01-06Added "IntelliSense for TypeScript" (#2326)astro@0.22.7Gravatar Morritz 1-0/+17 2022-01-06[ci] collect statsGravatar FredKSchott 1-0/+1 2022-01-06[ci] yarn formatGravatar FredKSchott 1-8/+8 2022-01-05[ci] release (#2320)Gravatar github-actions[bot] 31-54/+46 2022-01-05chore: update compiler (#2324)Gravatar Nate Moore 3-5/+10