From 878721f5d77bd37bba5a7b85c23612199033c831 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Fri, 13 Sep 2024 12:44:52 -0700 Subject: [PATCH] Refactor AddPlasma and AddPlasmaFlux (#5231) --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Axel Huebl --- Source/Particles/AddPlasmaUtilities.H | 210 +++++++++ Source/Particles/AddPlasmaUtilities.cpp | 158 +++++++ Source/Particles/CMakeLists.txt | 1 + Source/Particles/Make.package | 1 + Source/Particles/PhysicalParticleContainer.H | 8 + .../Particles/PhysicalParticleContainer.cpp | 436 +++++------------- 6 files changed, 481 insertions(+), 333 deletions(-) create mode 100644 Source/Particles/AddPlasmaUtilities.H create mode 100644 Source/Particles/AddPlasmaUtilities.cpp diff --git a/Source/Particles/AddPlasmaUtilities.H b/Source/Particles/AddPlasmaUtilities.H new file mode 100644 index 00000000000..8f0489e3921 --- /dev/null +++ b/Source/Particles/AddPlasmaUtilities.H @@ -0,0 +1,210 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + * Authors: Andrew Myers + */ +#ifndef WARPX_ADDPLASMAUTILITIES_H_ +#define WARPX_ADDPLASMAUTILITIES_H_ + +#include "Initialization/PlasmaInjector.H" + +#ifdef WARPX_QED +# include "Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H" +# include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H" +#endif + +#include +#include +#include +#include +#include +#include + +/* + Finds the overlap region between the given tile_realbox and part_realbox, returning true + if an overlap exists and false if otherwise. This also sets the parameters overlap_realbox, + overlap_box, and shifted to the appropriate values. + */ +bool find_overlap (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted); + +/* + Finds the overlap region between the given tile_realbox, part_realbox and the surface used for + flux injection, returning true if an overlap exists and false if otherwise. This also sets the + parameters overlap_realbox, overlap_box, and shifted to the appropriate values. + */ +bool find_overlap_flux (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + const PlasmaInjector& plasma_injector, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted); + +/* + This computes the scale_fac (used for setting the particle weights) on a volumetric basis. + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real compute_scale_fac_volume (const amrex::GpuArray& dx, + const amrex::Long pcount) { + using namespace amrex::literals; + return (pcount != 0) ? AMREX_D_TERM(dx[0],*dx[1],*dx[2])/pcount : 0.0_rt; +} + +/* + Given a refinement ratio, this computes the total increase in resolution for a plane + defined by the normal_axis. + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int compute_area_weights (const amrex::IntVect& iv, const int normal_axis) { + int r = AMREX_D_TERM(iv[0],*iv[1],*iv[2]); +#if defined(WARPX_DIM_3D) + r /= iv[normal_axis]; +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + if (normal_axis == 0) { r /= iv[0]; } + else if (normal_axis == 2) { r /= iv[1]; } +#elif defined(WARPX_DIM_1D_Z) + if (normal_axis == 2) { r /= iv[0]; } +#endif + return r; +} + +/* + This computes the scale_fac (used for setting the particle weights) on a on area basis + (used for flux injection). + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real compute_scale_fac_area (const amrex::GpuArray& dx, + const amrex::Real num_ppc_real, const int flux_normal_axis) { + using namespace amrex::literals; + amrex::Real scale_fac = AMREX_D_TERM(dx[0],*dx[1],*dx[2])/num_ppc_real; + // Scale particle weight by the area of the emitting surface, within one cell +#if defined(WARPX_DIM_3D) + scale_fac /= dx[flux_normal_axis]; +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + // When emission is in the r direction, the emitting surface is a cylinder. + // The factor 2*pi*r is added later below. + if (flux_normal_axis == 0) { scale_fac /= dx[0]; } + // When emission is in the z direction, the emitting surface is an annulus + // The factor 2*pi*r is added later below. + if (flux_normal_axis == 2) { scale_fac /= dx[1]; } + // When emission is in the theta direction (flux_normal_axis == 1), + // the emitting surface is a rectangle, within the plane of the simulation +#elif defined(WARPX_DIM_1D_Z) + if (flux_normal_axis == 2) { scale_fac /= dx[0]; } +#endif + return scale_fac; +} + +/* + These structs encapsulates several data structures needed for using the parser during plasma + injection. + */ +struct PlasmaParserWrapper +{ + PlasmaParserWrapper (std::size_t a_num_user_int_attribs, + std::size_t a_num_user_real_attribs, + const amrex::Vector< std::unique_ptr >& a_user_int_attrib_parser, + const amrex::Vector< std::unique_ptr >& a_user_real_attrib_parser); + + amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > m_user_int_attrib_parserexec_pinned; + amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > m_user_real_attrib_parserexec_pinned; +}; + +struct PlasmaParserHelper +{ + template + PlasmaParserHelper (SoAType& a_soa, std::size_t old_size, + const std::vector& a_user_int_attribs, + const std::vector& a_user_real_attribs, + std::map& a_particle_icomps, + std::map& a_particle_comps, + const PlasmaParserWrapper& wrapper) : + m_wrapper_ptr(&wrapper) { + m_pa_user_int_pinned.resize(a_user_int_attribs.size()); + m_pa_user_real_pinned.resize(a_user_real_attribs.size()); + +#ifdef AMREX_USE_GPU + m_d_pa_user_int.resize(a_user_int_attribs.size()); + m_d_pa_user_real.resize(a_user_real_attribs.size()); + m_d_user_int_attrib_parserexec.resize(a_user_int_attribs.size()); + m_d_user_real_attrib_parserexec.resize(a_user_real_attribs.size()); +#endif + + for (std::size_t ia = 0; ia < a_user_int_attribs.size(); ++ia) { + m_pa_user_int_pinned[ia] = a_soa.GetIntData(a_particle_icomps[a_user_int_attribs[ia]]).data() + old_size; + } + for (std::size_t ia = 0; ia < a_user_real_attribs.size(); ++ia) { + m_pa_user_real_pinned[ia] = a_soa.GetRealData(a_particle_comps[a_user_real_attribs[ia]]).data() + old_size; + } + +#ifdef AMREX_USE_GPU + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_pa_user_int_pinned.begin(), + m_pa_user_int_pinned.end(), m_d_pa_user_int.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_pa_user_real_pinned.begin(), + m_pa_user_real_pinned.end(), m_d_pa_user_real.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, wrapper.m_user_int_attrib_parserexec_pinned.begin(), + wrapper.m_user_int_attrib_parserexec_pinned.end(), m_d_user_int_attrib_parserexec.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, wrapper.m_user_real_attrib_parserexec_pinned.begin(), + wrapper.m_user_real_attrib_parserexec_pinned.end(), m_d_user_real_attrib_parserexec.begin()); +#endif + } + + int** getUserIntDataPtrs (); + amrex::ParticleReal** getUserRealDataPtrs (); + [[nodiscard]] amrex::ParserExecutor<7> const* getUserIntParserExecData () const; + [[nodiscard]] amrex::ParserExecutor<7> const* getUserRealParserExecData () const; + + amrex::Gpu::PinnedVector m_pa_user_int_pinned; + amrex::Gpu::PinnedVector m_pa_user_real_pinned; + +#ifdef AMREX_USE_GPU + // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, + // and them memcpy to device from host + amrex::Gpu::DeviceVector m_d_pa_user_int; + amrex::Gpu::DeviceVector m_d_pa_user_real; + amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > m_d_user_int_attrib_parserexec; + amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > m_d_user_real_attrib_parserexec; +#endif + const PlasmaParserWrapper* m_wrapper_ptr; +}; + +#ifdef WARPX_QED +struct QEDHelper +{ + template + QEDHelper (SoAType& a_soa, std::size_t old_size, + std::map& a_particle_comps, + bool a_has_quantum_sync, bool a_has_breit_wheeler, + const std::shared_ptr& a_shr_p_qs_engine, + const std::shared_ptr& a_shr_p_bw_engine) + : has_quantum_sync(a_has_quantum_sync), has_breit_wheeler(a_has_breit_wheeler) + { + if(has_quantum_sync){ + quantum_sync_get_opt = + a_shr_p_qs_engine->build_optical_depth_functor(); + p_optical_depth_QSR = a_soa.GetRealData( + a_particle_comps["opticalDepthQSR"]).data() + old_size; + } + if(has_breit_wheeler){ + breit_wheeler_get_opt = + a_shr_p_bw_engine->build_optical_depth_functor(); + p_optical_depth_BW = a_soa.GetRealData( + a_particle_comps["opticalDepthBW"]).data() + old_size; + } + } + + amrex::ParticleReal* p_optical_depth_QSR = nullptr; + amrex::ParticleReal* p_optical_depth_BW = nullptr; + + const bool has_quantum_sync; + const bool has_breit_wheeler; + + QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; + BreitWheelerGetOpticalDepth breit_wheeler_get_opt; +}; +#endif + +#endif /*WARPX_ADDPLASMAUTILITIES_H_*/ diff --git a/Source/Particles/AddPlasmaUtilities.cpp b/Source/Particles/AddPlasmaUtilities.cpp new file mode 100644 index 00000000000..31066516477 --- /dev/null +++ b/Source/Particles/AddPlasmaUtilities.cpp @@ -0,0 +1,158 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + * Authors: Andrew Myers + */ +#include "AddPlasmaUtilities.H" + +#include + +bool find_overlap (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted) +{ + using namespace amrex::literals; + + bool no_overlap = false; + for (int dir=0; dir= part_realbox.lo(dir) ) { + const amrex::Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); + } else { + no_overlap = true; break; + } + // Count the number of cells in this direction in overlap_realbox + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, + int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) + /dx[dir] )) - 1); + shifted[dir] = + static_cast(std::round((overlap_realbox.lo(dir)-prob_lo[dir])/dx[dir])); + // shifted is exact in non-moving-window direction. That's all we care. + } + return no_overlap; +} + +bool find_overlap_flux (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + const PlasmaInjector& plasma_injector, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted) +{ + using namespace amrex::literals; + + bool no_overlap = false; + for (int dir=0; dir 0) { + if (plasma_injector.surface_flux_pos < tile_realbox.lo(dir) || + plasma_injector.surface_flux_pos >= tile_realbox.hi(dir)) { + no_overlap = true; + break; + } + } else { + if (plasma_injector.surface_flux_pos <= tile_realbox.lo(dir) || + plasma_injector.surface_flux_pos > tile_realbox.hi(dir)) { + no_overlap = true; + break; + } + } + overlap_realbox.setLo( dir, plasma_injector.surface_flux_pos ); + overlap_realbox.setHi( dir, plasma_injector.surface_flux_pos ); + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, 0 ); + shifted[dir] = + static_cast(std::round((overlap_realbox.lo(dir)-prob_lo[dir])/dx[dir])); + } else { + if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { + const amrex::Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); + overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); + } else { + no_overlap = true; break; + } + if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { + const amrex::Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); + } else { + no_overlap = true; break; + } + // Count the number of cells in this direction in overlap_realbox + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, + int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) + /dx[dir] )) - 1); + shifted[dir] = + static_cast(std::round((overlap_realbox.lo(dir)-prob_lo[dir])/dx[dir])); + // shifted is exact in non-moving-window direction. That's all we care. + } + } + + return no_overlap; +} + +PlasmaParserWrapper::PlasmaParserWrapper (const std::size_t a_num_user_int_attribs, + const std::size_t a_num_user_real_attribs, + const amrex::Vector< std::unique_ptr >& a_user_int_attrib_parser, + const amrex::Vector< std::unique_ptr >& a_user_real_attrib_parser) + +{ + m_user_int_attrib_parserexec_pinned.resize(a_num_user_int_attribs); + m_user_real_attrib_parserexec_pinned.resize(a_num_user_real_attribs); + + for (std::size_t ia = 0; ia < a_num_user_int_attribs; ++ia) { + m_user_int_attrib_parserexec_pinned[ia] = a_user_int_attrib_parser[ia]->compile<7>(); + } + for (std::size_t ia = 0; ia < a_num_user_real_attribs; ++ia) { + m_user_real_attrib_parserexec_pinned[ia] = a_user_real_attrib_parser[ia]->compile<7>(); + } +} + +int** PlasmaParserHelper::getUserIntDataPtrs () { +#ifdef AMREX_USE_GPU + return m_d_pa_user_int.dataPtr(); +#else + return m_pa_user_int_pinned.dataPtr(); +#endif +} + +amrex::ParticleReal** PlasmaParserHelper::getUserRealDataPtrs () { +#ifdef AMREX_USE_GPU + return m_d_pa_user_real.dataPtr(); +#else + return m_pa_user_real_pinned.dataPtr(); +#endif +} + +amrex::ParserExecutor<7> const* PlasmaParserHelper::getUserIntParserExecData () const { +#ifdef AMREX_USE_GPU + return m_d_user_int_attrib_parserexec.dataPtr(); +#else + return m_wrapper_ptr->m_user_int_attrib_parserexec_pinned.dataPtr(); +#endif +} + +amrex::ParserExecutor<7> const* PlasmaParserHelper::getUserRealParserExecData () const { +#ifdef AMREX_USE_GPU + return m_d_user_real_attrib_parserexec.dataPtr(); +#else + return m_wrapper_ptr->m_user_real_attrib_parserexec_pinned.dataPtr(); +#endif +} diff --git a/Source/Particles/CMakeLists.txt b/Source/Particles/CMakeLists.txt index 67af14ef889..6b434c0a4e1 100644 --- a/Source/Particles/CMakeLists.txt +++ b/Source/Particles/CMakeLists.txt @@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS) warpx_set_suffix_dims(SD ${D}) target_sources(lib_${SD} PRIVATE + AddPlasmaUtilities.cpp MultiParticleContainer.cpp ParticleBoundaries.cpp PhotonParticleContainer.cpp diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 58cbe11a980..69918f69940 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -1,3 +1,4 @@ +CEXE_sources += AddPlasmaUtilities.cpp CEXE_sources += MultiParticleContainer.cpp CEXE_sources += WarpXParticleContainer.cpp CEXE_sources += RigidInjectedParticleContainer.cpp diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 5d9b41b8b75..8102fc96a91 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -392,6 +392,14 @@ public: } protected: + + /* + Finds the box defining the region where refine injection should be used, if that + option is enabled. Currently this only works for numLevels() == 2 and static mesh + refinement. + */ + bool findRefinedInjectionBox (amrex::Box& fine_injection_box, amrex::IntVect& rrfac); + std::string species_name; std::vector> plasma_injectors; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 94a65303cc5..d1a19f06993 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -15,6 +15,7 @@ #include "Initialization/InjectorMomentum.H" #include "Initialization/InjectorPosition.H" #include "MultiParticleContainer.H" +#include "Particles/AddPlasmaUtilities.H" #ifdef WARPX_QED # include "Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H" # include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H" @@ -217,8 +218,7 @@ namespace const GpuArray& pa, long& ip, const bool& do_field_ionization, int* pi #ifdef WARPX_QED - ,const bool& has_quantum_sync, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR - ,const bool& has_breit_wheeler, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_BW + ,const QEDHelper& qed_helper #endif ) noexcept { @@ -227,8 +227,8 @@ namespace } if (do_field_ionization) {pi[ip] = 0;} #ifdef WARPX_QED - if (has_quantum_sync) {p_optical_depth_QSR[ip] = 0._rt;} - if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;} + if (qed_helper.has_quantum_sync) {qed_helper.p_optical_depth_QSR[ip] = 0._rt;} + if (qed_helper.has_breit_wheeler) {qed_helper.p_optical_depth_BW[ip] = 0._rt;} #endif idcpu[ip] = amrex::ParticleIdCpus::Invalid; @@ -964,22 +964,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int amrex::LayoutData* cost = WarpX::getCosts(lev); - const int nlevs = numLevels(); - static bool refine_injection = false; - static Box fine_injection_box; - static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); - // This does not work if the mesh is dynamic. But in that case, we should - // not use refined injected either. We also assume there is only one fine level. - if (WarpX::moving_window_active(WarpX::GetInstance().getistep(0)+1) and WarpX::refine_plasma - and do_continuous_injection and nlevs == 2) - { - refine_injection = true; - fine_injection_box = ParticleBoxArray(1).minimalBox(); - fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits::lowest()/2); - fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits::max()/2); - rrfac = m_gdb->refRatio(0); - fine_injection_box.coarsen(rrfac); - } + Box fine_injection_box; + amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + const bool refine_injection = findRefinedInjectionBox(fine_injection_box, rrfac); InjectorPosition* inj_pos = plasma_injector.getInjectorPosition(); InjectorDensity* inj_rho = plasma_injector.getInjectorDensity(); @@ -995,18 +982,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const bool radially_weighted = plasma_injector.radially_weighted; #endif - - // User-defined integer and real attributes: prepare parsers - const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); - const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_int_attrib_parserexec_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_real_attrib_parserexec_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - user_int_attrib_parserexec_pinned[ia] = m_user_int_attrib_parser[ia]->compile<7>(); - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - user_real_attrib_parserexec_pinned[ia] = m_user_real_attrib_parser[ia]->compile<7>(); - } + auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); + auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); + const PlasmaParserWrapper plasma_parser_wrapper (m_user_int_attribs.size(), + m_user_real_attribs.size(), + m_user_int_attrib_parser, + m_user_real_attrib_parser); MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { @@ -1032,30 +1013,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int RealBox overlap_realbox; Box overlap_box; IntVect shifted; - bool no_overlap = false; - - for (int dir=0; dir= part_realbox.lo(dir) ) { - const Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - // Count the number of cells in this direction in overlap_realbox - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, - int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) - /dx[dir] )) - 1); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - // shifted is exact in non-moving-window direction. That's all we care. - } + const bool no_overlap = find_overlap(tile_realbox, part_realbox, dx, problo, overlap_realbox, overlap_box, shifted); if (no_overlap) { continue; // Go to the next tile } @@ -1110,19 +1068,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } return 0; }; - const int flag_pcount = checker(); - if (flag_pcount == 1) { - pcounts[index] = num_ppc*r; - } else { - pcounts[index] = 0; - } + pcounts[index] = checker() ? num_ppc*r : 0; } -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(k); -#endif -#if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); -#endif }); // Max number of new particles. All of them are created, @@ -1160,40 +1108,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int pa[ia] = soa.GetRealData(ia).data() + old_size; } uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; - // user-defined integer and real attributes - amrex::Gpu::PinnedVector pa_user_int_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector pa_user_real_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - pa_user_int_pinned[ia] = soa.GetIntData(particle_icomps[m_user_int_attribs[ia]]).data() + old_size; - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - pa_user_real_pinned[ia] = soa.GetRealData(particle_comps[m_user_real_attribs[ia]]).data() + old_size; - } -#ifdef AMREX_USE_GPU - // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, - // and them memcpy to device from host - amrex::Gpu::DeviceVector d_pa_user_int(n_user_int_attribs); - amrex::Gpu::DeviceVector d_pa_user_real(n_user_real_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_int_attrib_parserexec(n_user_int_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_real_attrib_parserexec(n_user_real_attribs); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_int_pinned.begin(), - pa_user_int_pinned.end(), d_pa_user_int.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_real_pinned.begin(), - pa_user_real_pinned.end(), d_pa_user_real.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_int_attrib_parserexec_pinned.begin(), - user_int_attrib_parserexec_pinned.end(), d_user_int_attrib_parserexec.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_real_attrib_parserexec_pinned.begin(), - user_real_attrib_parserexec_pinned.end(), d_user_real_attrib_parserexec.begin()); - int** pa_user_int_data = d_pa_user_int.dataPtr(); - ParticleReal** pa_user_real_data = d_pa_user_real.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = d_user_int_attrib_parserexec.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = d_user_real_attrib_parserexec.dataPtr(); -#else - int** pa_user_int_data = pa_user_int_pinned.dataPtr(); - ParticleReal** pa_user_real_data = pa_user_real_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = user_int_attrib_parserexec_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = user_real_attrib_parserexec_pinned.dataPtr(); -#endif + + PlasmaParserHelper plasma_parser_helper (soa, old_size, m_user_int_attribs, m_user_real_attribs, particle_icomps, particle_comps, plasma_parser_wrapper); + int** pa_user_int_data = plasma_parser_helper.getUserIntDataPtrs(); + ParticleReal** pa_user_real_data = plasma_parser_helper.getUserRealDataPtrs(); + amrex::ParserExecutor<7> const* user_int_parserexec_data = plasma_parser_helper.getUserIntParserExecData(); + amrex::ParserExecutor<7> const* user_real_parserexec_data = plasma_parser_helper.getUserRealParserExecData(); int* pi = nullptr; if (do_field_ionization) { @@ -1201,34 +1121,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } #ifdef WARPX_QED - //Pointer to the optical depth component - amrex::ParticleReal* p_optical_depth_QSR = nullptr; - amrex::ParticleReal* p_optical_depth_BW = nullptr; - - // If a QED effect is enabled, the corresponding optical depth - // has to be initialized - const bool loc_has_quantum_sync = has_quantum_sync(); - const bool loc_has_breit_wheeler = has_breit_wheeler(); - if (loc_has_quantum_sync) { - p_optical_depth_QSR = soa.GetRealData( - particle_comps["opticalDepthQSR"]).data() + old_size; - } - if(loc_has_breit_wheeler) { - p_optical_depth_BW = soa.GetRealData( - particle_comps["opticalDepthBW"]).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 = - m_shr_p_qs_engine->build_optical_depth_functor(); - } - if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = - m_shr_p_bw_engine->build_optical_depth_functor(); - } + const QEDHelper qed_helper(soa, old_size, particle_comps, + has_quantum_sync(), has_breit_wheeler(), + m_shr_p_qs_engine, m_shr_p_bw_engine); #endif const bool loc_do_field_ionization = do_field_ionization; @@ -1246,18 +1141,14 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); + amrex::ignore_unused(j,k); const auto index = overlap_box.index(iv); #ifdef WARPX_DIM_RZ Real theta_offset = 0._rt; if (rz_random_theta) { theta_offset = amrex::Random(engine) * 2._rt * MathConst::pi; } #endif - Real scale_fac = 0.0_rt; - if( pcounts[index] != 0) { - amrex::Real const dV = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); - scale_fac = dV/pcounts[index]; - } - + const Real scale_fac = compute_scale_fac_volume(dx, pcounts[index]); for (int i_part = 0; i_part < pcounts[index]; ++i_part) { long ip = poffset[index] + i_part; @@ -1281,8 +1172,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!box_contains) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1319,8 +1209,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!inj_pos->insideBounds(xb, yb, z0)) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1333,8 +1222,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if ( dens < density_min ){ ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1351,8 +1239,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!inj_pos->insideBounds(xb, yb, z0_lab)) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1363,8 +1250,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if ( dens < density_min ){ ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1388,12 +1274,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } #ifdef WARPX_QED - if(loc_has_quantum_sync){ - p_optical_depth_QSR[ip] = quantum_sync_get_opt(engine); + if(qed_helper.has_quantum_sync){ + qed_helper.p_optical_depth_QSR[ip] = qed_helper.quantum_sync_get_opt(engine); } - if(loc_has_breit_wheeler){ - p_optical_depth_BW[ip] = breit_wheeler_get_opt(engine); + if(qed_helper.has_breit_wheeler){ + qed_helper.p_optical_depth_BW[ip] = qed_helper.breit_wheeler_get_opt(engine); } #endif // Initialize user-defined integers with user-defined parser @@ -1481,25 +1367,6 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const auto dx = geom.CellSizeArray(); const auto problo = geom.ProbLoArray(); - Real scale_fac = 0._rt; - // Scale particle weight by the area of the emitting surface, within one cell -#if defined(WARPX_DIM_3D) - scale_fac = dx[0]*dx[1]*dx[2]/dx[plasma_injector.flux_normal_axis]/num_ppc_real; -#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - scale_fac = dx[0]*dx[1]/num_ppc_real; - // When emission is in the r direction, the emitting surface is a cylinder. - // The factor 2*pi*r is added later below. - if (plasma_injector.flux_normal_axis == 0) { scale_fac /= dx[0]; } - // When emission is in the z direction, the emitting surface is an annulus - // The factor 2*pi*r is added later below. - if (plasma_injector.flux_normal_axis == 2) { scale_fac /= dx[1]; } - // When emission is in the theta direction (flux_normal_axis == 1), - // the emitting surface is a rectangle, within the plane of the simulation -#elif defined(WARPX_DIM_1D_Z) - scale_fac = dx[0]/num_ppc_real; - if (plasma_injector.flux_normal_axis == 2) { scale_fac /= dx[0]; } -#endif - amrex::LayoutData* cost = WarpX::getCosts(0); // Create temporary particle container to which particles will be added; @@ -1510,19 +1377,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, for (int ic = 0; ic < NumRuntimeIntComps(); ++ic) { tmp_pc.AddIntComp(false); } tmp_pc.defineAllParticleTiles(); - const int nlevs = numLevels(); - static bool refine_injection = false; - static Box fine_injection_box; - static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); - // This does not work if the mesh is dynamic. But in that case, we should - // not use refined injected either. We also assume there is only one fine level. - if (WarpX::refine_plasma && nlevs == 2) - { - refine_injection = true; - fine_injection_box = ParticleBoxArray(1).minimalBox(); - rrfac = m_gdb->refRatio(0); - fine_injection_box.coarsen(rrfac); - } + Box fine_injection_box; + amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + const bool refine_injection = findRefinedInjectionBox(fine_injection_box, rrfac); InjectorPosition* flux_pos = plasma_injector.getInjectorFluxPosition(); InjectorFlux* inj_flux = plasma_injector.getInjectorFlux(); @@ -1536,6 +1393,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const bool radially_weighted = plasma_injector.radially_weighted; #endif + auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); + auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); + const PlasmaParserWrapper plasma_parser_wrapper (m_user_int_attribs.size(), + m_user_real_attribs.size(), + m_user_int_attrib_parser, + m_user_real_attrib_parser); + MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { info.EnableTiling(tile_size); @@ -1560,61 +1424,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, RealBox overlap_realbox; Box overlap_box; IntVect shifted; - bool no_overlap = false; - - for (int dir=0; dir 0) { - if (plasma_injector.surface_flux_pos < tile_realbox.lo(dir) || - plasma_injector.surface_flux_pos >= tile_realbox.hi(dir)) { - no_overlap = true; - break; - } - } else { - if (plasma_injector.surface_flux_pos <= tile_realbox.lo(dir) || - plasma_injector.surface_flux_pos > tile_realbox.hi(dir)) { - no_overlap = true; - break; - } - } - overlap_realbox.setLo( dir, plasma_injector.surface_flux_pos ); - overlap_realbox.setHi( dir, plasma_injector.surface_flux_pos ); - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, 0 ); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - } else { - if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { - const Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { - const Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - // Count the number of cells in this direction in overlap_realbox - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, - int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) - /dx[dir] )) - 1); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - // shifted is exact in non-moving-window direction. That's all we care. - } - } + const bool no_overlap = find_overlap_flux(tile_realbox, part_realbox, dx, problo, plasma_injector, overlap_realbox, overlap_box, shifted); if (no_overlap) { continue; // Go to the next tile } @@ -1632,6 +1442,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, Gpu::DeviceVector offset(overlap_box.numPts()); auto *pcounts = counts.data(); const amrex::IntVect lrrfac = rrfac; + const int flux_normal_axis = plasma_injector.flux_normal_axis; Box fine_overlap_box; // default Box is NOT ok(). if (refine_injection) { fine_overlap_box = overlap_box & amrex::shift(fine_injection_box, -shifted); @@ -1649,17 +1460,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, auto index = overlap_box.index(iv); int r; if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { - r = AMREX_D_TERM(lrrfac[0],*lrrfac[1],*lrrfac[2]); + r = compute_area_weights(lrrfac, flux_normal_axis); } else { r = 1; } pcounts[index] = num_ppc_int*r; } -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(k); -#elif defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); -#endif }); // Max number of new particles. All of them are created, @@ -1694,46 +1501,11 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; - // user-defined integer and real attributes - const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); - const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); - amrex::Gpu::PinnedVector pa_user_int_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector pa_user_real_pinned(n_user_real_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_int_attrib_parserexec_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_real_attrib_parserexec_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - pa_user_int_pinned[ia] = soa.GetIntData(particle_icomps[m_user_int_attribs[ia]]).data() + old_size; - user_int_attrib_parserexec_pinned[ia] = m_user_int_attrib_parser[ia]->compile<7>(); - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - pa_user_real_pinned[ia] = soa.GetRealData(particle_comps[m_user_real_attribs[ia]]).data() + old_size; - user_real_attrib_parserexec_pinned[ia] = m_user_real_attrib_parser[ia]->compile<7>(); - } -#ifdef AMREX_USE_GPU - // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, - // and them memcpy to device from host - amrex::Gpu::DeviceVector d_pa_user_int(n_user_int_attribs); - amrex::Gpu::DeviceVector d_pa_user_real(n_user_real_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_int_attrib_parserexec(n_user_int_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_real_attrib_parserexec(n_user_real_attribs); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_int_pinned.begin(), - pa_user_int_pinned.end(), d_pa_user_int.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_real_pinned.begin(), - pa_user_real_pinned.end(), d_pa_user_real.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_int_attrib_parserexec_pinned.begin(), - user_int_attrib_parserexec_pinned.end(), d_user_int_attrib_parserexec.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_real_attrib_parserexec_pinned.begin(), - user_real_attrib_parserexec_pinned.end(), d_user_real_attrib_parserexec.begin()); - int** pa_user_int_data = d_pa_user_int.dataPtr(); - ParticleReal** pa_user_real_data = d_pa_user_real.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = d_user_int_attrib_parserexec.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = d_user_real_attrib_parserexec.dataPtr(); -#else - int** pa_user_int_data = pa_user_int_pinned.dataPtr(); - ParticleReal** pa_user_real_data = pa_user_real_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = user_int_attrib_parserexec_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = user_real_attrib_parserexec_pinned.dataPtr(); -#endif + PlasmaParserHelper plasma_parser_helper (soa, old_size, m_user_int_attribs, m_user_real_attribs, particle_icomps, particle_comps, plasma_parser_wrapper); + int** pa_user_int_data = plasma_parser_helper.getUserIntDataPtrs(); + ParticleReal** pa_user_real_data = plasma_parser_helper.getUserRealDataPtrs(); + amrex::ParserExecutor<7> const* user_int_parserexec_data = plasma_parser_helper.getUserIntParserExecData(); + amrex::ParserExecutor<7> const* user_real_parserexec_data = plasma_parser_helper.getUserRealParserExecData(); int* p_ion_level = nullptr; if (do_field_ionization) { @@ -1741,34 +1513,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef WARPX_QED - //Pointer to the optical depth component - amrex::ParticleReal* p_optical_depth_QSR = nullptr; - amrex::ParticleReal* p_optical_depth_BW = nullptr; - - // If a QED effect is enabled, the corresponding optical depth - // has to be initialized - const bool loc_has_quantum_sync = has_quantum_sync(); - const bool loc_has_breit_wheeler = has_breit_wheeler(); - if (loc_has_quantum_sync) { - p_optical_depth_QSR = soa.GetRealData( - particle_comps["opticalDepthQSR"]).data() + old_size; - } - if(loc_has_breit_wheeler) { - p_optical_depth_BW = soa.GetRealData( - particle_comps["opticalDepthBW"]).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 = - m_shr_p_qs_engine->build_optical_depth_functor(); - } - if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = - m_shr_p_bw_engine->build_optical_depth_functor(); - } + const QEDHelper qed_helper(soa, old_size, particle_comps, + has_quantum_sync(), has_breit_wheeler(), + m_shr_p_qs_engine, m_shr_p_bw_engine); #endif const bool loc_do_field_ionization = do_field_ionization; @@ -1787,6 +1534,24 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); const auto index = overlap_box.index(iv); + + Real scale_fac = compute_scale_fac_area(dx, num_ppc_real, flux_normal_axis); + + auto lo = getCellCoords(overlap_corner, dx, {0._rt, 0._rt, 0._rt}, iv); + auto hi = getCellCoords(overlap_corner, dx, {1._rt, 1._rt, 1._rt}, iv); + + if (flux_pos->overlapsWith(lo, hi)) + { + int r; + if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { + r = compute_area_weights(lrrfac, flux_normal_axis); + } else { + r = 1; + } + scale_fac /= r; + } + amrex::ignore_unused(j,k); + for (int i_part = 0; i_part < pcounts[index]; ++i_part) { const long ip = poffset[index] + i_part; @@ -1878,14 +1643,15 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef WARPX_QED - if (loc_has_quantum_sync) { - p_optical_depth_QSR[ip] = quantum_sync_get_opt(engine); + if(qed_helper.has_quantum_sync){ + qed_helper.p_optical_depth_QSR[ip] = qed_helper.quantum_sync_get_opt(engine); } - if(loc_has_breit_wheeler){ - p_optical_depth_BW[ip] = breit_wheeler_get_opt(engine); + if(qed_helper.has_breit_wheeler){ + qed_helper.p_optical_depth_BW[ip] = qed_helper.breit_wheeler_get_opt(engine); } #endif + // Initialize user-defined integers with user-defined parser for (int ia = 0; ia < n_user_int_attribs; ++ia) { pa_user_int_data[ia][ip] = static_cast(user_int_parserexec_data[ia](pos.x, pos.y, pos.z, u.x, u.y, u.z, t)); @@ -1967,26 +1733,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // are in the right tile.) tmp_pc.Redistribute(); - // Add the particles to the current container, tile by tile - for (int lev=0; levaddParticles(tmp_pc, true); } void @@ -3412,6 +3160,28 @@ void PhysicalParticleContainer::resample (const int timestep, const bool verbose WARPX_PROFILE_VAR_STOP(blp_resample_actual); } +bool +PhysicalParticleContainer::findRefinedInjectionBox (amrex::Box& a_fine_injection_box, amrex::IntVect& a_rrfac) +{ + WARPX_PROFILE("PhysicalParticleContainer::findRefinedInjectionBox"); + + // This does not work if the mesh is dynamic. But in that case, we should + // not use refined injected either. We also assume there is only one fine level. + static bool refine_injection = false; + static Box fine_injection_box; + static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + if (!refine_injection and WarpX::moving_window_active(WarpX::GetInstance().getistep(0)+1) and WarpX::refine_plasma and do_continuous_injection and numLevels() == 2) { + refine_injection = true; + fine_injection_box = ParticleBoxArray(1).minimalBox(); + fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits::lowest()/2); + fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits::max()/2); + rrfac = m_gdb->refRatio(0); + fine_injection_box.coarsen(rrfac); + } + a_fine_injection_box = fine_injection_box; + a_rrfac = rrfac; + return refine_injection; +} #ifdef WARPX_QED