Content from Introduction to WarpX

Last updated on 2025-04-22 | Edit this page

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.

macroparticles in the cells of a grid
Some macroparticles traveling in space, across the cells of a grid.

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.

pic loop
The loop at the basis of standard PIC codes.

If you want to know more about PIC, here are a few references:

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.

Checklist

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!

Key Points

🔮 The particle-in-cell method is used to simulate the self-consistent dynamics of relativistic charged particles

🚀 WarpX is a open-source high-performance particle-in-cell code

✨ WarpX is used in a variety of scientif domains

Content from Install

Last updated on 2025-04-22 | Edit this page

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 Python and Paraview

Basic dependencies

Just a heads-up before we dive deeper.

Callout

📣 Everything you need to know to use WarpX is in the documentation, check it out!

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 😌!

Callout

BASH 

conda install -c conda-forge warpx

Ok, maybe two lines if you want to keep your system clean by creating a new environment.

BASH 

conda create -n warpx -c conda-forge warpx 
conda activate warpx

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:

BASH 

which warpx.1d warpx.2d warpx.3d warpx.rz

If you get 3 different paths that look something like:

BASH 

/home/<username>/anaconda3/envs/warpx/bin/warpx.xd

then you got this 🙌! You can also import pywarpx in Python.

Caution

Conda’s WarpX is serial! To get a parallel WarpX version, install it from source.

From source

Caution

Coming soon. For now, refer to the main documentation.

Key Points

🎯 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-04-22 | Edit this page

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

snapshot of a beam-beam simulation
A snapshot of a 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)

#################################
####### 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 = beam1 beam2 pho1 pho2 ele_nlbw1 pos_nlbw1 ele_nlbw2 pos_nlbw2
particles.photon_species = pho1 pho2

beam1.species_type = electron
beam1.injection_style = gaussian_beam
beam1.x_rms = sigmax
beam1.y_rms = sigmay
beam1.z_rms = sigmaz
beam1.x_m = - mux
beam1.y_m = - muy
beam1.z_m = - muz
beam1.npart = nmacropart
beam1.q_tot = -beam_charge
beam1.z_cut = 4
beam1.focal_distance = muz
beam1.momentum_distribution_type = gaussian
beam1.uz_m = beam_gamma
beam1.uy_m = 0.0
beam1.ux_m = 0.0
beam1.ux_th = dux
beam1.uy_th = duy
beam1.uz_th = beam_uth
beam1.initialize_self_fields = 1
beam1.do_qed_quantum_sync = 1
beam1.qed_quantum_sync_phot_product_species = pho1
beam1.do_classical_radiation_reaction = 0

beam2.species_type = positron
beam2.injection_style = gaussian_beam
beam2.x_rms = sigmax
beam2.y_rms = sigmay
beam2.z_rms = sigmaz
beam2.x_m = mux
beam2.y_m = muy
beam2.z_m = muz
beam2.npart = nmacropart
beam2.q_tot = beam_charge
beam2.z_cut = 4
beam2.focal_distance = muz
beam2.momentum_distribution_type = gaussian
beam2.uz_m = -beam_gamma
beam2.uy_m = 0.0
beam2.ux_m = 0.0
beam2.ux_th = dux
beam2.uy_th = duy
beam2.uz_th = beam_uth
beam2.initialize_self_fields = 1
beam2.do_qed_quantum_sync = 1
beam2.qed_quantum_sync_phot_product_species = pho2
beam2.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 = ele_nlbw1
pho1.qed_breit_wheeler_pos_product_species = pos_nlbw1

pho2.species_type = photon
pho2.injection_style = none
pho2.do_qed_breit_wheeler = 1
pho2.qed_breit_wheeler_ele_product_species = ele_nlbw2
pho2.qed_breit_wheeler_pos_product_species = pos_nlbw2

ele_nlbw1.species_type = electron
ele_nlbw1.injection_style = none
ele_nlbw1.do_qed_quantum_sync = 1
ele_nlbw1.qed_quantum_sync_phot_product_species = pho1
ele_nlbw1.do_classical_radiation_reaction = 0

pos_nlbw1.species_type = positron
pos_nlbw1.injection_style = none
pos_nlbw1.do_qed_quantum_sync = 1
pos_nlbw1.qed_quantum_sync_phot_product_species = pho1
pos_nlbw1.do_classical_radiation_reaction = 0

ele_nlbw2.species_type = electron
ele_nlbw2.injection_style = none
ele_nlbw2.do_qed_quantum_sync = 1
ele_nlbw2.qed_quantum_sync_phot_product_species = pho2
ele_nlbw2.do_classical_radiation_reaction = 0

pos_nlbw2.species_type = positron
pos_nlbw2.injection_style = none
pos_nlbw2.do_qed_quantum_sync = 1
pos_nlbw2.qed_quantum_sync_phot_product_species = pho2
pos_nlbw2.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 = bound trajs

