Remove persistent E+B (#1050)

* add functor for doing the tmp particles copy for the back-transformed diagnosti

* merge the particle push options into one kernel

* EOL

* fix assertion

* add a FieldGatherandPushPX method to PhysicalParticleContainer

* handle offset in copyAttribs

* allow this functor to be constructed even it we aren't doing the back transformed diagnostics

* EOL

* update the overloads of PushPX for the Photon and RigidInjected ParticleContainers

* function for dispatching the right field gather

* init this val to 0.0

* fix some typos

* handle scaling the fields for rigid injection

* EOL

* don't need to get pointers to E and B arrays in PushPX any more.

* actually I can't remove these yet

* EOL

* variable order bug

* move the QED stuff to the proper place

* EOL

* make sure we don't build these functors unless the runtime options are toggled

* EOL

* perform the field gather prior to the photon particle push

* remove E and B components and FieldGather methods. Reimplement PushP for rigid injected and physical particles

* update ionization to do field gather inline

* remove E and B from the particle diagnostics

* don't write E or B in these tests for particles

* add missing files

* remove EB from the Regtest ini file too

* no need to do this twice

* important typo

* also do the gather inline for the QED processes that need to

* move these sources inside ifdef for QED

* fix bug in RZ

* remove some fields from the Python tests.

* remove all particle E and B comps from json benchmarks

* don't assert that Ey is the langmuir output

* remove uy from this output

* update test

* restore the mesh fields I turned off by mistake

* turn off field IO for a few python tests I missed

* fix typo

* reset Langmuir_multi benchmark

* update Langmuir_multi_nodal benchmark

* update single precision langmuir bench

* update psatd single precision languir one too

* also do ionization_lab

* finally, ionizaiton_boost

* update benchmarks_json/Langmuir_multi_psatd.json

* update benchmarks_json/Langmuir_multi_psatd_current_correction.json

* update benchmarks_json/Langmuir_multi_psatd_momentum_conserving.json

* update benchmarks_json/Langmuir_multi_psatd_nodal.json

* remove the particle E and B from the choices in the docs

* fix offset bug

* also add the Gather subdirectory

* Update Source/WarpX.H

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

* add docstring for LowerCornerWithGalilean

* add new source files to CMakeLists.txt

* also need to update the GPU regression tests

* update the name of the output file for this python test

* remove field gather call from FieldDiagnostics

* fix typo in docstring

* init fields to 0

* add docstring to the CopyParticleAttribs constructor

* some explicit amrex::namepace

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

10 files changed, 415 insertions, 48 deletions

diff --git a/Source/Particles/ElementaryProcess/CMakeLists.txt b/Source/Particles/ElementaryProcess/CMakeLists.txt
index 1a4550a94..bc3bff031 100644
--- a/Source/Particles/ElementaryProcess/CMakeLists.txt
+++ b/Source/Particles/ElementaryProcess/CMakeLists.txt
@@ -1,3 +1,13 @@
+target_sources(WarpX
+  PRIVATE
+    Ionization.cpp
+)
+
 if(WarpX_HAVE_QED)
+    target_sources(WarpX
+      PRIVATE
+        QEDPairGeneration.cpp
+        QEDPhotonEmission.cpp
+    )
     add_subdirectory(QEDInternals)
 endif()
diff --git a/Source/Particles/ElementaryProcess/Ionization.H b/Source/Particles/ElementaryProcess/Ionization.H
index 6cf30bd4d..6275c7af4 100644
--- a/Source/Particles/ElementaryProcess/Ionization.H
+++ b/Source/Particles/ElementaryProcess/Ionization.H
@@ -10,38 +10,96 @@
 
 #include "Utils/WarpXConst.H"
 #include "Particles/WarpXParticleContainer.H"
+#include "Particles/Gather/GetExternalFields.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 
 struct IonizationFilterFunc
 {
-    const amrex::Real* const AMREX_RESTRICT m_ionization_energies;
-    const amrex::Real* const AMREX_RESTRICT m_adk_prefactor;
-    const amrex::Real* const AMREX_RESTRICT m_adk_exp_prefactor;
-    const amrex::Real* const AMREX_RESTRICT m_adk_power;
+    const amrex::Real* AMREX_RESTRICT m_ionization_energies;
+    const amrex::Real* AMREX_RESTRICT m_adk_prefactor;
+    const amrex::Real* AMREX_RESTRICT m_adk_exp_prefactor;
+    const amrex::Real* AMREX_RESTRICT m_adk_power;
 
     int comp;
     int m_atomic_number;
 
+    GetParticlePosition m_get_position;
+    GetExternalEField m_get_externalE;
+    GetExternalBField m_get_externalB;
+
+    amrex::Array4<const amrex::Real> m_ex_arr;
+    amrex::Array4<const amrex::Real> m_ey_arr;
+    amrex::Array4<const amrex::Real> m_ez_arr;
+    amrex::Array4<const amrex::Real> m_bx_arr;
+    amrex::Array4<const amrex::Real> m_by_arr;
+    amrex::Array4<const amrex::Real> m_bz_arr;
+
+    amrex::IndexType m_ex_type;
+    amrex::IndexType m_ey_type;
+    amrex::IndexType m_ez_type;
+    amrex::IndexType m_bx_type;
+    amrex::IndexType m_by_type;
+    amrex::IndexType m_bz_type;
+
+    amrex::GpuArray<amrex::Real, 3> m_dx_arr;
+    amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
+
+    int m_l_lower_order_in_v;
+    int m_nox;
+    int m_n_rz_azimuthal_modes;
+
+    amrex::Dim3 m_lo;
+
+    IonizationFilterFunc (const WarpXParIter& a_pti, int lev, int ngE,
+                          amrex::FArrayBox const& exfab,
+                          amrex::FArrayBox const& eyfab,
+                          amrex::FArrayBox const& ezfab,
+                          amrex::FArrayBox const& bxfab,
+                          amrex::FArrayBox const& byfab,
+                          amrex::FArrayBox const& bzfab,
+                          amrex::Array<amrex::Real,3> v_galilean,
+                          const amrex::Real* const AMREX_RESTRICT a_ionization_energies,
+                          const amrex::Real* const AMREX_RESTRICT a_adk_prefactor,
+                          const amrex::Real* const AMREX_RESTRICT a_adk_exp_prefactor,
+                          const amrex::Real* const AMREX_RESTRICT a_adk_power,
+                          int a_comp,
+                          int a_atomic_number,
+                          int a_offset = 0) noexcept;
+
     template <typename PData>
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool operator() (const PData& ptd, int i) const noexcept
     {
+        using namespace amrex::literals;
+
         const int ion_lev = ptd.m_runtime_idata[comp][i];
         if (ion_lev < m_atomic_number)
         {
             constexpr amrex::Real c = PhysConst::c;
             constexpr amrex::Real c2_inv = 1./c/c;
 
+            // gather E and B
+            amrex::ParticleReal xp, yp, zp;
+            m_get_position(i, xp, yp, zp);
+
+            amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
+            m_get_externalE(i, ex, ey, ez);
+
+            amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
+            m_get_externalB(i, bx, by, bz);
+
+            doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
+                           m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
+                           m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
+                           m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
+                           m_nox, m_l_lower_order_in_v);
+
             // Compute electric field amplitude in the particle's frame of
             // reference (particularly important when in boosted frame).
             amrex::ParticleReal ux = ptd.m_rdata[PIdx::ux][i];
             amrex::ParticleReal uy = ptd.m_rdata[PIdx::uy][i];
             amrex::ParticleReal uz = ptd.m_rdata[PIdx::uz][i];
-            amrex::ParticleReal ex = ptd.m_rdata[PIdx::Ex][i];
-            amrex::ParticleReal ey = ptd.m_rdata[PIdx::Ey][i];
-            amrex::ParticleReal ez = ptd.m_rdata[PIdx::Ez][i];
-            amrex::ParticleReal bx = ptd.m_rdata[PIdx::Bx][i];
-            amrex::ParticleReal by = ptd.m_rdata[PIdx::By][i];
-            amrex::ParticleReal bz = ptd.m_rdata[PIdx::Bz][i];
 
             amrex::Real ga = std::sqrt(1. + (ux*ux + uy*uy + uz*uz) * c2_inv);
             amrex::Real E = std::sqrt(
diff --git a/Source/Particles/ElementaryProcess/Ionization.cpp b/Source/Particles/ElementaryProcess/Ionization.cpp
new file mode 100644
index 000000000..4c30987f4
--- /dev/null
+++ b/Source/Particles/ElementaryProcess/Ionization.cpp
@@ -0,0 +1,72 @@
+/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
+ * Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpX.H"
+#include "Particles/ElementaryProcess/Ionization.H"
+
+IonizationFilterFunc::IonizationFilterFunc (const WarpXParIter& a_pti, int lev, int ngE,
+                                            amrex::FArrayBox const& exfab,
+                                            amrex::FArrayBox const& eyfab,
+                                            amrex::FArrayBox const& ezfab,
+                                            amrex::FArrayBox const& bxfab,
+                                            amrex::FArrayBox const& byfab,
+                                            amrex::FArrayBox const& bzfab,
+                                            amrex::Array<amrex::Real,3> v_galilean,
+                                            const amrex::Real* const AMREX_RESTRICT a_ionization_energies,
+                                            const amrex::Real* const AMREX_RESTRICT a_adk_prefactor,
+                                            const amrex::Real* const AMREX_RESTRICT a_adk_exp_prefactor,
+                                            const amrex::Real* const AMREX_RESTRICT a_adk_power,
+                                            int a_comp,
+                                            int a_atomic_number,
+                                            int a_offset) noexcept
+{
+    m_ionization_energies = a_ionization_energies;
+    m_adk_prefactor = a_adk_prefactor;
+    m_adk_exp_prefactor = a_adk_exp_prefactor;
+    m_adk_power = a_adk_power;
+    comp = a_comp;
+    m_atomic_number = a_atomic_number;
+
+    m_get_position  = GetParticlePosition(a_pti, a_offset);
+    m_get_externalE = GetExternalEField  (a_pti, a_offset);
+    m_get_externalB = GetExternalBField  (a_pti, a_offset);
+
+    m_ex_arr = exfab.array();
+    m_ey_arr = eyfab.array();
+    m_ez_arr = ezfab.array();
+    m_bx_arr = bxfab.array();
+    m_by_arr = byfab.array();
+    m_bz_arr = bzfab.array();
+
+    m_ex_type = exfab.box().ixType();
+    m_ey_type = eyfab.box().ixType();
+    m_ez_type = ezfab.box().ixType();
+    m_bx_type = bxfab.box().ixType();
+    m_by_type = byfab.box().ixType();
+    m_bz_type = bzfab.box().ixType();
+
+    amrex::Box box = a_pti.tilebox();
+    box.grow(ngE);
+
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev, 0));
+    m_dx_arr = {dx[0], dx[1], dx[2]};
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<amrex::Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, lev);
+    m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    m_l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    m_nox = WarpX::nox;
+    m_n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    m_lo = amrex::lbound(box);
+}
diff --git a/Source/Particles/ElementaryProcess/Make.package b/Source/Particles/ElementaryProcess/Make.package
index 588d017e6..cc0610731 100644
--- a/Source/Particles/ElementaryProcess/Make.package
+++ b/Source/Particles/ElementaryProcess/Make.package
@@ -1,4 +1,8 @@
+CEXE_sources += Ionization.cpp
+
 ifeq ($(QED),TRUE)
+     CEXE_sources += QEDPairGeneration.cpp
+     CEXE_sources += QEDPhotonEmission.cpp
     include $(WARPX_HOME)/Source/Particles/ElementaryProcess/QEDInternals/Make.package
     VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/ElementaryProcess/QEDInternals/
 endif
diff --git a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
index bed6337c3..7056ac884 100644
--- a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
+++ b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
@@ -84,6 +84,9 @@ public:
 class BreitWheelerEvolveOpticalDepth
 {
 public:
+
+    BreitWheelerEvolveOpticalDepth () = default;
+
     /**
      * Constructor acquires a reference to control parameters and
      * lookup tables data.
@@ -138,10 +141,10 @@ public:
     }
 
 private:
-    //laser wavelenght is not used with SI units
-    const amrex::Real m_dummy_lambda{1.0};
+    //laser wavelength is not used with SI units
+    amrex::Real m_dummy_lambda{1.0};
 
-    const PicsarBreitWheelerCtrl m_ctrl;
+    PicsarBreitWheelerCtrl m_ctrl;
 
     //lookup table data
     size_t m_TTfunc_size;
diff --git a/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
index 6d46c9019..fe7d560cd 100644
--- a/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
+++ b/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
@@ -84,6 +84,9 @@ public:
 class QuantumSynchrotronEvolveOpticalDepth
 {
 public:
+
+    QuantumSynchrotronEvolveOpticalDepth () = default;
+
     /**
      * Constructor acquires pointers to control parameters and
      * lookup tables data.
@@ -139,10 +142,10 @@ public:
     }
 
 private:
-    //laser wavelenght is not used with SI units
-    const amrex::Real m_dummy_lambda{1.0};
+    //laser wavelength is not used with SI units
+    amrex::Real m_dummy_lambda{1.0};
 
-    const PicsarQuantumSynchrotronCtrl m_ctrl;
+    PicsarQuantumSynchrotronCtrl m_ctrl;
 
     //lookup table data
     size_t m_KKfunc_size;
diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.H b/Source/Particles/ElementaryProcess/QEDPairGeneration.H
index 22f37a351..838afcaac 100644
--- a/Source/Particles/ElementaryProcess/QEDPairGeneration.H
+++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.H
@@ -10,7 +10,9 @@
 
 #include "Utils/WarpXConst.H"
 #include "Particles/WarpXParticleContainer.H"
-
+#include "Particles/Gather/GetExternalFields.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 #include "QEDInternals/BreitWheelerEngineWrapper.H"
 
 /** @file
@@ -32,9 +34,9 @@ public:
     *
     * @param[in] opt_depth_runtime_comp index of the optical depth component
     */
-    PairGenerationFilterFunc(int const opt_depth_runtime_comp):
-        m_opt_depth_runtime_comp{opt_depth_runtime_comp}
-    {}
+    PairGenerationFilterFunc(int const opt_depth_runtime_comp)
+        : m_opt_depth_runtime_comp(opt_depth_runtime_comp)
+        {}
 
     /**
     * \brief Functor call. This method determines if a given (photon) particle
@@ -75,9 +77,16 @@ public:
     *
     * @param[in] generate_functor functor to be called to determine the properties of the generated pairs
     */
-    PairGenerationTransformFunc(BreitWheelerGeneratePairs const generate_functor):
-        m_generate_functor{generate_functor}
-    {}
+    PairGenerationTransformFunc(BreitWheelerGeneratePairs const generate_functor,
+                                const WarpXParIter& a_pti, int lev, int ngE,
+                                amrex::FArrayBox const& exfab,
+                                amrex::FArrayBox const& eyfab,
+                                amrex::FArrayBox const& ezfab,
+                                amrex::FArrayBox const& bxfab,
+                                amrex::FArrayBox const& byfab,
+                                amrex::FArrayBox const& bzfab,
+                                amrex::Array<amrex::Real,3> v_galilean,
+                                int a_offset = 0);
 
     /**
     * \brief Functor call. It determines the properties of the generated pair
@@ -104,12 +113,22 @@ public:
         const ParticleReal ux = src.m_rdata[PIdx::ux][i_src];
         const ParticleReal uy = src.m_rdata[PIdx::uy][i_src];
         const ParticleReal uz = src.m_rdata[PIdx::uz][i_src];
-        const ParticleReal ex = src.m_rdata[PIdx::Ex][i_src];
-        const ParticleReal ey = src.m_rdata[PIdx::Ey][i_src];
-        const ParticleReal ez = src.m_rdata[PIdx::Ez][i_src];
-        const ParticleReal bx = src.m_rdata[PIdx::Bx][i_src];
-        const ParticleReal by = src.m_rdata[PIdx::By][i_src];
-        const ParticleReal bz = src.m_rdata[PIdx::Bz][i_src];
+
+        // gather E and B
+        amrex::ParticleReal xp, yp, zp;
+        m_get_position(i_src, xp, yp, zp);
+
+        amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
+        m_get_externalE(i_src, ex, ey, ez);
+
+        amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
+        m_get_externalB(i_src, bx, by, bz);
+
+        doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
+                       m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
+                       m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
+                       m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
+                       m_nox, m_l_lower_order_in_v);
 
         const auto px = ux*me;
         const auto py = uy*me;
@@ -153,7 +172,34 @@ public:
 private:
 
     const BreitWheelerGeneratePairs
-        m_generate_functor; /*!< A copy of the functor to generate pairs. It contains only pointers to the lookup tables.*/
+    m_generate_functor; /*!< A copy of the functor to generate pairs. It contains only pointers to the lookup tables.*/
+
+    GetParticlePosition m_get_position;
+    GetExternalEField m_get_externalE;
+    GetExternalBField m_get_externalB;
+
+    amrex::Array4<const amrex::Real> m_ex_arr;
+    amrex::Array4<const amrex::Real> m_ey_arr;
+    amrex::Array4<const amrex::Real> m_ez_arr;
+    amrex::Array4<const amrex::Real> m_bx_arr;
+    amrex::Array4<const amrex::Real> m_by_arr;
+    amrex::Array4<const amrex::Real> m_bz_arr;
+
+    amrex::IndexType m_ex_type;
+    amrex::IndexType m_ey_type;
+    amrex::IndexType m_ez_type;
+    amrex::IndexType m_bx_type;
+    amrex::IndexType m_by_type;
+    amrex::IndexType m_bz_type;
+
+    amrex::GpuArray<amrex::Real, 3> m_dx_arr;
+    amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
+
+    int m_l_lower_order_in_v;
+    int m_nox;
+    int m_n_rz_azimuthal_modes;
+
+    amrex::Dim3 m_lo;
 };
 
 #endif //QED_PAIR_GENERATION_H_
diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp b/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp
new file mode 100644
index 000000000..cdb88d281
--- /dev/null
+++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp
@@ -0,0 +1,62 @@
+/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
+ * Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpX.H"
+#include "Particles/ElementaryProcess/QEDPairGeneration.H"
+
+PairGenerationTransformFunc::
+PairGenerationTransformFunc (BreitWheelerGeneratePairs const generate_functor,
+                             const WarpXParIter& a_pti, int lev, int ngE,
+                             amrex::FArrayBox const& exfab,
+                             amrex::FArrayBox const& eyfab,
+                             amrex::FArrayBox const& ezfab,
+                             amrex::FArrayBox const& bxfab,
+                             amrex::FArrayBox const& byfab,
+                             amrex::FArrayBox const& bzfab,
+                             amrex::Array<amrex::Real,3> v_galilean,
+                             int a_offset)
+: m_generate_functor(generate_functor)
+{
+    m_get_position  = GetParticlePosition(a_pti, a_offset);
+    m_get_externalE = GetExternalEField  (a_pti, a_offset);
+    m_get_externalB = GetExternalBField  (a_pti, a_offset);
+
+    m_ex_arr = exfab.array();
+    m_ey_arr = eyfab.array();
+    m_ez_arr = ezfab.array();
+    m_bx_arr = bxfab.array();
+    m_by_arr = byfab.array();
+    m_bz_arr = bzfab.array();
+
+    m_ex_type = exfab.box().ixType();
+    m_ey_type = eyfab.box().ixType();
+    m_ez_type = ezfab.box().ixType();
+    m_bx_type = bxfab.box().ixType();
+    m_by_type = byfab.box().ixType();
+    m_bz_type = bzfab.box().ixType();
+
+    amrex::Box box = a_pti.tilebox();
+    box.grow(ngE);
+
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev, 0));
+    m_dx_arr = {dx[0], dx[1], dx[2]};
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<amrex::Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, lev);
+    m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    m_l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    m_nox = WarpX::nox;
+    m_n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    m_lo = amrex::lbound(box);
+}
diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
index c7c801eeb..4dc573f60 100644
--- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
+++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
@@ -10,7 +10,9 @@
 
 #include "Utils/WarpXConst.H"
 #include "Particles/WarpXParticleContainer.H"
