From 314acff8ddd7ca7af3af9a65c43908acc266533f Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 17 Mar 2020 15:25:42 +0100 Subject: Add QED particle creaction routines (#698) * Initial work to add back QED particle generation * Work in progress: port old QED routines * Add two distinct CopyFuncs * modified getMFItInfo and CopyFunc (not working) * bugfixing & work to add back QED particle creation routines * bugfixing * added back quantum photon emission * bugfixing * bugfixing * added back pair generation (still some bugs in photon emission) * removed unwanted check * bugfixing * bugfixing * bugfixing * Moved QED folder * added comments + some refactoring * added comments * remove some virtual functions to make lgtm happy * updated tests * added PhysicalParticleType * bugfixing * added copyright * improved comments * improved comments * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * moved inclusion of QEDInternals folder * moved some inclusion directives between Make files * moved some inclusion directives between Make files (forgot to add a file) * Update Source/Particles/PhysicalParticleContainer.H Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.H Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.H Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.H Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.H Co-Authored-By: Axel Huebl * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: Axel Huebl * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.H Co-Authored-By: Axel Huebl * corrected alignment * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: Axel Huebl * removed some unnecessary amrex:: * add missing comment * Replaced BL_PROFILE with WARPX_PROFILE * bugfixing and making some variables const * removed some moves * removed some moves * started to change tau into optical_depth_BW or optical_depth_QSR * Using initialization policy to initialize optical depth * bugfixing * forgot to add a file * fixed bug * Revert "fixed bug" This reverts commit a3fb98d10cc30327635aeaa71451a05ca2229ff4. * Update Source/Particles/ElementaryProcess/QEDPairGeneration.H Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Added doQEDEvents to OneStep_sub1 * add a bunch of const * add _rt suffix * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * added path to included files * Introduced a templated AmIA function using physical_species * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: MaxThevenet * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: MaxThevenet * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: MaxThevenet * added paths to included headers * updated documentation * updated examples * bugfixing * bugfixing * fixing examples * fixed example * fixed example * correct a misprint in error message * fixed issue related to 1./mass for photons * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/MultiParticleContainer.cpp Co-Authored-By: MaxThevenet * Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H Co-Authored-By: MaxThevenet * make the use of energy_threshold more transparent * remove unnecessary checks * bigfixing * added comment * separate checks for QED processes * added a CheckIonizationProductSpecies for consistency * bugfixing * now using a new variable for photon energy creation threshold * removed unwanted comment * added option to set a user-defined threshold for photon creation * bugfixing * updated documentation * updated example to include new option * updated doc * fixed merge conflict * correct bug in example * reorganized function Co-authored-by: Axel Huebl Co-authored-by: MaxThevenet --- Source/Particles/MultiParticleContainer.cpp | 232 +++++++++++++++++++++++++--- 1 file changed, 210 insertions(+), 22 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 5b6cfedc0..12b0a9d61 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -9,6 +9,7 @@ * License: BSD-3-Clause-LBNL */ #include "MultiParticleContainer.H" +#include "SpeciesPhysicalProperties.H" #include "Utils/WarpXUtil.H" #include "WarpX.H" @@ -267,7 +268,10 @@ MultiParticleContainer::InitData () // This is used for ionization and pair creation processes. mapSpeciesProduct(); + CheckIonizationProductSpecies(); + #ifdef WARPX_QED + CheckQEDProductSpecies(); InitQED(); #endif @@ -534,13 +538,30 @@ MultiParticleContainer::mapSpeciesProduct () // pc->ionization_product_name and store its ID into // pc->ionization_product. if (pc->do_field_ionization){ - int i_product = getSpeciesID(pc->ionization_product_name); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - i != i_product, - "ERROR: ionization product cannot be the same species"); + const int i_product = getSpeciesID(pc->ionization_product_name); pc->ionization_product = i_product; } + +#ifdef WARPX_QED + if (pc->has_breit_wheeler()){ + const int i_product_ele = getSpeciesID( + pc->m_qed_breit_wheeler_ele_product_name); + pc->m_qed_breit_wheeler_ele_product = i_product_ele; + + const int i_product_pos = getSpeciesID( + pc->m_qed_breit_wheeler_pos_product_name); + pc->m_qed_breit_wheeler_pos_product = i_product_pos; + } + + if(pc->has_quantum_sync()){ + const int i_product_phot = getSpeciesID( + pc->m_qed_quantum_sync_phot_product_name); + pc->m_qed_quantum_sync_phot_product = i_product_phot; + } +#endif + } + } /* \brief Given a species name, return its ID. @@ -594,7 +615,7 @@ MultiParticleContainer::doFieldIonization () for (int lev = 0; lev <= pc_source->finestLevel(); ++lev) { - auto info = getMFItInfo(*pc_source, *pc_product); + const auto info = getMFItInfo(*pc_source, *pc_product); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -604,8 +625,8 @@ MultiParticleContainer::doFieldIonization () auto& src_tile = pc_source ->ParticlesAt(lev, mfi); auto& dst_tile = pc_product->ParticlesAt(lev, mfi); - auto np_dst = dst_tile.numParticles(); - auto num_added = filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst, + const auto np_dst = dst_tile.numParticles(); + const auto num_added = filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst, Filter, Copy, Transform); setNewParticleIDs(dst_tile, np_dst, num_added); @@ -648,23 +669,15 @@ MultiParticleContainer::doCoulombCollisions () } } -MFItInfo MultiParticleContainer::getMFItInfo (const WarpXParticleContainer& pc_src, - const WarpXParticleContainer& pc_dst) const noexcept +void MultiParticleContainer::CheckIonizationProductSpecies() { - MFItInfo info; - - if (pc_src.do_tiling && Gpu::notInLaunchRegion()) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(pc_dst.do_tiling, - "For ionization, either all or none of the " - "particle species must use tiling."); - info.EnableTiling(pc_src.tile_size); + for (int i=0; ido_field_ionization){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + i != allcontainers[i]->ionization_product, + "ERROR: ionization product cannot be the same species"); + } } - -#ifdef _OPENMP - info.SetDynamic(true); -#endif - - return info; } #ifdef WARPX_QED @@ -701,6 +714,17 @@ void MultiParticleContainer::InitQuantumSync () { std::string lookup_table_mode; ParmParse pp("qed_qs"); + + //If specified, use a user-defined energy threshold for photon creaction + ParticleReal temp; + if(pp.query("photon_creation_energy_threshold", temp)){ + m_quantum_sync_photon_creation_energy_threshold = temp; + } + else{ + amrex::Print() << "Using default value (2*me*c^2)" << + " for photon energy creaction threshold \n" ; + } + pp.query("lookup_table_mode", lookup_table_mode); if(lookup_table_mode.empty()){ amrex::Abort("Quantum Synchrotron table mode should be provided"); @@ -952,4 +976,168 @@ MultiParticleContainer::BreitWheelerGenerateTable () m_shr_p_bw_engine->init_lookup_tables_from_raw_data(table_data); } } + +void MultiParticleContainer::doQedEvents() +{ + WARPX_PROFILE("MPC::doQedEvents"); + + doQedBreitWheeler(); + doQedQuantumSync(); +} + +void MultiParticleContainer::doQedBreitWheeler() +{ + WARPX_PROFILE("MPC::doQedBreitWheeler"); + + // Loop over all species. + // Photons undergoing Breit Wheeler process create electrons + // in pc_product_ele and positrons in pc_product_pos + + for (auto& pc_source : allcontainers){ + if(!pc_source->has_breit_wheeler()) continue; + + // Get product species + auto& pc_product_ele = + allcontainers[pc_source->m_qed_breit_wheeler_ele_product]; + auto& pc_product_pos = + allcontainers[pc_source->m_qed_breit_wheeler_pos_product]; + + SmartCopyFactory copy_factory_ele(*pc_source, *pc_product_ele); + SmartCopyFactory copy_factory_pos(*pc_source, *pc_product_pos); + auto phys_pc_ptr = static_cast(pc_source.get()); + + const auto Filter = phys_pc_ptr->getPairGenerationFilterFunc(); + const auto CopyEle = copy_factory_ele.getSmartCopy(); + const auto CopyPos = copy_factory_pos.getSmartCopy(); + + const auto pair_gen_functor = m_shr_p_bw_engine->build_pair_functor(); + auto Transform = PairGenerationTransformFunc(pair_gen_functor); + + pc_source ->defineAllParticleTiles(); + pc_product_pos->defineAllParticleTiles(); + pc_product_ele->defineAllParticleTiles(); + + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev) + { + const auto info = getMFItInfo(*pc_source, *pc_product_ele, *pc_product_pos); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + { + auto& src_tile = pc_source->ParticlesAt(lev, mfi); + auto& dst_ele_tile = pc_product_ele->ParticlesAt(lev, mfi); + auto& dst_pos_tile = pc_product_pos->ParticlesAt(lev, mfi); + + const auto np_dst_ele = dst_ele_tile.numParticles(); + const auto np_dst_pos = dst_pos_tile.numParticles(); + const auto num_added = filterCopyTransformParticles<1>( + dst_ele_tile, dst_pos_tile, + src_tile, np_dst_ele, np_dst_pos, + Filter, CopyEle, CopyPos, Transform); + + setNewParticleIDs(dst_ele_tile, np_dst_ele, num_added); + setNewParticleIDs(dst_pos_tile, np_dst_pos, num_added); + } + } + } +} + +void MultiParticleContainer::doQedQuantumSync() +{ + WARPX_PROFILE("MPC::doQedEvents::doQedQuantumSync"); + + // Loop over all species. + // Electrons or positrons undergoing Quantum photon emission process + // create photons in pc_product_phot + + for (auto& pc_source : allcontainers){ + if(!pc_source->has_quantum_sync()){ continue; } + + // Get product species + auto& pc_product_phot = + allcontainers[pc_source->m_qed_quantum_sync_phot_product]; + + SmartCopyFactory copy_factory_phot(*pc_source, *pc_product_phot); + auto phys_pc_ptr = + static_cast(pc_source.get()); + + const auto Filter = phys_pc_ptr->getPhotonEmissionFilterFunc(); + const auto CopyPhot = copy_factory_phot.getSmartCopy(); + + auto Transform = PhotonEmissionTransformFunc( + m_shr_p_qs_engine->build_optical_depth_functor(), + pc_source->particle_runtime_comps["optical_depth_QSR"], + m_shr_p_qs_engine->build_phot_em_functor()); + + pc_source ->defineAllParticleTiles(); + pc_product_phot->defineAllParticleTiles(); + + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev) + { + const auto info = getMFItInfo(*pc_source, *pc_product_phot); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + { + auto& src_tile = pc_source->ParticlesAt(lev, mfi); + auto& dst_tile = pc_product_phot->ParticlesAt(lev, mfi); + + const auto np_dst = dst_tile.numParticles(); + + const auto num_added = + filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst, + Filter, CopyPhot, Transform); + + setNewParticleIDs(dst_tile, np_dst, num_added); + + cleanLowEnergyPhotons( + dst_tile, np_dst, num_added, + m_quantum_sync_photon_creation_energy_threshold); + + } + } + } + +} + +void MultiParticleContainer::CheckQEDProductSpecies() +{ + for (int i=0; ihas_breit_wheeler()){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + i != pc->m_qed_breit_wheeler_ele_product, + "ERROR: Breit Wheeler product cannot be the same species"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + i != pc->m_qed_breit_wheeler_pos_product, + "ERROR: Breit Wheeler product cannot be the same species"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + allcontainers[pc->m_qed_breit_wheeler_ele_product]-> + AmIA() + && + allcontainers[pc->m_qed_breit_wheeler_pos_product]-> + AmIA(), + "ERROR: Breit Wheeler product species are of wrong type"); + } + + if(pc->has_quantum_sync()){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + i != pc->m_qed_quantum_sync_phot_product, + "ERROR: Quantum Synchrotron product cannot be the same species"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + allcontainers[pc->m_qed_quantum_sync_phot_product]-> + AmIA(), + "ERROR: Quantum Synchrotron product species is of wrong type"); + } + } +} + + #endif -- cgit v1.2.3