Content from Introduction to WarpX
Last updated on 2025-11-11 | Edit this page
Estimated time: 10 minutes
Overview
Questions
- 🤌 What is WarpX?
- 🤔 What is a PIC code?
- 🧐 What can I use WarpX for?
Objectives
- 💡 Understand the basics of PIC codes
- 🧑💻 Learn about the features of WarpX
- 🎯 Figure out if WarpX can be useful for you!
Overview of PIC
WarpX is a general purpose
open-source high-performance
Particle-In-Cell (PIC) code.
The PIC method is a very popular approach to simulate the
dynamics of physical systems governed by relativistic
electrodynamics. Plasmas ⭐ and beams ☄️,
which are often made of charged particles that may travel almost at the
speed of light, fall into this category. Additionally, other physical
effects can be integrated into the PIC algorithm, such as different
quantum 🎲 processes.
Here is a picture that condenses the core idea: there are particles and fields. The particles are approximated as macroparticles (usually each representative of many real particles), while the fields are approximated on a grid in space. The particles and the fields are updated, self-consistently, in a temporal loop.
Here is a more informative image that explains the core algorithmic steps. As the particles travel in space, they generate currents, which in turn generate an electromagnetic field. The electromagnetic field then push the particles via the Lorentz force. Therefore, the current density \(\textbf{J}\) and the force \(\textbf{F}_L\) are the quantitieis that connect the particles and the fields, or in PIC lingo, the macroparticles and the grid. Hence, the idea is the following. First, define a well-posed initial condition, and then iterate the following: * Interpolate the fields from the grid to the particles’ positions and compute the Lorentz force that acts on each macroparticle, * Advance the position and momenta of the macroparticles using Newton equations, * Project while cumulating the contribution to the current density of each macroparticles, * Solve Maxwell’s equations.
In some cases, one can choose to solve Poisson’s equation instead of Maxwell’s. In that case, the current \(\textbf{J}\) calculation is replaced with the charge density \(\rho\) calculation. Once \(\rho\) is known, the electrostic potential is computed to then find the electric field.
If you want to know more about PIC, here are a few references:
- The two bibles on PIC 📚
- An old review written by one of the pioneers: John M. Dawson,
Particle simulation of plasmas, Rev. Mod. Phys. 55, 403
- Browse our docs for many, many more references about advanced algorithms and methods.
Features and applications of WarpX
WarpX is developed and used by a wide range of researchers working in different fields, from beam physics to nuclear fusion and a lot more. To learn about what WarpX can be used for, check out our examples and the scientific publications that acknowledged WarpX.
Some cool features of WarpX:
📖 Open-source!
✈️ Runs on GPUs: NVIDIA, AMD, and Intel
🚀 Runs on multiple GPUs or CPUs, on systems ranging from laptops to supercomputers
🤓 Many many advanced algorithms and methods: mesh-refinement, embedded boundaries, electrostatic/electromagnetic/pseudospectral solvers, etc.
💾 Standards: openPMD for input/output data, PICMI for inputs
🤸 Active development and mainteinance: check out our GitHub repo
🗺️ International, cross-disciplinary community: plasma physics, fusion devices, laser-plasma interactions, beam physics, plasma-based acceleration, astrophysics, others?
We will add here more details soon!
Content from Install
Last updated on 2026-02-03 | Edit this page
Estimated time: 10 minutes
Overview
Questions
- 🔧 How can I install and run WarpX?
- 🕵️ How can I analyze the simulation results?
Objectives
- 💻 Install WarpX on your local machine either with Conda or from source
- 👀 Install the visualizations tools in
PythonandParaview
Basic dependencies
Just a heads-up before we dive deeper.
Via Conda-Forge
First, you need a Conda installation and we will
assume that you indeed have one.
If not, follow the instruction
at this link.
You can install Conda on most operative systems: Windows, macOS,
and Linux.
We will also assume you have some familiarity with the
terminal. Once you have Conda on your system,
WarpX is available as a package via Conda-Forge.
The installation is a one-liner 😌!
Ok, maybe two lines if you want to keep your system clean by creating a new environment.
The first warpx (after -n) is the name of
the conda environment; the second warpx is the name of the
package to be installed from the conda-forge channel.
Now you should have 4 different WarpX binaries in your
PATH called warpx.1d, warpx.2d,
warpx.3d, warpx.rz.
Each binary for a different dimensionality.
To check this, run:
If you get 3 different paths that look something like:
then you got this 🙌! You can also import pywarpx in
Python.
Conda’s WarpX is serial! To get a parallel WarpX version, install it from source.
From source
🎯 WarpX is easy to install via Conda:
conda -c conda-forge warpx
🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.
Content from A Beam-Beam Collision
Last updated on 2025-11-11 | Edit this page
Estimated time: 30 minutes
Overview
Questions
How to use WarpX for beam-beam simulations of colliders 💥?
Objectives
Learn how to setup, run, and visualize your own beam-beam ☄️☄️ simulation
Setup
In this example we will simulate a bunch of electrons colliding against a bunch of positrons. We have selected the parameters of the C\(^3\) linear collider.
Whenever you need to prepare an input file, this is where you want to go.
OUTPUT
#################################
########## MY CONSTANTS #########
#################################
my_constants.mc2 = m_e*clight*clight
my_constants.nano = 1.0e-9
my_constants.micro = 1.e-6
my_constants.milli = 1.e-3
my_constants.GeV = q_e*1.e9
# BEAMS
my_constants.beam_energy = 125*GeV
my_constants.beam_gamma = beam_energy/mc2
my_constants.beam_npart = 6.24e9
my_constants.nmacropart = 1e5
my_constants.beam_charge = q_e * beam_npart
#my_constants.sigmax = 210.0*nano
#my_constants.sigmay = 3.1*nano
my_constants.sigmaz = 100.*micro
my_constants.beam_uth = 0.3/100*beam_gamma
my_constants.mux = 0.*sigmax
my_constants.muy = 0.*sigmay
my_constants.muz = 4*sigmaz
my_constants.emitx = 900*nano
my_constants.emity = 20*nano
my_constants.dux = emitx / sigmax
my_constants.duy = emity / sigmay
my_constants.betax = 12*milli
my_constants.betay = 0.12*milli
my_constants.sigmax = sqrt( emitx * betax / beam_gamma )
my_constants.sigmay = sqrt( emity * betay / beam_gamma )
# BOX
my_constants.Lx = 20*sigmax
my_constants.Ly = 20*sigmay
my_constants.Lz = 16*sigmaz
my_constants.nx = 128
my_constants.ny = 128
my_constants.nz = 64
my_constants.dx = Lx/nx
my_constants.dy = Ly/ny
my_constants.dz = Lz/nz
# TIME
my_constants.T = 0.5*Lz/clight
my_constants.dt = T / nz
my_constants.nt = floor(T/dt)
# LUMI DIAG
my_constants.bin_num_1d = 2048
my_constants.bin_center_min_1d = 0.
my_constants.bin_center_max_1d = 2*beam_energy/q_e
my_constants.bin_size_1d = (bin_center_max_1d - bin_center_min_1d) / bin_num_1d
my_constants.bin_edge_min_1d = bin_center_min_1d - 0.5 * bin_size_1d
my_constants.bin_edge_max_1d = bin_center_max_1d + 0.5 * bin_size_1d
my_constants.bin_num_1d_eff = bin_num_1d + 1
#################################
####### GENERAL PARAMETERS ######
#################################
stop_time = T
amr.n_cell = nx ny nz
amr.max_level = 0
geometry.dims = 3
geometry.prob_lo = -0.5*Lx -0.5*Ly -0.5*Lz
geometry.prob_hi = 0.5*Lx 0.5*Ly 0.5*Lz
#################################
######## BOUNDARY CONDITION #####
#################################
boundary.field_lo = open open open
boundary.field_hi = open open open
boundary.particle_lo = Absorbing Absorbing Absorbing
boundary.particle_hi = Absorbing Absorbing Absorbing
#################################
############ NUMERICS ###########
#################################
warpx.do_electrostatic = relativistic
warpx.const_dt = dt
warpx.grid_type = collocated
algo.particle_shape = 3
algo.particle_pusher = vay
warpx.poisson_solver = fft
warpx.use_2d_slices_fft_solver = 1
#################################
########### PARTICLES ###########
#################################
particles.species_names = ele1 pos2 pho1 pho2 ele2 pos1
particles.photon_species = pho1 pho2
ele1.species_type = electron
ele1.injection_style = gaussian_beam
ele1.x_rms = sigmax
ele1.y_rms = sigmay
ele1.z_rms = sigmaz
ele1.x_m = - mux
ele1.y_m = - muy
ele1.z_m = - muz
ele1.npart = nmacropart
ele1.q_tot = -beam_charge
ele1.z_cut = 4
ele1.focal_distance = muz
ele1.momentum_distribution_type = gaussian
ele1.uz_m = beam_gamma
ele1.uy_m = 0.0
ele1.ux_m = 0.0
ele1.ux_th = dux
ele1.uy_th = duy
ele1.uz_th = beam_uth
ele1.initialize_self_fields = 1
ele1.do_qed_quantum_sync = 1
ele1.qed_quantum_sync_phot_product_species = pho1
ele1.do_classical_radiation_reaction = 0
pos2.species_type = positron
pos2.injection_style = gaussian_beam
pos2.x_rms = sigmax
pos2.y_rms = sigmay
pos2.z_rms = sigmaz
pos2.x_m = mux
pos2.y_m = muy
pos2.z_m = muz
pos2.npart = nmacropart
pos2.q_tot = beam_charge
pos2.z_cut = 4
pos2.focal_distance = muz
pos2.momentum_distribution_type = gaussian
pos2.uz_m = -beam_gamma
pos2.uy_m = 0.0
pos2.ux_m = 0.0
pos2.ux_th = dux
pos2.uy_th = duy
pos2.uz_th = beam_uth
pos2.initialize_self_fields = 1
pos2.do_qed_quantum_sync = 1
pos2.qed_quantum_sync_phot_product_species = pho2
pos2.do_classical_radiation_reaction = 0
pho1.species_type = photon
pho1.injection_style = none
pho1.do_qed_breit_wheeler = 1
pho1.qed_breit_wheeler_ele_product_species = ele1
pho1.qed_breit_wheeler_pos_product_species = pos1
pho2.species_type = photon
pho2.injection_style = none
pho2.do_qed_breit_wheeler = 1
pho2.qed_breit_wheeler_ele_product_species = ele2
pho2.qed_breit_wheeler_pos_product_species = pos2
ele2.species_type = electron
ele2.injection_style = none
ele2.do_qed_quantum_sync = 1
ele2.qed_quantum_sync_phot_product_species = pho2
ele2.do_classical_radiation_reaction = 0
pos1.species_type = positron
pos1.injection_style = none
pos1.do_qed_quantum_sync = 1
pos1.qed_quantum_sync_phot_product_species = pho1
pos1.do_classical_radiation_reaction = 0
#################################
############# QED ###############
#################################
qed_qs.photon_creation_energy_threshold = 0.
qed_qs.lookup_table_mode = builtin
qed_qs.chi_min = 1.e-7
qed_bw.lookup_table_mode = builtin
qed_bw.chi_min = 1.e-2
warpx.do_qed_schwinger = 0.
#################################
######### DIAGNOSTICS ###########
#################################
# FULL
diagnostics.diags_names = particles_out particles_in
# particles that exit the simulation domain at any given time
particles_out.dump_last_timestep = 1
particles_out.diag_type = BoundaryScraping
particles_out.format = openpmd
particles_out.openpmd_backend = h5
particles_out.intervals = -1
ele1.save_particles_at_xlo = 1
ele1.save_particles_at_ylo = 1
ele1.save_particles_at_zlo = 1
ele1.save_particles_at_xhi = 1
ele1.save_particles_at_yhi = 1
ele1.save_particles_at_zhi = 1
pos2.save_particles_at_xlo = 1
pos2.save_particles_at_ylo = 1
pos2.save_particles_at_zlo = 1
pos2.save_particles_at_xhi = 1
pos2.save_particles_at_yhi = 1
pos2.save_particles_at_zhi = 1
pho1.save_particles_at_xlo = 1
pho1.save_particles_at_ylo = 1
pho1.save_particles_at_zlo = 1
pho1.save_particles_at_xhi = 1
pho1.save_particles_at_yhi = 1
pho1.save_particles_at_zhi = 1
pho2.save_particles_at_xlo = 1
pho2.save_particles_at_ylo = 1
pho2.save_particles_at_zlo = 1
pho2.save_particles_at_xhi = 1
pho2.save_particles_at_yhi = 1
pho2.save_particles_at_zhi = 1
ele2.save_particles_at_xlo = 1
ele2.save_particles_at_ylo = 1
ele2.save_particles_at_zlo = 1
ele2.save_particles_at_xhi = 1
ele2.save_particles_at_yhi = 1
ele2.save_particles_at_zhi = 1
pos1.save_particles_at_xlo = 1
pos1.save_particles_at_ylo = 1
pos1.save_particles_at_zlo = 1
pos1.save_particles_at_xhi = 1
pos1.save_particles_at_yhi = 1
pos1.save_particles_at_zhi = 1
# particles inside the simulation domain
particles_in.intervals = 1
particles_in.diag_type = Full
particles_in.species = ele1 pos2
particles_in.fields_to_plot = none
particles_in.format = openpmd
particles_in.openpmd_backend = h5
particles_in.dump_last_timestep = 1
# REDUCED
warpx.reduced_diags_names = DiffLumi_ele1_pos2 DiffLumi_ele1_ele2 DiffLumi_pos1_pos2 DiffLumi_pos1_ele2 DiffLumi_pos1_pho2 DiffLumi_ele1_pho2 DiffLumi_pho1_pos2 DiffLumi_pho1_ele2 DiffLumi_pho1_pho2 CollRel_ele1_pos2
CollRel_ele1_pos2.type = ColliderRelevant
CollRel_ele1_pos2.species = ele1 pos2
CollRel_ele1_pos2.intervals = 1
DiffLumi_ele1_pos2.type = DifferentialLuminosity
DiffLumi_ele1_pos2.intervals = nt
DiffLumi_ele1_pos2.species = ele1 pos2
DiffLumi_ele1_pos2.bin_number = bin_num_1d_eff
DiffLumi_ele1_pos2.bin_max = bin_edge_max_1d
DiffLumi_ele1_pos2.bin_min = bin_edge_min_1d
DiffLumi_ele1_ele2.type = DifferentialLuminosity
DiffLumi_ele1_ele2.intervals = nt
DiffLumi_ele1_ele2.species = ele1 ele2
DiffLumi_ele1_ele2.bin_number = bin_num_1d_eff
DiffLumi_ele1_ele2.bin_max = bin_edge_max_1d
DiffLumi_ele1_ele2.bin_min = bin_edge_min_1d
DiffLumi_pos1_pos2.type = DifferentialLuminosity
DiffLumi_pos1_pos2.intervals = nt
DiffLumi_pos1_pos2.species = pos1 pos2
DiffLumi_pos1_pos2.bin_number = bin_num_1d_eff
DiffLumi_pos1_pos2.bin_max = bin_edge_max_1d
DiffLumi_pos1_pos2.bin_min = bin_edge_min_1d
DiffLumi_pos1_ele2.type = DifferentialLuminosity
DiffLumi_pos1_ele2.intervals = nt
DiffLumi_pos1_ele2.species = pos1 ele2
DiffLumi_pos1_ele2.bin_number = bin_num_1d_eff
DiffLumi_pos1_ele2.bin_max = bin_edge_max_1d
DiffLumi_pos1_ele2.bin_min = bin_edge_min_1d
DiffLumi_ele1_pho2.type = DifferentialLuminosity
DiffLumi_ele1_pho2.intervals = nt
DiffLumi_ele1_pho2.species = ele1 pho2
DiffLumi_ele1_pho2.bin_number = bin_num_1d_eff
DiffLumi_ele1_pho2.bin_max = bin_edge_max_1d
DiffLumi_ele1_pho2.bin_min = bin_edge_min_1d
DiffLumi_pos1_pho2.type = DifferentialLuminosity
DiffLumi_pos1_pho2.intervals = nt
DiffLumi_pos1_pho2.species = pos1 pho2
DiffLumi_pos1_pho2.bin_number = bin_num_1d_eff
DiffLumi_pos1_pho2.bin_max = bin_edge_max_1d
DiffLumi_pos1_pho2.bin_min = bin_edge_min_1d
DiffLumi_pho1_pos2.type = DifferentialLuminosity
DiffLumi_pho1_pos2.intervals = nt
DiffLumi_pho1_pos2.species = pho1 pos2
DiffLumi_pho1_pos2.bin_number = bin_num_1d_eff
DiffLumi_pho1_pos2.bin_max = bin_edge_max_1d
DiffLumi_pho1_pos2.bin_min = bin_edge_min_1d
DiffLumi_pho1_ele2.type = DifferentialLuminosity
DiffLumi_pho1_ele2.intervals = nt
DiffLumi_pho1_ele2.species = pho1 ele2
DiffLumi_pho1_ele2.bin_number = bin_num_1d_eff
DiffLumi_pho1_ele2.bin_max = bin_edge_max_1d
DiffLumi_pho1_ele2.bin_min = bin_edge_min_1d
DiffLumi_pho1_pho2.type = DifferentialLuminosity
DiffLumi_pho1_pho2.intervals = nt
DiffLumi_pho1_pho2.species = pho1 pho2
DiffLumi_pho1_pho2.bin_number = bin_num_1d_eff
DiffLumi_pho1_pho2.bin_max = bin_edge_max_1d
DiffLumi_pho1_pho2.bin_min = bin_edge_min_1dSome notable details:
Poisson solver: because the beams are ultra-relativistic (125 GeV electrons and positrons) and flat, we use an FFT-based electrostatic solver. Specifically, it solves many 2D Poisson equations in the \((x,y)\) plane for each \(z\). The full 3D version of this solver is also available with
warpx.use_2d_slices_fft_solver = 0.Resolution: the number of grid cells is reduced to fit in a laptop. For production simulations, make sure you increase the resolution.
Timestep: since the beams travel at the speed of light along \(z\) and the simulation frame is the center of mass frame, it makes sense to choose $dt = dz / ( 2 c ) $. However this is not strictly necessary. Sometimes it can be useful to save resources by choosing a larger timestep. Just make sure you resolve ‘’well-enough’’ the shortest timescale that you’re interested in.
QED lookup tables: here we use the built-in ones for convenience. For production runs, make sure to use tables with enough points and set the ranges of the \(\chi\) parameter to what you need.
-
Diagnostics:
- the trajectories of the beam particles. This diagnostic can easily take up too much memory. For simulations with many macroparticles, consider using a field diagnostic.
- all the particles that exit the domain
(
BoundaryScraping) - the differential luminosity of every pair of left-ward and right-ward moving species
Run
First things first. Create a new folder where you copy the input file. The simulation in small enough that you should be able to run it in serial with the Conda installation of WarpX.
If you want to make the simulation faster and/or bigger, then you should run with the parallel version of WarpX. The optimal setup to run the simulation depends on your hardware. This is an example that should work on many common laptops, even though it might not be ideal.
On my laptop’s CPU (12th Gen Intel® Core™ i9-12900H × 20) the serial simulation took ~195 seconds, while the parallel one ~44 seconds on 4 cores!
Visualize
With Python 🐍
Now that we have the results, we can analyze them using Python.
We will use the openPMD-viewer library
to grab the data that the simulation produced in openPMD
format. Here you can find a
few great tutorials on how to use the viewer. If you feel nerdy
and/or you need to deal with the data in parallel workflows, you can use
the openPMD-api.
As an example for the beam-beam simulation, we have developed simple Jupyter notebook where we retrieve the beams’ particles positions and project them on the \((z,x)\) and \((z,y)\) planes.
You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.
With Paraview
Coming soon!
💅 There are several details one needs to take care when setting up a beam-beam simulation
🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.
📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.
Content from A Magnetic Mirror
Last updated on 2025-11-11 | Edit this page
Estimated time: 12 minutes
Overview
Questions
🪞 How to simulate the dynamics of charged particles in an external field?
Objectives
🏃 Run and 👀 visualize some protons in a magnetic mirror!
Setup
In this example we will simulate a bunch of protons inside a
magnetic mirror machine. The protons are initialized with
random positions and velocities. The magnetic field is loaded from a
.h5 file. Make sure to download the input file.
Whenever you need to prepare an input file, this is where you want to go. By the way, analytics tell us that this is the most popular page of the documentation 👠!
OUTPUT
##########################
# USER-DEFINED CONSTANTS #
##########################
my_constants.Lx = 2 # [m]
my_constants.Ly = 2 # [m]
my_constants.Lz = 5 # [m]
my_constants.dt = 4.4e-7 # [s]
my_constants.Np = 1000
############
# NUMERICS #
############
geometry.dims = 3
geometry.prob_hi = 0.5*Lx 0.5*Ly Lz
geometry.prob_lo = -0.5*Lx -0.5*Ly 0
amr.n_cell = 40 40 40
max_step = 500
warpx.const_dt = dt
##############
# ALGORITHMS #
##############
algo.particle_shape = 1
amr.max_level = 0
warpx.do_electrostatic = labframe
warpx.grid_type = collocated
warpx.serialize_initial_conditions = 0
warpx.use_filter = 0
##############
# BOUNDARIES #
##############
boundary.field_hi = pec pec pec
boundary.field_lo = pec pec pec
boundary.particle_hi = absorbing absorbing absorbing
boundary.particle_lo = absorbing absorbing absorbing
#############
# PARTICLES #
#############
particles.species_names = protons
protons.charge = q_e
protons.mass = m_p
protons.do_not_deposit = 1 # test particles
protons.initialize_self_fields = 0
protons.injection_style = gaussian_beam
protons.x_rms = 0.1*Lx
protons.y_rms = 0.1*Ly
protons.z_rms = 0.1*Lz
protons.x_m = 0.
protons.y_m = 0.
protons.z_m = 0.5*Lz
protons.npart = Np
protons.q_tot = q_e*Np
protons.momentum_distribution_type = uniform
protons.ux_min = -9.5e-05
protons.uy_min = -9.5e-05
protons.uz_min = -0.000134
protons.ux_max = 9.5e-05
protons.uy_max = 9.5e-05
protons.uz_max = 0.000134
##########
# FIELDS #
##########
# field here is applied on directly the particles!
particles.B_ext_particle_init_style = read_from_file
particles.read_fields_from_path = example-femm-3d.h5
###############
# DIAGNOSTICS #
###############
diagnostics.diags_names = diag1
diag1.diag_type = Full
diag1.fields_to_plot = Bx By Bz
diag1.format = openpmd
diag1.intervals = 1
diag1.proton.variables = ux uy uz w x y z
diag1.species = protons
diag1.write_species = 1A few notable details:
The protons are test particles because of the parameter
protons.do_not_deposit=1. This means that the protons do not deposit their current density, therefore they do not contribute to the fields.The magnetic field is applied directly to the particles with the
particles.B_ext_particle_init_styleflag, so in principle the grid is not used at all. For technical reasons, we must define a grid nonetheless.
Now that we have an idea of what the input files looks like, let’s
set up our environment. Activate the warpx environment if
you need to. Create a new directory with your own copy of the input
file. Also, don’t forget to download the field file and place
it in the directory where you will run the input.
Run
Let’s run the code
How would you do it? 🤷
You should see a standard output flashing out a lot of info.
At the end, you should find in your folder:
- a subfolder called
diags: here is where the code stored the diagnostics
- a file called
warpx_used_inputs: this is a summary of the inputs that were used to run the simulation
If that’s the case, yey! 💯
If the run went wrong, you may find a Backtrace.0.0 file
which can be useful for debugging purposes. Let me know if the code
fails in any way!
Here we have loaded the field of hte magnetic bottle from a file. You can also you can define an external field analytically.
Visualize
With Python 🐍
Now that we have the results, we can analyze them using Python.
We will use the openPMD-viewer library
to grab the data that the simulation produced in openPMD
format. Here you can find a
few great tutorials on how to use the viewer. If you feel nerdy
and/or you need to deal with the data in parallel workflows, you can use
the openPMD-api.
As an example for the magnetic bottle simulation, we have developed simple Jupyter notebook where we retrieve the magnetic field and the particle attributes at the end of the simulation. With a little bit more work, we also plot the trajectories of the particles.
You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.
With Paraview
Now it’s time to produce some pretty cool images and videos! 😎 If
you don’t have it, you can download Paraview here. In
the diags/diag1 directory you should find a file named
paraview.pmd: Paraview can read .pmd files.
Just open Paraview and from there open the .pmd file. You
should see Meshes and Particles in your
pipeline browser (usually on the left). We can zhuzh up the pipeline so
that we can visualize the trajectories of the protons in time
This is the pipeline that I have used to produce the visualizations below.
If you make any other 3D visualization with this data, let me know! We can add it here 😉!
And that’s all for now! 👋
💡 The external B field is loaded from an openPMD file, while the protons are defined as test particles.
📷 To analyze and visualize the simulation results in openPMD format, you
can use the openPMD-viewer
library for Python or you can open .pmd files directly in
Paraview.
Content from A FODO Cell
Last updated on 2026-02-03 | Edit this page
Estimated time: 30 minutes
Overview
Questions
How to simulate a particle beam travelling through a FODO cell 🧲🛝🧲?
Objectives
Learn how to setup, run, and visualize a particle beam simulation through a FODO (Focusing-Defocusing) cell with WarpX 🎢
Setup
In this example we will simulate a particle beam travelling through a FODO cell. A FODO cell is a periodic focusing structure used in particle accelerators, consisting of alternating Focusing (F) and Defocusing (D) quadrupole magnets, with O (drift) sections in between. The FODO cell is one of the most fundamental building blocks in accelerator physics, providing transverse focusing to keep particle beams confined as they travel through the accelerator. The alternating focusing and defocusing quadrupoles create a net focusing effect in both transverse planes, allowing the beam to be transported over long distances while maintaining its size.
WarpX for Beam Dynamics
WarpX is a Particle-in-Cell (PIC) code that can track particles embedded in external fields. Two important features are:
Space charge effects: You can turn space charge effects on or off using the appropriate input parameters. When enabled, WarpX solves the electromagnetic (or electrostatic) fields self-consistently from the particle distributions, allowing you to study how the beam’s own fields affect its dynamics.
Arbitrary external fields: WarpX allows you to define arbitrary external electromagnetic fields, either analytically or by loading them from files. This makes it possible to model complex accelerator elements like quadrupoles, dipoles, and other magnets.
WarpX is not necessarily the best tool for all beam dynamics simulations. If space charge effects are negligible and you’re primarily interested in linear beam optics, there are more specialized and computationally efficient tools available, such as ImpactX.
WarpX vs ImpactX
ImpactX is another code from the same ecosystem, the Beam, Plasma & Accelerator Simulation Toolkit (BLAST), designed for beam dynamics. See the ImpactX documentation and repository. Key differences:
-
WarpX — t-based (time-based):
solves Maxwell/Poisson on a grid, advances particles and fields in time.
Best for:
- Strong space charge, wakefields, full EM solutions
- Plasma–beam and other collective effects
- When time evolution of fields matters
-
ImpactX — z-based
(position-based): tracks particles along the beamline through lattice
elements. Best for:
- Weak or approximated space charge
- Optics design, linear dynamics, envelope tracking
- Fast runs when detailed field evolution is not needed
For this tutorial we use WarpX to set up a FODO simulation with optional space charge.
Equivalent ImpactX Example
If you’re interested in seeing how a similar FODO cell simulation can be set up using ImpactX, check out the ImpactX FODO cell example with 2D space charge using envelope tracking. This example demonstrates the envelope tracking approach.
Whenever you need to prepare an input file, this is where you want to go.
OUTPUT
####################
### MY CONSTANTS ###
####################
# BEAM
my_constants.mass = m_p # kg
my_constants.current = 0.5 # A
my_constants.energy = 6.7e6 * q_e # J
my_constants.velocity = sqrt(2 * energy / mass) # m/s
my_constants.gamma = energy / (mass * clight**2) + 1. # -
my_constants.beta = velocity / clight # -
my_constants.beta_x = 0.737881 # m
my_constants.beta_y = 0.737881 # m
my_constants.beta_t = 0.25 # m
my_constants.alpha_x = +2.4685083 # -
my_constants.alpha_y = -2.4685083 # -
my_constants.alpha_t = 0. # -
my_constants.gamma_x = (1+alpha_x**2)/beta_x # m^-1
my_constants.gamma_y = (1+alpha_y**2)/beta_y # m^-1
my_constants.gamma_t = (1+alpha_t**2)/beta_t # m^-1
my_constants.emitt_x = 1.0e-6 # m
my_constants.emitt_y = 1.0e-6 # m
my_constants.emitt_t = 1.0e-6 # m
my_constants.sigma_x = sqrt(emitt_x * beta_x) # m
my_constants.sigma_y = sqrt(emitt_y * beta_y) # m
my_constants.sigma_z = sqrt(emitt_t * beta_t) # m
my_constants.sigma_ux = beta * sqrt(emitt_x * gamma_x) # -
my_constants.sigma_uy = beta * sqrt(emitt_y * gamma_y) # -
my_constants.sigma_uz = beta * sqrt(emitt_t * gamma_t) # -
my_constants.sigma_t = sigma_z / velocity # s
my_constants.charge = current *sigma_t # C
# BOX
my_constants.Lx = 20*sigma_x # m
my_constants.Ly = 20*sigma_y # m
my_constants.zmin = -20*sigma_z # m
my_constants.zmax = 0
# LATTICE
my_constants.L_drift1 = 7.44e-2 # m
my_constants.L_drift2 = 14.88e-2 # m
my_constants.L_drift3 = 7.44e-2 # m
my_constants.quad_gradient = 38.64 # T/m
my_constants.L_quad = 6.10e-2 # m
# TIME
my_constants.T = (L_drift1+L_drift2+L_drift3+2*L_quad)/velocity # s
my_constants.dt = sigma_t/4 # s
my_constants.nt = floor(T/dt) # -
##########################
### GENERAL PARAMETERS ###
##########################
stop_time = T
geometry.dims = 3
geometry.prob_lo = -0.5*Lx -0.5*Ly zmin
geometry.prob_hi = +0.5*Lx +0.5*Ly zmax
amr.n_cell = 64 64 128
amr.max_level = 0
warpx.limit_verbose_step = 1
warpx.const_dt = dt
warpx.do_electrostatic = labframe
warpx.poisson_solver = fft
warpx.do_moving_window = 1
warpx.moving_window_dir = z
warpx.moving_window_v = beta # in units of c
algo.particle_shape = 3
boundary.field_lo = open open open
boundary.field_hi = open open open
#################
### PARTICLES ###
#################
particles.species_names = beam
beam.species_type = proton
beam.injection_style = gaussian_beam
beam.x_rms = sigma_x
beam.y_rms = sigma_y
beam.z_rms = sigma_z
beam.x_m = 0
beam.y_m = 0
beam.z_m = 0.5*zmin
beam.npart = 1e3
beam.q_tot = charge
beam.momentum_distribution_type = gaussian
beam.uz_m = beta*gamma
beam.uy_m = 0.0
beam.ux_m = 0.0
beam.ux_th = sigma_ux
beam.uy_th = sigma_uy
beam.uz_th = sigma_uz
###############
### LATTICE ###
###############
lattice.elements = drift1 quad1 drift2 quad2 drift3
drift1.type = drift
drift1.ds = L_drift1
quad1.type = quad
quad1.ds = L_quad
quad1.dBdx = -quad_gradient
drift2.type = drift
drift2.ds = L_drift2
quad2.type = quad
quad2.ds = L_quad
quad2.dBdx = quad_gradient
drift3.type = drift
drift3.ds = L_drift3
###################
### DIAGNOSTICS ###
###################
warpx.reduced_diags_names = beam_stats
beam_stats.type = BeamRelevant
beam_stats.species = beam
beam_stats.intervals = 1
diagnostics.diags_names = particles_in
particles_in.intervals = floor(nt/20)
particles_in.diag_type = Full
particles_in.species = beam
particles_in.fields_to_plot = none
particles_in.format = openpmd
particles_in.openpmd_backend = bp
particles_in.dump_last_timestep = 1Some notable details:
Space charge effects: Space charge is computed self-consistently by solving the Poisson equation at each timestep using the updated charge density from the particle distribution. This allows the beam’s own fields to affect its dynamics.
-
Turning off space charge: To turn off space charge effects, set
algo.maxwell_solver = noneand comment out the following inputs: External fields: To initialize an arbitrary external field (such as the quadrupole magnets in a FODO cell), check out the external fields documentation for details on how to define fields analytically or load them from files.
Run
First things first. Create a new folder where you copy the input file. The simulation should be small enough that you can run it in serial with the Conda installation of WarpX.
If you want to make the simulation faster and/or bigger, then you should run with the parallel version of WarpX. The optimal setup to run the simulation depends on your hardware. This is an example that should work on many common laptops, even though it might not be ideal.
With the default parameters, it may take a while on a laptop. If you have access to a GPU, you can experience the speed-up yourself. My experiments: 35 seconds on an NVIDIA A100 GPU, more than 1 hour on 4 cores of my 12th Gen Intel(R) Core(TM) i9-12900H.
Visualize
With Python 🐍
Now that we have the results, we can analyze them using Python.
We will use the openPMD-viewer library
to grab the data that the simulation produced in openPMD
format. Here you can find a
few tutorials on how to use the viewer. If you feel nerdy and/or you
need to deal with the data in parallel workflows, you can use the openPMD-api.
As an example for the FODO cell simulation, we have developed a simple Jupyter notebook where we retrieve the beam variables and analyze the beam dynamics through the cell.
Here is a video of the beam travelling through the FODO cell:
Nbviewer does not load this notebook; you can view it on GitHub (GitHub renders it in the browser) or download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook. If the simulation is too expensive,
🎯 Particle tracking follows individual particles, while envelope tracking follows the beam envelope - particle tracking is more detailed but computationally expensive.
🔬 A FODO cell is a periodic focusing structure with alternating focusing and defocusing quadrupole magnets.
⚡ WarpX allows you to turn space charge effects on/off and define arbitrary external fields.
🚀 WarpX is a full PIC code best suited for strong space charge effects and complex electromagnetic interactions. For cases with negligible space charge, ImpactX may be a more efficient choice.
🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.
📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.