Read laser profile from parser

docs/source/run/parameters.rst

@@ -809,9 +809,6 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start
     The names of the laser pulses, separated by a space.
     To run without a laser, choose the name ``no_laser``.
 
-* ``lasers.lambda0`` (`float`)
-    Wavelength of the laser pulses. Currently, all pulses must have the same wavelength.
-
 * ``lasers.use_phase`` (`bool`) optional (default `true`)
     Whether the phase terms (:math:`\theta` in Eq. (6) of [C. Benedetti et al. Plasma Phys. Control. Fusion 60.1: 014002 (2017)]) are computed and used in the laser envelope advance. Keeping the phase should be more accurate, but can cause numerical issues in the presence of strong depletion/frequency shift.
 
@@ -837,22 +834,23 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start
 * ``lasers.MG_average_rhs`` (`0` or `1`) optional (default `1`)
     Whether to use the most stable discretization for the envelope solver.
 
-* ``lasers.input_file`` (`string`) optional (default `""`)
-    Path to an openPMD file containing a laser envelope.
-    The file should comply with the `LaserEnvelope extension of the openPMD-standard <https://github.com/openPMD/openPMD-standard/blob/upcoming-2.0.0/EXT_LaserEnvelope.md>`__, as generated by `LASY <https://github.com/LASY-org/LASY>`__.
-    Currently supported geometries: 3D or cylindrical profiles with azimuthal decomposition.
-    The laser pulse is injected in the HiPACE++ simulation so that the beginning of the temporal profile from the file corresponds to the head of the simulation box, and time (in the file) is converted to space (HiPACE++ longitudinal coordinate) with ``z = -c*t + const``.
-    If this parameter is set, then the file is used to initialize all lasers instead of using a gaussian profile.
+* ``<laser name>.init_type`` (list of `string`) optional (default `gaussian`)
+    The initializing method of laser. Possible options are:
 
-* ``lasers.openPMD_laser_name`` (`string`) optional (default `laserEnvelope`)
-    Name of the laser envelope field inside the openPMD file to be read in.
+      * ``gaussian``(default) the laser is iniliatized with an ideal gaussian pulse.
 
-* ``lasers.iteration`` (`int`) optional (default `0`)
-    Iteration of the openPMD file to be read in.
+      * ``from_file``, the laser is loaded from an openPMD file.
+
+      *  ``parser``, the laser is initialized with the expression of the complex envelope function.
+
+Option: ``gaussian``
 
 * ``<laser name>.a0`` (`float`) optional (default `0`)
     Peak normalized vector potential of the laser pulse.
 
+* ``lasers.lambda0`` (`float`)
+    Wavelength of the laser pulses. Currently, all pulses must have the same wavelength.
+
 * ``<laser name>.position_mean`` (3 `float`) optional (default `0 0 0`)
     The mean position of the laser in `x, y, z`.
 
@@ -876,6 +874,32 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start
     Pulse front tilt angle on yz plane - the angle between the pulse front (maximum intensity contour)and the propagation
     direction defined by [Selcuk Akturk Opt. Express 12 (2004)](pi/2 is no PFT)
 
+Option: ``from_file``
+
+* ``lasers.input_file`` (`string`) optional (default `""`)
+    Path to an openPMD file containing a laser envelope.
+    The file should comply with the `LaserEnvelope extension of the openPMD-standard <https://github.com/openPMD/openPMD-standard/blob/upcoming-2.0.0/EXT_LaserEnvelope.md>`__, as generated by `LASY <https://github.com/LASY-org/LASY>`__.
+    Currently supported geometries: 3D or cylindrical profiles with azimuthal decomposition.
+    The laser pulse is injected in the HiPACE++ simulation so that the beginning of the temporal profile from the file corresponds to the head of the simulation box, and time (in the file) is converted to space (HiPACE++ longitudinal coordinate) with ``z = -c*t + const``.
+    If this parameter is set, then the file is used to initialize all lasers instead of using a gaussian profile.
+
+* ``lasers.openPMD_laser_name`` (`string`) optional (default `laserEnvelope`)
+    Name of the laser envelope field inside the openPMD file to be read in.
+
+* ``lasers.iteration`` (`int`) optional (default `0`)
+    Iteration of the openPMD file to be read in.
+
+Option: ``parser``
+
+* ``<laser name>.laser_real(x,y,z)`` (`string`)
+    Expression for the real part of the laser evelope in `x, y, z`.
+
+* ``<laser name>.laser_imag(x,y,z)`` (`string`)
+    Expression for the imaginary part of the laser evelope `x, y, z`.
+
+* ``lasers.lambda0`` (`float`)
+    Wavelength of the laser pulses. Currently, all pulses must have the same wavelength.
+
 Diagnostic parameters
 ---------------------
 

src/laser/Laser.H

@@ -13,6 +13,7 @@
 #include <AMReX_Vector.H>
 #include <AMReX_RealVect.H>
 #include "utils/Constants.H"
+#include "utils/Parser.H"
 
 class Laser
 {
@@ -21,6 +22,8 @@ public:
     Laser (std::string name, bool laser_from_file);
 
     std::string m_name {""};
+    /** The way to initialize a laser (from_file/gaussian/parser)*/
+    std::string m_laser_init_type = "gaussian";
     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) */