-
+#include "Particles/Gather/GetExternalFields.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 #include "QEDInternals/QuantumSyncEngineWrapper.H"
 
 /** @file
@@ -32,9 +34,9 @@ public:
     *
     * @param[in] opt_depth_runtime_comp Index of the optical depth component
     */
-    PhotonEmissionFilterFunc(int const opt_depth_runtime_comp):
-        m_opt_depth_runtime_comp{opt_depth_runtime_comp}
-    {}
+    PhotonEmissionFilterFunc(int const opt_depth_runtime_comp)
+        : m_opt_depth_runtime_comp(opt_depth_runtime_comp)
+        {}
 
     /**
     * \brief Functor call. This method determines if a given (electron or positron)
@@ -56,7 +58,7 @@ public:
     }
 
 private:
-    int m_opt_depth_runtime_comp /*!< Index of the optical depth component of the source species*/;
+    int m_opt_depth_runtime_comp; /*!< Index of the optical depth component of the source species*/
 };
 
 /**
@@ -81,15 +83,19 @@ public:
     * @param[in] opt_depth_runtime_comp index of the optical depth component of the source species
     * @param[in] emission_functor functor to generate photons and update momentum of the source particles
     */
-    PhotonEmissionTransformFunc(
+    PhotonEmissionTransformFunc (
         QuantumSynchrotronGetOpticalDepth opt_depth_functor,
         int const opt_depth_runtime_comp,
-        QuantumSynchrotronGeneratePhotonAndUpdateMomentum const emission_functor
-        ):
-        m_opt_depth_functor{opt_depth_functor},
-        m_opt_depth_runtime_comp{opt_depth_runtime_comp},
-        m_emission_functor{emission_functor}
-        {}
+        QuantumSynchrotronGeneratePhotonAndUpdateMomentum const emission_functor,
+        const WarpXParIter& a_pti, int lev, int ngE,
+        amrex::FArrayBox const& exfab,
+        amrex::FArrayBox const& eyfab,
+        amrex::FArrayBox const& ezfab,
+        amrex::FArrayBox const& bxfab,
+        amrex::FArrayBox const& byfab,
+        amrex::FArrayBox const& bzfab,
+        amrex::Array<amrex::Real,3> v_galilean,
+        int a_offset = 0);
 
     /**
     * \brief Functor call. It determines the properties of the generated photon
@@ -113,12 +119,22 @@ public:
         const ParticleReal ux = src.m_rdata[PIdx::ux][i_src];
         const ParticleReal uy = src.m_rdata[PIdx::uy][i_src];
         const ParticleReal uz = src.m_rdata[PIdx::uz][i_src];
-        const ParticleReal ex = src.m_rdata[PIdx::Ex][i_src];
-        const ParticleReal ey = src.m_rdata[PIdx::Ey][i_src];
-        const ParticleReal ez = src.m_rdata[PIdx::Ez][i_src];
-        const ParticleReal bx = src.m_rdata[PIdx::Bx][i_src];
-        const ParticleReal by = src.m_rdata[PIdx::By][i_src];
-        const ParticleReal bz = src.m_rdata[PIdx::Bz][i_src];
+
+        // gather E and B
+        amrex::ParticleReal xp, yp, zp;
+        m_get_position(i_src, xp, yp, zp);
+
+        amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
+        m_get_externalE(i_src, ex, ey, ez);
+
+        amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
+        m_get_externalB(i_src, bx, by, bz);
+
+        doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
+                       m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
+                       m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
+                       m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
+                       m_nox, m_l_lower_order_in_v);
 
         // Particle momentum is stored as gamma * velocity.
         // Convert to m * gamma * velocity before applying the emission functor.
@@ -160,6 +176,33 @@ private:
     const QuantumSynchrotronGeneratePhotonAndUpdateMomentum
         m_emission_functor;  /*!< A copy of the functor to generate photons. It contains only pointers to the lookup tables.*/
     const int m_opt_depth_runtime_comp = 0;  /*!< Index of the optical depth component of source species*/
