Enable multiple laser pulses (#823)

* Rename laser multilaser

* option to have two lasers on 1 grid

* fix input files to specify lasers

* eol

* init lasers

* fix

* safer for GPU when laser class becomes bigger, although a few more kernel launches

* try fix Juwels compile issue

* update doc for multiple laser pulses

Co-authored-by: AlexanderSinn <alexander.sinn@desy.de>

docs/source/run/parameters.rst

@@ -508,57 +508,61 @@ Laser parameters
 ----------------
 
 The laser profile is defined by :math:`a(x,y,z) = a_0 * \mathrm{exp}[-(x^2/w0_x^2 + y^2/w0_y^2 + z^2/L0^2)]`.
-The laser pulse length :math:`L0 = \tau / c_0` can be specified via the pulse duration ``laser.tau``. The model implemented is the one from [C. Benedetti et al. Plasma Phys. Control. Fusion 60.1: 014002 (2017)].
+The model implemented is the one from [C. Benedetti et al. Plasma Phys. Control. Fusion 60.1: 014002 (2017)].
+Unlike for ``beams`` and ``plasmas``, all the laser pulses are currently stored on the same array,
+which you can find in the output openPMD file as `laser_real` (for the real part of the envelope) and `laser_imag` for its imaginary part.
+Parameters starting with ``lasers.`` apply to all laser pulses, parameters starting with ``<laser name>`` apply to a single laser pulse.
 
-* ``laser.use_laser`` (`0` or `1`) optional (default `0`)
-    Whether to activate the laser envelope solver.
+* ``lasers.names`` (list of `string`)
+    The names of the laser pulses, separated by a space.
+    To run without a laser, choose the name ``no_laser``.
 
-* ``laser.a0`` (`float`) optional (default `0`)
-    Peak normalized vector potential of the laser pulse.
-
-* ``laser.position_mean`` (3 `float`) optional (default `0 0 0`)
-    The mean position of the laser in `x, y, z`.
-
-* ``laser.w0`` (2 `float`) optional (default `0 0`)
-    The laser waist in `x, y`.
-
-* ``laser.L0`` (`float`) optional (default `0`)
-    The laser pulse length in `z`. Use either the pulse length or the pulse duration.
-
-* ``laser.tau`` (`float`) optional (default `0`)
-    The laser pulse duration. The pulse length will be set to `laser.tau`:math:`/c_0`.
-    Use either the pulse length or the pulse duration.
-
-* ``laser.lambda0`` (`float`)
-    The laser pulse wavelength.
-
-* ``laser.focal_distance`` (`float`)
-    Distance at which the laser pulse if focused (in the z direction, counted from laser initial position).
+* ``lasers.lambda0`` (`float`)
+    Wavelength of the laser pulses. Currently, all pulses must have the same wavelength.
 
-* ``laser.use_phase`` (`bool`) optional (default `true`)
+* ``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.
 
-* ``laser.solver_type`` (`string`) optional (default `multigrid`)
+* ``lasers.solver_type`` (`string`) optional (default `multigrid`)
     Type of solver for the laser envelope solver, either ``fft`` or ``multigrid``.
     Currently, the approximation that the phase is evaluated on-axis only is made with both solvers.
     With the multigrid solver, we could drop this assumption.
     For now, the fft solver should be faster, more accurate and more stable, so only use the multigrid one with care.
 
-* ``laser.MG_tolerance_rel`` (`float`) optional (default `1e-4`)
+* ``lasers.MG_tolerance_rel`` (`float`) optional (default `1e-4`)
     Relative error tolerance of the multigrid solver used for the laser pulse.
 
-* ``laser.MG_tolerance_abs`` (`float`) optional (default `0.`)
+* ``lasers.MG_tolerance_abs`` (`float`) optional (default `0.`)
     Absolute error tolerance of the multigrid solver used for the laser pulse.
 
-* ``laser.MG_verbose`` (`int`) optional (default `0`)
+* ``lasers.MG_verbose`` (`int`) optional (default `0`)
     Level of verbosity of the multigrid solver used for the laser pulse.
 
-* ``laser.MG_average_rhs`` (`0` or `1`) optional (default `1`)
+* ``lasers.MG_average_rhs`` (`0` or `1`) optional (default `1`)
     Whether to use the most stable discretization for the envelope solver.
 
-* ``laser.3d_on_host`` (`0` or `1`) optional (default `0`)
+* ``lasers.3d_on_host`` (`0` or `1`) optional (default `0`)
     When running on GPU: whether the 3D array containing the laser envelope is stored in host memory (CPU, slower but large memory available) or in device memory (GPU, faster but less memory available).
 
+* ``<laser name>.a0`` (`float`) optional (default `0`)
+    Peak normalized vector potential of the laser pulse.
+
+* ``<laser name>.position_mean`` (3 `float`) optional (default `0 0 0`)
+    The mean position of the laser in `x, y, z`.
+
+* ``<laser name>.w0`` (2 `float`) optional (default `0 0`)
+    The laser waist in `x, y`.
+
+* ``<laser name>.L0`` (`float`) optional (default `0`)
+    The laser pulse length in `z`. Use either the pulse length or the pulse duration ``<laser name>.tau``.
+
+* ``<laser name>.tau`` (`float`) optional (default `0`)
+    The laser pulse duration. The pulse length will be set to `laser.tau`:math:`/c_0`.
+    Use either the pulse length or the pulse duration.
+
+* ``<laser name>.focal_distance`` (`float`)
+    Distance at which the laser pulse if focused (in the z direction, counted from laser initial position).
+
 Diagnostic parameters
 ---------------------
 

examples/beam_in_vacuum/inputs_SI

@@ -39,4 +39,6 @@ beam2.ppc = 2 2 1
 
 plasmas.names = no_plasma
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/beam_in_vacuum/inputs_normalized

@@ -30,4 +30,6 @@ beam.ppc = 2 2 1
 
 plasmas.names = no_plasma
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/beam_in_vacuum/inputs_normalized_transverse

@@ -39,4 +39,7 @@ beam2.position_mean = 0. 0. 0
 beam2.position_std = 8. 0.3 1.41
 
 plasmas.names = no_plasma
+
+lasers.names = no_laser
+
 diagnostic.diag_type = xz

examples/blowout_wake/inputs_SI

@@ -44,4 +44,6 @@ plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
 plasma.element = electron
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/blowout_wake/inputs_ionization_SI

@@ -53,4 +53,6 @@ ion.mass_Da = 1.008
 ion.initial_ion_level = 0
 ion.ionization_product = elec
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/blowout_wake/inputs_normalized

@@ -40,4 +40,6 @@ plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
 plasma.element = electron
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/gaussian_weight/inputs_SI

@@ -29,4 +29,6 @@ beam.ppc = 2 2 1
 
 plasmas.names = no_plasma
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/gaussian_weight/inputs_normalized

@@ -29,4 +29,6 @@ beam.u_std = 3. 4. 0.
 
 plasmas.names = no_plasma
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/get_started/inputs_normalized

@@ -36,6 +36,10 @@ plasma.element = electron               # type of plasma: electron, proton, or a
 plasma.density(x,y,z) = 1.              # density function in the respective unit systems
 plasma.ppc = 1 1                        # particles per cell in x,y
 
+### Laser
+################################
+lasers.names = no_laser                 # name(s) of the laser pulses. Here, no laser pulse used.
+
 ### Diagnostics
 ################################
 diagnostic.diag_type = xz               # 2D xz slice output. Options: xyz, xz, yz

examples/laser/inputs_SI

@@ -12,12 +12,13 @@ hipace.do_tiling = 0
 geometry.prob_lo     = -6.*kp_inv   -6.*kp_inv   -8.*kp_inv
 geometry.prob_hi     =  6.*kp_inv    6.*kp_inv    6.*kp_inv
 
-laser.use_laser = 1
+lasers.names = laser
+lasers.lambda0 = .8e-6
+
 laser.a0 = 1
 laser.position_mean = 0. 0. 0
-laser.w0 = 2.*kp_inv 2.*kp_inv
+laser.w0 = 2.*kp_inv
 laser.L0 = 2.*kp_inv
-laser.lambda0 = .8e-6
 laser.focal_distance = 0.001
 
 hipace.bxby_solver=explicit

examples/linear_wake/inputs_SI

@@ -38,4 +38,6 @@ plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
 plasma.element = electron
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/linear_wake/inputs_ion_motion_SI

@@ -57,4 +57,6 @@ ions.neutralize_background = false
 "ions.density(x,y,z)" = ne
 ions.ppc = 1 1
 
+lasers.names = no_laser
+
 diagnostic.diag_type = xyz

examples/linear_wake/inputs_normalized

@@ -36,3 +36,5 @@ plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
 plasma.element = electron
 diagnostic.diag_type = xyz
+
+lasers.names = no_laser

src/Hipace.H

@@ -16,7 +16,7 @@
 #include "particles/beam/BeamParticleContainer.H"
 #include "utils/AdaptiveTimeStep.H"
 #include "utils/GridCurrent.H"
-#include "laser/Laser.H"
+#include "laser/MultiLaser.H"
 #include "utils/Constants.H"
 #include "utils/Parser.H"
 #include "diagnostics/Diagnostic.H"
@@ -361,7 +361,7 @@ public:
     /** Adaptive time step instance */
     AdaptiveTimeStep m_adaptive_time_step;
     /** Laser instance (soon to be multi laser container) */
-    Laser m_laser;
+    MultiLaser m_multi_laser;
     /** GridCurrent instance */
     GridCurrent m_grid_current;
 #ifdef HIPACE_USE_OPENPMD