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>

2 files changed, 174 insertions, 0 deletions

diff --git a/Examples/Modules/qed/schwinger/analysis_schwinger.py b/Examples/Modules/qed/schwinger/analysis_schwinger.py
new file mode 100755
index 000000000..7d7f4c852
--- /dev/null
+++ b/Examples/Modules/qed/schwinger/analysis_schwinger.py
@@ -0,0 +1,104 @@
+#! /usr/bin/env python
+
+# Copyright 2020 Luca Fedeli, Neil Zaim
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+## This module is imported by each of the scripts used to analyze the Schwinger tests.
+
+## The pair production rate is calculated using the formula described in
+## Bulanov, S. S., et al. Physical review letters 104.22 (2010): 220404.
+
+import yt
+import numpy as np
+import sys
+import re
+
+# define some parameters
+
+c = 299792458.
+m_e = 9.10938356e-31
+e = 1.6021766208e-19
+hbar = 1.054571800e-34
+E_S = m_e**2*c**3/e/hbar # Schwinger field
+
+dV = (1.e-6)**3 # total simulation volume
+dt = 2.407291502e-16
+filename = sys.argv[1]
+
+Ex_test = 0.
+Ey_test = 0.
+Ez_test = 0.
+Bx_test = 0.
+By_test = 0.
+Bz_test = 0.
+
+# Find which test we are doing
+test_number = re.search( 'qed_schwinger([1234])', filename ).group(1)
+if test_number == '1':
+    Ex_test = 1.e16
+    Bx_test = 16792888.570516706
+    By_test = 5256650.141557486
+    Bz_test = 18363530.799561853
+elif test_number == '2':
+    Ex_test = 1.e18
+    Bx_test = 1679288857.0516706
+    By_test = 525665014.1557486
+    Bz_test = 1836353079.9561853
+elif test_number == '3':
+    Ey_test = 1.0321239524474501e+17
+elif test_number == '4':
+    Ez_test = 2.5e+20
+    By_test = 833910154604.3563
+else:
+    assert(False)
+
+def calculate_rate(Ex,Ey,Ez,Bx,By,Bz):
+## Calculate theoretical pair production rate from EM field value
+
+    E_squared = Ex**2 + Ey**2 + Ez**2
+    H_squared = c**2*(Bx**2 + By**2 + Bz**2)
+
+    F = (E_squared - H_squared)/2.
+    G = c*(Ex*Bx + Ey*By + Ez*Bz)
+
+    epsilon = np.sqrt(np.sqrt(F**2+G**2)+F)/E_S
+    eta = np.sqrt(np.sqrt(F**2+G**2)-F)/E_S
+
+    if(epsilon != 0. and eta != 0.):
+        return  e**2*E_S**2/4./np.pi**2/c/hbar**2*epsilon*eta/np.tanh(np.pi*eta/epsilon)*np.exp(-np.pi/epsilon)
+    elif (epsilon == 0.):
+        return 0.
+    else:
+        return  e**2*E_S**2/4./np.pi**2/c/hbar**2*epsilon**2/np.pi*np.exp(-np.pi/epsilon)
+
+
+def do_analysis(Ex,Ey,Ez,Bx,By,Bz):
+
+    data_set = yt.load(filename)
+
+    expected_total_physical_pairs_created = dV*dt*calculate_rate(Ex,Ey,Ez,Bx,By,Bz)
+    if expected_total_physical_pairs_created < 0.01:
+        np_ele = data_set.particle_type_counts["ele_schwinger"] if \
+            "ele_schwinger" in data_set.particle_type_counts.keys() else 0
+        np_pos = data_set.particle_type_counts["pos_schwinger"] if \
+            "pos_schwinger" in data_set.particle_type_counts.keys() else 0
+        assert(np_ele == 0 and np_pos == 0)
+        ## Assert whether pairs are created or not.
+
+    else:
+        all_data = data_set.all_data()
+
+        ele_data = all_data["ele_schwinger",'particle_weight']
+        pos_data = all_data["pos_schwinger",'particle_weight']
+
+        std_total_physical_pairs_created = np.sqrt(expected_total_physical_pairs_created)
+
+        # This first assert only works if a single tile is used in the simulation
+        assert(np.array_equal(ele_data,pos_data))
+        # 5 sigma test that has an intrisic probability to fail of 1 over ~2 millions
+        assert(np.abs(np.sum(ele_data)-expected_total_physical_pairs_created)<5*std_total_physical_pairs_created)
+
+do_analysis(Ex_test, Ey_test, Ez_test, Bx_test, By_test, Bz_test)
diff --git a/Examples/Modules/qed/schwinger/inputs_3d_schwinger b/Examples/Modules/qed/schwinger/inputs_3d_schwinger
new file mode 100644
index 000000000..04fc6ac67
--- /dev/null
+++ b/Examples/Modules/qed/schwinger/inputs_3d_schwinger
@@ -0,0 +1,70 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  8 8 8
+amr.max_grid_size = 8   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 8 # minimum size of each AMReX box, used to decompose the domain
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 1     1     1  # Is periodic?
+geometry.prob_lo     =  -5.e-7  -5.e-7  -5.e-7   # physical domain
+geometry.prob_hi     =  5.e-7   5.e-7  5.e-7
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = momentum-conserving
+algo.particle_pusher = boris
+interpolation.nox = 1 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 1
+interpolation.noz = 1
+warpx.verbose = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 0 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+#################################
+###### EXTERNAL EM FIELD ########
+#################################
+
+warpx.B_ext_grid_init_style = "constant"
+warpx.E_ext_grid_init_style = "constant"
+warpx.B_external_grid = 0. 0. 0.
+warpx.E_external_grid = 0. 0. 0.
+
+#################################
+############ PLASMA #############
+#################################
+particles.nspecies = 2 # number of species
+particles.species_names =  ele_schwinger pos_schwinger
+
+ele_schwinger.species_type = "electron"
+pos_schwinger.species_type = "positron"
+ele_schwinger.injection_style = "NUniformPerCell"
+pos_schwinger.injection_style = "NUniformPerCell"
+ele_schwinger.num_particles_per_cell_each_dim = 0 0
+pos_schwinger.num_particles_per_cell_each_dim = 0 0
+ele_schwinger.profile = "constant"
+pos_schwinger.profile = "constant"
+ele_schwinger.density = 0
+pos_schwinger.density = 0
+ele_schwinger.momentum_distribution_type = "gaussian"
+pos_schwinger.momentum_distribution_type = "gaussian"
+
+#################################
+############## QED ##############
+#################################
+warpx.do_qed_schwinger = 1
+qed_schwinger.ele_product_species = ele_schwinger
+qed_schwinger.pos_product_species = pos_schwinger
+
+#################################
+########## DIAGNOSTICS ##########
+#################################
+diagnostics.diags_names = diag1
+diag1.diag_type = Full
+diag1.period = 1
+