Enable to initialize spatial/temporal/angular chirped laser + CI test

CMakeLists.txt

@@ -342,6 +342,12 @@ if(BUILD_TESTING)
                 WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
         )
 
+        add_test(NAME laser_STC.SI.1Rank
+                COMMAND bash ${HiPACE_SOURCE_DIR}/tests/laser_STC.SI.1Rank.sh
+                        $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
+                WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+        )
+
         if (NOT HiPACE_COMPUTE STREQUAL CUDA)
 
             # These tests only run on CPU

docs/source/run/parameters.rst

@@ -901,9 +901,36 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start
           Distance at which the laser pulse is focused (in the z direction, counted from laser initial position).
 
       * ``<laser name>.propagation_angle_yz`` (`float`) optional (default `0`)
-          Propagation angle of the pulse in the yz plane (0 is along the z axis)
+          Propagation angle of the pulse in the yz plane (0 is along z axis)
 
-      Option: ``from_file`` the laser is loaded from an openPMD file.
+      * ``<laser name>.STC_theta_xy`` (`float`) optional (default `pi/2`)
+          Direction of the linear spatial and angular chirp on xoy plane (0 is along x axis).
+
+      * ``<laser name>.beta`` (`float`) optional (default `0.`)
+          Angular dispersion (or angular chirp) at focus defined by `S. Akturk et al., Optics Express 12, 4399 (2004) <https://doi.org/10.1364/OPEX.12.004399>`__.
+
+      * ``<laser name>.zeta`` (`float`) optional (default `0.`)
+          Spatial chirp at focus defined by `S. Akturk et al., Optics Express 12, 4399 (2004) <https://doi.org/10.1364/OPEX.12.004399>`__.
+
+      * ``<laser name>.phi2`` (`float`) optional (default `pi/2`)
+          The amount of temporal chirp :math:`\phi^{(2)}` at focus (in the lab frame).
+          Namely, a wave packet centered on the frequency :math:`(\omega_0 + \delta \omega)` will reach its peak intensity at :math:`z(\delta \omega) = z_0 - c \phi^{(2)} \, \delta \omega`.
+          Thus, a positive :math:`\phi^{(2)}` corresponds to positive chirp, i.e., red part of the spectrum in the front of the pulse and blue part of the spectrum in the back.
+          More specifically, the electric field in the focal plane is of the form:
+
+          .. math::
+
+              E(\boldsymbol{x},t) \propto Re\left[ \exp\left(  -\frac{(t-t_{peak})^2}{\tau^2 + 2i\phi^{(2)}} + i\omega_0 (t-t_{peak}) + i\phi_0 \right) \right]
+
+          where :math:`\tau` is given by ``<laser_name>.tau`` and represents the Fourier-limited duration of the laser pulse. Thus, the actual duration of the chirped laser pulse is:
+
+          .. math::
+
+               \tau' = \sqrt{ \tau^2 + 4 (\phi^{(2)})^2/\tau^2 }
+
+         See also the definition in `S. Akturk et al., Optics Express 12, 4399 (2004) <https://doi.org/10.1364/OPEX.12.004399>`__.
+
+      Option: ``from_file`` the laser is loaded from an openPMD file
 
       * ``<laser name>.input_file`` (`string`) optional (default `""`)
           Path to an openPMD file containing a laser envelope.

examples/laser/analysis_laser_init_chirp.py

