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 07795799555d9843017269415195f2226755e4a7 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Mon, 7 Oct 2019 11:37:46 -0700 Subject: Add check 'if (WarpX::gamma_boost > 1.)' and removed unsused functions. --- Source/Particles/PhysicalParticleContainer.H | 3 -- Source/Particles/PhysicalParticleContainer.cpp | 30 ++------------------ Source/Particles/WarpXParticleContainer.H | 8 ------ Source/Particles/WarpXParticleContainer.cpp | 39 -------------------------- 4 files changed, 3 insertions(+), 77 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index eec516888..b4081e959 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -135,9 +135,6 @@ public: amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, amrex::Real q_tot, long npart, int do_symmetrize); - void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, - std::array u, amrex::Real weight); - void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, std::array u, amrex::Real weight, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7c8c81eed..75e438454 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -236,6 +236,9 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, Gpu::HostVector& particle_uz, Gpu::HostVector& particle_w) { + if (WarpX::gamma_boost > 1.) { + MapParticletoBoostedFrame(x, y, z, u); + } particle_x.push_back(x); particle_y.push_back(y); particle_z.push_back(z); @@ -245,33 +248,6 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, particle_w.push_back(weight); } -void -PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, - std::array u, - Real weight) -{ - std::array attribs; - attribs.fill(0.0); - - // update attribs with input arguments - if (WarpX::gamma_boost > 1.) { - MapParticletoBoostedFrame(x, y, z, u); - } - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - attribs[PIdx::w ] = weight; - - if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) - { - // need to create old values - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - } - - // add particle - AddOneParticle(0, 0, 0, x, y, z, attribs); -} - void PhysicalParticleContainer::AddParticles (int lev) { diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7b0d2d1d0..ceb88d024 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -236,14 +236,6 @@ public: const amrex::ParticleReal* vx, const amrex::ParticleReal* vy, const amrex::ParticleReal* vz, int nattr, const amrex::ParticleReal* attr, int uniqueparticles, int id=-1); - void AddOneParticle (int lev, int grid, int tile, - amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, - std::array& attribs); - - void AddOneParticle (ParticleTileType& particle_tile, - amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, - std::array& attribs); - virtual void ReadHeader (std::istream& is); virtual void WriteHeader (std::ostream& os) const; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 65a82f233..7fb57500d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -134,45 +134,6 @@ WarpXParticleContainer::AllocData () resizeData(); } -void -WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, - ParticleReal x, ParticleReal y, ParticleReal z, - std::array& attribs) -{ - auto& particle_tile = DefineAndReturnParticleTile(lev, grid, tile); - AddOneParticle(particle_tile, x, y, z, attribs); -} - -void -WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, - ParticleReal x, ParticleReal y, ParticleReal z, - std::array& attribs) -{ - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); -#if (AMREX_SPACEDIM == 3) - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; -#elif (AMREX_SPACEDIM == 2) -#ifdef WARPX_DIM_RZ - attribs[PIdx::theta] = std::atan2(y, x); - x = std::sqrt(x*x + y*y); -#endif - p.pos(0) = x; - p.pos(1) = z; -#endif - - particle_tile.push_back(p); - particle_tile.push_back_real(attribs); - - for (int i = PIdx::nattribs; i < NumRealComps(); ++i) - { - particle_tile.push_back_real(i, 0.0); - } -} - void WarpXParticleContainer::AddNParticles (int lev, int n, const ParticleReal* x, const ParticleReal* y, const ParticleReal* z, -- 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 From bdef308e6c5b4aeed8190b6ecdb25b00a51ca5f9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 6 Nov 2019 16:22:41 -0800 Subject: Implement space-charge initialization (non-relativistic, single-level) --- Docs/source/running_cpp/parameters.rst | 4 + .../space_charge_initialization/analysis.py | 85 +++++++++++++++ .../Modules/space_charge_initialization/inputs | 29 ++++++ Python/pywarpx/picmi.py | 12 +-- Regression/WarpX-tests.ini | 34 ++++++ Source/Initialization/InitSpaceChargeField.cpp | 116 +++++++++++++++++++++ Source/Initialization/Make.package | 2 + Source/Initialization/WarpXInitData.cpp | 83 ++------------- Source/Make.WarpX | 1 + Source/Particles/PhysicalParticleContainer.cpp | 1 + Source/Particles/WarpXParticleContainer.H | 3 +- Source/Particles/WarpXParticleContainer.cpp | 82 ++++++--------- Source/WarpX.H | 22 +--- 13 files changed, 326 insertions(+), 148 deletions(-) create mode 100755 Examples/Modules/space_charge_initialization/analysis.py create mode 100644 Examples/Modules/space_charge_initialization/inputs create mode 100644 Source/Initialization/InitSpaceChargeField.cpp (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index e5727e376..50803560d 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -224,6 +224,10 @@ Particle initialization initialization. This can be required whith a moving window and/or when running in a boosted frame. +* ``.initialize_self_fields`` (`0` or `1`) + Whether to calculate the space-charge fields associated with this species + at the beginning of the simulation. + * ``.profile`` (`string`) Density profile for this species. The options are: diff --git a/Examples/Modules/space_charge_initialization/analysis.py b/Examples/Modules/space_charge_initialization/analysis.py new file mode 100755 index 000000000..30954c41d --- /dev/null +++ b/Examples/Modules/space_charge_initialization/analysis.py @@ -0,0 +1,85 @@ +#! /usr/bin/env python +""" +This script checks the space-charge initialization routine, by +verifying that the space-charge field of a Gaussian beam corresponds to +the expected theoretical field. +""" +import sys +import yt +import numpy as np +import matplotlib.pyplot as plt +import scipy.constants as scc +from scipy.special import gammainc +yt.funcs.mylog.setLevel(0) + +# Parameters from the Simulation +Qtot = -1.e-20 +r0 = 2.e-6 + +# Open data file +filename = sys.argv[1] +ds = yt.load( filename ) +# Extract data +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray().squeeze() +if ds.dimensionality == 2: + # Rename the z dimension as y, so as to make this script work for 2d and 3d + Ey_array = ad0['Ez'].to_ndarray().squeeze() +elif ds.dimensionality == 3: + Ey_array = ad0['Ey'].to_ndarray() + Ez_array = ad0['Ez'].to_ndarray() + +# Extract grid coordinates +Nx, Ny, Nz = ds.domain_dimensions +xmin, ymin, zmin = ds.domain_left_edge.v +Lx, Ly, Lz = ds.domain_width.v +x = xmin + Lx/Nx*(0.5+np.arange(Nx)) +y = ymin + Ly/Ny*(0.5+np.arange(Ny)) +z = zmin + Lz/Nz*(0.5+np.arange(Nz)) + +# Compute theoretical field +if ds.dimensionality == 2: + x_2d, y_2d = np.meshgrid(x, y, indexing='ij') + r2 = x_2d**2 + y_2d**2 + factor = (Qtot/r0)/(2*np.pi*scc.epsilon_0*r2) * (1-np.exp(-r2/(2*r0**2))) + Ex_th = x_2d * factor + Ey_th = y_2d * factor +elif ds.dimensionality == 3: + x_2d, y_2d, z_2d = np.meshgrid(x, y, z, indexing='ij') + r2 = x_2d**2 + y_2d**2 + z_2d**2 + factor = Qtot/(4*np.pi*scc.epsilon_0*r2**1.5) * gammainc(3./2, r2/(2.*r0**2)) + Ex_th = factor*x_2d + Ey_th = factor*y_2d + Ez_th = factor*z_2d + +# Plot theory and data +def make_2d(arr): + if arr.ndim == 3: + return arr[:,:,Nz//2] + else: + return arr +plt.figure(figsize=(10,10)) +plt.subplot(221) +plt.title('Ex: Theory') +plt.imshow(make_2d(Ex_th)) +plt.colorbar() +plt.subplot(222) +plt.title('Ex: Simulation') +plt.imshow(make_2d(Ex_array)) +plt.colorbar() +plt.subplot(223) +plt.title('Ey: Theory') +plt.imshow(make_2d(Ey_th)) +plt.colorbar() +plt.subplot(224) +plt.title('Ey: Simulation') +plt.imshow(make_2d(Ey_array)) +plt.colorbar() +plt.savefig('Comparison.png') + +# Automatically check the results +print( 'Relative error in Ex: %.3f'%(abs(Ex_array-Ex_th).max()/Ex_th.max())) +assert np.allclose( Ex_array, Ex_th, atol=0.15*Ex_th.max() ) +assert np.allclose( Ey_array, Ey_th, atol=0.15*Ey_th.max() ) +if ds.dimensionality == 3: + assert np.allclose( Ez_array, Ez_th, atol=0.15*Ez_th.max() ) diff --git a/Examples/Modules/space_charge_initialization/inputs b/Examples/Modules/space_charge_initialization/inputs new file mode 100644 index 000000000..f26f0bf22 --- /dev/null +++ b/Examples/Modules/space_charge_initialization/inputs @@ -0,0 +1,29 @@ +max_step = 0 +amr.n_cell = 64 64 64 +amr.max_grid_size = 32 +amr.max_level = 0 + +amr.plot_int = 1 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 0 # Is periodic? +geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6 +geometry.prob_hi = 50.e-6 50.e-6 50.e-6 + +particles.nspecies = 1 +particles.species_names = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.initialize_self_fields = 1 +beam.x_rms = 2.e-6 +beam.y_rms = 2.e-6 +beam.z_rms = 2.e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = 0.e-6 +beam.npart = 2000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 0.0 diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index c4e6803d5..6e9c9153b 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -56,7 +56,7 @@ class Species(picmistandard.PICMI_Species): if self.mass is None: self.mass = element.mass*periodictable.constants.atomic_mass_constant - def initialize_inputs(self, layout): + def initialize_inputs(self, layout, initialize_self_fields=False): self.species_number = pywarpx.particles.nspecies pywarpx.particles.nspecies += 1 @@ -68,7 +68,8 @@ class Species(picmistandard.PICMI_Species): else: pywarpx.particles.species_names += ' ' + self.name - self.species = pywarpx.Bucket.Bucket(self.name, mass=self.mass, charge=self.charge, injection_style = 'python') + self.species = pywarpx.Bucket.Bucket(self.name, mass=self.mass, charge=self.charge, + injection_style = 'python', initialize_self_fields=int(initialize_self_fields)) pywarpx.Particles.particles_list.append(self.species) if self.initial_distribution is not None: @@ -77,9 +78,9 @@ class Species(picmistandard.PICMI_Species): picmistandard.PICMI_MultiSpecies.Species_class = Species class MultiSpecies(picmistandard.PICMI_MultiSpecies): - def initialize_inputs(self, layout): + def initialize_inputs(self, layout, initialize_self_fields=False): for species in self.species_instances_list: - species.initialize_inputs(layout) + species.initialize_inputs(layout, initialize_self_fields) class GaussianBunchDistribution(picmistandard.PICMI_GaussianBunchDistribution): @@ -544,8 +545,7 @@ class Simulation(picmistandard.PICMI_Simulation): self.solver.initialize_inputs() for i in range(len(self.species)): - assert not self.initialize_self_fields[i], Exception('WarpX does not support initializing self fields') - self.species[i].initialize_inputs(self.layouts[i]) + self.species[i].initialize_inputs(self.layouts[i], self.initialize_self_fields[i]) for i in range(len(self.lasers)): self.lasers[i].initialize_inputs() diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 19da4ba04..c23a751a8 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -616,6 +616,40 @@ compareParticles = 0 particleTypes = electrons tolerance = 1.e-14 +[space_charge_initialization_2d] +buildDir = . +inputFile = Examples/Modules/space_charge_initialization/inputs +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +runtime_params = warpx.do_dynamic_scheduling=0 +analysisRoutine = Examples/Modules/space_charge_initialization/analysis.py +analysisOutputImage = Comparison.png + +[space_charge_initialization] +buildDir = . +inputFile = Examples/Modules/space_charge_initialization/inputs +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +runtime_params = warpx.do_dynamic_scheduling=0 +analysisRoutine = Examples/Modules/space_charge_initialization/analysis.py +analysisOutputImage = Comparison.png + [particles_in_pml_2d] buildDir = . inputFile = Examples/Tests/particles_in_PML/inputs.2d diff --git a/Source/Initialization/InitSpaceChargeField.cpp b/Source/Initialization/InitSpaceChargeField.cpp new file mode 100644 index 000000000..771078af3 --- /dev/null +++ b/Source/Initialization/InitSpaceChargeField.cpp @@ -0,0 +1,116 @@ + +#include +#include +#include + +#include + +using namespace amrex; + +void +WarpX::InitSpaceChargeField (WarpXParticleContainer& pc) +{ + // Allocate fields for charge and potential + const int num_levels = max_level + 1; + Vector > rho(num_levels); + Vector > phi(num_levels); + const int ng = WarpX::nox; + for (int lev = 0; lev <= max_level; lev++) { + BoxArray nba = boxArray(lev); + nba.surroundingNodes(); + rho[lev].reset(new MultiFab(nba, dmap[lev], 1, ng)); // Make ng big enough/use rho from sim + phi[lev].reset(new MultiFab(nba, dmap[lev], 1, 0)); + phi[lev]->setVal(0.); + } + + // Deposit particle charge density (source of Poisson solver) + bool local = false; + bool reset = true; + pc.DepositCharge(rho, local, reset); + + // Call amrex's multigrid solver + // ----------------------------- + + // Define the boundary conditions + Array lobc, hibc; + for (int idim=0; idimarray(mfi); + const auto& Ex_arr = (*Efield_fp[lev][0])[mfi].array(); + const auto& Ey_arr = (*Efield_fp[lev][1])[mfi].array(); + const auto& Ez_arr = (*Efield_fp[lev][2])[mfi].array(); + +#if (AMREX_SPACEDIM == 3) + amrex::ParallelFor( tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ex_arr(i,j,k) += inv_dx*( phi_arr(i+1,j,k) - phi_arr(i,j,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ey_arr(i,j,k) += inv_dy*( phi_arr(i,j+1,k) - phi_arr(i,j,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ez_arr(i,j,k) += inv_dz*( phi_arr(i,j,k+1) - phi_arr(i,j,k) ); + } + ); +#else + amrex::ParallelFor( tbx, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ex_arr(i,j,k) += inv_dx*( phi_arr(i+1,j,k) - phi_arr(i,j,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ez_arr(i,j,k) += inv_dz*( phi_arr(i,j+1,k) - phi_arr(i,j,k) ); + } + ); +#endif + } + } +} diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 2c6458b6d..0b65c1ab1 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -14,5 +14,7 @@ CEXE_sources += InjectorMomentum.cpp CEXE_headers += CustomDensityProb.H CEXE_headers += CustomMomentumProb.H +CEXE_sources += InitSpaceChargeField.cpp + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Initialization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Initialization diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 0814f369b..4fd5d066d 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -111,9 +111,13 @@ WarpX::InitFromScratch () mypc->AllocData(); mypc->InitData(); -#ifdef USE_OPENBC_POISSON - InitOpenbc(); -#endif + // Loop through species and calculate their space-charge field + for (int ispecies=0; ispeciesnSpecies(); ispecies++){ + WarpXParticleContainer& species = mypc->GetParticleContainer(ispecies); + if (species.initialize_self_fields) { + InitSpaceChargeField(species); + } + } InitPML(); @@ -226,79 +230,6 @@ WarpX::PostRestart () mypc->PostRestart(); } -#ifdef USE_OPENBC_POISSON -void -WarpX::InitOpenbc () -{ -#ifndef BL_USE_MPI - static_assert(false, "must use MPI"); -#endif - - static_assert(AMREX_SPACEDIM == 3, "Openbc is 3D only"); - BL_ASSERT(finestLevel() == 0); - - const int lev = 0; - - const Geometry& gm = Geom(lev); - const Box& gbox = gm.Domain(); - int lohi[6]; - warpx_openbc_decompose(gbox.loVect(), gbox.hiVect(), lohi, lohi+3); - - int nprocs = ParallelDescriptor::NProcs(); - int myproc = ParallelDescriptor::MyProc(); - Vector alllohi(6*nprocs,100000); - - MPI_Allgather(lohi, 6, MPI_INT, alllohi.data(), 6, MPI_INT, ParallelDescriptor::Communicator()); - - BoxList bl{IndexType::TheNodeType()}; - for (int i = 0; i < nprocs; ++i) - { - bl.push_back(Box(IntVect(alllohi[6*i ],alllohi[6*i+1],alllohi[6*i+2]), - IntVect(alllohi[6*i+3],alllohi[6*i+4],alllohi[6*i+5]), - IndexType::TheNodeType())); - } - BoxArray ba{bl}; - - Vector iprocmap(nprocs+1); - std::iota(iprocmap.begin(), iprocmap.end(), 0); - iprocmap.back() = myproc; - - DistributionMapping dm{iprocmap}; - - MultiFab rho_openbc(ba, dm, 1, 0); - MultiFab phi_openbc(ba, dm, 1, 0); - - bool local = true; - const std::unique_ptr& rho = mypc->GetChargeDensity(lev, local); - - rho_openbc.setVal(0.0); - rho_openbc.copy(*rho, 0, 0, 1, rho->nGrow(), 0, gm.periodicity(), FabArrayBase::ADD); - - const Real* dx = gm.CellSize(); - - warpx_openbc_potential(rho_openbc[myproc].dataPtr(), phi_openbc[myproc].dataPtr(), dx); - - BoxArray nba = boxArray(lev); - nba.surroundingNodes(); - MultiFab phi(nba, DistributionMap(lev), 1, 0); - phi.copy(phi_openbc, gm.periodicity()); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(phi); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - warpx_compute_E(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_3D(phi[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][0])[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][1])[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][2])[mfi]), - dx); - } -} -#endif - void WarpX::InitLevelData (int lev, Real time) { diff --git a/Source/Make.WarpX b/Source/Make.WarpX index b81fed147..ec004f8fd 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -76,6 +76,7 @@ include $(AMREX_HOME)/Src/Base/Make.package include $(AMREX_HOME)/Src/Particle/Make.package include $(AMREX_HOME)/Src/Boundary/Make.package include $(AMREX_HOME)/Src/AmrCore/Make.package +include $(AMREX_HOME)/Src/LinearSolvers/MLMG/Make.package ifeq ($(USE_SENSEI_INSITU),TRUE) include $(AMREX_HOME)/Src/Amr/Make.package diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 51690d659..6c2aef8a6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -41,6 +41,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); + pp.query("initialize_self_fields", initialize_self_fields); // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index dbd913c5b..cdfe112af 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -185,7 +185,7 @@ public: const amrex::MultiFab& Bz) = 0; void DepositCharge(amrex::Vector >& rho, - bool local = false); + bool local = false, bool reset = true); std::unique_ptr GetChargeDensity(int lev, bool local = false); virtual void DepositCharge(WarpXParIter& pti, @@ -254,6 +254,7 @@ public: void setNextID(int next_id) { ParticleType::NextID(next_id); } bool do_splitting = false; + bool initialize_self_fields = false; // split along diagonals (0) or axes (1) int split_type = 0; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 7fb57500d..29a7d95c7 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -501,68 +501,56 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, } void -WarpXParticleContainer::DepositCharge (Vector >& rho, bool local) +WarpXParticleContainer::DepositCharge (Vector >& rho, + bool local, bool reset) { int num_levels = rho.size(); int finest_level = num_levels - 1; - // each level deposits it's own particles + // each level deposits its own particles const int ng = rho[0]->nGrow(); for (int lev = 0; lev < num_levels; ++lev) { - rho[lev]->setVal(0.0, ng); - - const auto& gm = m_gdb->Geom(lev); - const auto& ba = m_gdb->ParticleBoxArray(lev); - - const Real* dx = gm.CellSize(); - const Real* plo = gm.ProbLo(); - BoxArray nba = ba; - nba.surroundingNodes(); - - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = nba[pti]; + if (reset) rho[lev]->setVal(0.0, ng); +#ifdef _OPENMP +#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 long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - const auto& particles = pti.GetArrayOfStructs(); - int nstride = particles.dataShape().first; - const long np = pti.numParticles(); - FArrayBox& rhofab = (*rho[lev])[pti]; + 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(); + } else { + ion_lev = nullptr; + } - WRPX_DEPOSIT_CIC(particles.dataPtr(), nstride, np, - wp.dataPtr(), &this->charge, - rhofab.dataPtr(), box.loVect(), box.hiVect(), - plo, dx, &ng); + DepositCharge(pti, wp, ion_lev, rho[lev].get(), 0, 0, np, thread_num, lev, lev); } +#ifdef _OPENMP + } +#endif - if (!local) rho[lev]->SumBoundary(gm.periodicity()); - } +#ifdef WARPX_DIM_RZ + WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev); +#endif - // now we average down fine to crse - std::unique_ptr crse; - for (int lev = finest_level - 1; lev >= 0; --lev) { - const BoxArray& fine_BA = rho[lev+1]->boxArray(); - const DistributionMapping& fine_dm = rho[lev+1]->DistributionMap(); - BoxArray coarsened_fine_BA = fine_BA; - coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); - coarsened_fine_data.setVal(0.0); - - IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME - - for (MFIter mfi(coarsened_fine_data); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - const Box& crse_box = coarsened_fine_data[mfi].box(); - const Box& fine_box = (*rho[lev+1])[mfi].box(); - WRPX_SUM_FINE_TO_CRSE_NODAL(bx.loVect(), bx.hiVect(), ratio.getVect(), - coarsened_fine_data[mfi].dataPtr(), crse_box.loVect(), crse_box.hiVect(), - (*rho[lev+1])[mfi].dataPtr(), fine_box.loVect(), fine_box.hiVect()); + if (!local) { + const auto& gm = m_gdb->Geom(lev); + rho[lev]->SumBoundary(gm.periodicity()); } - - rho[lev]->copy(coarsened_fine_data, m_gdb->Geom(lev).periodicity(), FabArrayBase::ADD); } } @@ -841,5 +829,3 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, if (pld.m_lev == lev-1){ } } - - diff --git a/Source/WarpX.H b/Source/WarpX.H index f55670cfb..4275df3c5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -360,22 +360,6 @@ private: /// void EvolveES(int numsteps); - /// - /// Compute the gravitational potential from rho by solving Poisson's equation. - /// Both rho and phi are assumed to be node-centered. This method is only used - /// in electrostatic mode. - /// - void computePhi(const amrex::Vector >& rho, - amrex::Vector >& phi) const; - - /// - /// Compute the electric field in each direction by computing the gradient - /// the potential phi using 2nd order centered differences. Both rho and phi - /// are assumed to be node-centered. This method is only used in electrostatic mode. - /// - void computeE(amrex::Vector, 3> >& E, - const amrex::Vector >& phi) const; - // // This stuff is needed by the nodal multigrid solver when running in // electrostatic mode. @@ -411,7 +395,11 @@ private: void InitFromCheckpoint (); void PostRestart (); - void InitOpenbc (); + void InitSpaceChargeField (WarpXParticleContainer& pc); + void computePhi(const amrex::Vector >& rho, + amrex::Vector >& phi) const; + void computeE(amrex::Vector, 3> >& E, + const amrex::Vector >& phi) const; void InitPML (); void ComputePMLFactors (); -- cgit v1.2.3 From 0bf5b667a13e2e4a8cf58516ead3665a224209b9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 8 Nov 2019 08:54:48 -0800 Subject: Minor fixes --- Source/Particles/WarpXParticleContainer.cpp | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 29a7d95c7..9abeaf131 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -504,21 +504,16 @@ void WarpXParticleContainer::DepositCharge (Vector >& rho, bool local, bool reset) { + // Loop over the refinement levels + for (int lev = 0; lev < max_level; ++lev) { - int num_levels = rho.size(); - int finest_level = num_levels - 1; - - // each level deposits its own particles - const int ng = rho[0]->nGrow(); - for (int lev = 0; lev < num_levels; ++lev) { - - if (reset) rho[lev]->setVal(0.0, ng); + // Reset the `rho` array if `reset` is True + if (reset) rho[lev]->setVal(0.0, rho[lev]->nGrow()); + // Loop over particle tiles and deposit charge #ifdef _OPENMP -#pragma omp parallel + #pragma omp parallel { -#endif -#ifdef _OPENMP int thread_num = omp_get_thread_num(); #else int thread_num = 0; -- cgit v1.2.3 From 074741f7d5fb64f43c1ca692669e407c264c0b88 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 8 Nov 2019 12:18:48 -0800 Subject: Use mesh refinement averaging --- Source/Particles/WarpXParticleContainer.cpp | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9abeaf131..6f195fa87 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -510,7 +510,7 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, // Reset the `rho` array if `reset` is True if (reset) rho[lev]->setVal(0.0, rho[lev]->nGrow()); - // Loop over particle tiles and deposit charge + // Loop over particle tiles and deposit charge on each level #ifdef _OPENMP #pragma omp parallel { @@ -542,10 +542,23 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev); #endif - if (!local) { - const auto& gm = m_gdb->Geom(lev); - rho[lev]->SumBoundary(gm.periodicity()); - } + // Exchange guard cells + if (!local) rho[lev]->SumBoundary( m_gdb->Geom(lev).periodicity() ); + } + + // Now that the charge has been deposited at each level, + // we average down from fine to crse + for (int lev = finest_level - 1; lev >= 0; --lev) { + const DistributionMapping& fine_dm = rho[lev+1]->DistributionMap(); + BoxArray coarsened_fine_BA = rho[lev+1]->boxArray(); + coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); + MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); + coarsened_fine_data.setVal(0.0); + + int const refinement_ratio = 2; + + interpolateDensityFineToCoarse( rho[lev+1], coarsened_fine_data, refinement_ratio ); + rho[lev].ParallelAdd( coarsened_fine_data, m_gdb->Geom(lev).periodicity() ); } } -- cgit v1.2.3 From aea192d6ff5aabfedf09c9f139e831fb5d7bbdd6 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 11 Nov 2019 20:43:52 -0800 Subject: Update DepositCharge --- Source/Parallelization/Make.package | 1 + Source/Parallelization/WarpXComm.cpp | 14 +++++++------- Source/Particles/WarpXParticleContainer.cpp | 9 +++++---- Source/WarpX.H | 23 ----------------------- 4 files changed, 13 insertions(+), 34 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package index c74583522..7c3c38627 100644 --- a/Source/Parallelization/Make.package +++ b/Source/Parallelization/Make.package @@ -2,6 +2,7 @@ CEXE_sources += WarpXComm.cpp CEXE_sources += WarpXRegrid.cpp CEXE_headers += WarpXSumGuardCells.H CEXE_headers += WarpXComm_K.H +CEXE_headers += WarpXComm.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parallelization diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 52df3dc25..851b78748 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1,14 +1,14 @@ +#include "WarpXComm.H" #include #include #include #include -#include -#include +#include +#include #include #include - using namespace amrex; void @@ -503,9 +503,9 @@ WarpX::SyncCurrent () } void -WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, - std::array< amrex::MultiFab *, 3 > const & coarse, - int const refinement_ratio) +interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, + std::array< amrex::MultiFab *, 3 > const & coarse, + int const refinement_ratio) { BL_PROFILE("interpolateCurrentFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); @@ -563,7 +563,7 @@ WarpX::SyncRho () } void -WarpX::interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) +interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) { BL_PROFILE("interpolateDensityFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 6f195fa87..4be339707 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -7,7 +7,7 @@ #include #include #include - +#include // Import low-level single-particle kernels #include #include @@ -505,7 +505,8 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, bool local, bool reset) { // Loop over the refinement levels - for (int lev = 0; lev < max_level; ++lev) { + int const finest_level = rho.size() - 1; + for (int lev = 0; lev < finest_level; ++lev) { // Reset the `rho` array if `reset` is True if (reset) rho[lev]->setVal(0.0, rho[lev]->nGrow()); @@ -557,8 +558,8 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, int const refinement_ratio = 2; - interpolateDensityFineToCoarse( rho[lev+1], coarsened_fine_data, refinement_ratio ); - rho[lev].ParallelAdd( coarsened_fine_data, m_gdb->Geom(lev).periodicity() ); + interpolateDensityFineToCoarse( *rho[lev+1], coarsened_fine_data, refinement_ratio ); + rho[lev]->ParallelAdd( coarsened_fine_data, m_gdb->Geom(lev).periodicity() ); } } diff --git a/Source/WarpX.H b/Source/WarpX.H index 4275df3c5..d17787c57 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -419,29 +419,6 @@ private: std::array, 3> getInterpolatedB(int lev) const; - /** \brief Fills the values of the current on the coarse patch by - * averaging the values of the current of the fine patch (on the same level). - * Also fills the guards of the coarse patch. - * - * \param[in] fine fine patches to interpolate from - * \param[out] coarse coarse patches to interpolate to - * \param[in] refinement_ratio integer ratio between the two - */ - void interpolateCurrentFineToCoarse (std::array< amrex::MultiFab const *, 3 > const & fine, - std::array< amrex::MultiFab *, 3 > const & coarse, - int const refinement_ratio); - - /** \brief Fills the values of the charge density on the coarse patch by - * averaging the values of the charge density of the fine patch (on the same level). - * - * \param[in] fine fine patches to interpolate from - * \param[out] coarse coarse patches to interpolate to - * \param[in] refinement_ratio integer ratio between the two - */ - void interpolateDensityFineToCoarse (const amrex::MultiFab& fine, - amrex::MultiFab& coarse, - int const refinement_ratio); - void ExchangeWithPmlB (int lev); void ExchangeWithPmlE (int lev); void ExchangeWithPmlF (int lev); -- cgit v1.2.3 From 377eeead34efa5b77b666bd879756ed672b4e22f Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 12 Nov 2019 09:29:36 -0800 Subject: Modernize DepositCharge --- Source/Particles/WarpXParticleContainer.cpp | 82 ++++++++++++----------------- 1 file changed, 34 insertions(+), 48 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 7fb57500d..51ed84ea7 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -503,66 +503,54 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, void WarpXParticleContainer::DepositCharge (Vector >& rho, bool local) { + // Loop over the refinement levels + int const finest_level = rho.size() - 1; + for (int lev = 0; lev < finest_level; ++lev) { - int num_levels = rho.size(); - int finest_level = num_levels - 1; - - // each level deposits it's own particles - const int ng = rho[0]->nGrow(); - for (int lev = 0; lev < num_levels; ++lev) { - - rho[lev]->setVal(0.0, ng); - - const auto& gm = m_gdb->Geom(lev); - const auto& ba = m_gdb->ParticleBoxArray(lev); - - const Real* dx = gm.CellSize(); - const Real* plo = gm.ProbLo(); - BoxArray nba = ba; - nba.surroundingNodes(); - - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = nba[pti]; - + // Loop over particle tiles and deposit charge on each level +#ifdef _OPENMP + #pragma omp parallel + { + int thread_num = omp_get_thread_num(); +#else + int thread_num = 0; +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - const auto& particles = pti.GetArrayOfStructs(); - int nstride = particles.dataShape().first; - const long np = pti.numParticles(); - FArrayBox& rhofab = (*rho[lev])[pti]; + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - WRPX_DEPOSIT_CIC(particles.dataPtr(), nstride, np, - wp.dataPtr(), &this->charge, - rhofab.dataPtr(), box.loVect(), box.hiVect(), - plo, dx, &ng); - } + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } - if (!local) rho[lev]->SumBoundary(gm.periodicity()); + DepositCharge(pti, wp, ion_lev, rho[lev].get(), 0, 0, np, thread_num, lev, lev); + } +#ifdef _OPENMP + } +#endif + // Exchange guard cells + if (!local) rho[lev]->SumBoundary( m_gdb->Geom(lev).periodicity() ); } - // now we average down fine to crse - std::unique_ptr crse; + // Now that the charge has been deposited at each level, + // we average down from fine to crse for (int lev = finest_level - 1; lev >= 0; --lev) { - const BoxArray& fine_BA = rho[lev+1]->boxArray(); const DistributionMapping& fine_dm = rho[lev+1]->DistributionMap(); - BoxArray coarsened_fine_BA = fine_BA; + BoxArray coarsened_fine_BA = rho[lev+1]->boxArray(); coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); coarsened_fine_data.setVal(0.0); - IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME + int const refinement_ratio = 2; - for (MFIter mfi(coarsened_fine_data); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - const Box& crse_box = coarsened_fine_data[mfi].box(); - const Box& fine_box = (*rho[lev+1])[mfi].box(); - WRPX_SUM_FINE_TO_CRSE_NODAL(bx.loVect(), bx.hiVect(), ratio.getVect(), - coarsened_fine_data[mfi].dataPtr(), crse_box.loVect(), crse_box.hiVect(), - (*rho[lev+1])[mfi].dataPtr(), fine_box.loVect(), fine_box.hiVect()); - } - - rho[lev]->copy(coarsened_fine_data, m_gdb->Geom(lev).periodicity(), FabArrayBase::ADD); + interpolateDensityFineToCoarse( *rho[lev+1], coarsened_fine_data, refinement_ratio ); + rho[lev]->ParallelAdd( coarsened_fine_data, m_gdb->Geom(lev).periodicity() ); } } @@ -841,5 +829,3 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, if (pld.m_lev == lev-1){ } } - - -- cgit v1.2.3 From 38ffda500337187c6f51422ccdd0ac208b483691 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 12 Nov 2019 09:35:13 -0800 Subject: Remove Fortran files --- Source/Parallelization/WarpXComm.cpp | 9 +- Source/Particles/Make.package | 1 - Source/Particles/WarpXParticleContainer.cpp | 1 + Source/Particles/deposit_cic.F90 | 134 ---------------------------- 4 files changed, 6 insertions(+), 139 deletions(-) delete mode 100644 Source/Particles/deposit_cic.F90 (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 52df3dc25..529d52604 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -503,9 +504,9 @@ WarpX::SyncCurrent () } void -WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, - std::array< amrex::MultiFab *, 3 > const & coarse, - int const refinement_ratio) +interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, + std::array< amrex::MultiFab *, 3 > const & coarse, + int const refinement_ratio) { BL_PROFILE("interpolateCurrentFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); @@ -563,7 +564,7 @@ WarpX::SyncRho () } void -WarpX::interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) +interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) { BL_PROFILE("interpolateDensityFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index a6c124ddd..6d34fa0ef 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -1,4 +1,3 @@ -F90EXE_sources += deposit_cic.F90 F90EXE_sources += interpolate_cic.F90 F90EXE_sources += push_particles_ES.F90 diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 51ed84ea7..2a5be6f5d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include diff --git a/Source/Particles/deposit_cic.F90 b/Source/Particles/deposit_cic.F90 deleted file mode 100644 index 2831ce96b..000000000 --- a/Source/Particles/deposit_cic.F90 +++ /dev/null @@ -1,134 +0,0 @@ -module warpx_ES_deposit_cic - - use iso_c_binding - use amrex_fort_module, only : amrex_real, amrex_particle_real - - implicit none - -contains - -! This routine computes the charge density due to the particles using cloud-in-cell -! deposition. The Fab rho is assumed to be node-centered. -! -! Arguments: -! particles : a pointer to the particle array-of-structs -! ns : the stride length of particle struct (the size of the struct in number of reals) -! np : the number of particles -! weights : the particle weights -! charge : the charge of this particle species -! rho : a Fab that will contain the charge density on exit -! lo : a pointer to the lo corner of this valid box for rho, in index space -! hi : a pointer to the hi corner of this valid box for rho, in index space -! plo : the real position of the left-hand corner of the problem domain -! dx : the mesh spacing -! ng : the number of ghost cells in rho -! - subroutine warpx_deposit_cic_3d(particles, ns, np, & - weights, charge, rho, lo, hi, plo, dx, & - ng) & - bind(c,name='warpx_deposit_cic_3d') - integer, value, intent(in) :: ns, np - real(amrex_particle_real), intent(in) :: particles(ns,np) - real(amrex_particle_real), intent(in) :: weights(np) - real(amrex_real), intent(in) :: charge - integer, intent(in) :: lo(3) - integer, intent(in) :: hi(3) - integer, intent(in) :: ng - real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: plo(3) - real(amrex_real), intent(in) :: dx(3) - - integer i, j, k, n - real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi - real(amrex_real) lx, ly, lz - real(amrex_real) inv_dx(3) - real(amrex_real) qp, inv_vol - - inv_dx = 1.0d0/dx - - inv_vol = inv_dx(1) * inv_dx(2) * inv_dx(3) - - do n = 1, np - - qp = weights(n) * charge * inv_vol - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - lz = (particles(3, n) - plo(3))*inv_dx(3) - - i = floor(lx) - j = floor(ly) - k = floor(lz) - - wx_hi = lx - i - wy_hi = ly - j - wz_hi = lz - k - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - wz_lo = 1.0d0 - wz_hi - - rho(i, j, k ) = rho(i, j, k ) + wx_lo*wy_lo*wz_lo*qp - rho(i, j, k+1) = rho(i, j, k+1) + wx_lo*wy_lo*wz_hi*qp - rho(i, j+1, k ) = rho(i, j+1, k ) + wx_lo*wy_hi*wz_lo*qp - rho(i, j+1, k+1) = rho(i, j+1, k+1) + wx_lo*wy_hi*wz_hi*qp - rho(i+1, j, k ) = rho(i+1, j, k ) + wx_hi*wy_lo*wz_lo*qp - rho(i+1, j, k+1) = rho(i+1, j, k+1) + wx_hi*wy_lo*wz_hi*qp - rho(i+1, j+1, k ) = rho(i+1, j+1, k ) + wx_hi*wy_hi*wz_lo*qp - rho(i+1, j+1, k+1) = rho(i+1, j+1, k+1) + wx_hi*wy_hi*wz_hi*qp - - end do - - end subroutine warpx_deposit_cic_3d - - subroutine warpx_deposit_cic_2d(particles, ns, np, & - weights, charge, rho, lo, hi, plo, dx, & - ng) & - bind(c,name='warpx_deposit_cic_2d') - integer, value, intent(in) :: ns, np - real(amrex_particle_real), intent(in) :: particles(ns,np) - real(amrex_particle_real), intent(in) :: weights(np) - real(amrex_real), intent(in) :: charge - integer, intent(in) :: lo(2) - integer, intent(in) :: hi(2) - integer, intent(in) :: ng - real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: plo(2) - real(amrex_real), intent(in) :: dx(2) - - integer i, j, n - real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi - real(amrex_real) lx, ly - real(amrex_real) inv_dx(2) - real(amrex_real) qp, inv_vol - - inv_dx = 1.0d0/dx - - inv_vol = inv_dx(1) * inv_dx(2) - - do n = 1, np - - qp = weights(n) * charge * inv_vol - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - - i = floor(lx) - j = floor(ly) - - wx_hi = lx - i - wy_hi = ly - j - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - - rho(i, j ) = rho(i, j ) + wx_lo*wy_lo*qp - rho(i, j+1) = rho(i, j+1) + wx_lo*wy_hi*qp - rho(i+1, j ) = rho(i+1, j ) + wx_hi*wy_lo*qp - rho(i+1, j+1) = rho(i+1, j+1) + wx_hi*wy_hi*qp - - end do - - end subroutine warpx_deposit_cic_2d - -end module warpx_ES_deposit_cic -- cgit v1.2.3 From f4925a8a2c9ac1b88178dca21ac78d067cb35c37 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 13 Nov 2019 14:48:50 -0800 Subject: Fix error in DepositCharge --- Source/Particles/WarpXParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 01a2e5ccb..24b7e51a3 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -506,7 +506,7 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, { // Loop over the refinement levels int const finest_level = rho.size() - 1; - for (int lev = 0; lev < finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) { // Reset the `rho` array if `reset` is True if (reset) rho[lev]->setVal(0.0, rho[lev]->nGrow()); -- cgit v1.2.3 From 3983a4153a21164fe2a955a22c1824687e1148a8 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 16 Nov 2019 06:57:40 -0800 Subject: Minor fixes --- Source/Initialization/InitSpaceChargeField.cpp | 9 +++++---- Source/Particles/WarpXParticleContainer.H | 3 ++- Source/Particles/WarpXParticleContainer.cpp | 7 +++++-- 3 files changed, 12 insertions(+), 7 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Initialization/InitSpaceChargeField.cpp b/Source/Initialization/InitSpaceChargeField.cpp index 67cc46ccb..73ee0478e 100644 --- a/Source/Initialization/InitSpaceChargeField.cpp +++ b/Source/Initialization/InitSpaceChargeField.cpp @@ -29,9 +29,10 @@ WarpX::InitSpaceChargeField (WarpXParticleContainer& pc) } // Deposit particle charge density (source of Poisson solver) - bool local = false; - bool reset = true; - pc.DepositCharge(rho, local, reset); + bool const local = false; + bool const reset = true; + bool const do_rz_volume_scaling = true; + pc.DepositCharge(rho, local, reset, do_rz_volume_scaling); // Compute the potential phi, by solving the Poisson equation computePhi( rho, phi ); @@ -114,7 +115,7 @@ WarpX::computeE(amrex::Vector, 3> >& const Box& tby = mfi.tilebox(Ey_nodal_flag); const Box& tbz = mfi.tilebox(Ez_nodal_flag); - const auto& phi_arr = phi[0]->array(mfi); + const auto& phi_arr = phi[lev]->array(mfi); const auto& Ex_arr = (*E[lev][0])[mfi].array(); const auto& Ey_arr = (*E[lev][1])[mfi].array(); const auto& Ez_arr = (*E[lev][2])[mfi].array(); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7da25b452..89dbae7b5 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -185,7 +185,8 @@ public: const amrex::MultiFab& Bz) = 0; void DepositCharge(amrex::Vector >& rho, - bool local = false, bool reset = true); + bool local = false, bool reset = false, + bool do_rz_volume_scaling = false ); std::unique_ptr GetChargeDensity(int lev, bool local = false); virtual void DepositCharge(WarpXParIter& pti, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 24b7e51a3..09f51ef49 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -502,7 +502,8 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, void WarpXParticleContainer::DepositCharge (Vector >& rho, - bool local, bool reset) + bool local, bool reset, + bool do_rz_volume_scaling) { // Loop over the refinement levels int const finest_level = rho.size() - 1; @@ -540,7 +541,9 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, #endif #ifdef WARPX_DIM_RZ - if (reset) WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev); + if (do_rz_volume_scaling) { + WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev); + } #endif // Exchange guard cells -- cgit v1.2.3