[mini-PR] Fix bug in Breit-Wheeler engine (#852)

* fixed bug in BW engine

* fixed bug

* fixed bug

* fixed bug

* fixed bug

* fixed bug

* eliminate useless variable

* updated test

* updated inputfile

* Updated tests

* increase tolerance from .04 to .07 in QED 3D BW test

* do plot pos_bw and ele_bw

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

2 files changed, 77 insertions, 56 deletions

diff --git a/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py b/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py
index e5540ffa8..b11e4786d 100755
--- a/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py
+++ b/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py
@@ -14,6 +14,7 @@ import sys
 import math as m
 import scipy.special as spe
 import scipy.integrate as integ
+import scipy.stats as st
 
 # This script checks if the optical depth of photons undergoing the
 # Breit Wheeler process behaves as expected. Four populations of photons
@@ -31,7 +32,7 @@ import scipy.integrate as integ
 
 
 # Tolerance
-tol = 2e-2
+tol = 7.e-2
 
 # EM fields
 E_f = np.array([-2433321316961438, 973328526784575, 1459992790176863])
@@ -53,15 +54,14 @@ schwinger_field_norm = electron_mass*speed_of_light*lambda_ref/(2.0*reduced_plan
 
 # Initial parameters
 spec_names = ["p1", "p2", "p3", "p4"]
- #Assumes average = 1 at t = 0 (ok for a large number of particles)
-tau_begin_avg = np.array([1.0, 1.0, 1.0, 1.0])
 mec = electron_mass*speed_of_light
-p_begin = [
-    np.array([2000.0,0,0])*mec,
-    np.array([0.0,5000.0,0.0])*mec,
-    np.array([0.0,0.0,10000.0])*mec,
-    np.array([57735.02691896, 57735.02691896, 57735.02691896])*mec
-]
+p_begin = {
+    "p1": np.array([2000.0,0,0])*mec,
+    "p2": np.array([0.0,5000.0,0.0])*mec,
+    "p3": np.array([0.0,0.0,10000.0])*mec,
+    "p4": np.array([57735.02691896, 57735.02691896, 57735.02691896])*mec
+}
+initial_particle_number = 16384
 #______________
 
 def calc_chi_gamma(p, E, B):
@@ -100,23 +100,37 @@ def check():
     data_set_end = yt.load(filename_end)
 
     sim_time = data_set_end.current_time.to_value()
-
     all_data_end = data_set_end.all_data()
 
-    tau_end_avg = np.array([
-       np.average(all_data_end[name, 'particle_optical_depth_BW'])
-       for name in spec_names])
-
-    dNBW_dt_sim = (tau_begin_avg - tau_end_avg)/sim_time
-
-    dNBW_dt_theo = [
-        dNBW_dt(calc_chi_gamma(p, E_f, B_f), np.linalg.norm(p*speed_of_light))
-        for p in p_begin
-    ]
-
-    discrepancy = np.abs(dNBW_dt_sim-dNBW_dt_theo)/dNBW_dt_theo
-
-    assert(np.all(discrepancy < tol))
+    for name in spec_names:
+        opt_depth = all_data_end[name, 'particle_optical_depth_BW']
+
+        #check that the distribution is still exponential with scale 1 and loc 0
+        opt_depth_loc, opt_depth_scale = st.expon.fit(opt_depth)
+        exp_loc = 0.0
+        exp_scale = 1.0
+        loc_discrepancy = np.abs(opt_depth_loc-exp_loc)
+        scale_discrepancy = np.abs(opt_depth_scale-exp_scale)
+        print("species " + name)
+        print("exp distrib loc tol = " + str(tol))
+        print("exp distrib loc discrepancy = " + str(loc_discrepancy))
+        assert(loc_discrepancy < tol)
+        print("exp distrib scale tol = " + str(tol))
+        print("exp distrib scale discrepancy = " + str(scale_discrepancy/exp_scale))
+        assert(scale_discrepancy/exp_scale < tol)
+        ###
+
+        #check if number of lost photons is (n0* (1 - exp(-rate*t)) )
+        dNBW_dt_theo = dNBW_dt(
+            calc_chi_gamma(p_begin[name], E_f, B_f),
+                np.linalg.norm(p_begin[name]*speed_of_light))
+        exp_lost= initial_particle_number*(1.0 - np.exp(-dNBW_dt_theo*sim_time))
+        lost =  initial_particle_number-np.size(opt_depth)
+        discrepancy_lost = np.abs(exp_lost-lost)
+        print("lost fraction tol = " + str(tol))
+        print("lost fraction discrepancy = " + str(discrepancy_lost/exp_lost))
+        assert(discrepancy_lost/exp_lost < tol)
+        ###
 
 def main():
     check()
diff --git a/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution b/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution
index daceb5687..dd5098695 100644
--- a/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution
+++ b/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution
@@ -1,15 +1,15 @@
 #################################
 ####### GENERAL PARAMETERS ######
 #################################
-max_step = 10
-amr.n_cell =  64 64 64
+max_step = 5
+amr.n_cell =  32 32 16
 amr.max_grid_size = 32   # 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
 amr.plot_int = 10
 geometry.coord_sys   = 0                  # 0: Cartesian
 geometry.is_periodic = 1  1  1      # Is periodic?
-geometry.prob_lo     =  -1.e-6   -1.e-6  -1e-6  # physical domain
-geometry.prob_hi     =  1.e-6   1.e-6   1e-6
+geometry.prob_lo     =  -0.5e-6   -0.5e-6  -0.25e-6  # physical domain
+geometry.prob_hi     =  0.5e-6   0.5e-6   0.25e-6
 amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
 
 #################################
@@ -30,6 +30,7 @@ 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
+warpx.fields_to_plot = Ex
 
 #################################
 ############ PLASMA #############
@@ -42,13 +43,13 @@ particles.photon_species = p1 p2 p3 p4
 p1.species_type = "photon"
 p1.injection_style = "NUniformPerCell"
 p1.profile = "constant"
-p1.xmin = -0.5e-6
-p1.ymin = -0.5e-6
-p1.zmin = -0.5e-6
-p1.xmax = 0.5e6
-p1.ymax = 0.5e6
-p1.zmax = 0.5e6
-p1.num_particles_per_cell_each_dim = 2 2 2
+p1.xmin = -0.6e-6
+p1.ymin = -0.6e-6
+p1.zmin = -0.6e-6
+p1.xmax = 0.6e6
+p1.ymax = 0.6e6
+p1.zmax = 0.6e6
+p1.num_particles_per_cell_each_dim = 1 1 1
 p1.density = 1e19
 p1.profile = "constant"
 p1.momentum_distribution_type = "gaussian"
@@ -58,6 +59,7 @@ p1.uz_m = 0.0
 p1.ux_th = 0.
 p1.uy_th = 0.
 p1.uz_th = 0.
+p1.plot_vars = ux uy uz
 ##########QED####################
 p1.do_qed = 1
 p1.do_qed_breit_wheeler = 1
@@ -68,13 +70,13 @@ p1.qed_breit_wheeler_pos_product_species = pos_bw
 p2.species_type = "photon"
 p2.injection_style = "NUniformPerCell"
 p2.profile = "constant"
-p2.xmin = -0.5e-6
-p2.ymin = -0.5e-6
-p2.zmin = -0.5e-6
-p2.xmax = 0.5e6
-p2.ymax = 0.5e6
-p2.zmax = 0.5e6
-p2.num_particles_per_cell_each_dim = 2 2 2
+p2.xmin = -0.6e-6
+p2.ymin = -0.6e-6
+p2.zmin = -0.6e-6
+p2.xmax = 0.6e6
+p2.ymax = 0.6e6
+p2.zmax = 0.6e6
+p2.num_particles_per_cell_each_dim = 1 1 1
 p2.density = 1e19
 p2.profile = "constant"
 p2.momentum_distribution_type = "gaussian"
@@ -84,6 +86,7 @@ p2.uz_m = 0.0
 p2.ux_th = 0.
 p2.uy_th = 0.
 p2.uz_th = 0.
+p2.plot_vars = ux uy uz
 ##########QED####################
 p2.do_qed = 1
 p2.do_qed_breit_wheeler = 1
@@ -94,13 +97,13 @@ p2.qed_breit_wheeler_pos_product_species = pos_bw
 p3.species_type = "photon"
 p3.injection_style = "NUniformPerCell"
 p3.profile = "constant"
-p3.xmin = -0.5e-6
-p3.ymin = -0.5e-6
-p3.zmin = -0.5e-6
-p3.xmax = 0.5e6
-p3.ymax = 0.5e6
-p3.zmax = 0.5e6
-p3.num_particles_per_cell_each_dim = 2 2 2
+p3.xmin = -0.6e-6
+p3.ymin = -0.6e-6
+p3.zmin = -0.6e-6
+p3.xmax = 0.6e6
+p3.ymax = 0.6e6
+p3.zmax = 0.6e6
+p3.num_particles_per_cell_each_dim = 1 1 1
 p3.density = 1e19
 p3.profile = "constant"
 p3.momentum_distribution_type = "gaussian"
@@ -110,6 +113,7 @@ p3.uz_m = 10000.0
 p3.ux_th = 0.
 p3.uy_th = 0.
 p3.uz_th = 0.
+p3.plot_vars = ux uy uz
 ##########QED####################
 p3.do_qed = 1
 p3.do_qed_breit_wheeler = 1
@@ -120,13 +124,13 @@ p3.qed_breit_wheeler_pos_product_species = pos_bw
 p4.species_type = "photon"
 p4.injection_style = "NUniformPerCell"
 p4.profile = "constant"
-p4.xmin = -0.5e-6
-p4.ymin = -0.5e-6
-p4.zmin = -0.5e-6
-p4.xmax = 0.5e6
-p4.ymax = 0.5e6
-p4.zmax = 0.5e6
-p4.num_particles_per_cell_each_dim = 2 2 2
+p4.xmin = -0.6e-6
+p4.ymin = -0.6e-6
+p4.zmin = -0.6e-6
+p4.xmax = 0.6e6
+p4.ymax = 0.6e6
+p4.zmax = 0.6e6
+p4.num_particles_per_cell_each_dim = 1 1 1
 p4.density = 1e19
 p4.profile = "constant"
 p4.momentum_distribution_type = "gaussian"
@@ -136,6 +140,7 @@ p4.uz_m = 57735.02691896
 p4.ux_th = 0.
 p4.uy_th = 0.
 p4.uz_th = 0.
+p4.plot_vars = ux uy uz
 ##########QED####################
 p4.do_qed = 1
 p4.do_qed_breit_wheeler = 1
@@ -153,6 +158,7 @@ ele_bw.momentum_distribution_type = "gaussian"
 ele_bw.xmin = 1 ## Ugly trick to avoid electrons at T=0
 ele_bw.xmax = -1 ##  Ugly trick to avoid electrons at T=0
 ele_bw.do_qed = 0
+ele_bw.plot_species = 1
 
 pos_bw.species_type = "positron"
 pos_bw.injection_style = nuniformpercell
@@ -163,6 +169,7 @@ pos_bw.momentum_distribution_type = "gaussian"
 pos_bw.xmin = 1 ## Ugly trick to avoid positrons at T=0
 pos_bw.xmax = -1 ##  Ugly trick to avoid positrons at T=0
 pos_bw.do_qed = 0
+pos_bw.plot_species = 1
 #################################
 
 ##########QED TABLES####################