Skip to content

Commit 4abebc3

Browse files
RemiLehebnaraEZoni
authored and
dpgrote
committed
Add WarpX example for FEL simulation (BLAST-WarpX#5337)
This adds an example for how to run FEL simulations with the boosted-frame technique. https://warpx--5337.org.readthedocs.build/en/5337/usage/examples/free_electron_laser/README.html --------- Co-authored-by: Brian Naranjo <[email protected]> Co-authored-by: Edoardo Zoni <[email protected]>
1 parent f7168d7 commit 4abebc3

File tree

11 files changed

+417
-2
lines changed

11 files changed

+417
-2
lines changed

‎Docs/source/refs.bib‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -444,3 +444,38 @@ @article{Vranic2015
444444
issn = {0010-4655},
445445
doi = {https://doi.org/10.1016/j.cpc.2015.01.020},
446446
}
447+
448+
@misc{Fallahi2020,
449+
title={MITHRA 2.0: A Full-Wave Simulation Tool for Free Electron Lasers},
450+
author={Arya Fallahi},
451+
year={2020},
452+
eprint={2009.13645},
453+
archivePrefix={arXiv},
454+
primaryClass={physics.acc-ph},
455+
url={https://arxiv.org/abs/2009.13645},
456+
}
457+
458+
@article{VayFELA2009,
459+
title = {FULL ELECTROMAGNETIC SIMULATION OF FREE-ELECTRON LASER AMPLIFIER PHYSICS VIA THE LORENTZ-BOOSTED FRAME APPROACH},
460+
author = {Fawley, William M and Vay, Jean-Luc},
461+
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.},
462+
url = {https://www.osti.gov/biblio/964405},
463+
place = {United States},
464+
year = {2009},
465+
month = {4},
466+
}
467+
468+
@article{VayFELB2009,
469+
author = {Fawley, W. M. and Vay, J.‐L.},
470+
title = "{Use of the Lorentz‐Boosted Frame Transformation to Simulate Free‐Electron Laser Amplifier Physics}",
471+
journal = {AIP Conference Proceedings},
472+
volume = {1086},
473+
number = {1},
474+
pages = {346-350},
475+
year = {2009},
476+
month = {01},
477+
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.}",
478+
issn = {0094-243X},
479+
doi = {10.1063/1.3080930},
480+
url = {https://doi.org/10.1063/1.3080930},
481+
}

‎Docs/source/usage/examples.rst‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ Particle Accelerator & Beam Physics
4444

4545
examples/gaussian_beam/README.rst
4646
examples/beam_beam_collision/README.rst
47-
47+
examples/free_electron_laser/README.rst
4848

4949
High Energy Astrophysical Plasma Physics
5050
----------------------------------------

‎Docs/source/usage/examples/free_electron_laser‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
../../../../Examples/Physics_applications/free_electron_laser

‎Docs/source/usage/parameters.rst‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -257,7 +257,7 @@ Overall simulation parameters
257257
``warpx.self_fields_absolute_tolerance``).
258258

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

‎Examples/Physics_applications/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
add_subdirectory(beam_beam_collision)
55
add_subdirectory(capacitive_discharge)
6+
add_subdirectory(free_electron_laser)
67
add_subdirectory(laser_acceleration)
78
add_subdirectory(laser_ion)
89
add_subdirectory(plasma_acceleration)

‎Examples/Physics_applications/free_electron_laser/CMakeLists.txt‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
# Add tests (alphabetical order) ##############################################
2+
#
3+
4+
add_warpx_test(
5+
test_1d_fel # name
6+
1 # dims
7+
2 # nprocs
8+
inputs_test_1d_fel # inputs
9+
analysis_fel.py # analysis
10+
diags/diag_labframe # output
11+
OFF # dependency
12+
)

‎Examples/Physics_applications/free_electron_laser/README.rst‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
.. _examples-free-electron-laser:
2+
3+
Free-electron laser
4+
===================
5+
6+
This example shows how to simulate the physics of a free-electron laser (FEL) using WarpX.
7+
In this example, a relativistic electron beam is sent through an undulator (represented by an external,
8+
oscillating magnetic field). The radiation emitted by the beam grows exponentially
9+
as the beam travels through the undulator, due to the Free-Electron-Laser instability.
10+
11+
The parameters of the simulation are taken from section 5.1 of :cite:t:`ex-Fallahi2020`.
12+
13+
The simulation is performed in 1D, and uses the boosted-frame technique as described in
14+
: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).
15+
Even though the simulation is run in this boosted frame, the results are reconstructed in the
16+
laboratory frame, using WarpX's ``BackTransformed`` diagnostic.
17+
18+
The effect of space-charge is intentionally turned off in this example, as it may not be properly modeled in 1D.
19+
This is achieved by initializing two species of opposite charge (electrons and positrons) to
20+
represent the physical electron beam, as discussed in :cite:t:`ex-VayFELB2009`.
21+
22+
Run
23+
---
24+
25+
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.
26+
27+
.. literalinclude:: inputs_test_1d_fel
28+
:language: ini
29+
:caption: You can copy this file from ``Examples/Physics_applications/free_electron_laser/inputs_test_1d_fel``.
30+
31+
Visualize
32+
---------
33+
34+
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,
35+
1.6 m into the undulator. Microbunching of the beam is visible in the electron density (blue). One can also see the
36+
emitted FEL radiation (red) slipping ahead of the beam.
37+
38+
.. figure:: https://gist.githubusercontent.com/RemiLehe/871a1e24c69e353c5dbb4625cd636cd1/raw/7f4e3da7e0001cff6c592190fee8622580bbe37a/FEL.png
39+
:alt: Results of the WarpX FEL simulation.
40+
:width: 100%
41+
42+
This figure was obtained with the script below, which can be run with ``python3 plot_sim.py``.
43+
44+
.. literalinclude:: plot_sim.py
45+
:language: ini
46+
: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‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
#!/usr/bin/env python
2+
3+
"""
4+
This script tests that the FEL is correctly modelled in the simulation.
5+
6+
The physical parameters are the same as the ones from section 5.1
7+
of https://arxiv.org/pdf/2009.13645
8+
9+
The simulation uses the boosted-frame technique as described in
10+
https://www.osti.gov/servlets/purl/940581
11+
In particular, the effect of space-charge is effectively turned off
12+
by initializing an electron and positron beam on top of each other,
13+
each having half the current of the physical beam.
14+
15+
The script checks that the radiation wavelength and gain length
16+
are the expected ones. The check is performed both in the
17+
lab-frame diagnostics and boosted-frame diagnostics.
18+
"""
19+
20+
import os
21+
import sys
22+
23+
import numpy as np
24+
from openpmd_viewer import OpenPMDTimeSeries
25+
from scipy.constants import c, e, m_e
26+
27+
sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
28+
from checksumAPI import evaluate_checksum
29+
30+
# Physical parameters of the test
31+
gamma_bunch = 100.6
32+
Bu = 0.5
33+
lambda_u = 3e-2
34+
k_u = 2 * np.pi / lambda_u
35+
K = e * Bu / (m_e * c * k_u) # Undulator parameter
36+
gamma_boost = (
37+
gamma_bunch / (1 + K * K / 2) ** 0.5
38+
) # Lorentz factor of the ponderomotive frame
39+
beta_boost = (1 - 1.0 / gamma_boost**2) ** 0.5
40+
41+
42+
# Analyze the diagnostics showing quantities in the lab frame
43+
filename = sys.argv[1]
44+
ts_lab = OpenPMDTimeSeries(filename)
45+
46+
47+
# Extract the growth of the peak electric field
48+
def extract_peak_E_lab(iteration):
49+
"""
50+
Extract the position of the peak electric field
51+
"""
52+
Ex, info = ts_lab.get_field("E", "x", iteration=iteration)
53+
Ex_max = abs(Ex).max()
54+
z_max = info.z[abs(Ex).argmax()]
55+
return z_max, Ex_max
56+
57+
58+
# Loop through all iterations
59+
# Since the radiation power is proportional to the square of the peak electric field,
60+
# the log of the power is equal to the log of the square of the peak electric field,
61+
# up to an additive constant.
62+
z_lab_peak, E_lab_peak = ts_lab.iterate(extract_peak_E_lab)
63+
log_P_peak = np.log(E_lab_peak**2)
64+
65+
# Pick the iterations between which the growth of the log of the power is linear
66+
# (i.e. the growth of the power is exponential) and fit a line to extract the
67+
# gain length.
68+
i_start = 6
69+
i_end = 23
70+
# Perform linear fit
71+
p = np.polyfit(z_lab_peak[i_start:i_end], log_P_peak[i_start:i_end], 1)
72+
# Extract the gain length
73+
Lg = 1 / p[0]
74+
Lg_expected = 0.22 # Expected gain length from https://arxiv.org/pdf/2009.13645
75+
print(f"Gain length: {Lg}")
76+
assert abs(Lg - Lg_expected) / Lg_expected < 0.15
77+
78+
# Check that the radiation wavelength is the expected one
79+
iteration_check = 14
80+
Ex, info = ts_lab.get_field("E", "x", iteration=iteration_check)
81+
Nz = len(info.z)
82+
fft_E = abs(np.fft.fft(Ex))
83+
lambd = 1.0 / np.fft.fftfreq(Nz, d=info.dz)
84+
lambda_radiation_lab = lambd[fft_E[: Nz // 2].argmax()]
85+
lambda_expected = lambda_u / (2 * gamma_boost**2)
86+
print(f"lambda_radiation_lab: {lambda_radiation_lab}")
87+
print(f"lambda_expected: {lambda_expected}")
88+
assert abs(lambda_radiation_lab - lambda_expected) / lambda_expected < 0.01
89+
90+
# Analyze the diagnostics showing quantities in the boosted frame
91+
ts = OpenPMDTimeSeries("diags/diag_boostedframe")
92+
93+
94+
# Extract the growth of the peak electric field
95+
def extract_peak_E_boost(iteration):
96+
"""
97+
Extract the peak electric field in a *boosted-frame* snapshot.
98+
Also return the position of the peak in the lab frame.
99+
"""
100+
Ex, info = ts.get_field("E", "x", iteration=iteration)
101+
By, info = ts.get_field("B", "y", iteration=iteration)
102+
E_lab = gamma_boost * (Ex + c * beta_boost * By)
103+
E_lab_peak = abs(E_lab).max()
104+
z_boost_peak = info.z[abs(E_lab).argmax()]
105+
t_boost_peak = ts.current_t
106+
z_lab_peak = gamma_boost * (z_boost_peak + beta_boost * c * t_boost_peak)
107+
return z_lab_peak, E_lab_peak
108+
109+
110+
# Loop through all iterations
111+
z_lab_peak, E_lab_peak = ts.iterate(extract_peak_E_boost)
112+
log_P_peak = np.log(E_lab_peak**2)
113+
114+
# Pick the iterations between which the growth of the log of the power is linear
115+
# (i.e. the growth of the power is exponential) and fit a line to extract the
116+
# gain length.
117+
i_start = 16
118+
i_end = 25
119+
# Perform linear fit
120+
p = np.polyfit(z_lab_peak[i_start:i_end], log_P_peak[i_start:i_end], 1)
121+
# Extract the gain length
122+
Lg = 1 / p[0]
123+
Lg_expected = 0.22 # Expected gain length from https://arxiv.org/pdf/2009.13645
124+
print(f"Gain length: {Lg}")
125+
assert abs(Lg - Lg_expected) / Lg_expected < 0.15
126+
127+
# Check that the radiation wavelength is the expected one
128+
iteration_check = 2000
129+
Ex, info = ts.get_field("E", "x", iteration=iteration_check)
130+
By, info = ts.get_field("B", "y", iteration=iteration_check)
131+
E_lab = gamma_boost * (Ex + c * beta_boost * By)
132+
Nz = len(info.z)
133+
fft_E = abs(np.fft.fft(E_lab))
134+
lambd = 1.0 / np.fft.fftfreq(Nz, d=info.dz)
135+
lambda_radiation_boost = lambd[fft_E[: Nz // 2].argmax()]
136+
lambda_radiation_lab = lambda_radiation_boost / (2 * gamma_boost)
137+
lambda_expected = lambda_u / (2 * gamma_boost**2)
138+
assert abs(lambda_radiation_lab - lambda_expected) / lambda_expected < 0.01
139+
140+
# compare checksums
141+
evaluate_checksum(
142+
test_name=os.path.split(os.getcwd())[1],
143+
output_file=sys.argv[1],
144+
output_format="openpmd",
145+
)

‎Examples/Physics_applications/free_electron_laser/inputs_test_1d_fel‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
my_constants.gamma_bunch=100.6
2+
my_constants.Bu = 0.5
3+
my_constants.lambda_u = 3e-2
4+
my_constants.k_u= 2*pi/lambda_u
5+
my_constants.K = q_e*Bu/(m_e*clight*k_u) # Undulator parameter
6+
7+
warpx.gamma_boost = gamma_bunch/sqrt(1+K*K/2) # Lorentz factor of the ponderomotive frame
8+
warpx.boost_direction = z
9+
algo.maxwell_solver = yee
10+
algo.particle_shape = 2
11+
algo.particle_pusher = vay
12+
13+
# geometry
14+
geometry.dims = 1
15+
geometry.prob_hi = 0
16+
geometry.prob_lo = -192e-6
17+
18+
amr.max_grid_size = 1024
19+
amr.max_level = 0
20+
amr.n_cell = 1024
21+
22+
# boundary
23+
boundary.field_hi = absorbing_silver_mueller
24+
boundary.field_lo = absorbing_silver_mueller
25+
boundary.particle_hi = absorbing
26+
boundary.particle_lo = absorbing
27+
28+
# diagnostics
29+
diagnostics.diags_names = diag_labframe diag_boostedframe
30+
31+
# Diagnostic that show quantities in the frame
32+
# of the simulation (boosted-frame)
33+
diag_boostedframe.diag_type = Full
34+
diag_boostedframe.format = openpmd
35+
diag_boostedframe.intervals = 100
36+
37+
# Diagnostic that show quantities
38+
# reconstructed in the lab frame
39+
diag_labframe.diag_type = BackTransformed
40+
diag_labframe.num_snapshots_lab = 25
41+
diag_labframe.dz_snapshots_lab = 0.1
42+
diag_labframe.format = openpmd
43+
diag_labframe.buffer_size = 64
44+
45+
# Run the simulation long enough for
46+
# all backtransformed diagnostic to be complete
47+
warpx.compute_max_step_from_btd = 1
48+
49+
particles.species_names = electrons positrons
50+
particles.rigid_injected_species= electrons positrons
51+
52+
electrons.charge = -q_e
53+
electrons.injection_style = nuniformpercell
54+
electrons.mass = m_e
55+
electrons.momentum_distribution_type = constant
56+
electrons.num_particles_per_cell_each_dim = 8
57+
electrons.profile = constant
58+
electrons.density = 2.7e19/2
59+
electrons.ux = 0.0
60+
electrons.uy = 0.0
61+
electrons.uz = gamma_bunch
62+
electrons.zmax = -25e-6
63+
electrons.zmin = -125e-6
64+
electrons.zinject_plane=0.0
65+
electrons.rigid_advance=0
66+
67+
positrons.charge = q_e
68+
positrons.injection_style = nuniformpercell
69+
positrons.mass = m_e
70+
positrons.momentum_distribution_type = constant
71+
positrons.num_particles_per_cell_each_dim = 8
72+
positrons.profile = constant
73+
positrons.density = 2.7e19/2
74+
positrons.ux = 0.0
75+
positrons.uy = 0.0
76+
positrons.uz = gamma_bunch
77+
positrons.zmax = -25e-6
78+
positrons.zmin = -125e-6
79+
positrons.zinject_plane=0.0
80+
positrons.rigid_advance=0
81+
82+
warpx.do_moving_window = 1
83+
warpx.moving_window_dir = z
84+
warpx.moving_window_v = sqrt(1-(1+K*K/2)/(gamma_bunch*gamma_bunch))
85+
86+
# Undulator field
87+
particles.B_ext_particle_init_style = parse_B_ext_particle_function
88+
particles.Bx_external_particle_function(x,y,z,t) = 0
89+
particles.By_external_particle_function(x,y,z,t) = if( z>0, Bu*cos(k_u*z), 0 )
90+
particles.Bz_external_particle_function(x,y,z,t) =0.0
91+
92+
warpx.cfl = 0.99

0 commit comments

Comments
 (0)