bound.dump_last_timestep = 1
bound.diag_type = BoundaryScraping
bound.format = openpmd
bound.openpmd_backend = h5
bound.intervals = 1

beam1.save_particles_at_xlo = 1
beam1.save_particles_at_ylo = 1
beam1.save_particles_at_zlo = 1
beam1.save_particles_at_xhi = 1
beam1.save_particles_at_yhi = 1
beam1.save_particles_at_zhi = 1

beam2.save_particles_at_xlo = 1
beam2.save_particles_at_ylo = 1
beam2.save_particles_at_zlo = 1
beam2.save_particles_at_xhi = 1
beam2.save_particles_at_yhi = 1
beam2.save_particles_at_zhi = 1

ele_nlbw1.save_particles_at_xlo = 1
ele_nlbw1.save_particles_at_ylo = 1
ele_nlbw1.save_particles_at_zlo = 1
ele_nlbw1.save_particles_at_xhi = 1
ele_nlbw1.save_particles_at_yhi = 1
ele_nlbw1.save_particles_at_zhi = 1

ele_nlbw2.save_particles_at_xlo = 1
ele_nlbw2.save_particles_at_ylo = 1
ele_nlbw2.save_particles_at_zlo = 1
ele_nlbw2.save_particles_at_xhi = 1
ele_nlbw2.save_particles_at_yhi = 1
ele_nlbw2.save_particles_at_zhi = 1

pos_nlbw1.save_particles_at_xlo = 1
pos_nlbw1.save_particles_at_ylo = 1
pos_nlbw1.save_particles_at_zlo = 1
pos_nlbw1.save_particles_at_xhi = 1
pos_nlbw1.save_particles_at_yhi = 1
pos_nlbw1.save_particles_at_zhi = 1

pos_nlbw2.save_particles_at_xlo = 1
pos_nlbw2.save_particles_at_ylo = 1
pos_nlbw2.save_particles_at_zlo = 1
pos_nlbw2.save_particles_at_xhi = 1
pos_nlbw2.save_particles_at_yhi = 1
pos_nlbw2.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

trajs.intervals = 1
trajs.diag_type = Full
trajs.species = beam1 beam2
trajs.fields_to_plot = none
trajs.format = openpmd
trajs.openpmd_backend = h5
trajs.dump_last_timestep = 1

# REDUCED
warpx.reduced_diags_names = DiffLumi_beam1_beam2 DiffLumi_beam1_ele2 DiffLumi_beam1_pos2 DiffLumi_beam1_pho2 DiffLumi_ele1_beam2 DiffLumi_ele1_ele2 DiffLumi_ele1_pos2 DiffLumi_ele1_pho2 DiffLumi_pos1_beam2 DiffLumi_pos1_ele2 DiffLumi_pos1_pos2 DiffLumi_pos1_pho2 DiffLumi_pho1_beam2 DiffLumi_pho1_ele2 DiffLumi_pho1_pos2 DiffLumi_pho1_pho2

DiffLumi_beam1_beam2.type = DifferentialLuminosity
DiffLumi_beam1_ele2.type  = DifferentialLuminosity
DiffLumi_beam1_pos2.type  = DifferentialLuminosity
DiffLumi_beam1_pho2.type  = DifferentialLuminosity
DiffLumi_ele1_beam2.type  = DifferentialLuminosity
DiffLumi_ele1_ele2.type   = DifferentialLuminosity
DiffLumi_ele1_pos2.type   = DifferentialLuminosity
DiffLumi_ele1_pho2.type   = DifferentialLuminosity
DiffLumi_pos1_beam2.type  = DifferentialLuminosity
DiffLumi_pos1_ele2.type   = DifferentialLuminosity
DiffLumi_pos1_pos2.type   = DifferentialLuminosity
DiffLumi_pos1_pho2.type   = DifferentialLuminosity
DiffLumi_pho1_beam2.type  = DifferentialLuminosity
DiffLumi_pho1_ele2.type   = DifferentialLuminosity
DiffLumi_pho1_pos2.type   = DifferentialLuminosity
DiffLumi_pho1_pho2.type   = DifferentialLuminosity

DiffLumi_beam1_beam2.bin_max = 2.1*beam_energy/q_e
DiffLumi_beam1_ele2.bin_max  = 2.1*beam_energy/q_e
DiffLumi_beam1_pos2.bin_max  = 2.1*beam_energy/q_e
DiffLumi_beam1_pho2.bin_max  = 2.1*beam_energy/q_e
DiffLumi_ele1_beam2.bin_max  = 2.1*beam_energy/q_e
DiffLumi_ele1_ele2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_ele1_pos2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_ele1_pho2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_pos1_beam2.bin_max  = 2.1*beam_energy/q_e
DiffLumi_pos1_ele2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_pos1_pos2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_pos1_pho2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_pho1_beam2.bin_max  = 2.1*beam_energy/q_e
DiffLumi_pho1_ele2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_pho1_pos2.bin_max   = 2.1*beam_energy/q_e
DiffLumi_pho1_pho2.bin_max   = 2.1*beam_energy/q_e

