Merge pull request #344 from lucafedeli88/qed_phys_part_with_lambda

QED: now particles initialize a real component (optical depth) to be used for QED modules

4 files changed, 230 insertions, 0 deletions

diff --git a/Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py b/Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py
new file mode 100755
index 000000000..850ecc0fe
--- /dev/null
+++ b/Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py
@@ -0,0 +1,31 @@
+#! /usr/bin/env python3
+import yt
+import numpy as np
+import scipy.stats as st
+import sys
+
+# This script checks if photons initialized with Breit Wheeler process enabled
+# do actually have an exponentially distributed optical depth
+
+# Tolerance
+tol = 1e-2
+
+def check():
+    filename = sys.argv[1]
+    data_set = yt.load(filename)
+
+    all_data = data_set.all_data()
+    res_tau = all_data["photons", 'particle_tau']
+
+    loc, scale = st.expon.fit(res_tau)
+
+    # loc should be very close to 0, scale should be very close to 1
+    assert(np.abs(loc - 0) < tol)
+    assert(np.abs(scale - 1) < tol)
+
+def main():
+    check()
+
+if __name__ == "__main__":
+    main()
+
diff --git a/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init b/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init
new file mode 100644
index 000000000..06f87f92b
--- /dev/null
+++ b/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init
@@ -0,0 +1,69 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  128 128
+amr.max_grid_size = 128   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain
+amr.plot_int = 1
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 0     0      # Is periodic?
+geometry.prob_lo     =  -32.e-6   -32.e-6    # physical domain
+geometry.prob_hi     =  32.e-6   32.e-6
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+warpx.fine_tag_lo = -5.e-6   -35.e-6
+warpx.fine_tag_hi =  5.e-6   -25.e-6
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = standard
+algo.particle_pusher = boris
+interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.plot_raw_fields = 0
+warpx.plot_raw_fields_guards = 0
+warpx.plot_finepatch = 0
+warpx.plot_crsepatch = 0
+warpx.use_filter = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.nspecies = 1 # number of species
+particles.species_names = photons
+particles.photon_species = photons
+#################################
+
+photons.charge = -q_e
+photons.mass = m_e
+photons.injection_style = "NUniformPerCell"
+photons.profile = "constant"
+photons.xmin = -30e-6
+photons.ymin = -30e-6
+photons.zmin = -30e-6
+photons.xmax = 30e-6
+photons.ymax = 30e-6
+photons.zmax = 30e-6
+photons.num_particles_per_cell_each_dim = 2 2
+photons.density = 1e19
+photons.profile = "constant"
+photons.momentum_distribution_type = "gaussian"
+photons.ux_m = 0.0
+photons.uy_m = 0.0
+photons.uz_m = 0.0
+photons.ux_th = 100.
+photons.uy_th = 100.
+photons.uz_th = 100.
+##########QED####################
+photons.do_qed = 1
+photons.do_qed_breit_wheeler = 1
+#################################
diff --git a/Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py b/Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py
new file mode 100755
index 000000000..05b313ee6
--- /dev/null
+++ b/Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py
@@ -0,0 +1,36 @@
+#! /usr/bin/env python3
+import yt
+import numpy as np
+import scipy.stats as st
+import sys
+
+# This script checks if electrons and positrons initialized with
+# Quantum Synchrotron process enabled
+# do actually have an exponentially distributed optical depth
+
+# Tolerance
+tol = 1e-2
+
+def check():
+    filename = sys.argv[1]
+    data_set = yt.load(filename)
+
+    all_data = data_set.all_data()
+    res_ele_tau = all_data["electrons", 'particle_tau']
+    res_pos_tau = all_data["positrons", 'particle_tau']
+
+    loc_ele, scale_ele = st.expon.fit(res_ele_tau)
+    loc_pos, scale_pos = st.expon.fit(res_pos_tau)
+
+    # loc should be very close to 0, scale should be very close to 1
+    assert(np.abs(loc_ele - 0) < tol)
+    assert(np.abs(loc_pos - 0) < tol)
+    assert(np.abs(scale_ele - 1) < tol)
+    assert(np.abs(scale_pos - 1) < tol)
+
+def main():
+    check()
+
+if __name__ == "__main__":
+    main()
+
diff --git a/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init b/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init
new file mode 100644
index 000000000..cdc5be585
--- /dev/null
+++ b/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init
@@ -0,0 +1,94 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  128 128
+amr.max_grid_size = 128   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain
+amr.plot_int = 1
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 0     0      # Is periodic?
+geometry.prob_lo     =  -32.e-6   -32.e-6    # physical domain
+geometry.prob_hi     =  32.e-6   32.e-6
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+warpx.fine_tag_lo = -5.e-6   -35.e-6
+warpx.fine_tag_hi =  5.e-6   -25.e-6
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = standard
+algo.particle_pusher = boris
+interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.plot_raw_fields = 0
+warpx.plot_raw_fields_guards = 0
+warpx.plot_finepatch = 0
+warpx.plot_crsepatch = 0
+warpx.use_filter = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+
+#################################
+############ PLASMA #############
+#################################
+particles.nspecies = 2 # number of species
+particles.species_names = electrons positrons
+#################################
+
+electrons.charge = -q_e
+electrons.mass = m_e
+electrons.injection_style = "NUniformPerCell"
+electrons.profile = "constant"
+electrons.xmin = -30e-6
+electrons.ymin = -30e-6
+electrons.zmin = -30e-6
+electrons.xmax = 30e-6
+electrons.ymax = 30e-6
+electrons.zmax = 30e-6
+electrons.num_particles_per_cell_each_dim = 2 2
+electrons.density = 1e19
+electrons.profile = "constant"
+electrons.momentum_distribution_type = "gaussian"
+electrons.ux_m = 0.0
+electrons.uy_m = 0.0
+electrons.uz_m = 0.0
+electrons.ux_th = 100.
+electrons.uy_th = 100.
+electrons.uz_th = 100.
+##########QED####################
+electrons.do_qed = 1
+electrons.do_qed_quantum_sync = 1
+#################################
+
+positrons.charge = q_e
+positrons.mass = m_e
+positrons.injection_style = "NUniformPerCell"
+positrons.profile = "constant"
+positrons.xmin = -30e-6
+positrons.ymin = -30e-6
+positrons.zmin = -30e-6
+positrons.xmax = 30e-6
+positrons.ymax = 30e-6
+positrons.zmax = 30e-6
+positrons.num_particles_per_cell_each_dim = 2 2
+positrons.density = 1e19
+positrons.profile = "constant"
+positrons.momentum_distribution_type = "gaussian"
+positrons.ux_m = 0.0
+positrons.uy_m = 0.0
+positrons.uz_m = 0.0
+positrons.ux_th = 100.
+positrons.uy_th = 100.
+positrons.uz_th = 100.
+##########QED####################
+positrons.do_qed = 1
+positrons.do_qed_quantum_sync = 1
+#################################