Read wavelength and flag 'isLaserEnvelope' from lasy file

Docs/source/usage/parameters.rst

@@ -1209,17 +1209,17 @@ Laser initialization
       though ``<laser_name>.wavelength`` and ``<laser_name>.e_max`` should be included in the laser
       function, they still have to be specified as they are used for numerical purposes.
     - ``"from_file"``: the electric field of the laser is read from an external file. Currently both
-      the `lasy <https://lasydoc.readthedocs.io/en/latest/>` format as well as a custom binary format are supported. It requires to provide
+      the `lasy <https://lasydoc.readthedocs.io/en/latest/>`_ format as well as a custom binary format are supported. It requires to provide
       the name of the file to load setting the additional parameter ``<laser_name>.binary_file_name`` or ``<laser_name>.lasy_file_name`` (`string`).
-      It accepts an optional parameter ``<laser_name>.time_chunk_size`` (`int`) , only supported for a lasy file;
-      this allows to read only time_chunk_size timesteps from the lasy file. New timesteps are read as soon as they are needed.
+      It accepts an optional parameter ``<laser_name>.time_chunk_size`` (`int`), supported for both lasy and binary files;
+      this allows to read only time_chunk_size timesteps from the file. New timesteps are read as soon as they are needed.
 
-      The default value is automatically set to the number of timesteps contained in the lasy file
+      The default value is automatically set to the number of timesteps contained in the file
       (i.e. only one read is performed at the beginning of the simulation).
       It also accepts the optional parameter ``<laser_name>.delay`` (`float`; in seconds), which allows
       delaying (``delay > 0``) or anticipating (``delay < 0``) the laser by the specified amount of time.
 
-      Details about the usage of the lasy format: A lasy file is always 3D, but in the case where WarpX is compiled in 2D (or 1D), the laser antenna
+      Details about the usage of the lasy format: A lasy file can be read in 3D cartesian or in RZ geometry. In the cartesian geometry, in the case where WarpX is compiled in 2D (or 1D), the laser antenna
       will emit the field values that correspond to y=0 in the lasy file (and x=0 in the 1D case).
       One can generate a lasy file from Python, see an example at ``Examples/Tests/laser_injection_from_file``.
 
@@ -1229,16 +1229,14 @@ Laser initialization
       ``<laser_name>.e_max`` (i.e. in most cases the maximum of abs(E(x,y,t)) should be 1,
       so that the maximum field intensity can be set straightforwardly with ``<laser_name>.e_max``).
       The binary file has to respect the following format:
-
         * flag to indicate the grid is uniform(1 byte, 0 means non-uniform, !=0 means uniform) - only uniform is supported
-        * np, numbrer of timesteps (uint32_t, must be >=2)
+        * nt, number of timesteps (uint32_t, must be >=2)
         * nx, number of points along x (uint32_t, must be >=2)
         * ny, number of points along y (uint32_t, must be 1 for 2D simulations and >=2 for 3D simulations)
         * timesteps (double[2]=[t_min,t_max])
         * x_coords (double[2]=[x_min,x_max])
-        * y_coords (double[1] if 2D, double[2]=[y_min,y_max] if 3D)
+        * y_coords (double[1] in 2D, double[2]=[y_min,y_max] in 3D)
         * field_data (double[nt * nx * ny], with nt being the slowest coordinate).
-
       A binary file can be generated from Python, see an example at ``Examples/Tests/laser_injection_from_file``
 
 * ``<laser_name>.profile_t_peak`` (`float`; in seconds)

Examples/Tests/laser_injection_from_file/analysis_1d.py

@@ -25,6 +25,7 @@
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.constants import c, epsilon_0
 from scipy.signal import hilbert
 
 import yt ; yt.funcs.mylog.setLevel(50)
@@ -46,7 +47,8 @@
 tt = 10.*fs
 t_c = 20.*fs
 
-E_max= 16282454014843.37
+laser_energy = 1.0
+E_max = np.sqrt( 2*(2/np.pi)**(3/2)*laser_energy / (epsilon_0*w0**2*c*tt) )
 
 # Function for the envelope
 def gauss_env(T,Z):
@@ -120,7 +122,6 @@ def main() :
 
     # Create a laser using lasy
     pol = (1, 0)
-    laser_energy = 1.0  # J
     profile = GaussianProfile(wavelength, pol, laser_energy, w0, tt, t_peak=0)
     dim = "xyt"
     lo = (-25e-6, -25e-6, -20e-15)

Examples/Tests/laser_injection_from_file/analysis_2d.py

@@ -25,6 +25,7 @@
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.constants import c, epsilon_0
 from scipy.signal import hilbert
 
 import yt ; yt.funcs.mylog.setLevel(50)
@@ -46,7 +47,8 @@
 tt = 10.*fs
 t_c = 20.*fs
 
-E_max= 16282454014843.37
+laser_energy = 1.0
+E_max = np.sqrt( 2*(2/np.pi)**(3/2)*laser_energy / (epsilon_0*w0**2*c*tt) )
 
 # Function for the envelope
 def gauss_env(T, X, Y, Z):
