From 8bfa3488bae938d4c1c4ec9f71eababa51556324 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Fri, 8 Mar 2019 11:10:17 -0800 Subject: DeviceVector -> ManagedDeviceVector --- Source/Particles/PhysicalParticleContainer.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index afc1412d5..d98e79177 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -927,7 +927,7 @@ PhysicalParticleContainer::FieldGather (int lev, #pragma omp parallel #endif { - Cuda::DeviceVector xp, yp, zp; + Cuda::ManagedDeviceVector xp, yp, zp; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -1444,7 +1444,7 @@ PhysicalParticleContainer::SplitParticles(int lev) { auto& mypc = WarpX::GetInstance().GetPartContainer(); auto& pctmp_split = mypc.GetPCtmp(); - Cuda::DeviceVector xp, yp, zp; + Cuda::ManagedDeviceVector xp, yp, zp; RealVector psplit_x, psplit_y, psplit_z, psplit_w; RealVector psplit_ux, psplit_uy, psplit_uz; long np_split_to_add = 0; @@ -1593,10 +1593,10 @@ PhysicalParticleContainer::SplitParticles(int lev) void PhysicalParticleContainer::PushPX(WarpXParIter& pti, - Cuda::DeviceVector& xp, - Cuda::DeviceVector& yp, - Cuda::DeviceVector& zp, - Cuda::DeviceVector& giv, + Cuda::ManagedDeviceVector& xp, + Cuda::ManagedDeviceVector& yp, + Cuda::ManagedDeviceVector& zp, + Cuda::ManagedDeviceVector& giv, Real dt) { -- cgit v1.2.3 From d767bc8b4ac6c374928f98b9807cd361df5291d9 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Tue, 12 Mar 2019 14:03:56 -0700 Subject: Implementation of axisymmetric solver, mode 0 only --- Examples/Tests/Langmuir/inputs.2d.rz | 58 ++++++++++++++++ GNUmakefile | 1 + Source/Diagnostics/WarpXIO.cpp | 2 +- Source/Evolve/WarpXEvolveEM.cpp | 5 ++ Source/FieldSolver/WarpXPushFieldsEM.cpp | 6 +- Source/FortranInterface/WarpX_f.H | 15 ++++- Source/FortranInterface/WarpX_picsar.F90 | 91 ++++++++++++++++++++++++-- Source/Initialization/PlasmaInjector.H | 10 +-- Source/Initialization/PlasmaInjector.cpp | 3 + Source/Make.WarpX | 3 +- Source/Particles/PhysicalParticleContainer.cpp | 43 +++++++++++- Source/Particles/WarpXParticleContainer.H | 6 +- Source/Particles/WarpXParticleContainer.cpp | 33 ++++++++-- 13 files changed, 249 insertions(+), 27 deletions(-) create mode 100644 Examples/Tests/Langmuir/inputs.2d.rz (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Tests/Langmuir/inputs.2d.rz b/Examples/Tests/Langmuir/inputs.2d.rz new file mode 100644 index 000000000..f3aa30911 --- /dev/null +++ b/Examples/Tests/Langmuir/inputs.2d.rz @@ -0,0 +1,58 @@ +# Maximum number of time steps +max_step = 100 + +# number of grid points +amr.n_cell = 64 64 + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +amr.max_grid_size = 32 +amr.blocking_factor = 32 32 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +amr.plot_int = 1 # How often to write plotfiles. "<= 0" means no plotfiles. + +# Geometry +geometry.coord_sys = 1 # 1: RZ +geometry.is_periodic = 0 1 +geometry.prob_lo = 0.e0 -2e-05 # physical domain +geometry.prob_hi = 2e-05 2e-05 + +# Verbosity +warpx.verbose = 1 +warpx.do_moving_window = 0 +warpx.moving_window_dir = z +warpx.moving_window_v = 0.0 # in units of the speed of light + +# Algorithms +algo.current_deposition = 3 +algo.charge_deposition = 0 +algo.field_gathering = 0 +algo.particle_pusher = 0 +interpolation.nox = 1 +interpolation.noy = 1 +interpolation.noz = 1 + +# CFL +warpx.cfl = 1.0 + +particles.nspecies = 1 +particles.species_names = electrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = 2 2 + +electrons.zmax = 0.0 + +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 0.0 +electrons.uy = 0.0 +electrons.uz = 0.1 +electrons.radially_weighted = true # Only true is supported diff --git a/GNUmakefile b/GNUmakefile index 06f4c6229..02949555e 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -31,6 +31,7 @@ USE_SENSEI_INSITU = FALSE WarpxBinDir = Bin USE_PSATD = FALSE +USE_RZ = FALSE DO_ELECTROSTATIC = FALSE diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index e68b4a4e8..d43d74968 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -1253,7 +1253,7 @@ WarpX::WritePlotFile () const particle_varnames.push_back("By"); particle_varnames.push_back("Bz"); -#if (AMREX_SPACEDIM == 2) && WARPX_RZ +#ifdef WARPX_RZ particle_varnames.push_back("theta"); #endif diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 4c87d8321..60b0b5fa3 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -463,9 +463,14 @@ WarpX::ComputeDt () if (maxwell_fdtd_solver_id == 0) { // CFL time step Yee solver +#ifdef WARPX_RZ + // Derived semi-analytically by R. Lehe + deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); +#else deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), + 1./(dx[1]*dx[1]), + 1./(dx[2]*dx[2]))) * PhysConst::c ); +#endif } else { // CFL time step CKC solver #if (BL_SPACEDIM == 3) diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6d588ae86..97204cc17 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -76,6 +76,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) const Box& tbx = mfi.tilebox(Bx_nodal_flag); const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); + const Real xmin = mfi.tilebox().smallEnd(0)*dx[0]; if (do_nodal) { auto const& Bxfab = Bx->array(mfi); @@ -112,6 +113,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) BL_TO_FORTRAN_3D((*By)[mfi]), BL_TO_FORTRAN_3D((*Bz)[mfi]), &dtsdx, &dtsdy, &dtsdz, + &xmin, &dx[0], &WarpX::maxwell_fdtd_solver_id); } @@ -230,6 +232,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) const Box& tex = mfi.tilebox(Ex_nodal_flag); const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); + const Real xmin = mfi.tilebox().smallEnd(0)*dx[0]; if (do_nodal) { auto const& Exfab = Ex->array(mfi); @@ -272,7 +275,8 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) BL_TO_FORTRAN_3D((*jy)[mfi]), BL_TO_FORTRAN_3D((*jz)[mfi]), &mu_c2_dt, - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, + &xmin, &dx[0]); } if (F) diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 1f5256483..33f7f9639 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -120,6 +120,15 @@ extern "C" const long* nox, const long* noy,const long* noz, const long* lvect, const long* current_depo_algo); + // Current deposition finalize for RZ + void warpx_current_deposition_rz_volume_scaling( + amrex::Real* jx, const long* jx_ng, const int* jx_ntot, + amrex::Real* jy, const long* jy_ng, const int* jy_ntot, + amrex::Real* jz, const long* jz_ng, const int* jz_ntot, + const long* nmodes, + const amrex::Real* rmin, + const amrex::Real* dr); + // Field gathering void warpx_geteb_energy_conserving(const long* np, @@ -214,7 +223,9 @@ extern "C" const amrex::Real* mudt, const amrex::Real* dtsdx, const amrex::Real* dtsdy, - const amrex::Real* dtsdz); + const amrex::Real* dtsdz, + const amrex::Real* xmin, + const amrex::Real* dx); void warpx_push_bvec( const int* xlo, const int* xhi, @@ -229,6 +240,8 @@ extern "C" const amrex::Real* dtsdx, const amrex::Real* dtsdy, const amrex::Real* dtsdz, + const amrex::Real* xmin, + const amrex::Real* dx, const int* maxwell_fdtd_solver_id); void warpx_push_evec_f( diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 8e9d2db3a..346d43b8b 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -10,16 +10,29 @@ #elif (AMREX_SPACEDIM == 2) -#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2dxz_energy_conserving_generic -#define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic_2d -#define WRPX_PXR_PUSH_BVEC pxrpush_em2d_bvec #define WRPX_PXR_PUSH_BVEC_CKC pxrpush_em2d_bvec_ckc #define WRPX_PXR_PUSH_EVEC_F pxrpush_em2d_evec_f #define WRPX_PXR_PUSH_EVEC_F_CKC pxrpush_em2d_evec_f_ckc + +#ifdef WARPX_RZ + +#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2drz_energy_conserving_generic +#define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz +#define WRPX_PXR_RZ_VOLUME_SCALING apply_rz_volume_scaling +#define WRPX_PXR_PUSH_BVEC pxrpush_emrz_bvec +#define WRPX_PXR_PUSH_EVEC pxrpush_emrz_evec + +#else + +#define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2dxz_energy_conserving_generic +#define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic_2d +#define WRPX_PXR_PUSH_BVEC pxrpush_em2d_bvec #define WRPX_PXR_PUSH_EVEC pxrpush_em2d_evec #endif +#endif + #define LVECT_CURRDEPO 8_c_long #define LVECT_FIELDGATHE 64_c_long @@ -301,6 +314,58 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n end subroutine warpx_current_deposition + ! _________________________________________________________________ + !> + !> @brief + !> Applies the inverse volume scaling for RZ current deposition + !> + !> @details + !> The scaling is done for both single mode (FDTD) and + !> multi mode (spectral) + ! + !> @param[inout] jx,jy,jz current arrays + !> @param[in] nmodes number of spectral modes + !> @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_current_deposition_rz_volume_scaling( & + jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot,nmodes, & + rmin,dr) & + bind(C, name="warpx_current_deposition_rz_volume_scaling") + + integer, intent(in) :: jx_ntot(AMREX_SPACEDIM), jy_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) + integer(c_long), intent(in) :: jx_ng, jy_ng, jz_ng + integer(c_long), intent(IN) :: nmodes + real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) + real(amrex_real), intent(IN) :: rmin, dr + + integer(c_long) :: type_rz_depose = 1 + + ! Compute the number of valid cells and guard cells + integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & + jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) + jx_nvalid = jx_ntot - 2*jx_ng + jy_nvalid = jy_ntot - 2*jy_ng + jz_nvalid = jz_ntot - 2*jz_ng + jx_nguards = jx_ng + jy_nguards = jy_ng + jz_nguards = jz_ng + +#ifdef WARPX_RZ + CALL WRPX_PXR_RZ_VOLUME_SCALING( & + jx,jx_nguards,jx_nvalid, & + jy,jy_nguards,jy_nvalid, & + jz,jz_nguards,jz_nvalid, & +#ifdef WARPX_USE_PSATD + nmodes, & +#endif + rmin,dr,type_rz_depose) +#endif + + end subroutine warpx_current_deposition_rz_volume_scaling + ! _________________________________________________________________ !> !> @brief @@ -473,7 +538,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jx, jxlo, jxhi, & jy, jylo, jyhi, & jz, jzlo, jzhi, & - mudt, dtsdx, dtsdy, dtsdz) bind(C, name="warpx_push_evec") + mudt, dtsdx, dtsdy, dtsdz, xmin, dx) bind(C, name="warpx_push_evec") integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), & ylo(BL_SPACEDIM), yhi(BL_SPACEDIM), zlo(BL_SPACEDIM), zhi(BL_SPACEDIM), & @@ -489,6 +554,8 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN) :: mudt, dtsdx, dtsdy, dtsdz + real(amrex_real), intent(IN) :: xmin, dx + CALL WRPX_PXR_PUSH_EVEC(& xlo, xhi, ylo, yhi, zlo, zhi, & ex, exlo, exhi,& @@ -500,7 +567,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jx, jxlo, jxhi,& jy, jylo, jyhi,& jz, jzlo, jzhi,& - mudt, dtsdx, dtsdy, dtsdz); + mudt, dtsdx, dtsdy, dtsdz & +#ifdef WARPX_RZ + ,xmin,dx & +#endif + ) end subroutine warpx_push_evec @@ -526,7 +597,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n bx, bxlo, bxhi, & by, bylo, byhi, & bz, bzlo, bzhi, & - dtsdx, dtsdy, dtsdz, & + dtsdx, dtsdy, dtsdz, xmin, dx, & maxwell_fdtd_solver_id) bind(C, name="warpx_push_bvec") integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), & @@ -542,6 +613,8 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN) :: dtsdx, dtsdy, dtsdz + real(amrex_real), intent(IN) :: xmin, dx + IF (maxwell_fdtd_solver_id .eq. 0) THEN ! Yee FDTD solver CALL WRPX_PXR_PUSH_BVEC( & @@ -552,7 +625,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n bx, bxlo, bxhi, & by, bylo, byhi, & bz, bzlo, bzhi, & - dtsdx,dtsdy,dtsdz) + dtsdx,dtsdy,dtsdz & +#ifdef WARPX_RZ + ,xmin,dx & +#endif + ) ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN ! Cole-Karkkainen FDTD solver CALL WRPX_PXR_PUSH_BVEC_CKC( & diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index a25d1141b..bf5f10bba 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -275,19 +275,21 @@ public: amrex::Real q_tot; long npart; + bool radially_weighted = true; + std::string str_density_function; std::string str_momentum_function_ux; std::string str_momentum_function_uy; std::string str_momentum_function_uz; -protected: - - amrex::Real mass, charge; - amrex::Real xmin, xmax; amrex::Real ymin, ymax; amrex::Real zmin, zmax; +protected: + + amrex::Real mass, charge; + amrex::Real density; int species_id; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 3c120a7c1..2adfcbad9 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -278,6 +278,9 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) StringParseAbortMessage("Injection style", part_pos_s); } + pp.query("radially_weighted", radially_weighted); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(radially_weighted, "ERROR: Only radially_weighted=true is supported"); + // parse plasma boundaries xmin = std::numeric_limits::lowest(); ymin = std::numeric_limits::lowest(); diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 75383e633..06142deb8 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -83,7 +83,8 @@ ifeq ($(USE_PSATD),TRUE) endif endif -ifeq ($(WARPX_RZ),TRUE) +ifeq ($(USE_RZ),TRUE) + USERSuffix += .RZ DEFINES += -DWARPX_RZ endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b534de916..bd8a99359 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -257,6 +257,9 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) if (!part_realbox.ok()) part_realbox = geom.ProbDomain(); int num_ppc = plasma_injector->num_particles_per_cell; +#ifdef WARPX_RZ + Real rmax = std::min(plasma_injector->xmax, part_realbox.hi(0)); +#endif const Real* dx = geom.CellSize(); @@ -424,7 +427,21 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab ); u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab ); } - attribs[PIdx::w ] = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); + Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); +#ifdef WARPX_RZ + if (plasma_injector->radially_weighted) { + weight *= 2*MathConst::pi*x; + } else { + // This is not correct since it might shift the particle + // out of the local grid + x = std::sqrt(x*rmax); + weight *= dx[0]; + } + Real theta = 2.*MathConst::pi*amrex::Random(); + y = x*std::sin(theta); + x = x*std::cos(theta); +#endif + attribs[PIdx::w ] = weight; attribs[PIdx::ux] = u[0]; attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; @@ -466,6 +483,9 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) if (!part_realbox.ok()) part_realbox = geom.ProbDomain(); int num_ppc = plasma_injector->num_particles_per_cell; +#ifdef WARPX_RZ + Real rmax = std::min(plasma_injector->xmax, part_realbox.hi(0)); +#endif const Real* dx = geom.CellSize(); @@ -636,7 +656,22 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab ); u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab ); } - attribs[PIdx::w ] = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); + Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); +#ifdef WARPX_RZ + Real r = x; + if (plasma_injector->radially_weighted) { + weight *= 2*MathConst::pi*r; + } else { + // This is not correct since it might shift the particle + // out of the local grid + x = std::sqrt(r*rmax); + weight *= dx[0]; + } + Real theta = 2.*MathConst::pi*amrex::Random(); + x = r*std::cos(theta); + y = r*std::sin(theta); +#endif + attribs[PIdx::w ] = weight; attribs[PIdx::ux] = u[0]; attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; @@ -659,6 +694,10 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) p.pos(1) = y; p.pos(2) = z; #elif (AMREX_SPACEDIM == 2) +#ifdef WARPX_RZ + attribs[PIdx::theta] = theta; + x = r; +#endif p.pos(0) = x; p.pos(1) = z; #endif diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index d11ce49c9..924980f81 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -11,7 +11,7 @@ struct PIdx enum { // Particle Attributes stored in amrex::ParticleContainer's struct of array w = 0, // weight ux, uy, uz, Ex, Ey, Ez, Bx, By, Bz, -#if (BL_SPACEDIM == 2) && WARPX_RZ +#ifdef WARPX_RZ theta, // RZ needs all three position components #endif #ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS @@ -188,11 +188,11 @@ public: void AddOneParticle (int lev, int grid, int tile, amrex::Real x, amrex::Real y, amrex::Real z, - const std::array& attribs); + std::array& attribs); void AddOneParticle (ParticleTileType& particle_tile, amrex::Real x, amrex::Real y, amrex::Real z, - const std::array& attribs); + std::array& attribs); void ReadHeader (std::istream& is); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 1ba82e3d3..e307fbaf8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -39,14 +39,17 @@ void WarpXParIter::SetPosition (const Cuda::DeviceVector& x, const Cuda::DeviceVector& y, const Cuda::DeviceVector& z) { #ifdef WARPX_RZ - const auto& attribs = GetAttribs(); - const auto& theta = attribs[PIdx::theta]; + auto& attribs = GetAttribs(); + auto& theta = attribs[PIdx::theta]; + Cuda::DeviceVector r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { theta[i] = std::atan2(y[i], x[i]); - x[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); + r[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); } -#endif + amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(r, z); +#else amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(x, z); +#endif } #endif @@ -117,7 +120,7 @@ WarpXParticleContainer::AllocData () void WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, Real x, Real y, Real z, - const std::array& attribs) + std::array& attribs) { auto& particle_tile = GetParticles(lev)[std::make_pair(grid,tile)]; AddOneParticle(particle_tile, x, y, z, attribs); @@ -126,7 +129,7 @@ WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, void WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, Real x, Real y, Real z, - const std::array& attribs) + std::array& attribs) { ParticleType p; p.id() = ParticleType::NextID(); @@ -223,7 +226,7 @@ WarpXParticleContainer::AddNParticles (int lev, { #ifdef WARPX_RZ if (comp == PIdx::theta) { - particle_tile.push_back_real(comp, theta.data, theta.data+np); + particle_tile.push_back_real(comp, theta.front(), theta.back()); } else { particle_tile.push_back_real(comp, np, 0.0); @@ -383,6 +386,14 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } BL_PROFILE_VAR_STOP(blp_pxr_cd); +#ifdef WARPX_RZ + warpx_current_deposition_rz_volume_scaling( + jx_ptr, &ngJ, jxntot.getVect(), + jy_ptr, &ngJ, jyntot.getVect(), + jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &xyzmin[0], &dx[0]); +#endif + BL_PROFILE_VAR_START(blp_accumulate); FArrayBox const* local_jx_const_ptr = local_jx[thread_num].get(); @@ -535,6 +546,14 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } BL_PROFILE_VAR_STOP(blp_pxr_cd); +#ifdef WARPX_RZ + warpx_current_deposition_rz_volume_scaling( + jx_ptr, &ngJ, jxntot.getVect(), + jy_ptr, &ngJ, jyntot.getVect(), + jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &xyzmin[0], &dx[0]); +#endif + BL_PROFILE_VAR_START(blp_accumulate); FArrayBox const* local_jx_const_ptr = local_jx[thread_num].get(); -- cgit v1.2.3 From 0a0100123f1c2dddff50949c80fed4e5650c334c Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 8 Apr 2019 13:51:38 -0700 Subject: Bug fixes in AddPlasma for RZ --- Source/Particles/PhysicalParticleContainer.cpp | 52 +++++++++++++++++--------- 1 file changed, 35 insertions(+), 17 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index bd8a99359..07468a85d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -385,6 +385,19 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) if(!tile_realbox.contains( RealVect{x, z} )) continue; #endif + // Save the x and y values to use in the insideBounds checks. + // This is needed with WARPX_RZ since x and y are modified. + Real xb = x; + Real yb = y; + +#ifdef WARPX_RZ + // Replace the x and y, choosing the angle randomly. + // These x and y are used to get the momentum and density + Real theta = 2.*MathConst::pi*amrex::Random(); + y = x*std::sin(theta); + x = x*std::cos(theta); +#endif + Real dens; std::array u; if (WarpX::gamma_boost == 1.){ @@ -392,7 +405,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) // If the particle is not within the species's // xmin, xmax, ymin, ymax, zmin, zmax, go to // the next generated particle. - if (!plasma_injector->insideBounds(x, y, z)) continue; + if (!plasma_injector->insideBounds(xb, yb, z)) continue; plasma_injector->getMomentum(u, x, y, z); dens = plasma_injector->getDensity(x, y, z); } else { @@ -419,7 +432,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) ); // If the particle is not within the lab-frame zmin, zmax, etc. // go to the next generated particle. - if (!plasma_injector->insideBounds(x, y, z0_lab)) continue; + if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue; // call `getDensity` with lab-frame parameters dens = plasma_injector->getDensity(x, y, z0_lab); // At this point u and dens are the lab-frame quantities @@ -430,16 +443,13 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); #ifdef WARPX_RZ if (plasma_injector->radially_weighted) { - weight *= 2*MathConst::pi*x; + weight *= 2*MathConst::pi*xb; } else { // This is not correct since it might shift the particle // out of the local grid - x = std::sqrt(x*rmax); + x = std::sqrt(xb*rmax); weight *= dx[0]; } - Real theta = 2.*MathConst::pi*amrex::Random(); - y = x*std::sin(theta); - x = x*std::cos(theta); #endif attribs[PIdx::w ] = weight; attribs[PIdx::ux] = u[0]; @@ -614,6 +624,19 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) if(!tile_realbox.contains( RealVect{x, z} )) continue; #endif + // Save the x and y values to use in the insideBounds checks. + // This is needed with WARPX_RZ since x and y are modified. + Real xb = x; + Real yb = y; + +#ifdef WARPX_RZ + // Replace the x and y, choosing the angle randomly. + // These x and y are used to get the momentum and density + Real theta = 2.*MathConst::pi*amrex::Random(); + x = xb*std::cos(theta); + y = xb*std::sin(theta); +#endif + Real dens; std::array u; if (WarpX::gamma_boost == 1.){ @@ -621,7 +644,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) // If the particle is not within the species's // xmin, xmax, ymin, ymax, zmin, zmax, go to // the next generated particle. - if (!plasma_injector->insideBounds(x, y, z)) continue; + if (!plasma_injector->insideBounds(xb, yb, z)) continue; plasma_injector->getMomentum(u, x, y, z); dens = plasma_injector->getDensity(x, y, z); } else { @@ -648,7 +671,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) ); // If the particle is not within the lab-frame zmin, zmax, etc. // go to the next generated particle. - if (!plasma_injector->insideBounds(x, y, z0_lab)) continue; + if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue; // call `getDensity` with lab-frame parameters dens = plasma_injector->getDensity(x, y, z0_lab); // At this point u and dens are the lab-frame quantities @@ -658,18 +681,14 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) } Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); #ifdef WARPX_RZ - Real r = x; if (plasma_injector->radially_weighted) { - weight *= 2*MathConst::pi*r; + weight *= 2*MathConst::pi*xb; } else { // This is not correct since it might shift the particle // out of the local grid - x = std::sqrt(r*rmax); + x = std::sqrt(xb*rmax); weight *= dx[0]; } - Real theta = 2.*MathConst::pi*amrex::Random(); - x = r*std::cos(theta); - y = r*std::sin(theta); #endif attribs[PIdx::w ] = weight; attribs[PIdx::ux] = u[0]; @@ -696,9 +715,8 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) #elif (AMREX_SPACEDIM == 2) #ifdef WARPX_RZ attribs[PIdx::theta] = theta; - x = r; #endif - p.pos(0) = x; + p.pos(0) = xb; p.pos(1) = z; #endif -- cgit v1.2.3 From 2e8dda19ad45656a6e928bdc1a2fc5fd0a3fa5a9 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 11 Apr 2019 18:30:24 -0400 Subject: Switch over to using Array4 in PPC::Evolve, which didn't exist when we first ported it --- Source/Laser/LaserParticleContainer.cpp | 5 - Source/Particles/PhysicalParticleContainer.cpp | 4 - Source/Particles/WarpXParticleContainer.H | 10 +- Source/Particles/WarpXParticleContainer.cpp | 197 +++++++++++++------------ 4 files changed, 110 insertions(+), 106 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 3ef1be154..db5499b8e 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -342,11 +342,6 @@ LaserParticleContainer::Evolve (int lev, int thread_num = 0; #endif - if (local_rho[thread_num] == nullptr) local_rho[thread_num].reset( new amrex::FArrayBox()); - if (local_jx[thread_num] == nullptr) local_jx[thread_num].reset( new amrex::FArrayBox()); - if (local_jy[thread_num] == nullptr) local_jy[thread_num].reset( new amrex::FArrayBox()); - if (local_jz[thread_num] == nullptr) local_jz[thread_num].reset( new amrex::FArrayBox()); - Cuda::ManagedDeviceVector plane_Xp, plane_Yp, amplitude_E; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d98e79177..1031b488f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1048,10 +1048,6 @@ PhysicalParticleContainer::Evolve (int lev, #else int thread_num = 0; #endif - if (local_rho[thread_num] == nullptr) local_rho[thread_num].reset( new amrex::FArrayBox()); - if (local_jx[thread_num] == nullptr) local_jx[thread_num].reset( new amrex::FArrayBox()); - if (local_jy[thread_num] == nullptr) local_jy[thread_num].reset( new amrex::FArrayBox()); - if (local_jz[thread_num] == nullptr) local_jz[thread_num].reset( new amrex::FArrayBox()); FArrayBox filtered_Ex, filtered_Ey, filtered_Ez; FArrayBox filtered_Bx, filtered_By, filtered_Bz; diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 3b3fd81f4..6ac2ca621 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -215,12 +215,12 @@ protected: static int do_not_push; - amrex::Vector > local_rho; - amrex::Vector > local_jx; - amrex::Vector > local_jy; - amrex::Vector > local_jz; + amrex::Vector local_rho; + amrex::Vector local_jx; + amrex::Vector local_jy; + amrex::Vector local_jz; - amrex::Vector > m_xp, m_yp, m_zp, m_giv; + amrex::Vector > m_xp, m_yp, m_zp, m_giv; private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index d54dd261d..53ac9d3ff 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -56,14 +56,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) m_yp.resize(num_threads); m_zp.resize(num_threads); m_giv.resize(num_threads); - for (int i = 0; i < num_threads; ++i) - { - local_rho[i].reset(nullptr); - local_jx[i].reset(nullptr); - local_jy[i].reset(nullptr); - local_jz[i].reset(nullptr); - } - } void @@ -232,35 +224,38 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, tby.grow(ngJ); tbz.grow(ngJ); - local_jx[thread_num]->resize(tbx); - local_jy[thread_num]->resize(tby); - local_jz[thread_num]->resize(tbz); + local_jx[thread_num].resize(tbx); + local_jy[thread_num].resize(tby); + local_jz[thread_num].resize(tbz); - jx_ptr = local_jx[thread_num]->dataPtr(); - jy_ptr = local_jy[thread_num]->dataPtr(); - jz_ptr = local_jz[thread_num]->dataPtr(); + jx_ptr = local_jx[thread_num].dataPtr(); + jy_ptr = local_jy[thread_num].dataPtr(); + jz_ptr = local_jz[thread_num].dataPtr(); - FArrayBox* local_jx_ptr = local_jx[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, b, + auto jxarr = local_jx[thread_num].array(); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_jx_ptr->setVal(0.0, b, 0, 1); + jxarr(i,j,k) = 0.0; }); - FArrayBox* local_jy_ptr = local_jy[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, b, + auto jyarr = local_jy[thread_num].array(); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_jy_ptr->setVal(0.0, b, 0, 1); + jyarr(i,j,k) = 0.0; }); - FArrayBox* local_jz_ptr = local_jz[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, b, + auto jzarr = local_jz[thread_num].array(); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_jz_ptr->setVal(0.0, b, 0, 1); + jzarr(i,j,k) = 0.0; }); - auto jxntot = local_jx[thread_num]->length(); - auto jyntot = local_jy[thread_num]->length(); - auto jzntot = local_jz[thread_num]->length(); + auto jxntot = local_jx[thread_num].length(); + auto jyntot = local_jy[thread_num].length(); + auto jzntot = local_jz[thread_num].length(); BL_PROFILE_VAR_START(blp_pxr_cd); if (j_is_nodal) { @@ -341,28 +336,31 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, &lvect,&WarpX::current_deposition_algo); } BL_PROFILE_VAR_STOP(blp_pxr_cd); - + BL_PROFILE_VAR_START(blp_accumulate); - - FArrayBox const* local_jx_const_ptr = local_jx[thread_num].get(); - FArrayBox* global_jx_ptr = jx.fabPtr(pti); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, thread_bx, + + const auto local_jx_arr = local_jx[thread_num].array(); + auto global_jx_arr = jx.array(pti); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_jx_ptr->atomicAdd(*local_jx_const_ptr, thread_bx, thread_bx, 0, 0, 1); + Gpu::Atomic::Add(&global_jx_arr(i, j, k), local_jx_arr(i, j, k)); }); - FArrayBox const* local_jy_const_ptr = local_jy[thread_num].get(); - FArrayBox* global_jy_ptr = jy.fabPtr(pti); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, thread_bx, + const auto local_jy_arr = local_jy[thread_num].array(); + auto global_jy_arr = jy.array(pti); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_jy_ptr->atomicAdd(*local_jy_const_ptr, thread_bx, thread_bx, 0, 0, 1); + Gpu::Atomic::Add(&global_jy_arr(i, j, k), local_jy_arr(i, j, k)); }); - FArrayBox const* local_jz_const_ptr = local_jz[thread_num].get(); - FArrayBox* global_jz_ptr = jz.fabPtr(pti); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, thread_bx, + const auto local_jz_arr = local_jz[thread_num].array(); + auto global_jz_arr = jz.array(pti); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_jz_ptr->atomicAdd(*local_jz_const_ptr, thread_bx, thread_bx, 0, 0, 1); + Gpu::Atomic::Add(&global_jz_arr(i, j, k), local_jz_arr(i, j, k)); }); BL_PROFILE_VAR_STOP(blp_accumulate); @@ -382,34 +380,38 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, tby.grow(ngJ); tbz.grow(ngJ); - local_jx[thread_num]->resize(tbx); - local_jy[thread_num]->resize(tby); - local_jz[thread_num]->resize(tbz); + local_jx[thread_num].resize(tbx); + local_jy[thread_num].resize(tby); + local_jz[thread_num].resize(tbz); - jx_ptr = local_jx[thread_num]->dataPtr(); - jy_ptr = local_jy[thread_num]->dataPtr(); - jz_ptr = local_jz[thread_num]->dataPtr(); + jx_ptr = local_jx[thread_num].dataPtr(); + jy_ptr = local_jy[thread_num].dataPtr(); + jz_ptr = local_jz[thread_num].dataPtr(); - FArrayBox* local_jx_ptr = local_jx[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, b, + auto jxarr = local_jx[thread_num].array(); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_jx_ptr->setVal(0.0, b, 0, 1); + jxarr(i,j,k) = 0.0; }); - FArrayBox* local_jy_ptr = local_jy[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, b, + auto jyarr = local_jy[thread_num].array(); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_jy_ptr->setVal(0.0, b, 0, 1); + jyarr(i,j,k) = 0.0; }); - FArrayBox* local_jz_ptr = local_jz[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, b, + auto jzarr = local_jz[thread_num].array(); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_jz_ptr->setVal(0.0, b, 0, 1); + jzarr(i,j,k) = 0.0; }); - auto jxntot = local_jx[thread_num]->length(); - auto jyntot = local_jy[thread_num]->length(); - auto jzntot = local_jz[thread_num]->length(); + + auto jxntot = local_jx[thread_num].length(); + auto jyntot = local_jy[thread_num].length(); + auto jzntot = local_jz[thread_num].length(); long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); @@ -496,25 +498,28 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, BL_PROFILE_VAR_START(blp_accumulate); - FArrayBox const* local_jx_const_ptr = local_jx[thread_num].get(); - FArrayBox* global_jx_ptr = cjx->fabPtr(pti); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, thread_bx, + const auto local_jx_arr = local_jx[thread_num].array(); + auto global_jx_arr = cjx->array(pti); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_jx_ptr->atomicAdd(*local_jx_const_ptr, thread_bx, thread_bx, 0, 0, 1); + Gpu::Atomic::Add(&global_jx_arr(i, j, k), local_jx_arr(i, j, k)); }); - FArrayBox const* local_jy_const_ptr = local_jy[thread_num].get(); - FArrayBox* global_jy_ptr = cjy->fabPtr(pti); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, thread_bx, + const auto local_jy_arr = local_jy[thread_num].array(); + auto global_jy_arr = cjy->array(pti); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_jy_ptr->atomicAdd(*local_jy_const_ptr, thread_bx, thread_bx, 0, 0, 1); + Gpu::Atomic::Add(&global_jy_arr(i, j, k), local_jy_arr(i, j, k)); }); - FArrayBox const* local_jz_const_ptr = local_jz[thread_num].get(); - FArrayBox* global_jz_ptr = cjz->fabPtr(pti); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, thread_bx, + const auto local_jz_arr = local_jz[thread_num].array(); + auto global_jz_arr = cjz->array(pti); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_jz_ptr->atomicAdd(*local_jz_const_ptr, thread_bx, thread_bx, 0, 0, 1); + Gpu::Atomic::Add(&global_jz_arr(i, j, k), local_jz_arr(i, j, k)); }); BL_PROFILE_VAR_STOP(blp_accumulate); @@ -547,15 +552,17 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, { const std::array& xyzmin = xyzmin_tile; tile_box.grow(ngRho); - local_rho[thread_num]->resize(tile_box); - FArrayBox* local_rho_ptr = local_rho[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, b, + local_rho[thread_num].resize(tile_box); + + auto rhoarr = local_rho[thread_num].array(); + amrex::ParallelFor(tile_box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_rho_ptr->setVal(0.0, b, 0, 1); + rhoarr(i,j,k) = 0.0; }); - data_ptr = local_rho[thread_num]->dataPtr(); - auto rholen = local_rho[thread_num]->length(); + data_ptr = local_rho[thread_num].dataPtr(); + auto rholen = local_rho[thread_num].length(); #if (AMREX_SPACEDIM == 3) const long nx = rholen[0]-1-2*ngRho; const long ny = rholen[1]-1-2*ngRho; @@ -579,14 +586,16 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &lvect, &WarpX::charge_deposition_algo); BL_PROFILE_VAR_STOP(blp_pxr_chd); - const int ncomp = 1; - FArrayBox const* local_fab = local_rho[thread_num].get(); - FArrayBox* global_fab = rhomf->fabPtr(pti); BL_PROFILE_VAR_START(blp_accumulate); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, tbx, + + const auto local_rho_arr = local_rho[thread_num].array(); + auto global_rho_arr = rhomf->array(pti); + amrex::ParallelFor(tile_box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_fab->atomicAdd(*local_fab, tbx, tbx, 0, icomp, ncomp); + Gpu::Atomic::Add(&global_rho_arr(i, j, k), local_rho_arr(i, j, k)); }); + BL_PROFILE_VAR_STOP(blp_accumulate); } @@ -599,15 +608,17 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); tile_box.grow(ngRho); - local_rho[thread_num]->resize(tile_box); - FArrayBox* local_rho_ptr = local_rho[thread_num].get(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, b, + local_rho[thread_num].resize(tile_box); + + auto rhoarr = local_rho[thread_num].array(); + amrex::ParallelFor(tile_box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - local_rho_ptr->setVal(0.0, b, 0, 1); + rhoarr(i,j,k) = 0.0; }); - data_ptr = local_rho[thread_num]->dataPtr(); - auto rholen = local_rho[thread_num]->length(); + data_ptr = local_rho[thread_num].dataPtr(); + auto rholen = local_rho[thread_num].length(); #if (AMREX_SPACEDIM == 3) const long nx = rholen[0]-1-2*ngRho; const long ny = rholen[1]-1-2*ngRho; @@ -633,14 +644,16 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &lvect, &WarpX::charge_deposition_algo); BL_PROFILE_VAR_STOP(blp_pxr_chd); - const int ncomp = 1; - FArrayBox const* local_fab = local_rho[thread_num].get(); - FArrayBox* global_fab = crhomf->fabPtr(pti); BL_PROFILE_VAR_START(blp_accumulate); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, tbx, + + const auto local_rho_arr = local_rho[thread_num].array(); + auto global_rho_arr = crhomf->array(pti); + amrex::ParallelFor(tile_box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - global_fab->atomicAdd(*local_fab, tbx, tbx, 0, icomp, ncomp); + Gpu::Atomic::Add(&global_rho_arr(i, j, k), local_rho_arr(i, j, k)); }); + BL_PROFILE_VAR_STOP(blp_accumulate); } }; -- cgit v1.2.3 From d6d9c796040e14dc0fe50e96bc2b1231d953a01e Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 29 Apr 2019 10:16:37 -0700 Subject: remove the compilation flag for STORE_OLD_PARTICLE_ATTRIBS - this can be done at runtime now --- GNUmakefile | 2 -- Source/FortranInterface/WarpX_f.H | 11 +++++------ Source/Make.WarpX | 4 ---- Source/Particles/PhysicalParticleContainer.cpp | 17 ++++++----------- Tools/performance_tests/GNUmakefile_perftest | 1 - 5 files changed, 11 insertions(+), 24 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/GNUmakefile b/GNUmakefile index be8c80650..1acd53be7 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -17,8 +17,6 @@ TINY_PROFILE = TRUE #COMM_PROFILE = TRUE #TRACE_PROFILE = TRUE -STORE_OLD_PARTICLE_ATTRIBS = FALSE - USE_OMP = TRUE USE_GPU = FALSE diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 473ca645a..f246f9f54 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -81,17 +81,16 @@ extern "C" { #endif -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS void warpx_copy_attribs(const long* np, const amrex_real* xp, const amrex_real* yp, const amrex_real* zp, const amrex_real* uxp, const amrex_real* uyp, const amrex_real* uzp, amrex_real* xpold, amrex_real* ypold, amrex_real* zpold, amrex_real* uxpold, amrex_real* uypold, amrex_real* uzpold); -#endif - void WRPX_COPY_SLICE(const int* lo, const int* hi, - const amrex_real* tmp, const int* tlo, const int* thi, - amrex_real* buf, const int* blo, const int* bhi, - const int* ncomp, const int* i_boost, const int* i_lab); + + void WRPX_COPY_SLICE(const int* lo, const int* hi, + const amrex_real* tmp, const int* tlo, const int* thi, + amrex_real* buf, const int* blo, const int* bhi, + const int* ncomp, const int* i_boost, const int* i_lab); // Charge deposition void warpx_charge_deposition(amrex::Real* rho, diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 926523d82..95e2b4ec0 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -123,10 +123,6 @@ ifeq ($(USE_RZ),TRUE) DEFINES += -DWARPX_RZ endif -ifeq ($(STORE_OLD_PARTICLE_ATTRIBS),TRUE) - DEFINES += -DWARPX_STORE_OLD_PARTICLE_ATTRIBS -endif - ifeq ($(DO_ELECTROSTATIC),TRUE) include $(AMREX_HOME)/Src/LinearSolvers/C_to_F_MG/Make.package include $(AMREX_HOME)/Src/LinearSolvers/F_MG/FParallelMG.mak diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index e31d43204..74052f842 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1801,7 +1801,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real { BL_PROFILE("PhysicalParticleContainer::GetParticleSlice"); -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS // Assume that the boost in the positive z direction. #if (AMREX_SPACEDIM == 2) AMREX_ALWAYS_ASSERT(direction == 1); @@ -1864,12 +1863,12 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = attribs[PIdx::xold ]; - auto& yp_old = attribs[PIdx::yold ]; - auto& zp_old = attribs[PIdx::zold ]; - auto& uxp_old = attribs[PIdx::uxold]; - auto& uyp_old = attribs[PIdx::uyold]; - auto& uzp_old = attribs[PIdx::uzold]; + auto& xp_old = pti.GetAttribs(particle_comps["xold"]); + auto& yp_old = pti.GetAttribs(particle_comps["yold"]); + auto& zp_old = pti.GetAttribs(particle_comps["zold"]); + auto& uxp_old = pti.GetAttribs(particle_comps["uxold"]); + auto& uyp_old = pti.GetAttribs(particle_comps["uyold"]); + auto& uzp_old = pti.GetAttribs(particle_comps["uzold"]); const long np = pti.numParticles(); @@ -1919,10 +1918,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real } } } -#else - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false , -"ERROR: WarpX must be compiled with STORE_OLD_PARTICLE_ATTRIBS=TRUE to use the back-transformed diagnostics"); -#endif } int PhysicalParticleContainer::GetRefineFac(const Real x, const Real y, const Real z) diff --git a/Tools/performance_tests/GNUmakefile_perftest b/Tools/performance_tests/GNUmakefile_perftest index bbb2cc659..7740b0a64 100644 --- a/Tools/performance_tests/GNUmakefile_perftest +++ b/Tools/performance_tests/GNUmakefile_perftest @@ -6,7 +6,6 @@ DEBUG = FALSE DIM = 3 COMP=intel TINY_PROFILE = TRUE -STORE_OLD_PARTICLE_ATTRIBS=TRUE USE_OMP = TRUE EBASE = perf_tests USE_PYTHON_MAIN = FALSE -- cgit v1.2.3 From 7e58007c9eba2b51003ff813a8bfe65ae25e8392 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 29 Apr 2019 10:41:06 -0700 Subject: replace the compile time checks used for the old particle attributes with runtime ones --- Source/Particles/PhysicalParticleContainer.cpp | 71 ++++++++++++---------- .../Particles/RigidInjectedParticleContainer.cpp | 28 ++++----- 2 files changed, 52 insertions(+), 47 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 74052f842..ee219e355 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -455,16 +455,18 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::ux] = u[0]; attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - attribs[PIdx::xold] = x; - attribs[PIdx::yold] = y; - attribs[PIdx::zold] = z; - - attribs[PIdx::uxold] = u[0]; - attribs[PIdx::uyold] = u[1]; - attribs[PIdx::uzold] = u[2]; -#endif + + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); } @@ -695,15 +697,18 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - attribs[PIdx::xold] = x; - attribs[PIdx::yold] = y; - attribs[PIdx::zold] = z; - - attribs[PIdx::uxold] = u[0]; - attribs[PIdx::uyold] = u[1]; - attribs[PIdx::uzold] = u[2]; -#endif + // note - this will be slow on the GPU, need to revisit + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } ParticleType p; p.id() = ParticleType::NextID(); @@ -1666,20 +1671,20 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - auto& xpold = attribs[PIdx::xold]; - auto& ypold = attribs[PIdx::yold]; - auto& zpold = attribs[PIdx::zold]; - auto& uxpold = attribs[PIdx::uxold]; - auto& uypold = attribs[PIdx::uyold]; - auto& uzpold = attribs[PIdx::uzold]; - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - -#endif + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& xpold = pti.GetAttribs(particle_comps["xold"]); + auto& ypold = pti.GetAttribs(particle_comps["yold"]); + auto& zpold = pti.GetAttribs(particle_comps["zold"]); + auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); + auto& uypold = pti.GetAttribs(particle_comps["uyold"]); + auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); + + warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), + uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); + } warpx_particle_pusher(&np, xp.dataPtr(), diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 3ee4d87e5..a5acca281 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -225,20 +225,20 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - auto& xpold = attribs[PIdx::xold]; - auto& ypold = attribs[PIdx::yold]; - auto& zpold = attribs[PIdx::zold]; - auto& uxpold = attribs[PIdx::uxold]; - auto& uypold = attribs[PIdx::uyold]; - auto& uzpold = attribs[PIdx::uzold]; - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - -#endif + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& xpold = pti.GetAttribs(particle_comps["xold"]); + auto& ypold = pti.GetAttribs(particle_comps["yold"]); + auto& zpold = pti.GetAttribs(particle_comps["zold"]); + auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); + auto& uypold = pti.GetAttribs(particle_comps["uyold"]); + auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); + + warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), + uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); + } // Save the position and momenta, making copies Cuda::ManagedDeviceVector xp_save, yp_save, zp_save; -- cgit v1.2.3 From 9fec75fc1dcf8941a7a5012c50b2bd727ccfcdcf Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 29 Apr 2019 12:45:06 -0700 Subject: make sure we initialize the new particle components for all the different paths for particle initialization --- Source/Particles/PhysicalParticleContainer.cpp | 12 ++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 17 +++++++++++++++++ 2 files changed, 29 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ee219e355..17e6d98d9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -184,6 +184,18 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } + AddOneParticle(0, 0, 0, x, y, z, attribs); } } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index d80074af4..2edd3c636 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -230,6 +230,15 @@ WarpXParticleContainer::AddNParticles (int lev, #endif p.pos(1) = z[i]; #endif + + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x[i]); + particle_tile.push_back_real(particle_comps["yold"], y[i]); + particle_tile.push_back_real(particle_comps["zold"], z[i]); + } + particle_tile.push_back(p); } @@ -240,6 +249,14 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + } + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_RZ -- cgit v1.2.3