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
+from scipy.constants import e 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,19 @@
                     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()
+
+parser.add_argument('--fourth',
+                    dest='fourth',
+                    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 +52,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 +71,31 @@
 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_linear = f'./{args.third}/reduced_elec.0000.txt'
+all_data_linear = diag.read_file(insitu_path_linear)
+Tx2_l = all_data_linear['[ux^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+Ty2_l = all_data_linear['[uy^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+Tz2_l = all_data_linear['[uz^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+temp_eV_linear = 1./3*(Tx2_l+Ty2_l+Tz2_l)
+
+insitu_path_circular = f'./{args.fourth}/reduced_elec.0000.txt'
+all_data_circular = diag.read_file(insitu_path_circular)
+Tx2_c = all_data_circular['[ux^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+Ty2_c = all_data_circular['[uy^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+Tz2_c = all_data_circular['[uz^2]'][0,:]*scc.m_e*scc.c**2/scc.e
+temp_eV_circular = 1./3*(Tx2_c+Ty2_c+Tz2_c)
+
+temp_eV_warpx_linear = 1.00286009
+temp_eV_warpx_circular = 9.68224535
+
+print(f"temp_eV_warpx_linear = {temp_eV_warpx_linear}")
+print(f"temp_eV_hipace_linear = {temp_eV_linear}")
+
+print(f"temp_eV_warpx_circular = {temp_eV_warpx_circular}")
+print(f"temp_eV_hipace_circular = {temp_eV_circular}")
+
+relative_diff_temp_linear = np.abs( ( temp_eV_linear - temp_eV_warpx_linear ) / temp_eV_warpx_linear )
+relative_diff_temp_circular = 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_linear < tolerance) and  (relative_diff_temp_circular < tolerance)), 'Test laser_ionization did not pass'