@@ -138,7 +140,6 @@ def main() :
 
     # Create a laser using lasy
     pol = (1, 0)
-    laser_energy = 1.0  # J
     profile = GaussianProfile(wavelength, pol, laser_energy, w0, tt, t_peak=0)
     dim = "xyt"
     lo = (-25e-6, -25e-6, -20e-15)

Examples/Tests/laser_injection_from_file/analysis_3d.py

@@ -25,6 +25,7 @@
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.constants import c, epsilon_0
 from scipy.signal import hilbert
 
 import yt ; yt.funcs.mylog.setLevel(50)
@@ -46,7 +47,8 @@
 tt = 10.*fs
 t_c = 20.*fs
 
-E_max= 16282454014843.37
+laser_energy = 1.0
+E_max = np.sqrt( 2*(2/np.pi)**(3/2)*laser_energy / (epsilon_0*w0**2*c*tt) )
 
 # Function for the envelope
 def gauss_env(T, X, Y, Z):
@@ -142,7 +144,6 @@ def main() :
 
     # Create a laser using lasy
     pol = (1, 0)
-    laser_energy = 1.0  # J
     profile = GaussianProfile(wavelength, pol, laser_energy, w0, tt, t_peak=0)
     dim = "xyt"
     lo = (-25e-6, -25e-6, -20e-15)

Examples/Tests/laser_injection_from_file/analysis_RZ.py

@@ -25,6 +25,7 @@
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.constants import c, epsilon_0
 from scipy.signal import hilbert
 
 import yt ; yt.funcs.mylog.setLevel(50)
@@ -46,7 +47,8 @@
 tt = 10.*fs
 t_c = 20.*fs
 
-E_max= 16282454014843.37
+laser_energy = 1.0
+E_max = np.sqrt( 2*(2/np.pi)**(3/2)*laser_energy / (epsilon_0*w0**2*c*tt) )
 
 # Function for the envelope
 def gauss_env(T, X, Y, Z):
@@ -139,7 +141,6 @@ def main() :
 
     # Create a laser using lasy
     pol = (1, 0)
-    laser_energy = 1.0  # J
     profile = GaussianProfile(wavelength, pol, laser_energy, w0, tt, t_peak=0)
     dim = "xyt"
     lo = (-25e-6, -25e-6, -20e-15)

Examples/Tests/laser_injection_from_file/analysis_from_RZ_file.py

