Add WarpX example for FEL simulation

Docs/source/refs.bib

@@ -444,3 +444,38 @@ @article{Vranic2015
   issn = {0010-4655},
   doi = {https://doi.org/10.1016/j.cpc.2015.01.020},
 }
+
+@misc{Fallahi2020,
+      title={MITHRA 2.0: A Full-Wave Simulation Tool for Free Electron Lasers},
+      author={Arya Fallahi},
+      year={2020},
+      eprint={2009.13645},
+      archivePrefix={arXiv},
+      primaryClass={physics.acc-ph},
+      url={https://arxiv.org/abs/2009.13645},
+}
+
+@article{VayFELA2009,
+  title = {FULL ELECTROMAGNETIC SIMULATION OF FREE-ELECTRON LASER AMPLIFIER PHYSICS VIA THE LORENTZ-BOOSTED FRAME APPROACH},
+  author = {Fawley, William M and Vay, Jean-Luc},
+  abstractNote = {Numerical simulation of some systems containing charged particles with highly relativistic directed motion can by speeded up by orders of magnitude by choice of the proper Lorentz-boosted frame[1]. A particularly good example is that of short wavelength free-electron lasers (FELs) in which a high energy electron beam interacts with a static magnetic undulator. In the optimal boost frame with Lorentz factor gamma_F , the red-shifted FEL radiation and blue shifted undulator have identical wavelengths and the number of required time-steps (presuming the Courant condition applies) decreases by a factor of 2(gamma_F)**2 for fully electromagnetic simulation. We have adapted the WARP code [2]to apply this method to several FEL problems involving coherent spontaneous emission (CSE) from pre-bunched ebeams, including that in a biharmonic undulator.},
+  url = {https://www.osti.gov/biblio/964405},
+  place = {United States},
+  year = {2009},
+  month = {4},
+}
+
+@article{VayFELB2009,
+  author = {Fawley, W. M. and Vay, J.‐L.},
+  title = "{Use of the Lorentz‐Boosted Frame Transformation to Simulate Free‐Electron Laser Amplifier Physics}",
+  journal = {AIP Conference Proceedings},
+  volume = {1086},
+  number = {1},
+  pages = {346-350},
+  year = {2009},
+  month = {01},
+  abstract = "{Recently [1] it has been pointed out that numerical simulation of some systems containing charged particles with highly relativistic directed motion can by speeded up by orders of magnitude by choice of the proper Lorentz boosted frame. A particularly good example is that of short wavelength free‐electron lasers (FELs) in which a high energy (E0⩾250 MeV) electron beam interacts with a static magnetic undulator. In the optimal boost frame with Lorentz factor γF, the red‐shifted FEL radiation and blue shifted undulator have identical wavelengths and the number of required time‐steps (presuming the Courant condition applies) decreases by a factor of γF2 for fully electromagnetic simulation.We have adapted the WARP code [2] to apply this method to several FEL problems including coherent spontaneous emission (CSE) from pre‐bunched e‐beams, and strong exponential gain in a single pass amplifier configuration. We discuss our results and compare with those from the “standard” FEL simulation approach which adopts the eikonal approximation for propagation of the radiation field.}",
+  issn = {0094-243X},
+  doi = {10.1063/1.3080930},
+  url = {https://doi.org/10.1063/1.3080930},
+}

Docs/source/usage/examples.rst

@@ -44,7 +44,7 @@ Particle Accelerator & Beam Physics
 
    examples/gaussian_beam/README.rst
    examples/beam_beam_collision/README.rst
-
+   examples/free_electron_laser/README.rst
 
 High Energy Astrophysical Plasma Physics
 ----------------------------------------

Docs/source/usage/examples/free_electron_laser

@@ -0,0 +1 @@
+../../../../Examples/Physics_applications/free_electron_laser

Docs/source/usage/parameters.rst

@@ -257,7 +257,7 @@ Overall simulation parameters
         ``warpx.self_fields_absolute_tolerance``).
 
     * ``fft``: Poisson's equation is solved using an Integrated Green Function method (which requires FFT calculations).
-        See these references for more details :cite:t:`QiangPhysRevSTAB2006`, :cite:t:`QiangPhysRevSTAB2006err`.
+        See these references for more details :cite:t:`param-QiangPhysRevSTAB2006`, :cite:t:`param-QiangPhysRevSTAB2006err`.
         It only works in 3D and it requires the compilation flag ``-DWarpX_FFT=ON``.
         If mesh refinement is enabled, this solver only works on the coarsest level.
         On the refined patches, the Poisson equation is solved with the multigrid solver.

Examples/Physics_applications/CMakeLists.txt

@@ -3,6 +3,7 @@
 
 add_subdirectory(beam_beam_collision)
 add_subdirectory(capacitive_discharge)
+add_subdirectory(free_electron_laser)
 add_subdirectory(laser_acceleration)
 add_subdirectory(laser_ion)
 add_subdirectory(plasma_acceleration)

Examples/Physics_applications/free_electron_laser/CMakeLists.txt

@@ -0,0 +1,12 @@
+# Add tests (alphabetical order) ##############################################
+#
+
+add_warpx_test(
+    test_1d_fel  # name
+    1  # dims
+    2  # nprocs
+    inputs_test_1d_fel  # inputs
+    analysis_fel.py  # analysis
+    diags/diag_labframe  # output
+    OFF  # dependency
+)

Examples/Physics_applications/free_electron_laser/README.rst

@@ -0,0 +1,46 @@
+.. _examples-free-electron-laser:
+
+Free-electron laser
+===================
+
+This example shows how to simulate the physics of a free-electron laser (FEL) using WarpX.
+In this example, a relativistic electron beam is sent through an undulator (represented by an external,
+oscillating magnetic field). The radiation emitted by the beam grows exponentially
+as the beam travels through the undulator, due to the Free-Electron-Laser instability.
+
+The parameters of the simulation are taken from section 5.1 of :cite:t:`ex-Fallahi2020`.
+
+The simulation is performed in 1D, and uses the boosted-frame technique as described in
+:cite:t:`ex-VayFELA2009` and :cite:t:`ex-VayFELB2009` to reduce the computational cost (the Lorentz frame of the simulation is moving at the average speed of the beam in the undulator).
+Even though the simulation is run in this boosted frame, the results are reconstructed in the
+laboratory frame, using WarpX's ``BackTransformed`` diagnostic.
+
+The effect of space-charge is intentionally turned off in this example, as it may not be properly modeled in 1D.
+This is achieved by initializing two species of opposite charge (electrons and positrons) to
+represent the physical electron beam, as discussed in :cite:t:`ex-VayFELB2009`.
+
+Run
+---
+
+This example can be run with the WarpX executable using an input file: ``warpx.1d inputs_test_1d_fel``. For `MPI-parallel <https://www.mpi-forum.org>`__ runs, prefix these lines with ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
+
+.. literalinclude:: inputs_test_1d_fel
+   :language: ini
+   :caption: You can copy this file from ``Examples/Physics_applications/free_electron_laser/inputs_test_1d_fel``.
+
+Visualize
+---------
+
+The figure below shows the results of the simulation. The left panel shows the exponential growth of the radiation along the undulator (note that the vertical axis is plotted in log scale). The right panel shows a snapshot of the simulation,
+1.6 m into the undulator. Microbunching of the beam is visible in the electron density (blue). One can also see the
+emitted FEL radiation (red) slipping ahead of the beam.
+
+.. figure:: https://gist.githubusercontent.com/RemiLehe/871a1e24c69e353c5dbb4625cd636cd1/raw/7f4e3da7e0001cff6c592190fee8622580bbe37a/FEL.png
+   :alt: Results of the WarpX FEL simulation.
+   :width: 100%
+
+This figure was obtained with the script below, which can be run with ``python3 plot_sim.py``.
+
+.. literalinclude:: plot_sim.py
+   :language: ini
+   :caption: You can copy this file from ``Examples/Physics_applications/free_electron_laser/plot_sim.py``.

Examples/Physics_applications/free_electron_laser/analysis_fel.py

@@ -0,0 +1,145 @@
+#!/usr/bin/env python
+
+"""
+This script tests that the FEL is correctly modelled in the simulation.
+
+The physical parameters are the same as the ones from section 5.1
+of https://arxiv.org/pdf/2009.13645
+
+The simulation uses the boosted-frame technique as described in
+https://www.osti.gov/servlets/purl/940581
+In particular, the effect of space-charge is effectively turned off
+by initializing an electron and positron beam on top of each other,
+each having half the current of the physical beam.
+
+The script checks that the radiation wavelength and gain length
+are the expected ones. The check is performed both in the
+lab-frame diagnostics and boosted-frame diagnostics.
+"""
+
+import os
+import sys
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.constants import c, e, m_e
+
+sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
+from checksumAPI import evaluate_checksum
+
+# Physical parameters of the test
+gamma_bunch = 100.6
+Bu = 0.5
+lambda_u = 3e-2
+k_u = 2 * np.pi / lambda_u
+K = e * Bu / (m_e * c * k_u)  # Undulator parameter
+gamma_boost = (
+    gamma_bunch / (1 + K * K / 2) ** 0.5
+)  # Lorentz factor of the ponderomotive frame
+beta_boost = (1 - 1.0 / gamma_boost**2) ** 0.5
+
+
+# Analyze the diagnostics showing quantities in the lab frame
+filename = sys.argv[1]
+ts_lab = OpenPMDTimeSeries(filename)
+
+
+# Extract the growth of the peak electric field
+def extract_peak_E_lab(iteration):
+    """
+    Extract the position of the peak electric field
+    """
+    Ex, info = ts_lab.get_field("E", "x", iteration=iteration)
+    Ex_max = abs(Ex).max()
+    z_max = info.z[abs(Ex).argmax()]
+    return z_max, Ex_max
+
+
+# Loop through all iterations
+# Since the radiation power is proportional to the square of the peak electric field,
+# the log of the power is equal to the log of the square of the peak electric field,
+# up to an additive constant.
+z_lab_peak, E_lab_peak = ts_lab.iterate(extract_peak_E_lab)
+log_P_peak = np.log(E_lab_peak**2)
+
+# Pick the iterations between which the growth of the log of the power is linear
+# (i.e. the growth of the power is exponential) and fit a line to extract the
+# gain length.
+i_start = 6
+i_end = 23
+# Perform linear fit
+p = np.polyfit(z_lab_peak[i_start:i_end], log_P_peak[i_start:i_end], 1)
+# Extract the gain length
+Lg = 1 / p[0]
+Lg_expected = 0.22  # Expected gain length from https://arxiv.org/pdf/2009.13645
+print(f"Gain length: {Lg}")
+assert abs(Lg - Lg_expected) / Lg_expected < 0.15
+
+# Check that the radiation wavelength is the expected one
+iteration_check = 14
+Ex, info = ts_lab.get_field("E", "x", iteration=iteration_check)
+Nz = len(info.z)
+fft_E = abs(np.fft.fft(Ex))
+lambd = 1.0 / np.fft.fftfreq(Nz, d=info.dz)
+lambda_radiation_lab = lambd[fft_E[: Nz // 2].argmax()]
+lambda_expected = lambda_u / (2 * gamma_boost**2)
+print(f"lambda_radiation_lab: {lambda_radiation_lab}")
+print(f"lambda_expected: {lambda_expected}")
+assert abs(lambda_radiation_lab - lambda_expected) / lambda_expected < 0.01
+
+# Analyze the diagnostics showing quantities in the boosted frame
+ts = OpenPMDTimeSeries("diags/diag_boostedframe")
+
+
+# Extract the growth of the peak electric field
+def extract_peak_E_boost(iteration):
+    """
+    Extract the peak electric field in a *boosted-frame* snapshot.
+    Also return the position of the peak in the lab frame.
+    """
+    Ex, info = ts.get_field("E", "x", iteration=iteration)
+    By, info = ts.get_field("B", "y", iteration=iteration)
+    E_lab = gamma_boost * (Ex + c * beta_boost * By)
+    E_lab_peak = abs(E_lab).max()
+    z_boost_peak = info.z[abs(E_lab).argmax()]
+    t_boost_peak = ts.current_t
+    z_lab_peak = gamma_boost * (z_boost_peak + beta_boost * c * t_boost_peak)
+    return z_lab_peak, E_lab_peak
+
+
+# Loop through all iterations
+z_lab_peak, E_lab_peak = ts.iterate(extract_peak_E_boost)
+log_P_peak = np.log(E_lab_peak**2)
+
+# Pick the iterations between which the growth of the log of the power is linear
+# (i.e. the growth of the power is exponential) and fit a line to extract the
+# gain length.
+i_start = 16
+i_end = 25
+# Perform linear fit
+p = np.polyfit(z_lab_peak[i_start:i_end], log_P_peak[i_start:i_end], 1)
+# Extract the gain length
+Lg = 1 / p[0]
+Lg_expected = 0.22  # Expected gain length from https://arxiv.org/pdf/2009.13645
+print(f"Gain length: {Lg}")
+assert abs(Lg - Lg_expected) / Lg_expected < 0.15
+
+# Check that the radiation wavelength is the expected one
+iteration_check = 2000
+Ex, info = ts.get_field("E", "x", iteration=iteration_check)
+By, info = ts.get_field("B", "y", iteration=iteration_check)
+E_lab = gamma_boost * (Ex + c * beta_boost * By)
+Nz = len(info.z)
+fft_E = abs(np.fft.fft(E_lab))
+lambd = 1.0 / np.fft.fftfreq(Nz, d=info.dz)
+lambda_radiation_boost = lambd[fft_E[: Nz // 2].argmax()]
+lambda_radiation_lab = lambda_radiation_boost / (2 * gamma_boost)
+lambda_expected = lambda_u / (2 * gamma_boost**2)
+assert abs(lambda_radiation_lab - lambda_expected) / lambda_expected < 0.01
+
+# compare checksums
+evaluate_checksum(
+    test_name=os.path.split(os.getcwd())[1],
+    output_file=sys.argv[1],
+    output_format="openpmd",
+)

Examples/Physics_applications/free_electron_laser/inputs_test_1d_fel

@@ -0,0 +1,90 @@
+my_constants.gamma_bunch=100.6
+my_constants.Bu = 0.5
+my_constants.lambda_u = 3e-2
+my_constants.k_u= 2*pi/lambda_u
+my_constants.K = q_e*Bu/(m_e*clight*k_u) # Undulator parameter
+
+warpx.gamma_boost = gamma_bunch/sqrt(1+K*K/2) # Lorentz factor of the ponderomotive frame
+warpx.boost_direction = z
+algo.maxwell_solver = yee
+algo.particle_shape = 2
+algo.particle_pusher = vay
+
+# geometry
+geometry.dims = 1
+geometry.prob_hi = 0
+geometry.prob_lo = -192e-6
+
+amr.max_grid_size = 1024
+amr.max_level = 0
+amr.n_cell = 1024
+
+# boundary
+boundary.field_hi = absorbing_silver_mueller
+boundary.field_lo = absorbing_silver_mueller
+boundary.particle_hi = absorbing
+boundary.particle_lo = absorbing
+
+# diagnostics
+diagnostics.diags_names = diag_labframe diag_boostedframe
+
+# Diagnostic that show quantities in the frame
+# of the simulation (boosted-frame)
+diag_boostedframe.diag_type = Full
+diag_boostedframe.format = openpmd
+diag_boostedframe.intervals = 100
+
+# Diagnostic that show quantities
+# reconstructed in the lab frame
+diag_labframe.diag_type = BackTransformed
+diag_labframe.num_snapshots_lab = 25
+diag_labframe.dz_snapshots_lab = 0.1
+diag_labframe.format = openpmd
+diag_labframe.buffer_size = 64
+
+# Run the simulation long enough for
+# all backtransformed diagnostic to be complete
+warpx.compute_max_step_from_btd = 1
+
+particles.species_names = electrons positrons
+particles.rigid_injected_species= electrons positrons
+
+electrons.charge = -q_e
+electrons.injection_style = nuniformpercell
+electrons.mass = m_e
+electrons.momentum_distribution_type = constant
+electrons.num_particles_per_cell_each_dim = 8
+electrons.profile = constant
+electrons.density = 2.7e19/2
+electrons.ux = 0.0
+electrons.uy = 0.0
+electrons.uz = gamma_bunch
+electrons.zmax = -25e-6
+electrons.zmin = -125e-6
+electrons.zinject_plane=0.0
+
+positrons.charge = q_e
+positrons.injection_style = nuniformpercell
+positrons.mass = m_e
+positrons.momentum_distribution_type = constant
+positrons.num_particles_per_cell_each_dim = 8
+positrons.profile = constant
+positrons.density = 2.7e19/2
+positrons.ux = 0.0
+positrons.uy = 0.0
+positrons.uz = gamma_bunch
+positrons.zmax = -25e-6
+positrons.zmin = -125e-6
+positrons.zinject_plane=0.0
+
+warpx.do_moving_window = 1
+warpx.moving_window_dir = z
+warpx.moving_window_v = sqrt(1-(1+K*K/2)/(gamma_bunch*gamma_bunch))
+
+# Undulator field
+particles.B_ext_particle_init_style = parse_B_ext_particle_function
+particles.Bx_external_particle_function(x,y,z,t) = 0
+particles.By_external_particle_function(x,y,z,t) = Bu*cos(k_u*z)
+particles.Bz_external_particle_function(x,y,z,t) =0.0
+
+warpx.cfl = 0.99