Acceleration by RF Cavities
Beam accelerated through a sequence of 4 RF cavities (without space charge).
We use a 230 MeV electron beam with initial normalized rms emittance of 1 um.
The lattice and beam parameters are based on Example 2 of the IMPACT-Z examples folder:
https://github.com/impact-lbl/IMPACT-Z/tree/master/examples/Example2
The final target beam energy and beam moments are based on simulation in IMPACT-Z, without space charge.
In this test, the initial and final values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must agree with nominal values.
Run
This example can be run either as:
Python script:
python3 run_rfcavity.pyorImpactX executable using an input file:
impactx input_rfcavity.in
For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Marco Garten, Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-
from impactx import ImpactX, distribution, elements
sim = ImpactX()
# set numerical parameters and IO control
sim.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = False
# domain decomposition & space charge mesh
sim.init_grids()
# load a 230 MeV electron beam with an initial
# unnormalized rms emittance of 1 mm-mrad in all
# three phase planes
kin_energy_MeV = 230.0 # reference energy
bunch_charge_C = 1.0e-10 # used with space charge
npart = 10000 # number of macro particles (outside tests, use 1e5 or more)
# reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(-1.0).set_mass_MeV(0.510998950).set_kin_energy_MeV(kin_energy_MeV)
# particle bunch
distr = distribution.Waterbag(
lambdaX=0.352498964601e-3,
lambdaY=0.207443478142e-3,
lambdaT=0.70399950746e-4,
lambdaPx=5.161852770e-6,
lambdaPy=9.163582894e-6,
lambdaPt=0.260528852031e-3,
muxpx=0.5712386101751441,
muypy=-0.514495755427526,
mutpt=-5.05773e-10,
)
sim.add_particles(bunch_charge_C, distr, npart)
# design the accelerator lattice
# Drift elements
dr1 = elements.Drift(name="dr1", ds=0.4, nslice=1)
dr2 = elements.Drift(name="dr2", ds=0.032997, nslice=1)
# RF cavity element
rf = elements.RFCavity(
name="rf",
ds=1.31879807,
escale=62.0,
freq=1.3e9,
phase=85.5,
cos_coefficients=[
0.1644024074311037,
-0.1324009958969339,
4.3443060026047219e-002,
8.5602654094946495e-002,
-0.2433578169042885,
0.5297150596779437,
0.7164884680963959,
-5.2579522442877296e-003,
-5.5025369142193678e-002,
4.6845673335028933e-002,
-2.3279346335638568e-002,
4.0800777539657775e-003,
4.1378326533752169e-003,
-2.5040533340490805e-003,
-4.0654981400000964e-003,
9.6630592067498289e-003,
-8.5275895985990214e-003,
-5.8078747006425020e-002,
-2.4044337836660403e-002,
1.0968240064697212e-002,
-3.4461179858301418e-003,
-8.1201564869443749e-004,
2.1438992904959380e-003,
-1.4997753525697276e-003,
1.8685171825676386e-004,
],
sin_coefficients=[
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
],
mapsteps=100,
nslice=4,
)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
sim.lattice.extend(
[
monitor,
dr1,
dr2,
rf,
dr2,
dr2,
rf,
dr2,
dr2,
rf,
dr2,
dr2,
rf,
dr2,
monitor,
]
)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000 # outside tests, use 1e5 or more
beam.units = static
beam.kin_energy = 230
beam.charge = 1.0e-10
beam.particle = electron
beam.distribution = waterbag
beam.lambdaX = 0.352498964601e-3
beam.lambdaY = 0.207443478142e-3
beam.lambdaT = 0.70399950746e-4
beam.lambdaPx = 5.161852770e-6
beam.lambdaPy = 9.163582894e-6
beam.lambdaPt = 0.260528852031e-3
beam.muxpx = 0.5712386101751441
beam.muypy = -0.514495755427526
beam.mutpt = -5.05773e-10
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor dr1 dr2 rf dr2 dr2 rf dr2 dr2 rf dr2 dr2 rf dr2 monitor
monitor.type = beam_monitor
monitor.backend = h5
dr1.type = drift
dr1.ds = 0.4
dr1.nslice = 1
dr2.type = drift
dr2.ds = 0.032997
dr1.nslice = 1
rf.type = rfcavity
rf.ds = 1.31879807
rf.escale = 62.0
rf.freq = 1.3e9
rf.phase = 85.5
rf.mapsteps = 100
rf.nslice = 4
rf.cos_coefficients = \
0.1644024074311037 \
-0.1324009958969339 \
4.3443060026047219e-002 \
8.5602654094946495e-002 \
-0.2433578169042885 \
0.5297150596779437 \
0.7164884680963959 \
-5.2579522442877296e-003 \
-5.5025369142193678e-002 \
4.6845673335028933e-002 \
-2.3279346335638568e-002 \
4.0800777539657775e-003 \
4.1378326533752169e-003 \
-2.5040533340490805e-003 \
-4.0654981400000964e-003 \
9.6630592067498289e-003 \
-8.5275895985990214e-003 \
-5.8078747006425020e-002 \
-2.4044337836660403e-002 \
1.0968240064697212e-002 \
-3.4461179858301418e-003 \
-8.1201564869443749e-004 \
2.1438992904959380e-003 \
-1.4997753525697276e-003 \
1.8685171825676386e-004
rf.sin_coefficients = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 \
0 0 0 0 0 0 0
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = false
Analyze
We run the following script to analyze correctness:
Script analysis_rfcavity.py
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
import numpy as np
import openpmd_api as io
from scipy.stats import moment
def get_moments(beam):
"""Calculate standard deviations of beam position & momenta
and emittance values
Returns
-------
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
"""
sigx = moment(beam["position_x"], moment=2) ** 0.5 # variance -> std dev.
sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
sigy = moment(beam["position_y"], moment=2) ** 0.5
sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
sigt = moment(beam["position_t"], moment=2) ** 0.5
sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
epstrms = beam.cov(ddof=0)
emittance_x = (sigx**2 * sigpx**2 - epstrms["position_x"]["momentum_x"] ** 2) ** 0.5
emittance_y = (sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2) ** 0.5
emittance_t = (sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2) ** 0.5
return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
# initial/final beam
series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
last_step = list(series.iterations)[-1]
initial = series.iterations[1].particles["beam"].to_df()
final = series.iterations[last_step].particles["beam"].to_df()
# compare number of particles
num_particles = 10000
assert num_particles == len(initial)
assert num_particles == len(final)
print("Initial Beam:")
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(initial)
print(f" sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(
f" emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)
atol = 0.0 # ignored
rtol = 1.5 * num_particles**-0.5 # from random sampling of a smooth distribution
print(f" rtol={rtol} (ignored: atol~={atol})")
assert np.allclose(
[sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
[
4.29466150443e-4,
2.41918588389e-4,
7.0399951912e-5,
2.21684103818e-9,
2.21684103818e-9,
1.83412186547e-8,
],
rtol=rtol,
atol=atol,
)
print("")
print("Final Beam:")
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(final)
print(f" sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(
f" emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)
atol = 0.0 # ignored
rtol = 1.5 * num_particles**-0.5 # from random sampling of a smooth distribution
print(f" rtol={rtol} (ignored: atol~={atol})")
assert np.allclose(
[sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
[
3.52596000000e-4,
2.41775000000e-4,
7.0417917357e-5,
1.70893497973e-9,
1.70893497973e-9,
1.413901564889e-8,
],
rtol=rtol,
atol=atol,
)
Acceleration by RF Cavities (Reference Particle Tracking)
This test is identical to the test examples-rfcavity above, except that the code is run in reference orbit tracking mode.
It is used to validate the numerical integration of the reference energy gain in a chain of 4 RF cavities, and to demonstrate this tracking mode.
In this test, the initial and final reference values of \(s\) and \(\gamma\) must agree with nominal values.
Run
This example can be run either as:
Python script:
python3 run_rfcavity_ref_part.pyorImpactX executable using an input file:
impactx input_rfcavity_ref_part.in
For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Marco Garten, Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-
from impactx import ImpactX, elements
sim = ImpactX()
# set numerical parameters and IO control
sim.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
# domain decomposition & space charge mesh
sim.init_grids()
# load a 230 MeV electron beam with an initial
# unnormalized rms emittance of 1 mm-mrad in all
# three phase planes
kin_energy_MeV = 230.0 # reference energy
# reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(-1.0).set_mass_MeV(0.510998950).set_kin_energy_MeV(kin_energy_MeV)
# design the accelerator lattice
# Drift elements
dr1 = elements.Drift(name="dr1", ds=0.4, nslice=1)
dr2 = elements.Drift(name="dr2", ds=0.032997, nslice=1)
# RF cavity element
rf = elements.RFCavity(
name="rf",
ds=1.31879807,
escale=62.0,
freq=1.3e9,
phase=85.5,
cos_coefficients=[
0.1644024074311037,
-0.1324009958969339,
4.3443060026047219e-002,
8.5602654094946495e-002,
-0.2433578169042885,
0.5297150596779437,
0.7164884680963959,
-5.2579522442877296e-003,
-5.5025369142193678e-002,
4.6845673335028933e-002,
-2.3279346335638568e-002,
4.0800777539657775e-003,
4.1378326533752169e-003,
-2.5040533340490805e-003,
-4.0654981400000964e-003,
9.6630592067498289e-003,
-8.5275895985990214e-003,
-5.8078747006425020e-002,
-2.4044337836660403e-002,
1.0968240064697212e-002,
-3.4461179858301418e-003,
-8.1201564869443749e-004,
2.1438992904959380e-003,
-1.4997753525697276e-003,
1.8685171825676386e-004,
],
sin_coefficients=[
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
],
mapsteps=100,
nslice=4,
)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
sim.lattice.extend(
[
monitor,
dr1,
dr2,
rf,
dr2,
dr2,
rf,
dr2,
dr2,
rf,
dr2,
dr2,
rf,
dr2,
monitor,
]
)
# run simulation
sim.track_reference(ref)
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.kin_energy = 230
beam.particle = electron
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor dr1 dr2 rf dr2 dr2 rf dr2 dr2 rf dr2 dr2 rf dr2 monitor
monitor.type = beam_monitor
monitor.backend = h5
dr1.type = drift
dr1.ds = 0.4
dr1.nslice = 1
dr2.type = drift
dr2.ds = 0.032997
dr1.nslice = 1
rf.type = rfcavity
rf.ds = 1.31879807
rf.escale = 62.0
rf.freq = 1.3e9
rf.phase = 85.5
rf.mapsteps = 100
rf.nslice = 4
rf.cos_coefficients = \
0.1644024074311037 \
-0.1324009958969339 \
4.3443060026047219e-002 \
8.5602654094946495e-002 \
-0.2433578169042885 \
0.5297150596779437 \
0.7164884680963959 \
-5.2579522442877296e-003 \
-5.5025369142193678e-002 \
4.6845673335028933e-002 \
-2.3279346335638568e-002 \
4.0800777539657775e-003 \
4.1378326533752169e-003 \
-2.5040533340490805e-003 \
-4.0654981400000964e-003 \
9.6630592067498289e-003 \
-8.5275895985990214e-003 \
-5.8078747006425020e-002 \
-2.4044337836660403e-002 \
1.0968240064697212e-002 \
-3.4461179858301418e-003 \
-8.1201564869443749e-004 \
2.1438992904959380e-003 \
-1.4997753525697276e-003 \
1.8685171825676386e-004
rf.sin_coefficients = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 \
0 0 0 0 0 0 0
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
algo.track = "reference_orbit"
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true