@@ -0,0 +1,97 @@
+#! /usr/bin/env python3
+
+# Copyright 2024
+#
+# This file is part of HiPACE++.
+#
+# Authors: Xingjian Hui
+# License: BSD-3-Clause-LBNL
+
+import argparse
+import numpy as np
+import scipy.constants as scc
+from openpmd_viewer.addons import LpaDiagnostics
+
+def get_zeta(Ar, m, w0, L):
+    nu = 0
+    summ = 0
+    laser_module = np.abs(Ar)
+    phi_envelop = np.array(np.arctan2(Ar.imag, Ar.real))
+    # unwrap phi_envelop
+    phi_envelop = np.unwrap(np.unwrap(phi_envelop, axis=0), axis=1)
+    # calculate pphi_pz
+    z_diff = np.diff(m.z)
+    x_diff = np.diff(m.x)
+    pphi_pz = (np.diff(phi_envelop, axis=0)).T / (z_diff / scc.c)
+    pphi_pzpy = (np.diff(pphi_pz, axis=0)).T / x_diff
+    for i in range(len(m.z) - 2):
+        for j in range(len(m.x) - 2):
+            nu = nu + pphi_pzpy[i, j] * laser_module[i, j]
+            summ = summ + laser_module[i, j]
+    nu = nu / scc.c / summ
+    a = 4 * nu * w0**2 * L**4
+    b = -4 * scc.c
+    c = nu * w0**2 * L**2
+    zeta_roots = np.roots([a, b, c])
+    return np.min(zeta_roots)
+
+def get_phi2 (Ar, m, tau):
+    # get temporal chirp phi2
+    temp_chirp = 0
+    summ = 0
+    laser_module1 = np.abs(Ar)
+    phi_envelop = np.unwrap(np.array(np.arctan2(Ar.imag, Ar.real)), axis=0)
+    # calculate pphi_pz
+    z_diff = np.diff(m.z)
+    pphi_pz = (np.diff(phi_envelop, axis=0)).T/ (z_diff/scc.c)
+    pphi_pz2 = ((np.diff(pphi_pz, axis=1)) / (z_diff[:len(z_diff)-1]) / scc.c).T
+    for i in range(len(m.z)-2):
+        for j in range(len(m.x)-2):
+            temp_chirp = temp_chirp + pphi_pz2[i,j] * laser_module1[i,j]
+            summ = summ + laser_module1[i,j]
+    x = temp_chirp * scc.c**2 / summ
+    a = 4 * x
+    b = -4
+    c = tau**4 * x
+    zeta_roots = np.roots([a, b, c])
+    return np.max(zeta_roots)
+
+def get_centroids(F, x, z):
+    index_array = np.mgrid[0:F.shape[0], 0:F.shape[1]][1]
+    centroids = np.sum(index_array * np.abs(F**2), axis=1) / np.sum(np.abs(F**2), axis=1)
+    return z[centroids.astype(int)]
+
+def get_beta(F, m, k0):
+    z_centroids = get_centroids(F.T, m.x, m.z)
+    weight = np.mean(np.abs(F.T)**2, axis = np.ndim(F) - 1)
+    derivative = np.gradient(z_centroids) / (m.x[1] - m.x[0])
+    return (np.sum(derivative * weight) / np.sum(weight)) / k0 / scc.c
+
+parser = argparse.ArgumentParser(description = 'Verify the chirp initialization')
+parser.add_argument('--output-dir',
+                    dest='output_dir',
+                    default='diags/hdf5',
+                    help='Path to the directory containing output files')
+parser.add_argument('--chirp_type',
+                    dest='chirp_type',
+                    default='phi2',
+                    help='the type of the initialized chirp')
+args = parser.parse_args()
+
+ts = LpaDiagnostics(args.output_dir)
+
+Ar, m = ts.get_field(field='laserEnvelope', iteration=0)
+lambda0 = .8e-6          # Laser wavelength
+w0 = 30.e-6              # Laser waist
+L0 = 5e-6
+tau = L0 / scc.c         # Laser duration
+k0 = 2 * scc.pi / lambda0
+if args.chirp_type == 'phi2':
+    phi2 = get_phi2(Ar, m, tau)
+    assert(np.abs(phi2 - 2.4e-26) / 2.4e-26 < 1e-2)
+elif args.chirp_type == 'zeta':
+    zeta = get_zeta(Ar, m, w0, L0)
+    assert(np.abs(zeta - 2.4e-19) / 2.4e-19 < 1e-2)
+elif args.chirp_type == 'beta':
+    beta = get_beta(Ar, m, k0)
+    assert(np.abs(beta - 2e-17) / 2e-17 < 1e-2)

