Update analysis_laser_ionization.py

examples/laser_ionization/analysis_laser_ionization.py

@@ -10,7 +10,11 @@
 import math
 from openpmd_viewer import OpenPMDTimeSeries
 import statistics
-from scipy.constants import e as qe
+import scipy.constants as scc
+
+import read_insitu_diagnostics as diag
+from read_insitu_diagnostics import temperature_in_eV
+
 parser = argparse.ArgumentParser(
     description='Script to analyze the equality of two simulations')
 parser.add_argument('--first',
@@ -21,8 +25,13 @@
                     dest='second',
                     required=True,
                     help='Path to the directory containing output files')
+parser.add_argument('--third',
+                    dest='third',
+                    required=True,
+                    help='Path to the directory containing output files')
 args = parser.parse_args()
 
+# diagnostics for calculation of the fraction of ionization in linear and circular polarization
 ts_linear = OpenPMDTimeSeries(args.first)
 ts_circular = OpenPMDTimeSeries(args.second)
 
@@ -37,12 +46,12 @@
 rho_elec_linear, _ = ts_linear.get_field(field='rho_elec', coord='z', iteration=iteration)
 rho_elec_mean_linear = np.mean(rho_elec_linear, axis=(1, 2))
 rho_average_linear = statistics.mean(rho_elec_mean_linear[0:10])
-fraction_linear = -rho_average_linear / qe / n0
+fraction_linear = -rho_average_linear / scc.e / n0
 
 rho_elec_circular, _ = ts_circular.get_field(field='rho_elec', coord='z', iteration=iteration)
 rho_elec_mean_circular = np.mean(rho_elec_circular, axis=(1, 2))
 rho_average_circular = statistics.mean(rho_elec_mean_circular[0:10]) #average over a thickness in the ionized region
-fraction_circular = -rho_average_circular / qe / n0
+fraction_circular = -rho_average_circular / scc.e / n0
 
 fraction_warpx_linear = 0.41014984 # result from WarpX simulation
 fraction_warpx_circular = 0.502250841 # result from WarpX simulation
@@ -56,4 +65,19 @@
 print(f"fraction_warpx_circular = {fraction_warpx_circular}")
 print(f"fraction_hipace_circular = {fraction_circular}")
 
-assert ( (relative_diff_linear < tolerance) and (relative_diff_circular < tolerance) ), 'Test laser_ionization did not pass'
+# in-situ diagnostics for calculation of the temperature in all directions in circular polarization
+insitu_path = f'./{args.first}/reduced_elec.0000.txt'
+all_data = diag.read_file(insitu_path)
+Tx2 = all_data['[ux^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+Ty2 = all_data['[uy^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+Tz2 = all_data['[uz^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+temp_eV_circular = 1./3*(Tx2+Ty2+Tz2)
+
+temp_eV_warpx_circular = 9.68224535
+
+print(f"temp_eV_warpx_circular = {temp_eV_warpx_circular}")
+print(f"temp_eV_hipace_circular = {temp_eV_circular}")
+
+relative_diff_temp = np.abs( ( temp_eV_circular - temp_eV_warpx_circular ) / temp_eV_warpx_circular )
+
+assert ( (relative_diff_linear < tolerance) and (relative_diff_circular < tolerance) and  (relative_diff_temp < tolerance)), 'Test laser_ionization did not pass'