diff --git a/examples/laser_ionization/analysis_laser_ionization.py b/examples/laser_ionization/analysis_laser_ionization.py index 46b562a946..798b60571e 100755 --- a/examples/laser_ionization/analysis_laser_ionization.py +++ b/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'