examples/laser/inputs_SI

@@ -2,7 +2,7 @@ max_step = 30
 hipace.dt = 70.e-6/clight
 hipace.verbose=3
 
-amr.n_cell = 128 128 50 # 128 128 100
+amr.n_cell = 128 128 50
 
 my_constants.kp_inv = 10.e-6
 

examples/laser/inputs_chirp

@@ -0,0 +1,32 @@
+max_step = 1
+hipace.dt = 70.e-6/clight
+hipace.verbose = 3
+
+amr.n_cell = 255 127 200
+
+my_constants.kp_inv = 10.e-6
+
+hipace.file_prefix = new
+hipace.do_tiling = 0
+
+geometry.prob_lo = -15.*kp_inv   -15.*kp_inv   -2.*kp_inv
+geometry.prob_hi =  15.*kp_inv    15.*kp_inv    2.*kp_inv
+
+lasers.names = laser
+lasers.lambda0 = .8e-6
+
+laser.a0 = 1
+laser.position_mean = 0. 0. 0
+laser.focal_distance = 0.0
+laser.w0 = 30e-6
+laser.L0 = 5e-6
+amr.max_level = 0
+
+diagnostic.output_period = 1
+
+hipace.depos_order_xy = 0
+
+boundary.field = Dirichlet
+boundary.particle = Absorbing
+
+diagnostic.diag_type = xz

src/laser/Laser.H

@@ -2,8 +2,8 @@
  *
  * This file is part of HiPACE++.
  *
- * Authors: MaxThevenet, AlexanderSinn
- * Severin Diederichs, atmyers, Angel Ferran Pousa
+ * Authors: MaxThevenet, AlexanderSinn, Severin Diederichs,
+ *          atmyers, Angel Ferran Pousa, Xingjian Hui
  * License: BSD-3-Clause-LBNL
  */
 
@@ -23,6 +23,7 @@ public:
      *\param[in] laser_geom_3D 3D laser geometry for level 0
      */
     void GetEnvelopeFromFileHelper (amrex::Geometry laser_geom_3D);
+
     template<typename input_type>
     void GetEnvelopeFromFile (amrex::Geometry laser_geom_3D);
     Laser (std::string name, amrex::Geometry laser_geom_3D);
@@ -33,12 +34,15 @@ public:
     amrex::Real m_a0 {0.}; /**< Laser peak normalized amplitude */
     amrex::Real m_w0 {0.}; /**< Laser waist */
     amrex::Real m_CEP {0.}; /**< Laser carrier-envelope phase (CEP) */
-    /** Propagation angle of the pulse in the yz plane (0 is the along the z axis) */
+    /** Propagation angle of the pulse in the yz plane (0 is along z axis) */
     amrex::Real m_propagation_angle_yz {0.};
-    /**Pulse front tilt angle of the pulse in yz plane (pi/2 is no PFT)*/
-    amrex::Real m_PFT_yz {MathConst::pi/2.0};
     amrex::Real m_L0 {0.}; /**< Laser length (HW 1/e in amplitude) */
     amrex::Real m_tau {0.}; /**< Laser duration (HW 1/e in amplitude) */
+    /** Time stretching factors */
+    amrex::Real m_STC_theta_xy {MathConst::pi/2.0};
+    amrex::Real m_zeta {0.};
+    amrex::Real m_beta {0.};
+    amrex::Real m_phi2 {0.};
     /** Focal distance of the laser pulse */
     amrex::Real m_focal_distance {0.};
     /** Average position of the Gaussian laser pulse */

src/laser/Laser.cpp

@@ -35,14 +35,18 @@ Laser::Laser (std::string name,  amrex::Geometry laser_geom_3D)
         queryWithParser(pp, "w0", m_w0);
         queryWithParser(pp, "CEP", m_CEP);
         queryWithParser(pp, "propagation_angle_yz", m_propagation_angle_yz);