+
+    GetParticlePosition m_get_position;
+    GetExternalEField m_get_externalE;
+    GetExternalBField m_get_externalB;
+
+    amrex::Array4<const amrex::Real> m_ex_arr;
+    amrex::Array4<const amrex::Real> m_ey_arr;
+    amrex::Array4<const amrex::Real> m_ez_arr;
+    amrex::Array4<const amrex::Real> m_bx_arr;
+    amrex::Array4<const amrex::Real> m_by_arr;
+    amrex::Array4<const amrex::Real> m_bz_arr;
+
+    amrex::IndexType m_ex_type;
+    amrex::IndexType m_ey_type;
+    amrex::IndexType m_ez_type;
+    amrex::IndexType m_bx_type;
+    amrex::IndexType m_by_type;
+    amrex::IndexType m_bz_type;
+
+    amrex::GpuArray<amrex::Real, 3> m_dx_arr;
+    amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
+
+    int m_l_lower_order_in_v;
+    int m_nox;
+    int m_n_rz_azimuthal_modes;
+
+    amrex::Dim3 m_lo;
 };
 
 
diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp b/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp
new file mode 100644
index 000000000..6134ff542
--- /dev/null
+++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp
@@ -0,0 +1,66 @@
+/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
+ * Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpX.H"
+#include "Particles/ElementaryProcess/QEDPhotonEmission.H"
+
+PhotonEmissionTransformFunc::
+PhotonEmissionTransformFunc (QuantumSynchrotronGetOpticalDepth opt_depth_functor,
+                             int const opt_depth_runtime_comp,
+                             QuantumSynchrotronGeneratePhotonAndUpdateMomentum const emission_functor,
+                             const WarpXParIter& a_pti, int lev, int ngE,
+                             amrex::FArrayBox const& exfab,
+                             amrex::FArrayBox const& eyfab,
+                             amrex::FArrayBox const& ezfab,
+                             amrex::FArrayBox const& bxfab,
+                             amrex::FArrayBox const& byfab,
+                             amrex::FArrayBox const& bzfab,
+                             amrex::Array<amrex::Real,3> v_galilean,
+                             int a_offset)
+:m_opt_depth_functor{opt_depth_functor},
+ m_opt_depth_runtime_comp{opt_depth_runtime_comp},
+ m_emission_functor{emission_functor}
+{
+    m_get_position  = GetParticlePosition(a_pti, a_offset);
+    m_get_externalE = GetExternalEField  (a_pti, a_offset);
+    m_get_externalB = GetExternalBField  (a_pti, a_offset);
+
+    m_ex_arr = exfab.array();
+    m_ey_arr = eyfab.array();
+    m_ez_arr = ezfab.array();
+    m_bx_arr = bxfab.array();
+    m_by_arr = byfab.array();
+    m_bz_arr = bzfab.array();
+
+    m_ex_type = exfab.box().ixType();
+    m_ey_type = eyfab.box().ixType();
+    m_ez_type = ezfab.box().ixType();
+    m_bx_type = bxfab.box().ixType();
+    m_by_type = byfab.box().ixType();
+    m_bz_type = bzfab.box().ixType();
+
+    amrex::Box box = a_pti.tilebox();
+    box.grow(ngE);
+
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev, 0));
+    m_dx_arr = {dx[0], dx[1], dx[2]};
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<amrex::Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, lev);
+    m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    m_l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    m_nox = WarpX::nox;
+    m_n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    m_lo = amrex::lbound(box);
+}