diff options
26 files changed, 360 insertions, 375 deletions
diff --git a/.gitignore b/.gitignore index 13147b8da..1e9fbc336 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,3 @@ -.DS_Store -plt* -chk* main?d.*.ex d/ f/ @@ -11,3 +8,32 @@ Backtrace.* *.vth *.sw[opqrs] +#################### +# AMReX data files # +#################### +plt* +chk* + +########## +# Python # +########## +*.pyc +__pycache__ + +####### +# IDE # +####### +.idea/ +cmake-build-*/ +.kdev?/ +*.kdev? + +# File-based project format: +*.iws + +###### +# OS # +###### +.DS_Store +.AppleDouble +.LSOverride diff --git a/Docs/source/conf.py b/Docs/source/conf.py index 54e533469..7f8ad0785 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -58,7 +58,7 @@ author = 'WarpX collaboration' # built documents. # # The short X.Y version. -version = '19.05' +version = '19.08' # The full version, including alpha/beta/rc tags. release = '' diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 518aae76a..e0d115e7f 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -668,6 +668,12 @@ Boundary conditions The characteristic depth, in number of cells, over which the absorption coefficients of the PML increases. +* ``warpx.do_pml_Lo`` (`2 ints in 2D`, `3 ints in 3D`; default: `1 1 1`) + The directions along which one wants a pml boundary condition for lower boundaries on mother grid. + +* ``warpx.do_pml_Hi`` (`2 floats in 2D`, `3 floats in 3D`; default: `1 1 1`) + The directions along which one wants a pml boundary condition for upper boundaries on mother grid. + Diagnostics and output ---------------------- diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost index a6d45426a..d90c75ada 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost +++ b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost @@ -19,7 +19,7 @@ geometry.prob_hi = 128.e-6 0.96e-6 ################################# warpx.verbose = 1 amrex.v = 1 -algo.current_deposition = direct +algo.current_deposition = esirkepov algo.charge_deposition = standard algo.field_gathering = standard algo.particle_pusher = vay diff --git a/LICENSE.txt b/LICENSE.txt index aeeded3a7..346d92dc3 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,4 +1,4 @@ -WarpX v19.05 Copyright (c) 2018, The Regents of the University of California, through Lawrence Berkeley National Laboratory, and Lawrence Livermore National Security, LLC, for the operation of Lawrence Livermore National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. +WarpX v19.08 Copyright (c) 2018, The Regents of the University of California, through Lawrence Berkeley National Laboratory, and Lawrence Livermore National Security, LLC, for the operation of Lawrence Livermore National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: diff --git a/Python/setup.py b/Python/setup.py index 8429cf95e..76d99fbaa 100644 --- a/Python/setup.py +++ b/Python/setup.py @@ -20,7 +20,7 @@ else: package_data = {} setup (name = 'pywarpx', - version = '19.05', + version = '19.08', packages = ['pywarpx'], package_dir = {'pywarpx':'pywarpx'}, description = """Wrapper of WarpX""", diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index c9ccb910a..23b9c0acf 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -40,7 +40,7 @@ goUpLink = 1 # email sendEmailWhenFail = 1 -emailTo = weiqunzhang@lbl.gov, jlvay@lbl.gov, rlehe@lbl.gov, atmyers@lbl.gov, mthevenet@lbl.gov, jaehongpark@lbl.gov, oshapoval@lbl.gov, henri.vincenti@cea.fr, ldianaamorim@lbl.gov, rjambunathan@lbl.gov +emailTo = weiqunzhang@lbl.gov, jlvay@lbl.gov, rlehe@lbl.gov, atmyers@lbl.gov, mthevenet@lbl.gov, jaehongpark@lbl.gov, oshapoval@lbl.gov, henri.vincenti@cea.fr, ldianaamorim@lbl.gov, rjambunathan@lbl.gov, axelhuebl@lbl.gov emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more details. [AMReX] 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"); |