DiffLumi_beam1_beam2.bin_min = 0.
DiffLumi_beam1_ele2.bin_min = 0.
DiffLumi_beam1_pos2.bin_min = 0.
DiffLumi_beam1_pho2.bin_min = 0.
DiffLumi_ele1_beam2.bin_min = 0.
DiffLumi_ele1_ele2.bin_min = 0.
DiffLumi_ele1_pos2.bin_min = 0.
DiffLumi_ele1_pho2.bin_min = 0.
DiffLumi_pos1_beam2.bin_min = 0.
DiffLumi_pos1_ele2.bin_min = 0.
DiffLumi_pos1_pos2.bin_min = 0.
DiffLumi_pos1_pho2.bin_min = 0.
DiffLumi_pho1_beam2.bin_min = 0.
DiffLumi_pho1_ele2.bin_min = 0.
DiffLumi_pho1_pos2.bin_min = 0.
DiffLumi_pho1_pho2.bin_min = 0.

DiffLumi_beam1_beam2.species = beam1 beam2
DiffLumi_beam1_ele2.species = beam1 ele_nlbw2
DiffLumi_beam1_pos2.species = beam1 pos_nlbw2
DiffLumi_beam1_pho2.species = beam1 pho2
DiffLumi_ele1_beam2.species = ele_nlbw1 beam2
DiffLumi_ele1_ele2.species= ele_nlbw1 ele_nlbw2
DiffLumi_ele1_pos2.species= ele_nlbw1 pos_nlbw2
DiffLumi_ele1_pho2.species= ele_nlbw1 pho2
DiffLumi_pos1_beam2.species = pos_nlbw1 beam2
DiffLumi_pos1_ele2.species= pos_nlbw1 ele_nlbw2
DiffLumi_pos1_pos2.species= pos_nlbw1 pos_nlbw2
DiffLumi_pos1_pho2.species= pos_nlbw1 pho2
DiffLumi_pho1_beam2.species = pho1 beam2
DiffLumi_pho1_ele2.species= pho1 ele_nlbw2
DiffLumi_pho1_pos2.species= pho1 pos_nlbw2
DiffLumi_pho1_pho2.species= pho1 pho2


DiffLumi_beam1_beam2.bin_number = 256
DiffLumi_beam1_ele2.bin_number = 256
DiffLumi_beam1_pos2.bin_number = 256
DiffLumi_beam1_pho2.bin_number = 256
DiffLumi_ele1_beam2.bin_number = 256
DiffLumi_ele1_ele2.bin_number = 256
DiffLumi_ele1_pos2.bin_number = 256
DiffLumi_ele1_pho2.bin_number = 256
DiffLumi_pos1_beam2.bin_number = 256
DiffLumi_pos1_ele2.bin_number = 256
DiffLumi_pos1_pos2.bin_number = 256
DiffLumi_pos1_pho2.bin_number = 256
DiffLumi_pho1_beam2.bin_number = 256
DiffLumi_pho1_ele2.bin_number = 256
DiffLumi_pho1_pos2.bin_number = 256
DiffLumi_pho1_pho2.bin_number = 256


DiffLumi_beam1_beam2.intervals = nt
DiffLumi_beam1_ele2.intervals = nt
DiffLumi_beam1_pos2.intervals = nt
DiffLumi_beam1_pho2.intervals = nt
DiffLumi_ele1_beam2.intervals = nt
DiffLumi_ele1_ele2.intervals = nt
DiffLumi_ele1_pos2.intervals = nt
DiffLumi_ele1_pho2.intervals = nt
DiffLumi_pos1_beam2.intervals = nt
DiffLumi_pos1_ele2.intervals = nt
DiffLumi_pos1_pos2.intervals = nt
DiffLumi_pos1_pho2.intervals = nt
DiffLumi_pho1_beam2.intervals = nt
DiffLumi_pho1_ele2.intervals = nt
DiffLumi_pho1_pos2.intervals = nt
DiffLumi_pho1_pho2.intervals = nt

Some 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.

BASH 

warpx.3d inputs_3d_beambeam_C3.txt

Just like that! 💃 Note that with Conda’s WarpX you can run this anywhere in your filesystem (provided that you copied there the input of course) because WarpX’s executable are in your $PATH.

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.

This is just one way of doing it!

BASH 

export OMP_NUM_THREADS=2
mpirun -np 4 <path/to/your/build/bin/warpx.3d> inputs_3d_beambeam_C3.txt

Testimonial

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

Caution

Coming soon!

Key Points

💅 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-04-22 | Edit this page

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 = 1

A 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_style flag, 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? 🤷

BASH 

warpx.3d inputs_3d_magnetic_mirror.txt

As simple as that! 😉

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.

paraview pipeline
simulation of proton trajectories inside a magnetic mirror
Protons trajectories in a magnetic mirror

If you make any other 3D visualization with this data, let me know! We can add it here 😉!

And that’s all for now! 👋

Key Points

💡 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.