-        queryWithParser(pp, "PFT_yz", m_PFT_yz);
+        queryWithParser(pp, "STC_theta_xy", m_STC_theta_xy);
         bool length_is_specified = queryWithParser(pp, "L0", m_L0);
         bool duration_is_specified = queryWithParser(pp, "tau", m_tau);
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE( length_is_specified + duration_is_specified == 1,
-        "Please specify exlusively either the pulse length L0 or the duration tau of gaussian lasers");
-        if (duration_is_specified) m_L0 = m_tau*get_phys_const().c;
+        "Please specify exlusively either the pulse length L0 or the duration tau of the laser");
+        if (duration_is_specified) m_L0 = m_tau * get_phys_const().c;
+        if (length_is_specified) m_tau = m_L0 / get_phys_const().c;
         queryWithParser(pp, "focal_distance", m_focal_distance);
         queryWithParser(pp, "position_mean",  m_position_mean);
+        queryWithParser(pp, "zeta",  m_zeta);
+        queryWithParser(pp, "beta",  m_beta);
+        queryWithParser(pp, "phi2",  m_phi2);
         return;
     }
     else if (m_laser_init_type == "parser") {

src/laser/MultiLaser.cpp

@@ -874,15 +874,20 @@ MultiLaser::InitLaserSlice (const int islice, const int comp)
             }
             else if (laser.m_laser_init_type == "gaussian") {
                 const amrex::Real a0 = laser.m_a0;
-                const amrex::Real w0 = laser.m_w0;
+                const amrex::Real w0_2 = laser.m_w0 * laser.m_w0;
                 const amrex::Real cep = laser.m_CEP;
                 const amrex::Real propagation_angle_yz = laser.m_propagation_angle_yz;
-                const amrex::Real PFT_yz = laser.m_PFT_yz - MathConst::pi/2.0;
                 const amrex::Real x0 = laser.m_position_mean[0];
                 const amrex::Real y0 = laser.m_position_mean[1];
                 const amrex::Real z0 = laser.m_position_mean[2];
                 const amrex::Real L0 = laser.m_L0;
+                const amrex::Real inv_tau2 = 1/(laser.m_tau*laser.m_tau);
                 const amrex::Real zfoc = laser.m_focal_distance;
+                const amrex::Real zeta = laser.m_zeta;
+                const amrex::Real beta = laser.m_beta;
+                const amrex::Real phi2 = laser.m_phi2;
+                const amrex::Real clight = get_phys_const().c;
+                const amrex::Real theta_xy = laser.m_STC_theta_xy;
                 amrex::ParallelFor(
                 bx,
                 [=] AMREX_GPU_DEVICE(int i, int j, int k)
@@ -891,21 +896,28 @@ MultiLaser::InitLaserSlice (const int islice, const int comp)
                     const amrex::Real y = j * dx_arr[1] + poff_y - y0;
                     const amrex::Real z = islice * dx_arr[2] + poff_z - z0;
                     // Coordinate rotation in yz plane for a laser propagating at an angle.
-                    const amrex::Real yp = std::cos( propagation_angle_yz + PFT_yz ) * y \
-                        - std::sin( propagation_angle_yz + PFT_yz ) * z;
-                    const amrex::Real zp = std::sin( propagation_angle_yz + PFT_yz ) * y \
-                        + std::cos( propagation_angle_yz + PFT_yz ) * z;
+                    const amrex::Real yp = std::cos(propagation_angle_yz) * y \
+                        - std::sin(propagation_angle_yz) * z;
+                    const amrex::Real zp = std::sin(propagation_angle_yz) * y \
+                        + std::cos(propagation_angle_yz) * z;
                     // For first laser, setval to 0.
                     if (ilaser == 0) {
                         arr(i, j, k, comp ) = 0._rt;
                         arr(i, j, k, comp + 1 ) = 0._rt;
                     }
                     // Compute envelope for time step 0
-                    Complex diffract_factor = 1._rt + I * ( zp - zfoc + z0 * std::cos( propagation_angle_yz ) ) \
-                       * 2._rt/( k0 * w0 * w0 );
-                    Complex inv_complex_waist_2 = 1._rt /( w0 * w0 * diffract_factor );
+                    Complex diffract_factor = 1._rt + I * (zp - zfoc + z0 * std::cos(propagation_angle_yz)) \
+                       * 2._rt/(k0 * w0_2);
+                    Complex inv_complex_waist_2 = 1._rt /(w0_2 * diffract_factor);
+                    // Time stretching due to STCs and phi2 complex envelope
+                    // (1 if zeta=0, beta=0, phi2=0)
+                    Complex stretch_factor = 1._rt \
+                        + 4._rt * (zeta + beta * zfoc * inv_tau2) * (zeta + beta * zfoc * inv_complex_waist_2) \
+                        + 2._rt * I * (phi2 - beta * beta * k0 * zfoc) * inv_tau2;
                     Complex prefactor = a0 / diffract_factor;
-                    Complex time_exponent = zp * zp / ( L0 * L0 );
+                    Complex time_exponent = 1._rt / ( stretch_factor * L0 * L0 ) *
+                        amrex::pow(zp - beta * k0 * (x * std::cos(theta_xy) + yp * std::sin(theta_xy)) * clight  - 2._rt * I * (x * std::cos(theta_xy) + yp * std::sin(theta_xy))
+                        * (zeta - beta * zfoc) * clight * inv_complex_waist_2, 2);
                     Complex stcfactor = prefactor * amrex::exp( - time_exponent );
                     Complex exp_argument = - ( x * x + yp * yp ) * inv_complex_waist_2;
                     Complex envelope = stcfactor * amrex::exp( exp_argument ) * \

tests/laser_STC.SI.1Rank.sh

@@ -0,0 +1,60 @@
+#! /usr/bin/env bash
+
+# Copyright 2024
+#
+# This file is part of HiPACE++.
+#
+# Authors: Xingjian Hui
+# License: BSD-3-Clause-LBNL
+
+
+# This file is part of the HiPACE++ test suite.
+# It initializes a Hipace simulation of a gaussian laser with initial STC parameters
+# and test the correctness of the STC factors
+
+# Abort on first encountered error
+set -eu -o pipefail
+
+# Read input parameters
+HIPACE_EXECUTABLE=$1
+HIPACE_SOURCE_DIR=$2
+
+HIPACE_EXAMPLE_DIR=${HIPACE_SOURCE_DIR}/examples/laser
+HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests
+
+FILE_NAME=`basename "$0"`
+TEST_NAME="${FILE_NAME%.*}"
+
+# Run the simulation with initial phi2
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_chirp \
+        laser.phi2 = 2.4e-26 \
+        laser.STC_theta_xy = 0 \
+        hipace.file_prefix = $TEST_NAME
+# Compare the result with theory
+$HIPACE_EXAMPLE_DIR/analysis_laser_init_chirp.py --output-dir=$TEST_NAME \
+        --chirp_type="phi2"
+
+rm -rf $TEST_NAME
+
+# Run the simulation with initial zeta
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_chirp \
+        laser.zeta = 2.4e-19 \
+        laser.STC_theta_xy = 0 \
+        hipace.file_prefix = $TEST_NAME
+
+# Compare the result with theory
+$HIPACE_EXAMPLE_DIR/analysis_laser_init_chirp.py --output-dir=$TEST_NAME \
+        --chirp_type="zeta"
+
+rm -rf $TEST_NAME
+
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_chirp \
+        laser.beta = 2e-17 \
+        laser.STC_theta_xy = 0 \
+        hipace.file_prefix = $TEST_NAME
+
+# Compare the result with theory
+$HIPACE_EXAMPLE_DIR/analysis_laser_init_chirp.py --output-dir=$TEST_NAME \
+        --chirp_type="beta"
+
+rm -rf $TEST_NAME

tests/laser_evolution.SI.2Rank.sh

@@ -50,3 +50,5 @@ $HIPACE_TEST_DIR/checksum/checksumAPI.py \
     --rtol $RTOL \
     --file_name $TEST_NAME \
     --test-name $TEST_NAME
+
+rm -rf $TEST_NAME