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 \(\lambda_x\), \(\lambda_y\), \(\lambda_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.particle_shape = 2 # B-spline order
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(ds=0.4, nslice=1)
dr2 = elements.Drift(ds=0.032997, nslice=1)
# RF cavity element
rf = elements.RFCavity(
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.evolve()
# 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.particle_shape = 2
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,
)