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

.github/workflows/cudaLocal.yml

@@ -46,4 +46,5 @@ jobs:
       run: |
         source $HOME/profile.hipace
         conda activate openpmd
+        python -m pip install git+https://github.com/huixingjian/lasy-.git@fix_bug_STC
         ctest --test-dir build --output-on-failure

.github/workflows/setup/ubuntu.sh

@@ -32,3 +32,4 @@ sudo update-alternatives --set python /usr/bin/python3
 
 python -m pip install --upgrade pip
 python -m pip install --upgrade matplotlib numpy scipy openpmd-viewer openpmd-api
+python -m pip install git+https://github.com/huixingjian/lasy-.git@fix_bug_STC

.github/workflows/setup/ubuntu_ompi.sh

@@ -34,3 +34,4 @@ sudo update-alternatives --set python /usr/bin/python3
 
 python -m pip install --upgrade pip
 python -m pip install --upgrade matplotlib numpy scipy openpmd-viewer openpmd-api
+python -m pip install git+https://github.com/huixingjian/lasy-.git@fix_bug_STC

CMakeLists.txt

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

docs/source/run/parameters.rst

@@ -911,6 +911,29 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start
       * ``<laser name>.propagation_angle_yz`` (`float`) optional (default `0`)
           Propagation angle of the pulse in the yz plane (0 is along the z axis)
 
+      * ``<laser name>.STC_theta_xy`` (`float`) optional (default `0`)
+          Direction of the linear spatial and angular chirps in the xy plane (`0` is along x, `pi/2` along y).
+          In what follows, all chirps are given as defined in `S. Akturk et al., Optics Express 12, 4399 (2004) <https://doi.org/10.1364/OPEX.12.004399>`__.
+
+      * ``<laser name>.beta`` (`float`) optional (default `0.`)
+          Angular dispersion (or angular chirp) at focus in :math:`second`.
+
+      * ``<laser name>.zeta`` (`float`) optional (default `0.`)
+          Spatial chirp at focus in :math:`second \cdot meter`.
+
+      * ``<laser name>.phi2`` (`float`) optional (default `pi/2`)
+          Temporal chirp :math:`\phi^{(2)}` at focus in :math:`second^2`.
+          Namely, a wave packet centered on frequency :math:`(\omega_0 + \delta \omega)` reaches 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 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 }
+
       Option: ``from_file`` the laser is loaded from an openPMD file.
 
       * ``<laser name>.input_file`` (`string`) optional (default `""`)

examples/laser_STC/analysis_laser_STC.py

@@ -0,0 +1,43 @@
+#! /usr/bin/env python3
+
+# Copyright 2025
+#
+# 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 lasy.utils.laser_utils import get_phi2, get_zeta, get_beta
+from lasy.laser import Laser
+from lasy.profiles import FromOpenPMDProfile
+
+lambda0 = .6e-6          # Laser wavelength
+k0 = 2 * scc.pi / lambda0
+
+parser = argparse.ArgumentParser(description='Compare laser propagation in vacuum with theory')
+parser.add_argument('--output-dir',
+                    dest='output_dir',
+                    default='diags/hdf5',
+                    help='Path to the directory containing output files')
+args = parser.parse_args()
+
+profile = FromOpenPMDProfile (path = args.output_dir, iteration = 0, pol=[1,0], field='laserEnvelope'\
+                              ,coord='', is_envelope=True, prefix='openpmd')
+laser = Laser(
+        dim="xyt",
+        lo=(-15e-6, -15e-6, -30e-15),
+        hi=(15e-6, 15e-6, +30e-15),
+        npoints=(127, 127, 200),
+        profile=profile,
+     )
+
+Phi2,phi2 = get_phi2(laser.dim, laser.grid)
+[beta_x, beta_y] = get_beta(laser.dim, laser.grid, k0)
+[zeta_x, zeta_y], [nu_x, nu_y] = get_zeta(laser.dim, laser.grid, k0)
+
+assert ((phi2-2.4e-24)/2.4e-24<0.01), 'Test phi2 did not pass'
+assert ((zeta_y-2.4e-22)/2.4e-22<0.01), 'Test zeta did not pass'
+assert ((beta_y-3e-18)/3e-18<0.01), 'Test zeta did not pass'

examples/laser_STC/inputs_STC

@@ -0,0 +1,27 @@
+max_step = 1
+amr.n_cell = 127 127 200
+
+hipace.verbose = 1
+hipace.dt = 1e-13
+diagnostic.output_period = 1
+
+amr.max_level = 0
+boundary.field = Dirichlet
+boundary.particle = Absorbing
+geometry.prob_lo     = -10e-6 -10e-6 -30e-6
+geometry.prob_hi     =  10e-6  10e-6  30e-6
+
+lasers.names = laser
+lasers.lambda0 = 600e-9
+lasers.solver_type = fft
+laser.a0 = 10
+laser.position_mean = 0. 0. 0
+laser.w0 = 5e-6
+laser.L0 = 15e-6
+laser.focal_distance = 0
+laser.phi2 = 2.4e-24
+laser.beta = 3e-18
+laser.zeta = 2.4e-22
+laser.STC_theta_xy = 3.14150265358/2
+diagnostic.field_data = laserEnvelope
+diagnostic.diag_type = xyz

src/laser/Laser.H

@@ -33,12 +33,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 the 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 {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,19 @@ 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");
+       "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;
         if (duration_is_specified) m_L0 = m_tau*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

@@ -882,15 +882,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 inv_tau2 = 1/(laser.m_tau*laser.m_tau);
                 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 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)
@@ -899,21 +904,29 @@ 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.1Rank.sh

@@ -0,0 +1,29 @@
+#! /usr/bin/env bash
+#
+# This file is part of HiPACE++.
+#
+# Authors: Xingjian Hui
+# This file is part of the HiPACE++ test suite.
+# It tests the STC initialisation with LASY functions.
+
+# abort on first encounted error
+set -eu -o pipefail
+
+# Read input parameters
+HIPACE_EXECUTABLE=$1
+HIPACE_SOURCE_DIR=$2
+HIPACE_COMPUTE=$3
+
+HIPACE_EXAMPLE_DIR=${HIPACE_SOURCE_DIR}/examples/laser_STC
+HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests
+
+FILE_NAME=`basename "$0"`
+TEST_NAME="${FILE_NAME%.*}"
+
+rm -rf $TEST_NAME
+
+# Run the simulation
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_STC \
+        hipace.file_prefix = $TEST_NAME
+# Compare the result with theory
+$HIPACE_EXAMPLE_DIR/analysis_laser_STC.py --output-dir=$TEST_NAME