diff options
Diffstat (limited to 'Examples/Tests/laser_injection_from_file/analysis.py')
| -rwxr-xr-x | Examples/Tests/laser_injection_from_file/analysis.py | 113 |
1 files changed, 89 insertions, 24 deletions
diff --git a/Examples/Tests/laser_injection_from_file/analysis.py b/Examples/Tests/laser_injection_from_file/analysis.py index 818e7aaf0..4a61a6df0 100755 --- a/Examples/Tests/laser_injection_from_file/analysis.py +++ b/Examples/Tests/laser_injection_from_file/analysis.py @@ -3,13 +3,17 @@ #ADD COMMENT #ADD COMMENT #ADD COMMENT - +import yt ; yt.funcs.mylog.setLevel(50) +import numpy as np +import scipy.constants as scc +import matplotlib.pyplot as plt +from scipy.signal import hilbert import sys import glob import os -import numpy as np -import matplotlib.pyplot as plt -import yt + +#Maximum acceptable error for this test +relative_error_threshold = 0.06 #Physical parameters um = 1.e-6 @@ -24,6 +28,7 @@ x_c = 0.*um t_c = 20.*fs foc_dist = 10*um E_max = 1e12 +rot_angle = -np.pi/4.0 #Parameters of the tx grid x_l = -12.0*um @@ -54,6 +59,18 @@ def gauss(T,X,Y,opt): return np.real(pre_fact * np.exp(exp_arg)) +# Function for the envelope +def gauss_env(T,XX,ZZ): + + X = np.cos(rot_angle)*XX + np.sin(rot_angle)*ZZ + Z = -np.sin(rot_angle)*XX + np.cos(rot_angle)*ZZ + + k0 = 2.0*np.pi/wavelength + inv_tau2 = 1./tt/tt + inv_w_2 = 1.0/(w0*w0) + exp_arg = - (X*X)*inv_w_2 - inv_tau2 / c/c * (Z-T*c)*(Z-T*c) + return E_max * np.real(np.exp(exp_arg)) + def write_file(fname, x, y, t, E): with open(fname, 'wb') as file: flag_unif = 0 @@ -90,32 +107,80 @@ def create_gaussian_2d(): write_file_unf("gauss_2d_unf.txye", xcoords, np.array([0.0]), tcoords, E_t) -def do_analysis(fname): - data_set_end = yt.load(fname) - sim_time = data_set_end.current_time.to_value() - ray0 = data_set_end.ray((0*um,0*um,0), (15*um, 15*um,0)) - - xx0 = np.array(ray0["t"])*np.sqrt(2)*15*um - EE0 = np.array(ray0["Ey"])/E_max +def do_analysis(fname, compname, steps): + ds = yt.load(fname) + + dt = ds.current_time.to_value()/steps + + # Define 2D meshes + x = np.linspace( + ds.domain_left_edge[0], + ds.domain_right_edge[0], + ds.domain_dimensions[0]).v + z = np.linspace( + ds.domain_left_edge[ds.dimensionality-1], + ds.domain_right_edge[ds.dimensionality-1], + ds.domain_dimensions[ds.dimensionality-1]).v + X, Z = np.meshgrid(x, z, sparse=False, indexing='ij') + + # Compute the theory for envelope + env_theory = gauss_env(+t_c-ds.current_time.to_value(), X,Z)+gauss_env(-t_c+ds.current_time.to_value(), X,Z) + + # Read laser field in PIC simulation, and compute envelope + all_data_level_0 = ds.covering_grid(level=0,left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) + F_laser = all_data_level_0['boxlib', 'Ey'].v.squeeze() + env = abs(hilbert(F_laser)) + extent = [ds.domain_left_edge[ds.dimensionality-1], ds.domain_right_edge[ds.dimensionality-1], + ds.domain_left_edge[0], ds.domain_right_edge[0] ] + + # Plot results + plt.figure(figsize=(8,6)) + plt.subplot(221) + plt.title('PIC field') + plt.imshow(F_laser, extent=extent) + plt.colorbar() + plt.subplot(222) + plt.title('PIC envelope') + plt.imshow(env, extent=extent) + plt.colorbar() + plt.subplot(223) + plt.title('Theory envelope') + plt.imshow(env_theory, extent=extent) + plt.colorbar() + plt.subplot(224) + plt.title('Difference') + plt.imshow(env-env_theory, extent=extent) + plt.colorbar() + plt.tight_layout() + plt.savefig(compname, bbox_inches='tight') + + relative_error_env = np.sum(np.abs(env-env_theory)) / np.sum(np.abs(env)) + print("Relative error envelope: ", relative_error_env) + assert(relative_error_env < relative_error_threshold) + + fft_F_laser = np.fft.fft2(F_laser) + + freq_rows = np.fft.fftfreq(F_laser.shape[0],dt) + freq_cols = np.fft.fftfreq(F_laser.shape[1],dt) + + pos_max = np.unravel_index(np.abs(fft_F_laser).argmax(), fft_F_laser.shape) + + freq = np.sqrt((freq_rows[pos_max[0]])**2 + (freq_cols[pos_max[1]]**2)) + exp_freq = c/wavelength + + relative_error_freq = np.abs(freq-exp_freq)/exp_freq + print("Relative error frequency: ", relative_error_freq) + assert(relative_error_freq < relative_error_threshold) - expected0 = [-gauss((sim_time)-x/c , 0, 0, '2d') for x in xx0] - - #DEBUG - plt.plot(xx0,EE0,'bo') - plt.plot(xx0,expected0,'ro') - plt.savefig('graph.png') - #__DEBUG__ - - return True def launch_analysis(executable): create_gaussian_2d() os.system("./" + executable + " inputs.2d_test_txye") - assert(do_analysis("diags/plotfiles/plt00507/")) - os.system("sed -i 's/gauss_2d_unf.txye/gauss_2d.txye/g' inputs.2d_test_txye") - os.system("./" + executable + " inputs.2d_test_txye") - assert(do_analysis("diags/plotfiles/plt00507/")) + do_analysis("diags/plotfiles/plt00250/", "comp_unf.pdf", 250) + os.system("sed 's/gauss_2d_unf.txye/gauss_2d.txye/g' inputs.2d_test_txye > inputs.2d_test_txye_non_unf") + os.system("./" + executable + " inputs.2d_test_txye_non_unf") + do_analysis("diags/plotfiles/plt00250/", "comp_non_unf.pdf", 250) def main() : |