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 <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* 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 <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* corrected alignment

* Update Source/Particles/ElementaryProcess/QEDPairGeneration.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* 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 <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Added doQEDEvents to OneStep_sub1

* add a bunch of const

* add _rt suffix

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* added path to included files

* Introduced a templated AmIA<something> function using physical_species

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* 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 <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/MultiParticleContainer.cpp

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* Update Source/Particles/ElementaryProcess/QEDPhotonEmission.H

Co-Authored-By: MaxThevenet <mthevenet@lbl.gov>

* 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 <axel.huebl@plasma.ninja>
Co-authored-by: MaxThevenet <mthevenet@lbl.gov>

1 files changed, 210 insertions, 22 deletions

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; i<nspecies; i++){
+        if (allcontainers[i]->do_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<PhysicalParticleContainer*>(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<PhysicalParticleContainer*>(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; i<nspecies; i++){
+        const auto& pc = allcontainers[i];
+        if (pc->has_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<PhysicalSpecies::electron>()
+                &&
+                allcontainers[pc->m_qed_breit_wheeler_pos_product]->
+                    AmIA<PhysicalSpecies::positron>(),
+                "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<PhysicalSpecies::photon>(),
+                "ERROR: Quantum Synchrotron product species is of wrong type");
+        }
+    }
+}
+
+
 #endif