diff options
Diffstat (limited to 'Source')
25 files changed, 695 insertions, 207 deletions
diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 2a9c16aa8..d55660b39 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -112,6 +112,12 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const int_flags.resize(1, 1); } +#ifdef WARPX_QED + if(pc->do_qed){ + real_names.push_back("tau"); + } +#endif + // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index a517a91f5..7a3262703 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -562,6 +562,18 @@ WarpX::ComputeDt () if (do_electrostatic) { dt[0] = const_dt; } + + for (int lev=0; lev <= max_level; lev++) { + const Real* dx_lev = geom[lev].CellSize(); + Print()<<"Level "<<lev<<": dt = "<<dt[lev] + <<" ; dx = "<<dx_lev[0] +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + <<" ; dz = "<<dx_lev[1]<<'\n'; +#elif (defined WARPX_DIM_3D) + <<" ; dy = "<<dx_lev[1] + <<" ; dz = "<<dx_lev[2]<<'\n'; +#endif + } } /* \brief computes max_step for wakefield simulation in boosted frame. @@ -578,13 +590,14 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ WarpX::moving_window_dir == AMREX_SPACEDIM-1, "Can use zmax_plasma_to_compute_max_step only if " + "moving window along z. TODO: all directions."); - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + - (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + - (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, - "Can use zmax_plasma_to_compute_max_step only if " + - "warpx.boost_direction = z. TODO: all directions."); + if (gamma_boost > 1){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + + (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + + (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, + "Can use zmax_plasma_to_compute_max_step in boosted frame only if " + + "warpx.boost_direction = z. TODO: all directions."); + } // Lower end of the simulation domain. All quantities are given in boosted // frame except zmax_plasma_to_compute_max_step. diff --git a/Source/FortranInterface/WarpX_f.F90 b/Source/FortranInterface/WarpX_f.F90 index 2ddf9409d..762ed2cd8 100644 --- a/Source/FortranInterface/WarpX_f.F90 +++ b/Source/FortranInterface/WarpX_f.F90 @@ -51,65 +51,6 @@ contains end subroutine warpx_compute_E - subroutine warpx_sync_rho_2d (lo, hi, crse, clo, chi, fine, flo, fhi, nc) & - bind(c, name='warpx_sync_rho_2d') - integer, intent(in) :: lo(2), hi(2), flo(2), fhi(2), clo(2), chi(2), nc - real(amrex_real), intent(in ) :: fine(flo(1):fhi(1),flo(2):fhi(2),nc) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),nc) - - integer :: i,j,ii,jj,m - - do m = 1, nc - do j = lo(2), hi(2) - jj = j*2 - do i = lo(1), hi(1) - ii = i*2 - crse(i,j,m) = 0.25d0 * & - ( fine(ii,jj,m) + 0.5d0 *(fine(ii-1,jj ,m)+fine(ii+1,jj ,m) & - & + fine(ii ,jj-1,m)+fine(ii ,jj+1,m)) & - & + 0.25d0*(fine(ii-1,jj-1,m)+fine(ii+1,jj-1,m) & - & + fine(ii-1,jj+1,m)+fine(ii+1,jj+1,m)) ) - end do - end do - end do - end subroutine warpx_sync_rho_2d - - subroutine warpx_sync_rho_3d (lo, hi, crse, clo, chi, fine, flo, fhi, nc) & - bind(c, name='warpx_sync_rho_3d') - integer, intent(in) :: lo(3), hi(3), flo(3), fhi(3), clo(3), chi(3), nc - real(amrex_real), intent(in ) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3),nc) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3),nc) - - integer :: i,j,k,ii,jj,kk,m - - do m = 1, nc - do k = lo(3), hi(3) - kk = k*2 - do j = lo(2), hi(2) - jj = j*2 - do i = lo(1), hi(1) - ii = i*2 - crse(i,j,k,m) = 0.125d0 * & - (fine(ii,jj,kk,m) + 0.5d0 *(fine(ii-1,jj ,kk ,m)+fine(ii+1,jj ,kk ,m) & - & + fine(ii ,jj-1,kk ,m)+fine(ii ,jj+1,kk ,m) & - & + fine(ii ,jj ,kk-1,m)+fine(ii ,jj ,kk+1,m)) & - & + 0.25d0 *(fine(ii-1,jj-1,kk ,m)+fine(ii+1,jj-1,kk ,m) & - & + fine(ii-1,jj+1,kk ,m)+fine(ii+1,jj+1,kk ,m) & - & + fine(ii-1,jj ,kk-1,m)+fine(ii+1,jj ,kk-1,m) & - & + fine(ii-1,jj ,kk+1,m)+fine(ii+1,jj ,kk+1,m) & - & + fine(ii ,jj-1,kk-1,m)+fine(ii ,jj+1,kk-1,m) & - & + fine(ii ,jj-1,kk+1,m)+fine(ii ,jj+1,kk+1,m)) & - & + 0.125d0*(fine(ii-1,jj-1,kk-1,m)+fine(ii-1,jj-1,kk+1,m) & - & + fine(ii-1,jj+1,kk-1,m)+fine(ii-1,jj+1,kk+1,m) & - & + fine(ii+1,jj-1,kk-1,m)+fine(ii+1,jj-1,kk+1,m) & - & + fine(ii+1,jj+1,kk-1,m)+fine(ii+1,jj+1,kk+1,m))) - end do - end do - end do - end do - - end subroutine warpx_sync_rho_3d - subroutine warpx_build_buffer_masks (lo, hi, msk, mlo, mhi, gmsk, glo, ghi, ng) & bind(c, name='warpx_build_buffer_masks') integer, dimension(3), intent(in) :: lo, hi, mlo, mhi, glo, ghi diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 2d5d346e2..b51a83387 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -16,8 +16,6 @@ #if (AMREX_SPACEDIM == 3) -#define WRPX_SYNC_RHO warpx_sync_rho_3d - #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_3d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_3d #define WRPX_BUILD_MASK warpx_build_mask_3d @@ -30,8 +28,6 @@ #elif (AMREX_SPACEDIM == 2) -#define WRPX_SYNC_RHO warpx_sync_rho_2d - #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_2d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_2d #define WRPX_BUILD_MASK warpx_build_mask_2d @@ -229,11 +225,6 @@ extern "C" #endif const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); - void WRPX_SYNC_RHO (const int* lo, const int* hi, - BL_FORT_FAB_ARG_ANYD(crse), - const BL_FORT_FAB_ARG_ANYD(fine), - const int* ncomp); - #ifdef WARPX_USE_PSATD void warpx_fft_mpi_init (int fcomm); void warpx_fft_domain_decomp (int* warpx_local_nz, int* warpx_local_z0, diff --git a/Source/Make.WarpX b/Source/Make.WarpX index afd70ad46..54baecbf6 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -4,6 +4,7 @@ OPENBC_HOME ?= ../../../openbc_poisson USE_MPI = TRUE USE_PARTICLES = TRUE +USE_RPATH = TRUE ifeq ($(USE_GPU),TRUE) USE_OMP = FALSE @@ -27,24 +28,24 @@ ifeq ($(USE_ASCENT_INSITU),TRUE) USE_ASCENT = TRUE endif -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - ifndef USE_PYTHON_MAIN USE_PYTHON_MAIN = FALSE endif ifeq ($(USE_PYTHON_MAIN),TRUE) - CXXFLAGS += -fPIC + XTRA_CXXFLAGS += -fPIC ifeq ($(USE_OMP),TRUE) LDFLAGS += -fopenmp endif - CFLAGS += -fPIC - FFLAGS += -fPIC - F90FLAGS += -fPIC + XTRA_CFLAGS += -fPIC + XTRA_FFLAGS += -fPIC + XTRA_F90FLAGS += -fPIC USERSuffix := .Lib DEFINES += -DWARPX_USE_PY endif +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + ifeq ($(DIM),3) DEFINES += -DWARPX_DIM_3D else @@ -93,6 +94,7 @@ ifeq ($(QED),TRUE) FFLAGS += -DWARPX_QED F90FLAGS += -DWARPX_QED include $(WARPX_HOME)/Source/QED/Make.package + USERSuffix := $(USERSuffix).QED endif @@ -115,8 +117,8 @@ ifeq ($(USE_OPENPMD), TRUE) # try pkg-config query ifeq (0, $(shell pkg-config "openPMD >= 0.9.0"; echo $$?)) CXXFLAGS += $(shell pkg-config --cflags openPMD) - LDFLAGS += $(shell pkg-config --libs openPMD) - LDFLAGS += -Xlinker -rpath -Xlinker $(shell pkg-config --variable=libdir openPMD) + LIBRARY_LOCATIONS += $(shell pkg-config --variable=libdir openPMD) + libraries += $(shell pkg-config --libs-only-l openPMD) # fallback to manual settings else OPENPMD_LIB_PATH ?= NOT_SET @@ -209,7 +211,11 @@ libwarpx$(PYDIM).a: $(objForExecs) libwarpx$(PYDIM).so: $(objForExecs) @echo Making dynamic library $@ ... +ifeq ($(USE_CUDA),TRUE) + $(SILENT) $(CXX) $(LINKFLAGS) $(SHARED_OPTION) -Xlinker=$(SHARED_OPTION) -Xlinker=-fPIC -o $@ $^ $(LDFLAGS) $(libraries) +else $(SILENT) $(CXX) $(SHARED_OPTION) -fPIC -o $@ $^ $(LDFLAGS) $(libraries) +endif $(SILENT) $(RM) AMReX_buildInfo.cpp @echo SUCCESS diff --git a/Source/Parallelization/InterpolateDensityFineToCoarse.H b/Source/Parallelization/InterpolateDensityFineToCoarse.H new file mode 100644 index 000000000..947db2aac --- /dev/null +++ b/Source/Parallelization/InterpolateDensityFineToCoarse.H @@ -0,0 +1,120 @@ +/* Copyright 2019 Axel Huebl, Weiqun Zhang + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef INTERPOLATE_DENSITY_FINE_TO_COARSE_H +#define INTERPOLATE_DENSITY_FINE_TO_COARSE_H + +#include <AMReX_Array4.H> +#include <AMReX_REAL.H> +#include <AMReX_BLassert.H> +#include <AMReX_Extension.H> +#include <AMReX_GpuQualifiers.H> + +#include <utility> // std::move + + +/** Fill a charge density (rho) coarse patch with averaged values from a fine patch + * + * Fills the values of the charge density on the coarse patch + * by averaging the values of the charge density of the fine patch. + */ +class InterpolateDensityFineToCoarse +{ +public: + /** Construct with fine and coarse patch and their refinement ratio + * + * @param[in] fine read-only fine patch + * @param[in,out] coarse overwritten coarse patch + * @param[in] refinement_ratio ratio between coarse and fine patch granularity + * (currently, only a value of is implemented) + * @param[in] number_of_components the number of components + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + InterpolateDensityFineToCoarse( + amrex::Array4<amrex::Real const> const fine, + amrex::Array4<amrex::Real > const coarse, + int const refinement_ratio, + int const number_of_components + ) : m_fine(std::move(fine)), + m_coarse(std::move(coarse)), + m_refinement_ratio(std::move(refinement_ratio)), + m_number_of_components(std::move(number_of_components)) + { + //! @note constants and stencils in operator() implementation assume 2x refinement + BL_ASSERT(refinement_ratio == 2); + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void + operator()( + int const i, + int const j, + int const k + ) const noexcept + { + using namespace amrex; + + auto const & fine_unsafe = m_fine; // out-of-bounds access not secured with zero-values yet + auto const & coarse = m_coarse; // out-of-bounds access not secured but will also not occur + + // return zero for out-of-bounds accesses during interpolation + // this is efficiently used as a method to add neutral elements beyond guards in the average below + auto const fine = [fine_unsafe] AMREX_GPU_DEVICE (int const j, int const k, int const l, int const m) noexcept + { + return fine_unsafe.contains(j, k, l) ? fine_unsafe(j, k, l, m) : amrex::Real{0.}; + }; + + int const ii = i * m_refinement_ratio; + int const jj = j * m_refinement_ratio; + int const kk = k * m_refinement_ratio; + for( int m = 0; m < m_number_of_components; ++m ) { +#if AMREX_SPACEDIM == 2 + coarse(i,j,kk,m) = 0.25_rt * ( + fine(ii,jj,kk,m) + + 0.5_rt * ( + fine(ii-1,jj ,kk,m) + fine(ii+1,jj ,kk,m) + + fine(ii ,jj-1,kk,m) + fine(ii ,jj+1,kk,m) + ) + + 0.25_rt * ( + fine(ii-1,jj-1,kk,m) + fine(ii+1,jj-1,kk,m) + + fine(ii-1,jj+1,kk,m) + fine(ii+1,jj+1,kk,m) + ) + ); +#elif AMREX_SPACEDIM == 3 + coarse(i,j,k,m) = 0.125_rt * ( + fine(ii,jj,kk,m) + + 0.5_rt * ( + fine(ii-1,jj ,kk ,m) + fine(ii+1,jj ,kk ,m) + + fine(ii ,jj-1,kk ,m) + fine(ii ,jj+1,kk ,m) + + fine(ii ,jj ,kk-1,m) + fine(ii ,jj ,kk+1,m) + ) + + 0.25_rt * ( + fine(ii-1,jj-1,kk ,m) + fine(ii+1,jj-1,kk ,m) + + fine(ii-1,jj+1,kk ,m) + fine(ii+1,jj+1,kk ,m) + + fine(ii-1,jj ,kk-1,m) + fine(ii+1,jj ,kk-1,m) + + fine(ii-1,jj ,kk+1,m) + fine(ii+1,jj ,kk+1,m) + + fine(ii ,jj-1,kk-1,m) + fine(ii ,jj+1,kk-1,m) + + fine(ii ,jj-1,kk+1,m) + fine(ii ,jj+1,kk+1,m) + ) + + 0.125_rt * ( + fine(ii-1,jj-1,kk-1,m) + fine(ii-1,jj-1,kk+1,m) + + fine(ii-1,jj+1,kk-1,m) + fine(ii-1,jj+1,kk+1,m) + + fine(ii+1,jj-1,kk-1,m) + fine(ii+1,jj-1,kk+1,m) + + fine(ii+1,jj+1,kk-1,m) + fine(ii+1,jj+1,kk+1,m) + ) + ); +#endif + } + } +private: + amrex::Array4< amrex::Real const > m_fine; + amrex::Array4< amrex::Real > m_coarse; + int m_refinement_ratio; + int m_number_of_components; +}; + +#endif // INTERPOLATE_DENSITY_FINE_TO_COARSE_H diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index c5a8877d3..52df3dc25 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -3,10 +3,12 @@ #include <WarpX_f.H> #include <WarpXSumGuardCells.H> #include <Parallelization/InterpolateCurrentFineToCoarse.H> +#include <Parallelization/InterpolateDensityFineToCoarse.H> #include <algorithm> #include <cstdlib> + using namespace amrex; void @@ -505,7 +507,7 @@ WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > std::array< amrex::MultiFab *, 3 > const & coarse, int const refinement_ratio) { - BL_PROFILE("InterpolateCurrentFineToCoarse()"); + BL_PROFILE("interpolateCurrentFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); const IntVect& ng = (fine[0]->nGrowVect() + 1) / refinement_ratio; // add equivalent no. of guards to coarse patch @@ -537,6 +539,8 @@ WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > void WarpX::SyncRho () { + BL_PROFILE("SyncRho()"); + if (!rho_fp[0]) return; const int ncomp = rho_fp[0]->nComp(); @@ -546,7 +550,7 @@ WarpX::SyncRho () { rho_cp[lev]->setVal(0.0); const IntVect& refinement_ratio = refRatio(lev-1); - SyncRho(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]); + interpolateDensityFineToCoarse(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]); } // For each level @@ -559,29 +563,26 @@ WarpX::SyncRho () } void -WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int refinement_ratio) +WarpX::interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) { + BL_PROFILE("interpolateDensityFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); - const IntVect& ng = (fine.nGrowVect()+1)/refinement_ratio; + const IntVect& ng = (fine.nGrowVect() + 1) / refinement_ratio; // add equivalent no. of guards to coarse patch const int nc = fine.nComp(); #ifdef _OPEMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif { - FArrayBox ffab; - for (MFIter mfi(crse,true); mfi.isValid(); ++mfi) + // OMP in-box decomposition of coarse into tilebox + for (MFIter mfi(coarse, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,refinement_ratio),1); - ffab.resize(fbx, nc); - fbx &= fine[mfi].box(); - ffab.setVal(0.0); - ffab.copy(fine[mfi], fbx, 0, fbx, 0, nc); - WRPX_SYNC_RHO(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_ANYD(crse[mfi]), - BL_TO_FORTRAN_ANYD(ffab), - &nc); + const Box& bx = mfi.growntilebox(ng); // only grow to outer directions of tileboxes for filling guards + + amrex::ParallelFor( + bx, + InterpolateDensityFineToCoarse(fine.const_array(mfi), coarse.array(mfi), refinement_ratio, nc) + ); } } } @@ -713,7 +714,7 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev) if (rho_fp[lev]) { rho_cp[lev]->setVal(0.0); const IntVect& refinement_ratio = refRatio(lev-1); - SyncRho(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]); + interpolateDensityFineToCoarse(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]); } } diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index c5ec6fb5b..bdeab16e8 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -90,13 +90,6 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, const int l0 = compute_shape_factor<depos_order - lower_in_v>( sz0, z-stagger_shift); - // Set fields on particle to zero - Exp[ip] = 0; - Eyp[ip] = 0; - Ezp[ip] = 0; - Bxp[ip] = 0; - Byp[ip] = 0; - Bzp[ip] = 0; // Each field is gathered in a separate block of // AMREX_SPACEDIM nested loops because the deposition // order can differ for each component of each field diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 4ce83685d..58546a106 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -9,6 +9,11 @@ #include <PhotonParticleContainer.H> #include <LaserParticleContainer.H> +#ifdef WARPX_QED + #include <BreitWheelerEngineWrapper.H> + #include <QuantumSyncEngineWrapper.H> +#endif + #include <memory> #include <map> #include <string> @@ -208,6 +213,17 @@ protected: std::vector<PCTypes> species_types; +#ifdef WARPX_QED + // The QED engines + BreitWheelerEngine bw_engine; + QuantumSynchrotronEngine qs_engine; + //_______________________________ + + //Initialize QED engines and provides smart pointers + //to species who need QED processes + void InitQED (); +#endif + private: // physical particles (+ laser) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 715c97b99..c860d21f5 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,11 +1,12 @@ -#include <limits> -#include <algorithm> -#include <string> - #include <MultiParticleContainer.H> #include <WarpX_f.H> #include <WarpX.H> +#include <limits> +#include <algorithm> +#include <string> +#include <memory> + using namespace amrex; MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) @@ -149,6 +150,11 @@ MultiParticleContainer::InitData () // For each species, get the ID of its product species. // This is used for ionization and pair creation processes. mapSpeciesProduct(); + +#ifdef WARPX_QED + InitQED(); +#endif + } @@ -726,3 +732,19 @@ MultiParticleContainer::doFieldIonization () } // lev } // pc_source } + +#ifdef WARPX_QED +void MultiParticleContainer::InitQED () +{ + for (auto& pc : allcontainers) { + if(pc->has_quantum_sync()){ + pc->set_quantum_sync_engine_ptr + (std::make_shared<QuantumSynchrotronEngine>(qs_engine)); + } + if(pc->has_breit_wheeler()){ + pc->set_breit_wheeler_engine_ptr + (std::make_shared<BreitWheelerEngine>(bw_engine)); + } + } +} +#endif diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index a6ffd1d76..407cf26f3 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -46,6 +46,15 @@ public: amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; + // Don't push momenta for photons + virtual void PushP (int lev, + amrex::Real dt, + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) override {}; // DepositCurrent should do nothing for photons @@ -63,8 +72,7 @@ public: int thread_num, int lev, int depos_lev, - amrex::Real dt) {}; - + amrex::Real dt) override {}; }; #endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 4a75ec9f3..3c70a957f 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -21,7 +21,25 @@ using namespace amrex; PhotonParticleContainer::PhotonParticleContainer (AmrCore* amr_core, int ispecies, const std::string& name) : PhysicalParticleContainer(amr_core, ispecies, name) -{} +{ + + ParmParse pp(species_name); + +#ifdef WARPX_QED + //IF do_qed is enabled, find out if Breit Wheeler process is enabled + if(do_qed) + pp.query("do_qed_breit_wheeler", do_qed_breit_wheeler); + + //Check for processes which do not make sense for photons + bool test_quantum_sync; + pp.query("do_qed_quantum_sync", test_quantum_sync); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + test_quantum_sync == 0, + "ERROR: do_qed_quantum_sync can be 1 for species NOT listed in particles.photon_species only!"); + //_________________________________________________________ +#endif + +} void PhotonParticleContainer::InitData() { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b558323a3..c953aa2d7 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -7,6 +7,11 @@ #include <AMReX_IArrayBox.H> +#ifdef WARPX_QED + #include <QuantumSyncEngineWrapper.H> + #include <BreitWheelerEngineWrapper.H> +#endif + #include <map> class PhysicalParticleContainer @@ -136,7 +141,15 @@ public: amrex::Real q_tot, long npart, int do_symmetrize); void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, - std::array<amrex::Real, 3> u, amrex::Real weight); + std::array<amrex::Real, 3> u, + amrex::Real weight, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_x, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_y, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_z, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_ux, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_uy, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_uz, + amrex::Gpu::HostVector<amrex::ParticleReal>& particle_w); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, const amrex::Real z_new, const amrex::Real t_boost, @@ -172,6 +185,16 @@ public: amrex::FArrayBox const * & ezfab, amrex::FArrayBox const * & bxfab, amrex::FArrayBox const * & byfab, amrex::FArrayBox const * & bzfab); +#ifdef WARPX_QED + bool has_quantum_sync() override; + bool has_breit_wheeler() override; + + void set_breit_wheeler_engine_ptr + (std::shared_ptr<BreitWheelerEngine>) override; + void set_quantum_sync_engine_ptr + (std::shared_ptr<QuantumSynchrotronEngine>) override; +#endif + protected: std::string species_name; @@ -186,6 +209,20 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; +#ifdef WARPX_QED + // A flag to enable quantum_synchrotron process for leptons + bool do_qed_quantum_sync = false; + + // A flag to enable breit_wheeler process [photons only!!] + bool do_qed_breit_wheeler = false; + + // A smart pointer to an instance of a Quantum Synchrotron engine + std::shared_ptr<QuantumSynchrotronEngine> shr_ptr_qs_engine; + + // A smart pointer to an instance of a Breit Wheeler engine [photons only!] + std::shared_ptr<BreitWheelerEngine> shr_ptr_bw_engine; +#endif + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 02dee1967..d10803320 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1,6 +1,8 @@ #include <limits> #include <sstream> +#include <PhysicalParticleContainer.H> + #include <MultiParticleContainer.H> #include <WarpX_f.H> #include <WarpX.H> @@ -49,6 +51,33 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp "version of Redistribute in AMReX does not work with runtime parameters"); #endif + +#ifdef WARPX_QED + //Add real component if QED is enabled + pp.query("do_qed", do_qed); + if(do_qed) + AddRealComp("tau"); + + //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled + if(do_qed) + pp.query("do_qed_quantum_sync", do_qed_quantum_sync); + + //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!! +#endif + + //variable to set plot_flags size + int plot_flag_size = PIdx::nattribs; + if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + plot_flag_size += 6; + +#ifdef WARPX_QED + if(do_qed){ + // plot_flag will have an entry for the optical depth + plot_flag_size++; + } +#endif + //_______________________________ + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); @@ -56,18 +85,11 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // By default, all particle variables are dumped to plotfiles, // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ - plot_flags.resize(PIdx::nattribs + 6, 1); - } else { - plot_flags.resize(PIdx::nattribs, 1); - } + plot_flags.resize(plot_flag_size, 1); } else { // Set plot_flag to 0 for all attribs - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ - plot_flags.resize(PIdx::nattribs + 6, 0); - } else { - plot_flags.resize(PIdx::nattribs, 0); - } + plot_flags.resize(plot_flag_size, 0); + // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ @@ -79,6 +101,13 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } } } + + #ifdef WARPX_QED + if(do_qed){ + //Optical depths is always plotted if QED is on + plot_flags[plot_flag_size-1] = 1; + } + #endif } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) @@ -92,7 +121,9 @@ void PhysicalParticleContainer::InitData() // Init ionization module here instead of in the PhysicalParticleContainer // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} + AddParticles(0); // Note - add on level 0 + Redistribute(); // We then redistribute } @@ -156,6 +187,16 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution<double> disty(y_m, y_rms); std::normal_distribution<double> distz(z_m, z_rms); + // Allocate temporary vectors on the CPU + Gpu::HostVector<ParticleReal> particle_x; + Gpu::HostVector<ParticleReal> particle_y; + Gpu::HostVector<ParticleReal> particle_z; + Gpu::HostVector<ParticleReal> particle_ux; + Gpu::HostVector<ParticleReal> particle_uy; + Gpu::HostVector<ParticleReal> particle_uz; + Gpu::HostVector<ParticleReal> particle_w; + int np = 0; + if (ParallelDescriptor::IOProcessor()) { // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) @@ -181,44 +222,61 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, u.z *= PhysConst::c; if (do_symmetrize){ // Add four particles to the beam: - CheckAndAddParticle( x, y, z, { u.x, u.y, u.z}, weight/4. ); - CheckAndAddParticle( x,-y, z, { u.x,-u.y, u.z}, weight/4. ); - CheckAndAddParticle(-x, y, z, {-u.x, u.y, u.z}, weight/4. ); - CheckAndAddParticle(-x,-y, z, {-u.x,-u.y, u.z}, weight/4. ); + CheckAndAddParticle(x, y, z, { u.x, u.y, u.z}, weight/4., + particle_x, particle_y, particle_z, + particle_ux, particle_uy, particle_uz, + particle_w); + CheckAndAddParticle(x, -y, z, { u.x, -u.y, u.z}, weight/4., + particle_x, particle_y, particle_z, + particle_ux, particle_uy, particle_uz, + particle_w); + CheckAndAddParticle(-x, y, z, { -u.x, u.y, u.z}, weight/4., + particle_x, particle_y, particle_z, + particle_ux, particle_uy, particle_uz, + particle_w); + CheckAndAddParticle(-x, -y, z, { -u.x, -u.y, u.z}, weight/4., + particle_x, particle_y, particle_z, + particle_ux, particle_uy, particle_uz, + particle_w); } else { - CheckAndAddParticle(x, y, z, {u.x,u.y,u.z}, weight); + CheckAndAddParticle(x, y, z, { u.x, u.y, u.z}, weight, + particle_x, particle_y, particle_z, + particle_ux, particle_uy, particle_uz, + particle_w); } } } } - Redistribute(); + // Add the temporary CPU vectors to the particle structure + np = particle_z.size(); + AddNParticles(0,np, + particle_x.dataPtr(), particle_y.dataPtr(), particle_z.dataPtr(), + particle_ux.dataPtr(), particle_uy.dataPtr(), particle_uz.dataPtr(), + 1, particle_w.dataPtr(),1); } void PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, std::array<Real, 3> u, - Real weight) + Real weight, + Gpu::HostVector<ParticleReal>& particle_x, + Gpu::HostVector<ParticleReal>& particle_y, + Gpu::HostVector<ParticleReal>& particle_z, + Gpu::HostVector<ParticleReal>& particle_ux, + Gpu::HostVector<ParticleReal>& particle_uy, + Gpu::HostVector<ParticleReal>& particle_uz, + Gpu::HostVector<ParticleReal>& particle_w) { - std::array<ParticleReal,PIdx::nattribs> attribs; - attribs.fill(0.0); - - // update attribs with input arguments if (WarpX::gamma_boost > 1.) { MapParticletoBoostedFrame(x, y, z, u); } - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - attribs[PIdx::w ] = weight; - - if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) - { - // need to create old values - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - } - - // add particle - AddOneParticle(0, 0, 0, x, y, z, attribs); + particle_x.push_back(x); + particle_y.push_back(y); + particle_z.push_back(z); + particle_ux.push_back(u[0]); + particle_uy.push_back(u[1]); + particle_uz.push_back(u[2]); + particle_w.push_back(weight); } void @@ -450,6 +508,30 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } +#ifdef WARPX_QED + //Pointer to the optical depth component + amrex::Real* p_tau; + + //If a QED effect is enabled, tau has to be initialized + bool loc_has_quantum_sync = has_quantum_sync(); + bool loc_has_breit_wheeler = has_breit_wheeler(); + if(loc_has_quantum_sync || loc_has_breit_wheeler){ + p_tau = soa.GetRealData(particle_comps["tau"]).data() + old_size; + } + + //If needed, get the appropriate functors from the engines + QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; + BreitWheelerGetOpticalDepth breit_wheeler_get_opt; + if(loc_has_quantum_sync){ + quantum_sync_get_opt = + shr_ptr_qs_engine->build_optical_depth_functor(); + } + if(loc_has_breit_wheeler){ + breit_wheeler_get_opt = + shr_ptr_bw_engine->build_optical_depth_functor(); + } +#endif + const GpuArray<Real,AMREX_SPACEDIM> overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -597,6 +679,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pi[ip] = loc_ionization_initial_level; } +#ifdef WARPX_QED + if(loc_has_quantum_sync){ + p_tau[ip] = quantum_sync_get_opt(); + } + + if(loc_has_breit_wheeler){ + p_tau[ip] = breit_wheeler_get_opt(); + } +#endif + u.x *= PhysConst::c; u.y *= PhysConst::c; u.z *= PhysConst::c; @@ -2047,3 +2139,30 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } ); } + +#ifdef WARPX_QED + +bool PhysicalParticleContainer::has_quantum_sync() +{ + return do_qed_quantum_sync; +} + +bool PhysicalParticleContainer::has_breit_wheeler() +{ + return do_qed_breit_wheeler; +} + +void +PhysicalParticleContainer:: +set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine> ptr) +{ + shr_ptr_bw_engine = ptr; +} + +void +PhysicalParticleContainer:: +set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine> ptr) +{ + shr_ptr_qs_engine = ptr; +} +#endif diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp index d1455249a..e88af017f 100644 --- a/Source/Particles/Sorting/Partition.cpp +++ b/Source/Particles/Sorting/Partition.cpp @@ -56,7 +56,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( gather_masks : current_masks; // - For each particle, find whether it is in the larger buffer, // by looking up the mask. Store the answer in `inexflag`. - amrex::ParallelFor( np, fillBufferFlag(pti, bmasks, inexflag, Geom(lev), 0) ); + amrex::ParallelFor( np, fillBufferFlag(pti, bmasks, inexflag, Geom(lev)) ); // - Find the indices that reorder particles so that the last particles // are in the larger buffer fillWithConsecutiveIntegers( pid ); @@ -93,7 +93,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // - For each particle in the large buffer, find whether it is in // the smaller buffer, by looking up the mask. Store the answer in `inexflag`. amrex::ParallelFor( np - n_fine, - fillBufferFlag(pti, bmasks, inexflag, Geom(lev), n_fine) ); + fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, n_fine) ); auto const sep2 = stablePartition( sep, pid.end(), inexflag ); if (bmasks == gather_masks) { diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index 852b8a73f..28a1b73b7 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -78,27 +78,84 @@ int iteratorDistance(ForwardIterator const first, /* \brief Functor that fills the elements of the particle array `inexflag` * with the value of the spatial array `bmasks`, at the corresponding particle position. * - * This is done only for the elements from `start_index` to the end of `inexflag` - * * \param[in] pti Contains information on the particle positions * \param[in] bmasks Spatial array, that contains a flag indicating * whether each cell is part of the gathering/deposition buffers * \param[out] inexflag Vector to be filled with the value of `bmasks` * \param[in] geom Geometry object, necessary to locate particles within the array `bmasks` - * \param[in] start_index Index that which elements start to be modified + * */ class fillBufferFlag { public: fillBufferFlag( WarpXParIter const& pti, amrex::iMultiFab const* bmasks, amrex::Gpu::DeviceVector<int>& inexflag, - amrex::Geometry const& geom, long const start_index=0 ) : + amrex::Geometry const& geom ) { + + // Extract simple structure that can be used directly on the GPU + m_particles = &(pti.GetArrayOfStructs()[0]); + m_buffer_mask = (*bmasks)[pti].array(); + m_inexflag_ptr = inexflag.dataPtr(); + m_domain = geom.Domain(); + for (int idim=0; idim<AMREX_SPACEDIM; idim++) { + m_prob_lo[idim] = geom.ProbLo(idim); + m_inv_cell_size[idim] = geom.InvCellSize(idim); + } + }; + + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator()( const long i ) const { + // Select a particle + auto const& p = m_particles[i]; + // Find the index of the cell where this particle is located + amrex::IntVect const iv = amrex::getParticleCell( p, + m_prob_lo, m_inv_cell_size, m_domain ); + // Find the value of the buffer flag in this cell and + // store it at the corresponding particle position in the array `inexflag` + m_inexflag_ptr[i] = m_buffer_mask(iv); + }; + + private: + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_prob_lo; + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_inv_cell_size; + amrex::Box m_domain; + int* m_inexflag_ptr; + WarpXParticleContainer::ParticleType const* m_particles; + amrex::Array4<int const> m_buffer_mask; +}; + +/* \brief Functor that fills the elements of the particle array `inexflag` + * with the value of the spatial array `bmasks`, at the corresponding particle position. + * + * Contrary to `fillBufferFlag`, here this is done only for the particles that + * the last elements of `particle_indices` point to (from the element at + * index `start_index` in `particle_indices`, to the last element of `particle_indices`) + * + * \param[in] pti Contains information on the particle positions + * \param[in] bmasks Spatial array, that contains a flag indicating + * whether each cell is part of the gathering/deposition buffers + * \param[out] inexflag Vector to be filled with the value of `bmasks` + * \param[in] geom Geometry object, necessary to locate particles within the array `bmasks` + * \param[in] start_index Index that which elements start to be modified + */ +class fillBufferFlagRemainingParticles +{ + public: + fillBufferFlagRemainingParticles( + WarpXParIter const& pti, + amrex::iMultiFab const* bmasks, + amrex::Gpu::DeviceVector<int>& inexflag, + amrex::Geometry const& geom, + amrex::Gpu::DeviceVector<long> const& particle_indices, + long const start_index ) : m_start_index(start_index) { // Extract simple structure that can be used directly on the GPU m_particles = &(pti.GetArrayOfStructs()[0]); m_buffer_mask = (*bmasks)[pti].array(); m_inexflag_ptr = inexflag.dataPtr(); + m_indices_ptr = particle_indices.dataPtr(); m_domain = geom.Domain(); for (int idim=0; idim<AMREX_SPACEDIM; idim++) { m_prob_lo[idim] = geom.ProbLo(idim); @@ -110,13 +167,13 @@ class fillBufferFlag AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()( const long i ) const { // Select a particle - auto const& p = m_particles[i+m_start_index]; + auto const& p = m_particles[m_indices_ptr[i+m_start_index]]; // Find the index of the cell where this particle is located amrex::IntVect const iv = amrex::getParticleCell( p, m_prob_lo, m_inv_cell_size, m_domain ); // Find the value of the buffer flag in this cell and // store it at the corresponding particle position in the array `inexflag` - m_inexflag_ptr[i+m_start_index] = m_buffer_mask(iv); + m_inexflag_ptr[m_indices_ptr[i+m_start_index]] = m_buffer_mask(iv); }; private: @@ -127,6 +184,7 @@ class fillBufferFlag WarpXParticleContainer::ParticleType const* m_particles; amrex::Array4<int const> m_buffer_mask; long const m_start_index; + long const* m_indices_ptr; }; /* \brief Functor that copies the elements of `src` into `dst`, diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7b0d2d1d0..879f4782e 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -6,6 +6,11 @@ #include <AMReX_Particles.H> #include <AMReX_AmrCore.H> +#ifdef WARPX_QED + #include <QuantumSyncEngineWrapper.H> + #include <BreitWheelerEngineWrapper.H> +#endif + #include <memory> enum struct ConvertDirection{WarpX_to_SI, SI_to_WarpX}; @@ -236,14 +241,6 @@ public: const amrex::ParticleReal* vx, const amrex::ParticleReal* vy, const amrex::ParticleReal* vz, int nattr, const amrex::ParticleReal* attr, int uniqueparticles, int id=-1); - void AddOneParticle (int lev, int grid, int tile, - amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, - std::array<amrex::ParticleReal,PIdx::nattribs>& attribs); - - void AddOneParticle (ParticleTileType& particle_tile, - amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, - std::array<amrex::ParticleReal,PIdx::nattribs>& attribs); - virtual void ReadHeader (std::istream& is); virtual void WriteHeader (std::ostream& os) const; @@ -321,6 +318,18 @@ protected: int do_boosted_frame_diags = 1; +#ifdef WARPX_QED + bool do_qed = false; + + virtual bool has_quantum_sync(){return false;}; + virtual bool has_breit_wheeler(){return false;}; + + virtual void + set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine>){}; + virtual void + set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine>){}; +#endif + amrex::Vector<amrex::FArrayBox> local_rho; amrex::Vector<amrex::FArrayBox> local_jx; amrex::Vector<amrex::FArrayBox> local_jy; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 65a82f233..7fb57500d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -135,45 +135,6 @@ WarpXParticleContainer::AllocData () } void -WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, - ParticleReal x, ParticleReal y, ParticleReal z, - std::array<ParticleReal,PIdx::nattribs>& attribs) -{ - auto& particle_tile = DefineAndReturnParticleTile(lev, grid, tile); - AddOneParticle(particle_tile, x, y, z, attribs); -} - -void -WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, - ParticleReal x, ParticleReal y, ParticleReal z, - std::array<ParticleReal,PIdx::nattribs>& attribs) -{ - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); -#if (AMREX_SPACEDIM == 3) - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; -#elif (AMREX_SPACEDIM == 2) -#ifdef WARPX_DIM_RZ - attribs[PIdx::theta] = std::atan2(y, x); - x = std::sqrt(x*x + y*y); -#endif - p.pos(0) = x; - p.pos(1) = z; -#endif - - particle_tile.push_back(p); - particle_tile.push_back_real(attribs); - - for (int i = PIdx::nattribs; i < NumRealComps(); ++i) - { - particle_tile.push_back_real(i, 0.0); - } -} - -void WarpXParticleContainer::AddNParticles (int lev, int n, const ParticleReal* x, const ParticleReal* y, const ParticleReal* z, const ParticleReal* vx, const ParticleReal* vy, const ParticleReal* vz, diff --git a/Source/QED/BreitWheelerEngineWrapper.H b/Source/QED/BreitWheelerEngineWrapper.H new file mode 100644 index 000000000..7aeb7bca2 --- /dev/null +++ b/Source/QED/BreitWheelerEngineWrapper.H @@ -0,0 +1,55 @@ +#ifndef WARPX_breit_wheeler_engine_wrapper_h_ +#define WARPX_breit_wheeler_engine_wrapper_h_ + +#include "QedWrapperCommons.H" + +//BW ENGINE from PICSAR +#include "breit_wheeler_engine.hpp" + +using WarpXBreitWheelerWrapper = + picsar::multi_physics::breit_wheeler_engine<amrex::Real, DummyStruct>; + +// Functors ================================== + +// These functors provide the core elementary functions of the library +// Can be included in GPU kernels + +/** + * \brief Functor to initialize the optical depth of photons for the + * Breit-Wheeler process + */ +class BreitWheelerGetOpticalDepth +{ +public: + BreitWheelerGetOpticalDepth () + {}; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real operator() () const + { + return WarpXBreitWheelerWrapper:: + internal_get_optical_depth(amrex::Random()); + } +}; +//____________________________________________ + +// Factory class ============================= + +/** + * \brief Wrapper for the Breit Wheeler engine of the PICSAR library + */ +class BreitWheelerEngine +{ +public: + BreitWheelerEngine (); + + /** + * \brief Builds the functor to initialize the optical depth + */ + BreitWheelerGetOpticalDepth build_optical_depth_functor (); +}; + +//============================================ + +#endif //WARPX_breit_wheeler_engine_wrapper_H_ diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp new file mode 100644 index 000000000..33e6d48d7 --- /dev/null +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -0,0 +1,16 @@ +#include "BreitWheelerEngineWrapper.H" +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the PICSAR library + +using namespace picsar::multi_physics; + +// Factory class ============================= + +BreitWheelerEngine::BreitWheelerEngine (){} + +BreitWheelerGetOpticalDepth BreitWheelerEngine::build_optical_depth_functor () +{ + return BreitWheelerGetOpticalDepth(); +} + +//============================================ diff --git a/Source/QED/Make.package b/Source/QED/Make.package new file mode 100644 index 000000000..d4bad3f18 --- /dev/null +++ b/Source/QED/Make.package @@ -0,0 +1,8 @@ +CEXE_headers += QedWrapperCommons.H +CEXE_headers += BreitWheelerEngineWrapper.H +CEXE_headers += QuantumSyncsEngineWrapper.H +CEXE_sources += BreitWheelerEngineWrapper.cpp +CEXE_sources += QuantumSyncEngineWrapper.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED +VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/QedWrapperCommons.H b/Source/QED/QedWrapperCommons.H new file mode 100644 index 000000000..821034c06 --- /dev/null +++ b/Source/QED/QedWrapperCommons.H @@ -0,0 +1,16 @@ +#ifndef WARPX_amrex_qed_wrapper_commons_h_ +#define WARPX_amrex_qed_wrapper_commons_h_ + +//Common definitions for the QED library wrappers + +#include <AMReX_AmrCore.H> + +//Sets the decorator for GPU +#define PXRMP_GPU AMREX_GPU_DEVICE +//Sets SI units in the library +#define PXRMP_WITH_SI_UNITS + +//An empty data type +struct DummyStruct{}; + +#endif //WARPX_amrex_qed_wrapper_commons_h_ diff --git a/Source/QED/QuantumSyncEngineWrapper.H b/Source/QED/QuantumSyncEngineWrapper.H new file mode 100644 index 000000000..b0e90995d --- /dev/null +++ b/Source/QED/QuantumSyncEngineWrapper.H @@ -0,0 +1,55 @@ +#ifndef WARPX_quantum_sync_engine_wrapper_h_ +#define WARPX_quantum_sync_engine_wrapper_h_ + +#include "QedWrapperCommons.H" + +//QS ENGINE from PICSAR +#include "quantum_sync_engine.hpp" + +using WarpXQuantumSynchrotronWrapper = + picsar::multi_physics::quantum_synchrotron_engine<amrex::Real, DummyStruct>; + +// Functors ================================== + +// These functors provide the core elementary functions of the library +// Can be included in GPU kernels + +/** +* Functor to initialize the optical depth of leptons for the +* Quantum Synchrotron process +*/ +class QuantumSynchrotronGetOpticalDepth +{ +public: + QuantumSynchrotronGetOpticalDepth () + {}; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real operator() () const + { + return WarpXQuantumSynchrotronWrapper:: + internal_get_optical_depth(amrex::Random()); + } +}; +//____________________________________________ + +// Factory class ============================= + +/** + * \brief Wrapper for the Quantum Synchrotron engine of the PICSAR library + */ +class QuantumSynchrotronEngine +{ +public: + QuantumSynchrotronEngine (); + + /** + * \brief Builds the functor to initialize the optical depth + */ + QuantumSynchrotronGetOpticalDepth build_optical_depth_functor (); +}; + +//============================================ + +#endif //WARPX_quantum_sync_engine_wrapper_h_ diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp new file mode 100644 index 000000000..b55149187 --- /dev/null +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -0,0 +1,17 @@ +#include "QuantumSyncEngineWrapper.H" +//This file provides a wrapper aroud the quantum_sync engine +//provided by the PICSAR library + +using namespace picsar::multi_physics; + +// Factory class ============================= + +QuantumSynchrotronEngine::QuantumSynchrotronEngine (){} + +QuantumSynchrotronGetOpticalDepth +QuantumSynchrotronEngine::build_optical_depth_functor () +{ + return QuantumSynchrotronGetOpticalDepth(); +} + +//============================================ diff --git a/Source/WarpX.H b/Source/WarpX.H index dd86dfc14..55fd6e83d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -438,10 +438,12 @@ private: * averaging the values of the charge density of the fine patch (on the same level). * * \param[in] fine fine patches to interpolate from - * \param[out] crse coarse patches to interpolate to - * \param[in] ref_ratio integer ratio between the two + * \param[out] coarse coarse patches to interpolate to + * \param[in] refinement_ratio integer ratio between the two */ - void SyncRho (const amrex::MultiFab& fine, amrex::MultiFab& crse, int ref_ratio); + void interpolateDensityFineToCoarse (const amrex::MultiFab& fine, + amrex::MultiFab& coarse, + int const refinement_ratio); void ExchangeWithPmlB (int lev); void ExchangeWithPmlE (int lev); |