diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/BoundaryConditions/PML.H | 7 | ||||
| -rw-r--r-- | Source/BoundaryConditions/PML.cpp | 17 | ||||
| -rw-r--r-- | Source/Evolve/WarpXEvolveEM.cpp | 8 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 49 | ||||
| -rw-r--r-- | Source/FortranInterface/WarpX_f.H | 16 | ||||
| -rw-r--r-- | Source/FortranInterface/WarpX_picsar.F90 | 144 | ||||
| -rw-r--r-- | Source/Initialization/WarpXInitData.cpp | 18 | ||||
| -rw-r--r-- | Source/Laser/LaserParticleContainer.cpp | 14 | ||||
| -rw-r--r-- | Source/Make.WarpX | 2 | ||||
| -rwxr-xr-x | Source/Particles/Deposition/ChargeDeposition.H | 97 | ||||
| -rw-r--r-- | Source/Particles/Deposition/Make.package | 1 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 14 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 4 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 10 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 241 | ||||
| -rw-r--r-- | Source/Utils/WarpXAlgorithmSelection.H | 6 | ||||
| -rw-r--r-- | Source/Utils/WarpXAlgorithmSelection.cpp | 5 | ||||
| -rw-r--r-- | Source/WarpX.H | 7 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 27 |
19 files changed, 320 insertions, 367 deletions
diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index b04938cd8..b34cbe88b 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -101,7 +101,8 @@ public: #ifdef WARPX_USE_PSATD amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif - int do_dive_cleaning, int do_moving_window); + int do_dive_cleaning, int do_moving_window, + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi); void ComputePMLFactors (amrex::Real dt); @@ -172,7 +173,9 @@ private: #endif static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, - const amrex::BoxArray& grid_ba, int ncell); + const amrex::BoxArray& grid_ba, int ncell, + const amrex::IntVect do_pml_Lo, + const amrex::IntVect do_pml_Hi); static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 96bc08af9..21d348482 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -319,11 +319,12 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, #ifdef WARPX_USE_PSATD Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif - int do_dive_cleaning, int do_moving_window) + int do_dive_cleaning, int do_moving_window, + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) : m_geom(geom), m_cgeom(cgeom) { - const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell); + const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi); if (ba.size() == 0) { m_ok = false; return; @@ -399,7 +400,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, BoxArray grid_cba = grid_ba; grid_cba.coarsen(ref_ratio); - const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell); + const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi); DistributionMapping cdm{cba}; @@ -438,12 +439,18 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } BoxArray -PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell) +PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell, + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) { Box domain = geom.Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if ( ! geom.isPeriodic(idim) ) { - domain.grow(idim, ncell); + if (do_pml_Lo[idim]){ + domain.growLo(idim, ncell); + } + if (do_pml_Hi[idim]){ + domain.growHi(idim, ncell); + } } } diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 3f6e5aab4..a532a22fd 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -483,11 +483,17 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time) Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(), cur_time, dt[lev]); #ifdef WARPX_DIM_RZ - // This is called after all particles have deposited their current. + // This is called after all particles have deposited their current and charge. ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev); if (current_buf[lev][0].get()) { ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); } + if (rho_fp[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), lev); + if (charge_buf[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), lev-1); + } + } #endif } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 21000424d..1df05bc0f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -671,4 +671,53 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu }); } } + +void +WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev) +{ + const long ngRho = Rho->nGrow(); + const std::array<Real,3>& dx = WarpX::CellSize(lev); + const Real dr = dx[0]; + + Box tilebox; + + for ( MFIter mfi(*Rho, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + + Array4<Real> const& Rho_arr = Rho->array(mfi); + + tilebox = mfi.tilebox(); + Box tb = convert(tilebox, IntVect::TheUnitVector()); + + // Lower corner of tile box physical domain + // Note that this is done before the tilebox.grow so that + // these do not include the guard cells. + const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, lev); + const Dim3 lo = lbound(tilebox); + const Real rmin = xyzmin[0]; + const int irmin = lo.x; + + // Rescale charge in r-z mode since the inverse volume factor was not + // included in the charge deposition. + amrex::ParallelFor(tb, Rho->nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + { + // Wrap the charge density deposited in the guard cells around + // to the cells above the axis. + // Rho is located on the boundary + if (rmin == 0. && 0 < i && i <= ngRho) { + Rho_arr(i,j,0,icomp) += Rho_arr(-i,j,0,icomp); + } + + // Apply the inverse volume scaling + const amrex::Real r = std::abs(rmin + (i - irmin)*dr); + if (r == 0.) { + // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis + Rho_arr(i,j,0,icomp) /= (MathConst::pi*dr/3.); + } else { + Rho_arr(i,j,0,icomp) /= (2.*MathConst::pi*r); + } + }); + } +} #endif diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 010ca0f4a..6dff3469d 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -75,22 +75,6 @@ extern "C" { #endif - // Charge deposition - void warpx_charge_deposition(amrex::Real* rho, - const long* np, const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, const amrex::Real* w, - const amrex::Real* q, const amrex::Real* xmin, const amrex::Real* ymin, const amrex::Real* zmin, - const amrex::Real* dx, const amrex::Real* dy, const amrex::Real* dz, - const long* nx, const long* ny, const long* nz, - const long* nxguard, const long* nyguard, const long* nzguard, - const long* nox, const long* noy,const long* noz, - const long* lvect, const long* charge_depo_algo); - - // Charge deposition finalize for RZ - void warpx_charge_deposition_rz_volume_scaling( - amrex::Real* rho, const long* rho_ng, const int* rho_ntot, - const amrex::Real* rmin, - const amrex::Real* dr); - // Current deposition void warpx_current_deposition( amrex::Real* jx, const long* jx_ng, const int* jx_ntot, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 2eeaaf5cf..3a9f8f41e 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -7,7 +7,6 @@ #ifdef WARPX_DIM_RZ #define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz -#define WRPX_PXR_RZ_VOLUME_SCALING_RHO apply_rz_volume_scaling_rho #else @@ -49,149 +48,6 @@ module warpx_to_pxr_module contains -! _________________________________________________________________ -!> -!> @brief -!> Main subroutine for the charge deposition -!> -!> @details -!> This subroutines enable to controle the interpolation order -!> via the parameters nox,noy,noz and the type of algorithm via -!> the parameter charge_depo_algo -! -!> @param[inout] rho charge array -!> @param[in] np number of particles -!> @param[in] xp,yp,zp particle position arrays -!> @param[in] w particle weight arrays -!> @param[in] q particle species charge -!> @param[in] xmin,ymin,zmin tile grid minimum position -!> @param[in] dx,dy,dz space discretization steps -!> @param[in] nx,ny,nz number of cells -!> @param[in] nxguard,nyguard,nzguard number of guard cells -!> @param[in] nox,noy,noz interpolation order -!> @param[in] lvect vector length -!> @param[in] charge_depo_algo algorithm choice for the charge deposition -!> -subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,nox,noy,noz,lvect,charge_depo_algo) & - bind(C, name="warpx_charge_deposition") - - integer(c_long), intent(IN) :: np - integer(c_long), intent(IN) :: nx,ny,nz - integer(c_long), intent(IN) :: nxguard,nyguard,nzguard - integer(c_long), intent(IN) :: nox,noy,noz - real(amrex_real), intent(IN OUT) :: rho(*) - real(amrex_real), intent(IN) :: q - real(amrex_real), intent(IN) :: dx,dy,dz - real(amrex_real), intent(IN) :: xmin,ymin,zmin - real(amrex_real), intent(IN), dimension(np) :: xp,yp,zp,w - integer(c_long), intent(IN) :: lvect - integer(c_long), intent(IN) :: charge_depo_algo - - - ! Dimension 3 -#if (AMREX_SPACEDIM==3) - - SELECT CASE(charge_depo_algo) - - ! Scalar classical charge deposition subroutines - CASE(1) - IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN - - CALL depose_rho_scalar_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,lvect) - - ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN - - CALL depose_rho_scalar_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,lvect) - - ELSE IF ((nox.eq.3).and.(noy.eq.3).and.(noz.eq.3)) THEN - - CALL depose_rho_scalar_3_3_3(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,lvect) - - ELSE - CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,nox,noy,noz, & - .TRUE._c_long,.FALSE._c_long) - ENDIF - - ! Optimized subroutines - CASE DEFAULT - - IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN - CALL depose_rho_vecHVv2_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,lvect) - - ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN - CALL depose_rho_vecHVv2_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,lvect) - - ELSE - CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& - nxguard,nyguard,nzguard,nox,noy,noz, & - .TRUE._c_long,.FALSE._c_long) - ENDIF - END SELECT - - ! Dimension 2 -#elif (AMREX_SPACEDIM==2) - -#ifdef WARPX_DIM_RZ - logical(pxr_logical) :: l_2drz = .TRUE._c_long -#else - logical(pxr_logical) :: l_2drz = .FALSE._c_long -#endif - - CALL pxr_depose_rho_n_2dxz(rho,np,xp,yp,zp,w,q,xmin,zmin,dx,dz,nx,nz,& - nxguard,nzguard,nox,noz, & - .TRUE._c_long, .FALSE._c_long, l_2drz, 0_c_long) - -#endif - - end subroutine warpx_charge_deposition - - ! _________________________________________________________________ - !> - !> @brief - !> Applies the inverse volume scaling for RZ charge deposition - !> - !> @details - !> The scaling is done for both single mode (FDTD) and - !> multi mode (spectral) (todo) - ! - !> @param[inout] rho charge array - !> @param[in] rmin tile grid minimum radius - !> @param[in] dr radial space discretization steps - !> @param[in] nx,ny,nz number of cells - !> @param[in] nxguard,nyguard,nzguard number of guard cells - !> - subroutine warpx_charge_deposition_rz_volume_scaling(rho,rho_ng,rho_ntot,rmin,dr) & - bind(C, name="warpx_charge_deposition_rz_volume_scaling") - - integer, intent(in) :: rho_ntot(AMREX_SPACEDIM) - integer(c_long), intent(in) :: rho_ng - real(amrex_real), intent(IN OUT):: rho(*) - real(amrex_real), intent(IN) :: rmin, dr - -#ifdef WARPX_DIM_RZ - integer(c_long) :: type_rz_depose = 1 - - ! Compute the number of valid cells and guard cells - integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) - rho_nvalid = rho_ntot - 2*rho_ng - rho_nguards = rho_ng - -#ifdef WARPX_DIM_RZ - CALL WRPX_PXR_RZ_VOLUME_SCALING_RHO( & - rho,rho_nguards,rho_nvalid, & - rmin,dr,type_rz_depose) - -#endif - - end subroutine warpx_charge_deposition_rz_volume_scaling - ! _________________________________________________________________ !> !> @brief diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 332acb473..590c11b84 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -131,21 +131,35 @@ WarpX::InitPML () { if (do_pml) { + amrex::IntVect do_pml_Lo_corrected = do_pml_Lo; + +#ifdef WARPX_DIM_RZ + do_pml_Lo_corrected[0] = 0; // no PML at r=0, in cylindrical geometry +#endif pml[0].reset(new PML(boxArray(0), DistributionMap(0), &Geom(0), nullptr, pml_ncell, pml_delta, 0, #ifdef WARPX_USE_PSATD dt[0], nox_fft, noy_fft, noz_fft, do_nodal, #endif - do_dive_cleaning, do_moving_window)); + do_dive_cleaning, do_moving_window, + do_pml_Lo_corrected, do_pml_Hi)); for (int lev = 1; lev <= finest_level; ++lev) { + amrex::IntVect do_pml_Lo_MR = amrex::IntVect::TheUnitVector(); +#ifdef WARPX_DIM_RZ + //In cylindrical geometry, if the edge of the patch is at r=0, do not add PML + if ((max_level > 0) && (fine_tag_lo[0]==0.)) { + do_pml_Lo_MR[0] = 0; + } +#endif pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), &Geom(lev), &Geom(lev-1), pml_ncell, pml_delta, refRatio(lev-1)[0], #ifdef WARPX_USE_PSATD dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, #endif - do_dive_cleaning, do_moving_window)); + do_dive_cleaning, do_moving_window, + do_pml_Lo_MR, amrex::IntVect::TheUnitVector())); } } } diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 515aa1f5d..786ebc622 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -453,7 +453,12 @@ LaserParticleContainer::Evolve (int lev, pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + if (crho) { + DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + } + } // // Particle Push @@ -522,7 +527,12 @@ LaserParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + if (crho) { + DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + } + } if (cost) { const Box& tbx = pti.tilebox(); diff --git a/Source/Make.WarpX b/Source/Make.WarpX index a52bd9cfd..e3a33a00f 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -101,7 +101,7 @@ ifeq ($(USE_OPENPMD), TRUE) ifeq (0, $(shell pkg-config "openPMD >= 0.9.0"; echo $$?)) CXXFLAGS += $(shell pkg-config --cflags openPMD) LDFLAGS += $(shell pkg-config --libs openPMD) - LDFLAGS += -Wl,-rpath,$(shell pkg-config --variable=libdir openPMD) + LDFLAGS += -Xlinker -rpath -Xlinker $(shell pkg-config --variable=libdir openPMD) # fallback to manual settings else OPENPMD_LIB_PATH ?= NOT_SET diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H new file mode 100755 index 000000000..a6573b7ab --- /dev/null +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -0,0 +1,97 @@ +#ifndef CHARGEDEPOSITION_H_ +#define CHARGEDEPOSITION_H_ + +#include "ShapeFactors.H" + +/* \brief Charge Deposition for thread thread_num + * /param xp, yp, zp : Pointer to arrays of particle positions. + * \param wp : Pointer to array of particle weights. + * \param rho_arr : Array4 of charge density, either full array or tile. + * \param np_to_depose : Number of particles for which current is deposited. + * \param dx : 3D cell size + * \param xyzmin : Physical lower bounds of domain. + * \param lo : Index lower bounds of domain. + * /param q : species charge. + */ +template <int depos_order> +void doChargeDepositionShapeN(const amrex::Real * const xp, + const amrex::Real * const yp, + const amrex::Real * const zp, + const amrex::Real * const wp, + const amrex::Array4<amrex::Real>& rho_arr, + const long np_to_depose, + const std::array<amrex::Real,3>& dx, + const std::array<amrex::Real, 3> xyzmin, + const amrex::Dim3 lo, + const amrex::Real q) +{ + const amrex::Real dxi = 1.0/dx[0]; + const amrex::Real dzi = 1.0/dx[2]; +#if (AMREX_SPACEDIM == 2) + const amrex::Real invvol = dxi*dzi; +#elif (defined WARPX_DIM_3D) + const amrex::Real dyi = 1.0/dx[1]; + const amrex::Real invvol = dxi*dyi*dzi; +#endif + + const amrex::Real xmin = xyzmin[0]; + const amrex::Real ymin = xyzmin[1]; + const amrex::Real zmin = xyzmin[2]; + + // Loop over particles and deposit into rho_arr + amrex::ParallelFor( + np_to_depose, + [=] AMREX_GPU_DEVICE (long ip) { + // --- Get particle quantities + const amrex::Real wq = q*wp[ip]*invvol; + + // --- Compute shape factors + // x direction + // Get particle position in grid coordinates +#if (defined WARPX_DIM_RZ) + const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real x = (r - xmin)*dxi; +#else + const amrex::Real x = (xp[ip] - xmin)*dxi; +#endif + // Compute shape factors for node-centered quantities + amrex::Real AMREX_RESTRICT sx[depos_order + 1]; + // i: leftmost grid point (node-centered) that the particle touches + const int i = compute_shape_factor<depos_order>(sx, x); + +#if (defined WARPX_DIM_3D) + // y direction + const amrex::Real y = (yp[ip] - ymin)*dyi; + amrex::Real AMREX_RESTRICT sy[depos_order + 1]; + const int j = compute_shape_factor<depos_order>(sy, y); +#endif + // z direction + const amrex::Real z = (zp[ip] - zmin)*dzi; + amrex::Real AMREX_RESTRICT sz[depos_order + 1]; + const int k = compute_shape_factor<depos_order>(sz, z); + + // Deposit charge into rho_arr +#if (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) + for (int iz=0; iz<=depos_order; iz++){ + for (int ix=0; ix<=depos_order; ix++){ + amrex::Gpu::Atomic::Add( + &rho_arr(lo.x+i+ix, lo.y+k+iz, 0), + sx[ix]*sz[iz]*wq); + } + } +#elif (defined WARPX_DIM_3D) + for (int iz=0; iz<=depos_order; iz++){ + for (int iy=0; iy<=depos_order; iy++){ + for (int ix=0; ix<=depos_order; ix++){ + amrex::Gpu::Atomic::Add( + &rho_arr(lo.x+i+ix, lo.y+j+iy, lo.z+k+iz), + sx[ix]*sy[iy]*sz[iz]*wq); + } + } + } +#endif + } + ); +} + +#endif // CHARGEDEPOSITION_H_ diff --git a/Source/Particles/Deposition/Make.package b/Source/Particles/Deposition/Make.package index 0d5ebe2a7..e1aace998 100644 --- a/Source/Particles/Deposition/Make.package +++ b/Source/Particles/Deposition/Make.package @@ -1,3 +1,4 @@ CEXE_headers += CurrentDeposition.H +CEXE_headers += ChargeDeposition.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d4cb71059..c44a9a8b7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1157,7 +1157,12 @@ PhysicalParticleContainer::Evolve (int lev, pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + if (has_buffer){ + DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + } + } if (! do_not_push) { @@ -1295,7 +1300,12 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); } - if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + if (has_buffer){ + DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + } + } if (cost) { const Box& tbx = pti.tilebox(); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index f89efe100..695955e15 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -347,7 +347,9 @@ RigidInjectedParticleContainer::Evolve (int lev, // particles have crossed the inject plane. const Real* plo = Geom(lev).ProbLo(); const Real* phi = Geom(lev).ProbHi(); - done_injecting[lev] = (zinject_plane_levels[lev] < plo[2] || zinject_plane_levels[lev] > phi[2]); + const int zdir = AMREX_SPACEDIM-1; + done_injecting[lev] = ((zinject_plane_levels[lev] < plo[zdir] && WarpX::moving_window_v + WarpX::beta_boost*PhysConst::c >= 0.) || + (zinject_plane_levels[lev] > phi[zdir] && WarpX::moving_window_v + WarpX::beta_boost*PhysConst::c <= 0.)); done_injecting_lev = done_injecting[lev]; PhysicalParticleContainer::Evolve (lev, diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 08aca752d..ee263c27b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -170,13 +170,13 @@ public: virtual void DepositCharge(WarpXParIter& pti, RealVector& wp, - amrex::MultiFab* rhomf, - amrex::MultiFab* crhomf, + amrex::MultiFab* rho, int icomp, - const long np_current, - const long np, + const long offset, + const long np_to_depose, int thread_num, - int lev ); + int lev, + int depos_lev); virtual void DepositCurrent(WarpXParIter& pti, RealVector& wp, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2107df35f..0ff80941d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -12,6 +12,7 @@ #include <GetAndSetPosition.H> #include <UpdatePosition.H> #include <CurrentDeposition.H> +#include <ChargeDeposition.H> using namespace amrex; @@ -553,140 +554,87 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } void -WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, - MultiFab* rhomf, MultiFab* crhomf, int icomp, - const long np_current, - const long np, int thread_num, int lev ) +WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + MultiFab* rho, int icomp, + const long offset, const long np_to_depose, + int thread_num, int lev, int depos_lev) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || + (depos_lev==(lev )), + "Deposition buffers only work for lev-1"); - BL_PROFILE_VAR_NS("PICSAR::ChargeDeposition", blp_pxr_chd); - BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); + // If no particles, do not do anything + if (np_to_depose == 0) return; - const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - const long lvect = 8; + const long ngRho = rho->nGrow(); + const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0)); + const Real q = this->charge; - long ngRho = rhomf->nGrow(); - Real* data_ptr; - Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); + BL_PROFILE_VAR_NS("PPC::ChargeDeposition", blp_ppc_chd); + BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); - const std::array<Real,3>& dx = WarpX::CellSize(lev); - const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); + // Get tile box where charge is deposited. + // The tile box is different when depositing in the buffers (depos_lev<lev) + // or when depositing inside the level (depos_lev=lev) + Box tilebox; + if (lev == depos_lev) { + tilebox = pti.tilebox(); + } else { + const IntVect& ref_ratio = WarpX::RefRatio(depos_lev); + tilebox = amrex::coarsen(pti.tilebox(),ref_ratio); + } + + tilebox.grow(ngRho); - // Deposit charge for particles that are not in the current buffers - if (np_current > 0) - { - const std::array<Real, 3>& xyzmin = xyzmin_tile; + const int nc = (rhomf->nComp() == 1 ? 1 : rhomf->nComp()/2); #ifdef AMREX_USE_GPU - data_ptr = (*rhomf)[pti].dataPtr(icomp*(rhomf->nComp()/2)); - auto rholen = (*rhomf)[pti].length(); + // No tiling on GPU: rho_arr points to the full rho array. + MultiFab rhoi(*rho, amrex::make_alias, icomp*nc, nc); + Array4<Real> const& rho_arr = rhoi.array(pti); #else - tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box, rhomf->nComp()); + // Tiling is on: rho_arr points to local_rho[thread_num] + const Box tb = amrex::convert(tilebox, IntVect::TheUnitVector()); - data_ptr = local_rho[thread_num].dataPtr(); - auto rholen = local_rho[thread_num].length(); + local_rho[thread_num].resize(tb, nc); - local_rho[thread_num].setVal(0.0); -#endif + // local_rho[thread_num] is set to zero + local_rho[thread_num].setVal(0.0); -#if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ngRho; - const long ny = rholen[1]-1-2*ngRho; - const long nz = rholen[2]-1-2*ngRho; -#else - const long nx = rholen[0]-1-2*ngRho; - const long ny = 0; - const long nz = rholen[1]-1-2*ngRho; + Array4<Real> const& rho_arr = local_rho[thread_num].array(); #endif - BL_PROFILE_VAR_START(blp_pxr_chd); - warpx_charge_deposition(data_ptr, &np_current, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - wp.dataPtr(), - &this->charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, - &ngRho, &ngRho, &ngRho, - &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_DIM_RZ - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &xyzmin[0], &dx[0]); -#endif - BL_PROFILE_VAR_STOP(blp_pxr_chd); - -#ifndef AMREX_USE_GPU - BL_PROFILE_VAR_START(blp_accumulate); - - (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp*(rhomf->nComp()/2), (rhomf->nComp()/2)); - BL_PROFILE_VAR_STOP(blp_accumulate); -#endif - } - - // Deposit charge for particles that are in the current buffers - if (np_current < np) - { - const IntVect& ref_ratio = WarpX::RefRatio(lev-1); - const Box& ctilebox = amrex::coarsen(pti.tilebox(), ref_ratio); - const std::array<Real,3>& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); - -#ifdef AMREX_USE_GPU - data_ptr = (*crhomf)[pti].dataPtr(icomp); - auto rholen = (*crhomf)[pti].length(); -#else - tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); - tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box, crhomf->nComp()); - - data_ptr = local_rho[thread_num].dataPtr(); - auto rholen = local_rho[thread_num].length(); - - local_rho[thread_num].setVal(0.0); -#endif + Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; -#if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ngRho; - const long ny = rholen[1]-1-2*ngRho; - const long nz = rholen[2]-1-2*ngRho; -#else - const long nx = rholen[0]-1-2*ngRho; - const long ny = 0; - const long nz = rholen[1]-1-2*ngRho; -#endif + // Lower corner of tile box physical domain + // Note that this includes guard cells since it is after tilebox.ngrow + const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); + // Indices of the lower bound + const Dim3 lo = lbound(tilebox); - long ncrse = np - np_current; - BL_PROFILE_VAR_START(blp_pxr_chd); - warpx_charge_deposition(data_ptr, &ncrse, - m_xp[thread_num].dataPtr() + np_current, - m_yp[thread_num].dataPtr() + np_current, - m_zp[thread_num].dataPtr() + np_current, - wp.dataPtr() + np_current, - &this->charge, - &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], - &cdx[0], &cdx[1], &cdx[2], &nx, &ny, &nz, - &ngRho, &ngRho, &ngRho, - &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_DIM_RZ - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &cxyzmin_tile[0], &cdx[0]); -#endif - BL_PROFILE_VAR_STOP(blp_pxr_chd); + BL_PROFILE_VAR_START(blp_ppc_chd); + if (WarpX::nox == 1){ + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, + np_to_depose, dx, xyzmin, lo, q); + } else if (WarpX::nox == 2){ + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, + np_to_depose, dx, xyzmin, lo, q); + } else if (WarpX::nox == 3){ + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, + np_to_depose, dx, xyzmin, lo, q); + } + BL_PROFILE_VAR_STOP(blp_ppc_chd); #ifndef AMREX_USE_GPU - BL_PROFILE_VAR_START(blp_accumulate); + BL_PROFILE_VAR_START(blp_accumulate); - (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp*(crhomf->nComp()/2), (crhomf->nComp()/2)); + (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp*nc, nc); - BL_PROFILE_VAR_STOP(blp_accumulate); + BL_PROFILE_VAR_STOP(blp_accumulate); #endif - } -}; +} void WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, bool local) @@ -763,8 +711,6 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) BoxArray nba = ba; nba.surroundingNodes(); - const std::array<Real,3>& dx = WarpX::CellSize(lev); - const int ng = WarpX::nox; auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,WarpX::ncomps,ng)); @@ -774,75 +720,28 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #pragma omp parallel { #endif - Cuda::ManagedDeviceVector<Real> xp, yp, zp; #ifdef _OPENMP - FArrayBox rho_loc; + 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 long np = pti.numParticles(); - - pti.GetPosition(xp, yp, zp); + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - // Data on the grid - Real* data_ptr; - FArrayBox& rhofab = (*rho)[pti]; + DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); + } #ifdef _OPENMP - const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); - const std::array<Real, 3>& xyzmin = xyzmin_tile; - tile_box.grow(ng); - rho_loc.resize(tile_box, rho->nComp()); - rho_loc = 0.0; - data_ptr = rho_loc.dataPtr(); - auto rholen = rho_loc.length(); -#else - const Box& box = pti.validbox(); - const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); - const std::array<Real, 3>& xyzmin = xyzmin_grid; - data_ptr = rhofab.dataPtr(); - auto rholen = rhofab.length(); -#endif - -#if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ng; - const long ny = rholen[1]-1-2*ng; - const long nz = rholen[2]-1-2*ng; -#else - const long nx = rholen[0]-1-2*ng; - const long ny = 0; - const long nz = rholen[1]-1-2*ng; + } #endif - long nxg = ng; - long nyg = ng; - long nzg = ng; - long lvect = 8; - - warpx_charge_deposition(data_ptr, - &np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), wp.dataPtr(), - &this->charge, &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, - &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); #ifdef WARPX_DIM_RZ - long ngRho = WarpX::nox; - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &xyzmin[0], &dx[0]); -#endif - -#ifdef _OPENMP - rhofab.atomicAdd(rho_loc); - } + WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho.get(), lev); #endif - } if (!local) rho->SumBoundary(gm.periodicity()); diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 3fb23698a..6a32513b7 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -34,11 +34,9 @@ struct CurrentDepositionAlgo { }; struct ChargeDepositionAlgo { - // These numbers corresponds to the algorithm code in WarpX's - // `warpx_charge_deposition` function + // Only the Standard algorithm is implemented enum { - Vectorized = 0, - Standard = 1 + Standard = 0 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 3aa4eb7b7..842085a36 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -31,12 +31,7 @@ const std::map<std::string, int> current_deposition_algo_to_int = { const std::map<std::string, int> charge_deposition_algo_to_int = { {"standard", ChargeDepositionAlgo::Standard }, -#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D - {"vectorized", ChargeDepositionAlgo::Vectorized }, - {"default", ChargeDepositionAlgo::Vectorized } -#else {"default", ChargeDepositionAlgo::Standard } -#endif }; const std::map<std::string, int> gathering_algo_to_int = { diff --git a/Source/WarpX.H b/Source/WarpX.H index c4bcca4bb..cc44541ee 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -187,6 +187,9 @@ public: amrex::MultiFab* Jy, amrex::MultiFab* Jz, int lev); + + void ApplyInverseVolumeScalingToChargeDensity(amrex::MultiFab* Rho, + int lev); #endif void DampPML (); @@ -258,6 +261,7 @@ public: static int do_moving_window; static int moving_window_dir; + static amrex::Real moving_window_v; // slice generation // void InitializeSliceMultiFabs (); @@ -500,10 +504,11 @@ private: int do_pml = 1; int pml_ncell = 10; int pml_delta = 10; + amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(); + amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector(); amrex::Vector<std::unique_ptr<PML> > pml; amrex::Real moving_window_x = std::numeric_limits<amrex::Real>::max(); - amrex::Real moving_window_v = std::numeric_limits<amrex::Real>::max(); amrex::Real current_injection_position = 0; // Plasma injection parameters diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c26d83253..2d6df625f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -30,6 +30,7 @@ Vector<Real> WarpX::B_external(3, 0.0); int WarpX::do_moving_window = 0; int WarpX::moving_window_dir = -1; +Real WarpX::moving_window_v = std::numeric_limits<amrex::Real>::max(); Real WarpX::gamma_boost = 1.; Real WarpX::beta_boost = 0.; @@ -386,6 +387,22 @@ WarpX::ReadParameters () pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); + Vector<int> parse_do_pml_Lo(AMREX_SPACEDIM,1); + pp.queryarr("do_pml_Lo", parse_do_pml_Lo); + do_pml_Lo[0] = parse_do_pml_Lo[0]; + do_pml_Lo[1] = parse_do_pml_Lo[1]; +#if (AMREX_SPACEDIM == 3) + do_pml_Lo[2] = parse_do_pml_Lo[2]; +#endif + Vector<int> parse_do_pml_Hi(AMREX_SPACEDIM,1); + pp.queryarr("do_pml_Hi", parse_do_pml_Hi); + do_pml_Hi[0] = parse_do_pml_Hi[0]; + do_pml_Hi[1] = parse_do_pml_Hi[1]; +#if (AMREX_SPACEDIM == 3) + do_pml_Hi[2] = parse_do_pml_Hi[2]; +#endif + + pp.query("dump_openpmd", dump_openpmd); pp.query("dump_plotfiles", dump_plotfiles); pp.query("plot_raw_fields", plot_raw_fields); @@ -396,7 +413,7 @@ WarpX::ReadParameters () if (not user_fields_to_plot){ // If not specified, set default values fields_to_plot = {"Ex", "Ey", "Ez", "Bx", "By", - "Bz", "jx", "jy", "jz", + "Bz", "jx", "jy", "jz", "part_per_cell"}; } // set plot_rho to true of the users requests it, so that @@ -414,9 +431,9 @@ WarpX::ReadParameters () // If user requests to plot proc_number for a serial run, // delete proc_number from fields_to_plot if (ParallelDescriptor::NProcs() == 1){ - fields_to_plot.erase(std::remove(fields_to_plot.begin(), - fields_to_plot.end(), - "proc_number"), + fields_to_plot.erase(std::remove(fields_to_plot.begin(), + fields_to_plot.end(), + "proc_number"), fields_to_plot.end()); } @@ -504,7 +521,7 @@ WarpX::ReadParameters () { ParmParse pp("algo"); // If not in RZ mode, read use_picsar_deposition - // In RZ mode, use_picsar_deposition is on, as the C++ version + // In RZ mode, use_picsar_deposition is on, as the C++ version // of the deposition does not support RZ pp.query("use_picsar_deposition", use_picsar_deposition); current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); |