Merge branch 'dev' of https://github.com/ECP-WarpX/WarpX into ParticlesEBParser

2 files changed, 262 insertions, 0 deletions

diff --git a/Examples/Modules/laser_injection_from_file/analysis.py b/Examples/Modules/laser_injection_from_file/analysis.py
new file mode 100755
index 000000000..d1ef8a001
--- /dev/null
+++ b/Examples/Modules/laser_injection_from_file/analysis.py
@@ -0,0 +1,216 @@
+#!/usr/bin/python3
+
+# This file is part of the WarpX automated test suite. It is used to test the
+# injection of a laser pulse from an external binary file.
+#
+# - Generate an input binary file with a gaussian laser pulse.
+# - Run the WarpX simulation for time T, when the pulse is fully injected
+# - Compute the theory for laser envelope at time T
+# - Compare theory and simulation, for both envelope and central frequency
+
+import yt ; yt.funcs.mylog.setLevel(50)
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import hilbert
+import glob
+import os
+
+#Maximum acceptable error for this test
+relative_error_threshold = 0.06
+
+#Physical parameters
+um = 1.e-6
+fs = 1.e-15
+c = 299792458
+
+#Parameters of the gaussian beam
+wavelength = 1.*um
+w0 = 6.*um
+tt = 10.*fs
+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
+x_r = 12.0*um
+x_points = 480
+t_l = 0.0*fs
+t_r = 40.0*fs
+t_points = 400
+tcoords = np.linspace(t_l, t_r, t_points)
+xcoords = np.linspace(x_l, x_r, x_points)
+
+def gauss(T,X,Y,opt):
+    """Compute the electric field for a Gaussian laser pulse.
+       This is used to write the binary input file.
+    """
+
+    k0 = 2.0*np.pi/wavelength
+    inv_tau2 = 1./tt/tt
+    osc_phase = k0*c*(T-t_c)
+
+    diff_factor = 1.0 + 1.0j* foc_dist * 2/(k0*w0*w0)
+    inv_w_2 = 1.0/(w0*w0*diff_factor)
+
+    pre_fact = np.exp(1.0j * osc_phase)
+
+    if opt == '3d':
+        pre_fact = pre_fact/diff_factor
+    else:
+        pre_fact = pre_fact/np.sqrt(diff_factor)
+
+    exp_arg = - (X*X + Y*Y)*inv_w_2 - inv_tau2 * (T-t_c)*(T-t_c)
+
+    return np.real(pre_fact * np.exp(exp_arg))
+
+# Function for the envelope
+def gauss_env(T,XX,ZZ):
+    '''Function to compute the theory for the envelope
+    '''
+
+    X = np.cos(rot_angle)*XX + np.sin(rot_angle)*ZZ
+    Z = -np.sin(rot_angle)*XX + np.cos(rot_angle)*ZZ
+
+    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):
+    """ For a given filename fname, space coordinates x and y, time coordinate t
+    and field E, write a WarpX-compatible input binary file containing the
+    profile of the laser pulse
+    """
+
+    with open(fname, 'wb') as file:
+        flag_unif = 0
+        file.write(flag_unif.to_bytes(1, byteorder='little'))
+        file.write((len(t)).to_bytes(4, byteorder='little', signed=False))
+        file.write((len(x)).to_bytes(4, byteorder='little', signed=False))
+        file.write((len(y)).to_bytes(4, byteorder='little', signed=False))
+        file.write(t.tobytes())
+        file.write(x.tobytes())
+        file.write(y.tobytes())
+        file.write(E.tobytes())
+
+
+def write_file_unf(fname, x, y, t, E):
+    """ For a given filename fname, space coordinates x and y, time coordinate t
+    and field E, write a WarpX-compatible input binary file containing the
+    profile of the laser pulse. This function should be used in the case
+    of a uniform spatio-temporal mesh
+    """
+
+    with open(fname, 'wb') as file:
+        flag_unif = 1
+        file.write(flag_unif.to_bytes(1, byteorder='little'))
+        file.write((len(t)).to_bytes(4, byteorder='little', signed=False))
+        file.write((len(x)).to_bytes(4, byteorder='little', signed=False))
+        file.write((len(y)).to_bytes(4, byteorder='little', signed=False))
+        file.write(t[0].tobytes())
+        file.write(t[-1].tobytes())
+        file.write(x[0].tobytes())
+        file.write(x[-1].tobytes())
+        if len(y) == 1 :
+            file.write(y[0].tobytes())
+        else :
+            file.write(y[0].tobytes())
+            file.write(y[-1].tobytes())
+        file.write(E.tobytes())
+
+def create_gaussian_2d():
+   T, X, Y = np.meshgrid(tcoords, xcoords, np.array([0.0]), indexing='ij')
+   E_t = gauss(T,X,Y,'2d')
+   write_file("gauss_2d.txye", xcoords, np.array([0.0]), tcoords, E_t)
+   write_file_unf("gauss_2d_unf.txye", xcoords, np.array([0.0]), tcoords, E_t)
+
+
+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)
+
+
+
+def launch_analysis(executable):
+    create_gaussian_2d()
+    os.system("./" + executable + " inputs.2d_test_txye")
+    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() :
+    executables = glob.glob("main2d*")
+    if len(executables) == 1 :
+        launch_analysis(executables[0])
+    else :
+        assert(False)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/Examples/Modules/laser_injection_from_file/inputs.2d_test_txye b/Examples/Modules/laser_injection_from_file/inputs.2d_test_txye
new file mode 100644
index 000000000..63ef9813a
--- /dev/null
+++ b/Examples/Modules/laser_injection_from_file/inputs.2d_test_txye
@@ -0,0 +1,46 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 400
+amr.n_cell =  672 672
+amr.plot_int = 50
+amr.max_grid_size = 512
+amr.blocking_factor = 32
+amr.max_level = 0
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 1 1 1     # Is periodic?
+geometry.prob_lo     = -25.e-6  -25.0e-6     # physical domain
+geometry.prob_hi     =  25.e-6   25.e-6
+warpx.verbose = 1
+warpx.serialize_ics = 1
+
+#################################
+############ NUMERICS ###########
+#################################
+interpolation.nox = 3
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.cfl = 0.98
+warpx.do_dynamic_scheduling = 0
+warpx.load_balance_int = -1
+warpx.use_filter = 0
+algo.maxwell_fdtd_solver = ckc
+
+#################################
+############ PLASMA #############
+#################################
+particles.nspecies = 0
+
+#################################
+############# LASER #############
+#################################
+lasers.nlasers      = 1
+lasers.names        = txye_laser
+txye_laser.position     = 0. 0. 0. # This point is on the laser plane
+txye_laser.direction    = 1. 0. 1.     # The plane normal direction
+txye_laser.polarization = 0. 1. 0.     # The main polarization vector
+txye_laser.e_max        = 1.e12        # Maximum amplitude of the laser field (in V/m)
+txye_laser.wavelength = 1.0e-6         # The wavelength of the laser (in meters)
+txye_laser.profile      = from_txye_file
+txye_laser.txye_file_name = "gauss_2d_unf.txye"
+txye_laser.time_chunk_size = 50