Add Schwinger process (#784)

* 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.

* Define doQEDSchwinger function

* Read input parameters for Schwinger process

* 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

* Add skeleton of doQEDSchwinger function

* Linked to PICSAR to calculate pair production rate

* Added missing header file

* fix conflicts in QEDPhotonEmission.H

* Written first version of FilterCreateTransformFromFAB.H

* Cleanup some useless comments

* Update Source/Particles/ElementaryProcess/QEDInternals/SchwingerProcessWrapper.H

Co-Authored-By: Luca Fedeli <luca.fedeli.88@gmail.com>

* Update Source/Particles/ParticleCreation/FilterCreateTransformFromFAB.H

Co-Authored-By: Luca Fedeli <luca.fedeli.88@gmail.com>

* Update Source/Particles/ParticleCreation/FilterCreateTransformFromFAB.H

Co-Authored-By: Luca Fedeli <luca.fedeli.88@gmail.com>

* Minor revisions and improvements

* Remove trailing white spaces

* Write filter function

* Remove print

* Some debugging and cleaning

* Write tranform function

* remove transform_dummy

* Added some tests

* Remove EOL whitespaces

* update prepare_file_travis.py

* Should fix error in automated tests

* Make the tests run in parallel

* Put path relative to Source/ in includes

* update include path

* Actually resolve conflicts

* put dVdt and weight_index as members of filter and transform functions

* a bit of debugging on GPU

* Add comments and documentation

* Fix typos

* Add assert for single precision

* Update Docs/source/running_cpp/parameters.rst

Co-Authored-By: Luca Fedeli <luca.fedeli.88@gmail.com>

* Update Source/Particles/MultiParticleContainer.H

Co-Authored-By: Luca Fedeli <luca.fedeli.88@gmail.com>

* Update Source/Particles/ParticleCreation/FilterCreateTransformFromFAB.H

Co-Authored-By: Luca Fedeli <luca.fedeli.88@gmail.com>

* Add profiler and minor modifications

* Fix typo

* Update asserts so that module works with momentum conserving algo

* update assert if particles are created in analysis script

* Apply suggestions from code review

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

* Add comments to FilterCopyTransform and FilterCreateTransform functions

* Add const and comment

* Update Source/Particles/ParticleCreation/FilterCreateTransformFromFAB.H

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>

* Update test for new diags and remove AllContainerType

* Update Regression/WarpX-tests.ini

* update test

* update tests

* Split Schwinger test into 4 separate tests

* Change name of analysis script

* update tests

* Combine all analysis scripts into a single one

Co-authored-by: Luca Fedeli <luca.fedeli@cea.fr>
Co-authored-by: Luca Fedeli <luca.fedeli.88@gmail.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: MaxThevenet <mthevenet@lbl.gov>
Co-authored-by: Neil <nzaim@iram-fe-003981.extra.cea.fr>

1 files changed, 148 insertions, 2 deletions

diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index a3b5d256e..123771c84 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -1,7 +1,8 @@
 /* Copyright 2019-2020 Andrew Myers, Ann Almgren, Axel Huebl
  * David Grote, Jean-Luc Vay, Luca Fedeli
- * Mathieu Lobet, Maxence Thevenet, Remi Lehe
- * Revathi Jambunathan, Weiqun Zhang, Yinjian Zhao
+ * Mathieu Lobet, Maxence Thevenet, Neil Zaim
+ * Remi Lehe, Revathi Jambunathan, Weiqun Zhang
+ * Yinjian Zhao
  *
  *
  * This file is part of WarpX.
@@ -12,6 +13,11 @@
 #include "SpeciesPhysicalProperties.H"
 #include "Utils/WarpXUtil.H"
 #include "WarpX.H"
+#ifdef WARPX_QED
+    #include "Particles/ElementaryProcess/QEDInternals/SchwingerProcessWrapper.H"
+    #include "Particles/ElementaryProcess/QEDSchwingerProcess.H"
+    #include "Particles/ParticleCreation/FilterCreateTransformFromFAB.H"
+#endif
 
 #include <AMReX_Vector.H>
 
@@ -258,6 +264,21 @@ MultiParticleContainer::ReadParameters ()
             AMREX_ALWAYS_ASSERT(lasers_names.size() == nlasers);
         }
 
+#ifdef WARPX_QED
+        ParmParse ppw("warpx");
+        ppw.query("do_qed_schwinger", m_do_qed_schwinger);
+
+        if (m_do_qed_schwinger) {
+            ParmParse ppq("qed_schwinger");
+            ppq.get("ele_product_species", m_qed_schwinger_ele_product_name);
+            ppq.get("pos_product_species", m_qed_schwinger_pos_product_name);
+#if (AMREX_SPACEDIM == 2)
+            ppq.get("y_size",m_qed_schwinger_y_size);
+#endif
+            ppq.query("threshold_poisson_gaussian",
+                      m_qed_schwinger_threshold_poisson_gaussian);
+        }
+#endif
         initialized = true;
     }
 }
@@ -576,6 +597,14 @@ MultiParticleContainer::mapSpeciesProduct ()
 
     }
 
+#ifdef WARPX_QED
+    if (m_do_qed_schwinger) {
+    m_qed_schwinger_ele_product =
+        getSpeciesID(m_qed_schwinger_ele_product_name);
+    m_qed_schwinger_pos_product =
+        getSpeciesID(m_qed_schwinger_pos_product_name);
+    }
+#endif
 }
 
 /* \brief Given a species name, return its ID.
@@ -991,6 +1020,112 @@ MultiParticleContainer::BreitWheelerGenerateTable ()
     }
 }
 
+void
+MultiParticleContainer::doQEDSchwinger ()
+{
+    WARPX_PROFILE("MPC::doQEDSchwinger");
+
+    if (!m_do_qed_schwinger) {return;}
+
+    auto & warpx = WarpX::GetInstance();
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(warpx.do_nodal ||
+       warpx.field_gathering_algo == GatheringAlgo::MomentumConserving,
+          "ERROR: Schwinger process only implemented for warpx.do_nodal = 1"
+                                 "or algo.field_gathering = momentum-conserving");
+
+    const int level_0 = 0;
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(warpx.maxLevel() == level_0,
+        "ERROR: Schwinger process not implemented with mesh refinement");
+
+#ifdef WARPX_DIM_RZ
+    amrex::Abort("Schwinger process not implemented in rz geometry");
+#endif
+
+#ifdef AMREX_USE_FLOAT
+    amrex::Abort("Schwinger process not implemented in single precision");
+#endif
+
+// Get cell volume multiplied by temporal step. In 2D the transverse size is
+// chosen by the user in the input file.
+    amrex::Geometry const & geom = warpx.Geom(level_0);
+    auto domain_box = geom.Domain();
+#if (AMREX_SPACEDIM == 2)
+    const auto dVdt = geom.CellSize(0) * geom.CellSize(1)
+        * m_qed_schwinger_y_size * warpx.getdt(level_0);
+#elif (AMREX_SPACEDIM == 3)
+    const auto dVdt = geom.CellSize(0) * geom.CellSize(1)
+        * geom.CellSize(2) * warpx.getdt(level_0);
+#endif
+
+    auto& pc_product_ele =
+            allcontainers[m_qed_schwinger_ele_product];
+    auto& pc_product_pos =
+            allcontainers[m_qed_schwinger_pos_product];
+
+    const MultiFab & Ex = warpx.getEfield(level_0,0);
+    const MultiFab & Ey = warpx.getEfield(level_0,1);
+    const MultiFab & Ez = warpx.getEfield(level_0,2);
+    const MultiFab & Bx = warpx.getBfield(level_0,0);
+    const MultiFab & By = warpx.getBfield(level_0,1);
+    const MultiFab & Bz = warpx.getBfield(level_0,2);
+
+    MFItInfo info;
+    if (TilingIfNotGPU()) {
+        info.EnableTiling();
+    }
+#ifdef _OPENMP
+    info.SetDynamic(WarpX::do_dynamic_scheduling);
+#pragma omp parallel
+#endif
+
+    for (MFIter mfi(Ex, info); mfi.isValid(); ++mfi )
+    {
+        // Make the box cell centered to avoid creating particles twice on the tile edges
+        const Box& box = enclosedCells(mfi.tilebox());
+
+        const auto& arrEx = Ex[mfi].array();
+        const auto& arrEy = Ey[mfi].array();
+        const auto& arrEz = Ez[mfi].array();
+        const auto& arrBx = Bx[mfi].array();
+        const auto& arrBy = By[mfi].array();
+        const auto& arrBz = Bz[mfi].array();
+
+        const Array4<const amrex::Real> array_EMFAB [] = {arrEx,arrEy,arrEz,
+                                           arrBx,arrBy,arrBz};
+
+        pc_product_ele->defineAllParticleTiles();
+        pc_product_pos->defineAllParticleTiles();
+
+        auto& dst_ele_tile = pc_product_ele->ParticlesAt(level_0, mfi);
+        auto& dst_pos_tile = pc_product_pos->ParticlesAt(level_0, mfi);
+
+        const auto np_ele_dst = dst_ele_tile.numParticles();
+        const auto np_pos_dst = dst_pos_tile.numParticles();
+
+        const auto Filter  = SchwingerFilterFunc{
+                              m_qed_schwinger_threshold_poisson_gaussian,dVdt};
+
+        const SmartCreateFactory create_factory_ele(*pc_product_ele);
+        const SmartCreateFactory create_factory_pos(*pc_product_pos);
+        const auto CreateEle = create_factory_ele.getSmartCreate();
+        const auto CreatePos = create_factory_pos.getSmartCreate();
+
+        const auto Transform = SchwingerTransformFunc{m_qed_schwinger_y_size,
+                            ParticleStringNames::to_index.find("w")->second};
+
+        const auto num_added = filterCreateTransformFromFAB<1>( dst_ele_tile,
+                              dst_pos_tile, box, array_EMFAB, np_ele_dst,
+                               np_pos_dst,Filter, CreateEle, CreatePos,
+                                Transform);
+
+        setNewParticleIDs(dst_ele_tile, np_ele_dst, num_added);
+        setNewParticleIDs(dst_pos_tile, np_pos_dst, num_added);
+
+    }
+}
+
 void MultiParticleContainer::doQedEvents()
 {
     WARPX_PROFILE("MPC::doQedEvents");
@@ -1151,6 +1286,17 @@ void MultiParticleContainer::CheckQEDProductSpecies()
                 "ERROR: Quantum Synchrotron product species is of wrong type");
         }
     }
+
+    if (m_do_qed_schwinger) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                allcontainers[m_qed_schwinger_ele_product]->
+                    AmIA<PhysicalSpecies::electron>()
+                &&
+                allcontainers[m_qed_schwinger_pos_product]->
+                    AmIA<PhysicalSpecies::positron>(),
+                "ERROR: Schwinger process product species are of wrong type");
+    }
+
 }
 
 #endif