Integrate picsar_qed in ablastr : move in the Schwinger process

Docs/source/refs.bib

@@ -2791,3 +2791,14 @@ @article{Zhang2017
 numpages = {28},
 pages = {011304},
 }
+
+@article{NarozhnyPRA2004,
+  title={On e+ e- pair production by colliding electromagnetic pulses},
+  author={Narozhny, Nikolay Borisovich and Bulanov, SS and Mur, Vadim Davydovich and Popov, Vladimir Stepanovich},
+  journal={Journal of Experimental and Theoretical Physics Letters},
+  volume={80},
+  number={6},
+  pages={382--385},
+  year={2004},
+  publisher={Springer}
+}

Docs/source/theory/multiphysics/qed.rst

@@ -20,7 +20,31 @@ Breit-Wheeler
 Schwinger process
 -----------------
 
-If the code is compiled with QED and the user activates the Schwinger process in the input file,
+In the Schwinger process, electron-positron pairs are created in vacuum by a sufficiently strong electromagnetic field.
+An expression for the Schwinger pair production rate can be found in :cite:p:`mqed-NarozhnyPRA2004`.
+Rewritten in SI units, the pair production rate per unit volume is:
+
+.. math::
+
+    \dfrac{d^2N}{dt dV} = \dfrac{e^2 E_s^2}{4 \pi^2 \hbar^2 c} \epsilon \eta \coth{\left({\dfrac{\pi \eta}{\epsilon}}\right)}\exp{\left({-\dfrac{\pi}{\epsilon}}\right)}
+
+
+where $e$ is the elementary charge, :math:`E_s`` is the Schwinger field, :math:`\epsilon = \mathcal{E}/E_s` and :math:`\eta = \mathcal{H}/E_s`.
+:math:`\mathcal{E}` and :math:`\mathcal{H}$` are given by:
+
+.. math::
+
+    \mathcal{E} = \sqrt{\sqrt{\mathcal{F}^2 + \mathcal{G}^2} + \mathcal{F}} \\
+    \mathcal{H} = \sqrt{\sqrt{\mathcal{F}^2 + \mathcal{G}^2} - \mathcal{F}} \\
+
+:math:`\mathcal{F}` and :math:`\mathcal{G}$` are the invariants of the electromagnetic field and are equal to:
+
+.. math::
+
+    \mathcal{F} = (\mathbf{E}^2 - c^2 \mathbf{B}^2)/2 \\
+    \mathcal{G} = c \mathbf{E} \cdot \mathbf{B} \\
+
+If the user activates the Schwinger process in the input file,
 electron-positron pairs can be created in vacuum in the function
 ``MultiParticleContainer::doQEDSchwinger``:
 
@@ -44,3 +68,6 @@ the total number of physical pairs created.
 So far the Schwinger module requires using ``warpx.grid_type = collocated`` or
 ``algo.field_gathering = momentum-conserving`` (so that the auxiliary fields are calculated on the nodes)
 and is not compatible with either mesh refinement, RZ, RCYLINDER, and RSPHERE coordinates or single precision.
+
+.. bibliography::
+    :keyprefix: mqed-

Docs/source/usage/parameters.rst

@@ -3902,7 +3902,7 @@ Lookup tables and other settings for QED modules
 ------------------------------------------------
 
 Lookup tables store pre-computed values for functions used by the QED modules.
-**This feature requires to compile with QED=TRUE (and also with QED_TABLE_GEN=TRUE for table generation)**
+Unless otherwise specified, **these feature requires to compile with QED=TRUE (and also with QED_TABLE_GEN=TRUE for table generation)**
 
 * ``qed_bw.lookup_table_mode`` (`string`)
     There are three options to prepare the lookup table required by the Breit-Wheeler module:
@@ -3991,7 +3991,7 @@ Lookup tables store pre-computed values for functions used by the QED modules.
 
 * ``warpx.do_qed_schwinger`` (`bool`) optional (default `0`)
     If this is 1, Schwinger electron-positron pairs can be generated in vacuum in the cells where the EM field is high enough.
-    Activating the Schwinger process requires the code to be compiled with ``QED=TRUE`` and ``PICSAR``.
+    Activating the Schwinger process does not require the code to be compiled with ``QED=TRUE``.
     If ``warpx.do_qed_schwinger = 1``, Schwinger product species must be specified with
     ``qed_schwinger.ele_product_species`` and ``qed_schwinger.pos_product_species``.
     Schwinger process requires either ``warpx.grid_type = collocated`` or

Regression/Checksum/benchmarks_json/test_3d_qed_schwinger_4.json

@@ -1,31 +1,31 @@
 {
-  "lev=0": {
-    "Bx": 10.282332458788613,
-    "By": 1707847966720000.0,
-    "Bz": 0.5093466319155224,
-    "Ex": 3460290602.1197195,
-    "Ey": 0.0,
-    "Ez": 5.119999999999012e+23,
-    "jx": 354749197616543.3,
-    "jy": 0.0,
-    "jz": 8899048833947625.0
-  },
   "ele_schwinger": {
     "particle_momentum_x": 3.695301335661867e-12,
-    "particle_momentum_y": 1.0868696459770461e-26,
+    "particle_momentum_y": 1.0868696454896924e-26,
     "particle_momentum_z": 2.8810249733308265e-12,
     "particle_position_x": 0.00028800000000000006,
     "particle_position_y": 0.000216,
     "particle_position_z": 0.0002540000004700524,
-    "particle_weight": 63968.698937700814
+    "particle_weight": 63968.69690641211
+  },
+  "lev=0": {
+    "Bx": 10.282332458788613,
+    "By": 1707847966720000.0,
+    "Bz": 0.5093466195466327,
+    "Ex": 3460290503.614944,
+    "Ey": 0.0,
+    "Ez": 5.119999999999012e+23,
+    "jx": 354749188059060.5,
+    "jy": 0.0,
+    "jz": 8899048551363524.0
   },
   "pos_schwinger": {
     "particle_momentum_x": 3.695301335661867e-12,
-    "particle_momentum_y": 1.0881422563719172e-26,
+    "particle_momentum_y": 1.0881422558277477e-26,
     "particle_momentum_z": 2.881024973330826e-12,
     "particle_position_x": 0.000288,
     "particle_position_y": 0.000216,
     "particle_position_z": 0.00023799999952994755,
-    "particle_weight": 63968.698937700814
+    "particle_weight": 63968.696906412115
   }
 }

Source/Evolve/WarpXEvolve.cpp

@@ -218,11 +218,12 @@ WarpX::Evolve (int numsteps)
         // multi-physics: field ionization
         doFieldIonization();
 
-#ifdef WARPX_QED
         // multi-physics: QED effects
+#ifdef WARPX_QED
         doQEDEvents();
-        mypc->doQEDSchwinger();
 #endif
+        mypc->doQEDSchwinger();
+
 
         // perform particle injection
         ExecutePythonCallback("particleinjection");

Source/Particles/ElementaryProcess/QEDSchwingerProcess.H

@@ -8,11 +8,12 @@
 #ifndef WARPX_QED_SCHWINGER_PROCESS_H_
 #define WARPX_QED_SCHWINGER_PROCESS_H_
 
-#include "Particles/ElementaryProcess/QEDInternals/SchwingerProcessWrapper.H"
 #include "Utils/TextMsg.H"
 
+#include "ablastr/qed/SchwingerProcess.H"
+
 /**
- * This structure is a functor which calls getSchwingerProductionNumber to
+ * This structure is a functor which calls ablastr::qed::schwinger::get_schwinger_pair_number to
  * calculate the number of pairs created during a given timestep at a given cell.
  */
 struct SchwingerFilterFunc
@@ -45,10 +46,11 @@ struct SchwingerFilterFunc
         const auto& arrBy = src_FABs.By;
         const auto& arrBz = src_FABs.Bz;
 
-        return getSchwingerProductionNumber( m_dV, m_dt,
-                    arrEx(i,j,k),arrEy(i,j,k),arrEz(i,j,k),
-                    arrBx(i,j,k),arrBy(i,j,k),arrBz(i,j,k),
-                    m_threshold_poisson_gaussian,engine);
+        return ablastr::qed::schwinger::get_schwinger_pair_number(
+            m_dV, m_dt,
+            arrEx(i,j,k), arrEy(i,j,k), arrEz(i,j,k),
+            arrBx(i,j,k), arrBy(i,j,k), arrBz(i,j,k),
+            m_threshold_poisson_gaussian, engine);
     }
 };
 

Source/Particles/MultiParticleContainer.H

@@ -205,7 +205,6 @@ public:
     */
     void doResampling (const amrex::Vector<amrex::Geometry>& geom, int timestep, bool verbose);
 
-#ifdef WARPX_QED
     /** If Schwinger process is activated, this function is called at every
      * timestep in Evolve and is used to create Schwinger electron-positron pairs.
      * Within this function we loop over all cells to calculate the number of
@@ -220,7 +219,6 @@ public:
      * object (to make the link between Real and int quantities).
      */
     [[nodiscard]] amrex::Box ComputeSchwingerGlobalBox () const;
-#endif
 
     void Restart (const std::string& dir);
 
@@ -498,33 +496,6 @@ protected:
      */
     void BreitWheelerGenerateTable();
 
-    /** Whether or not to activate Schwinger process */
-    bool m_do_qed_schwinger = false;
-    /** Name of Schwinger electron product species */
-    std::string m_qed_schwinger_ele_product_name;
-    /** Name of Schwinger positron product species */
-    std::string m_qed_schwinger_pos_product_name;
-    /** Index for Schwinger electron product species in allcontainers */
-    int m_qed_schwinger_ele_product;
-    /** Index for Schwinger positron product species in allcontainers */
-    int m_qed_schwinger_pos_product;
-    /** Transverse size used in 2D Schwinger pair production rate calculations */
-    amrex::Real m_qed_schwinger_y_size;
-    /** If the number of physical Schwinger pairs created within a cell is
-     * higher than this threshold we use a Gaussian distribution rather than
-     * a Poisson distribution for the pair production rate calculations
-     */
-    int m_qed_schwinger_threshold_poisson_gaussian = 25;
-    /** The 6 following variables are spatial boundaries beyond which Schwinger process is
-     *  deactivated
-     */
-    amrex::Real m_qed_schwinger_xmin = std::numeric_limits<amrex::Real>::lowest();
-    amrex::Real m_qed_schwinger_xmax = std::numeric_limits<amrex::Real>::max();
-    amrex::Real m_qed_schwinger_ymin = std::numeric_limits<amrex::Real>::lowest();
-    amrex::Real m_qed_schwinger_ymax = std::numeric_limits<amrex::Real>::max();
-    amrex::Real m_qed_schwinger_zmin = std::numeric_limits<amrex::Real>::lowest();
-    amrex::Real m_qed_schwinger_zmax = std::numeric_limits<amrex::Real>::max();
-
 #endif
 
 private:
@@ -561,14 +532,40 @@ private:
     */
     void CheckIonizationProductSpecies();
 
-#ifdef WARPX_QED
     /**
     * Should be called right after mapSpeciesProduct in InitData.
     * It checks the physical correctness of product particle species
     * selected by the user for QED processes.
     */
     void CheckQEDProductSpecies();
-#endif
+
+    /** Whether or not to activate Schwinger process */
+    bool m_do_qed_schwinger = false;
+    /** Name of Schwinger electron product species */
+    std::string m_qed_schwinger_ele_product_name;
+    /** Name of Schwinger positron product species */
+    std::string m_qed_schwinger_pos_product_name;
+    /** Index for Schwinger electron product species in allcontainers */
+    int m_qed_schwinger_ele_product;
+    /** Index for Schwinger positron product species in allcontainers */
+    int m_qed_schwinger_pos_product;
+    /** Transverse size used in 2D Schwinger pair production rate calculations */
+    amrex::Real m_qed_schwinger_y_size;
+    /** If the number of physical Schwinger pairs created within a cell is
+     * higher than this threshold we use a Gaussian distribution rather than
+     * a Poisson distribution for the pair production rate calculations
+     */
+    int m_qed_schwinger_threshold_poisson_gaussian = 25;
+    /** The 6 following variables are spatial boundaries beyond which Schwinger process is
+     *  deactivated
+     */
+    amrex::Real m_qed_schwinger_xmin = std::numeric_limits<amrex::Real>::lowest();
+    amrex::Real m_qed_schwinger_xmax = std::numeric_limits<amrex::Real>::max();
+    amrex::Real m_qed_schwinger_ymin = std::numeric_limits<amrex::Real>::lowest();
+    amrex::Real m_qed_schwinger_ymax = std::numeric_limits<amrex::Real>::max();
+    amrex::Real m_qed_schwinger_zmin = std::numeric_limits<amrex::Real>::lowest();
+    amrex::Real m_qed_schwinger_zmax = std::numeric_limits<amrex::Real>::max();
+
 
 
 };