@@ -0,0 +1,183 @@
+#!/usr/bin/env python3
+
+# Copyright 2020 Andrew Myers, Axel Huebl, Luca Fedeli
+# Remi Lehe, Ilian Kara-Mostefa
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+
+# This file is part of the WarpX automated test suite. It is used to test the
+# injection of a laser pulse from an external lasy file.
+#
+# - Generate an input lasy 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 in RZ, for both envelope and central frequency
+
+import glob
+import os
+import sys
+
+import matplotlib
+
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.constants import c, epsilon_0
+from scipy.signal import hilbert
+from scipy.special import genlaguerre
+
+import yt ; yt.funcs.mylog.setLevel(50)
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+#Maximum acceptable error for this test
+relative_error_threshold = 0.065
+
+#Physical parameters
+um = 1.e-6
+fs = 1.e-15
+c = 299792458
+
+#Parameters of the gaussian beam
+wavelength = 1.*um
+w0 = 12.*um
+tt = 10.*fs
+t_c = 20.*fs
+
+laser_energy = 1.0
+E_max = np.sqrt( 2*(2/np.pi)**(3/2)*laser_energy / (epsilon_0*w0**2*c*tt) )
+
+# Function for the envelope
+def laguerre_env(T, X, Y, Z, p, m):
+    if m>0:
+        complex_position= X -1j * Y
+    else:
+        complex_position= X +1j * Y
+    inv_w0_2 = 1.0/(w0**2)
+    inv_tau2 = 1.0/(tt**2)
+    radius = abs(complex_position)
+    scaled_rad_squared = (radius**2)*inv_w0_2
+    envelope = (
+            ( np.sqrt(2) * complex_position / w0 )** m
+            *  genlaguerre(p, m)(2 * scaled_rad_squared)
+            * np.exp(-scaled_rad_squared)
+            * np.exp(-( inv_tau2 / (c**2) ) * (Z-T*c)**2)
+        )
+    return E_max * np.real(envelope)
+
+def do_analysis(fname, compname, steps):
+    ds = yt.load(fname)
+    dt = ds.current_time.to_value()/steps
+
+    # Define 3D meshes
+    x = np.linspace(
+        ds.domain_left_edge[0],
+        ds.domain_right_edge[0],
+        ds.domain_dimensions[0]).v
+    y = np.linspace(
+        ds.domain_left_edge[1],
+        ds.domain_right_edge[1],
+        ds.domain_dimensions[1]).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, Y, Z = np.meshgrid(x, y, z, sparse=False, indexing='ij')
+
+    # Compute the theory for envelope
+    env_theory = laguerre_env(+t_c-ds.current_time.to_value(), X,Y,Z,p=0,m=1)+laguerre_env(-t_c+ds.current_time.to_value(), X,Y,Z,p=0,m=1)
+
+    # 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', 'Et'].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] ]
+
+    env_theory_slice= env_theory[:,env_theory.shape[1]//2, :]
+
+    # 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_slice, extent=extent)
+    plt.colorbar()
+    plt.subplot(224)
+    plt.title('Difference')
+    plt.imshow(env-env_theory_slice, extent=extent)
+    plt.colorbar()
+    plt.tight_layout()
+    plt.savefig(compname, bbox_inches='tight')
+
+    relative_error_env = np.sum(np.abs(env-env_theory_slice)) / np.sum(np.abs(env))
+    print("Relative error envelope: ", relative_error_env)
+    assert(relative_error_env < relative_error_threshold)
+
+    fft_F_laser = np.fft.fftn(F_laser)
+
+    freq_x = np.fft.fftfreq(F_laser.shape[0],dt)
+    freq_z = 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_x[pos_max[0]])**2 + (freq_z[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):
+    os.system("./" + executable + " inputs.from_RZ_file_test diag1.file_prefix=diags/plotfiles/plt")
+    do_analysis("diags/plotfiles/plt000628/", "comp_unf.pdf", 628)
+
+
+def main() :
+
+    from lasy.laser import Laser
+    from lasy.profiles import (CombinedLongitudinalTransverseProfile,
+                               GaussianProfile)
+    from lasy.profiles.longitudinal import GaussianLongitudinalProfile
+    from lasy.profiles.transverse import LaguerreGaussianTransverseProfile
+
+    # Create a Laguerre Gaussian laser in RZ geometry using lasy
+    pol = (1, 0)
+    profile = CombinedLongitudinalTransverseProfile(
+    wavelength,pol,laser_energy,
+    GaussianLongitudinalProfile(wavelength, tt, t_peak=0),
+    LaguerreGaussianTransverseProfile(w0, p=0, m=1),
+    )
+    dim = "rt"
+    lo = (0e-6, -20e-15)
+    hi = (+25e-6, +20e-15)
+    npoints = (100,100)
+    laser = Laser(dim, lo, hi, npoints, profile, n_azimuthal_modes=2)
+    laser.normalize(laser_energy, kind="energy")
+    laser.write_to_file("laguerrelaserRZ")
+    executables = glob.glob("*.ex")
+    if len(executables) == 1 :
+        launch_analysis(executables[0])
+    else :
+        assert(False)
+
+    # Do the checksum test
+    filename_end = "diags/plotfiles/plt000628/"
+    test_name = "LaserInjectionFromRZLASYFile"
+    checksumAPI.evaluate_checksum(test_name, filename_end)
+    print('Passed')
+
+if __name__ == "__main__":
+    main()

Examples/Tests/laser_injection_from_file/inputs.from_RZ_file_test

@@ -0,0 +1,50 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+stop_time = 40.e-15
+amr.n_cell = 32 1024
+amr.max_grid_size = 512
+amr.blocking_factor = 32
+amr.max_level = 0
+geometry.dims = RZ
+geometry.prob_lo     = 0.e-6  -10.e-6   # physical domain
+geometry.prob_hi     =  25.e-6   10.e-6
+warpx.verbose = 1
+warpx.serialize_initial_conditions = 1
+
+#################################
+####### Boundary condition ######
+#################################
+boundary.field_lo = none periodic
+boundary.field_hi = pec periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.cfl = 0.98
+warpx.use_filter = 0
+algo.load_balance_intervals = -1
+warpx.n_rz_azimuthal_modes=3
+
+# Order of particle shape factors
+algo.particle_shape = 3
+
+#################################
+############# LASER #############
+#################################
+lasers.names        = lasy_RZ_laser
+lasy_RZ_laser.position     = 0. 0. 0.     # This point is on the laser plane
+lasy_RZ_laser.direction    = 0. 0. 1.     # The plane normal direction
+lasy_RZ_laser.polarization = 0. 1. 0.     # The main polarization vector
+lasy_RZ_laser.e_max        = 1.e14        # Maximum amplitude of the laser field (in V/m)
+lasy_RZ_laser.wavelength = 1.0e-6         # The wavelength of the laser (in meters)
+lasy_RZ_laser.profile      = from_file
+lasy_RZ_laser.time_chunk_size = 50
+lasy_RZ_laser.lasy_file_name = "laguerrelaserRZ_00000.h5"
+lasy_RZ_laser.delay = 0.0
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = -1
+diag1.fields_to_plot = Er Et Ez Br Bz jr jt jz
+diag1.diag_type = Full

Regression/Checksum/benchmarks_json/LaserInjectionFromRZLASYFile.json

@@ -0,0 +1,12 @@
+{
+  "lev=0": {
+    "Br": 200017483.11124495,
+    "Bz": 4853774.363010781,
+    "Er": 195991045571624.53,
+    "Et": 5.673443453613198e+16,
+    "Ez": 1434607603497291.2,
+    "jr": 0.0,
+    "jt": 9.338788430554544,
+    "jz": 0.0
+  }
+}