@@ -34,6 +37,13 @@ public:
     amrex::Real m_focal_distance {0.};
     /** Average position of the Gaussian laser pulse */
     amrex::RealVect m_position_mean {0., 0., 0.};
+    /** Custom profile*/
+    amrex::Parser m_parser_lr;/**< owns data for profile_real_str */
+    amrex::Parser m_parser_li;/**< owns data for profile_imag_str */
+    /** stores the output function of makeFunctionWithParser for profile_real_str */
+    amrex::ParserExecutor<3> m_profile_real;
+    /** stores the output function of makeFunctionWithParser for profile_imag_str */
+    amrex::ParserExecutor<3> m_profile_imag;
 };
 
 #endif // LASER_H_

src/laser/Laser.cpp

@@ -17,19 +17,31 @@
 Laser::Laser (std::string name, bool laser_from_file)
 {
     m_name = name;
-    if (laser_from_file) return;
     amrex::ParmParse pp(m_name);
-    queryWithParser(pp, "a0", m_a0);
-    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);
-    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,
+    queryWithParser(pp, "init_type", m_laser_init_type);
+    if (m_laser_init_type == "from_file") return;
+    else if (m_laser_init_type == "gaussian"){
+        queryWithParser(pp, "a0", m_a0);
+        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);
+        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 the laser");
-    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);
+        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);
+        return;
+    }
+    else if (m_laser_init_type == "parser"){
+        std::string profile_real_str = "";
+        std::string profile_imag_str = "";
+        getWithParser(pp, "laser_real(x,y,z)", profile_real_str);
+        getWithParser(pp, "laser_imag(x,y,z)", profile_imag_str);
+        m_profile_real = makeFunctionWithParser<3>( profile_real_str, m_parser_lr, {"x", "y", "z"});
+        m_profile_imag = makeFunctionWithParser<3>( profile_imag_str, m_parser_li, {"x", "y", "z"});
+        return;
+    }
 }

src/laser/MultiLaser.cpp

@@ -21,7 +21,6 @@
 #ifdef HIPACE_USE_OPENPMD
 #   include <openPMD/auxiliary/Filesystem.hpp>
 #endif
-
 #include <AMReX_GpuComplex.H>
 
 void
@@ -1139,47 +1138,71 @@ MultiLaser::InitLaserSlice (const int islice, const int comp)
         const amrex::Box& bx = mfi.tilebox();
         amrex::Array4<amrex::Real> const & arr = m_slices.array(mfi);
         // Initialize a Gaussian laser envelope on slice islice
-        for (int ilaser=0; ilaser<m_nlasers; ilaser++) {
-            const auto& laser = m_all_lasers[ilaser];
-            const amrex::Real a0 = laser.m_a0;
-            const amrex::Real 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 zfoc = laser.m_focal_distance;
-            amrex::ParallelFor(
+        //check point
+        for (int ilaser=0; ilaser < m_nlasers; ilaser++) {
+            auto& laser = m_all_lasers[ilaser];
+            if (laser.m_laser_init_type == "parser") {
+                auto profile_real = laser.m_profile_real;
+                auto profile_imag = laser.m_profile_imag;
+                amrex::ParallelFor(
+                bx,
+                [=] AMREX_GPU_DEVICE(int i, int j, int k)
+                {
+                    const amrex::Real x = i * dx_arr[0] + poff_x;
+                    const amrex::Real y = j * dx_arr[1] + poff_y;
+                    const amrex::Real z = islice * dx_arr[2] + poff_z;
+                    if (ilaser == 0) {
+                        arr(i, j, k, comp ) = 0._rt;
+                        arr(i, j, k, comp + 1 ) = 0._rt;
+                    }
+                    arr(i, j, k, comp ) += profile_real(x,y,z);
+                    arr(i, j, k, comp + 1 ) += profile_imag(x,y,z);
+                }
+                );
+            }
+            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 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;
+                amrex::ParallelFor(
                 bx,
                 [=] AMREX_GPU_DEVICE(int i, int j, int k)
                 {
                     const amrex::Real x = i * dx_arr[0] + poff_x - x0;
                     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 + 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;
                     // 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 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 prefactor = a0/diffract_factor;
-                    Complex time_exponent = zp*zp/(L0*L0);
+                    Complex prefactor = a0 / diffract_factor;
+                    Complex time_exponent = zp * zp / ( L0 * L0 );
                     Complex stcfactor = prefactor * amrex::exp( - time_exponent );
-                    Complex exp_argument = - ( x*x + yp*yp ) * inv_complex_waist_2;
+                    Complex exp_argument = - ( x * x + yp * yp ) * inv_complex_waist_2;
                     Complex envelope = stcfactor * amrex::exp( exp_argument ) * \
-                        amrex::exp(I * yp * k0 * propagation_angle_yz + cep);
+                       amrex::exp(I * yp * k0 * propagation_angle_yz + cep );
                     arr(i, j, k, comp ) += envelope.real();
                     arr(i, j, k, comp + 1 ) += envelope.imag();
-                }
+                    }
                 );
+            }
         }
     }
 }