Expanding Beam Scraping Against a Vacuum Pipe
This example describes a coasting bunch, expanding transversely and encountering the aperture defined by the vacuum pipe. Space charge is neglected, making the problem analytically soluble.
We use a cold (zero emittance) 250 MeV proton bunch whose initial distribution is a uniformly-populated cylinder of transverse radius \(r_b = 2 \mathrm{mm}\) with zero momentum spread.
The beam propagates in a drift with a vacuum chamber radius of \(R = 3.5 \mathrm{mm}\).
To generate an expanding beam, a linear map is first applied. This map applies a radial kick to each particle that is proportional to the particle’s initial distance from the axis. This induces a phase space correlation within the beam, such that \(p_x = k \cdot x\) and \(p_y = k \cdot y\), similar to what would be induced by a space charge kick.
The beam remains cylindrical with zero emittance during its evolution in a 6 m drift. In the absence of an aperture, the beam radius evolves as \(r_b(s) = r_b(1 + k\cdot s)\). In the presence of an aperture, particles are lost during the transverse expansion. The fraction of charge remaining after a distance s is given by:
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.
In addition, the initial and final values of \(\sigma_{p_x}\), \(\sigma_{p_y}\), and \(\sigma_{p_t}\) must agree with nominal values.
Finally, the fraction of charge lost against the aperture at the exit of the drift must agree with nominal values.
The physical problem is defined by four relevant parameters, defined within run_scraping.py, that can be modified by the user:
# problem parameters
beam_radius = r_b
aperture_radius = R
correlation_k = k
drift_distance = s
These parameters should also be modified inside analysis_scraping.py for testing.
Run
This example can be run as a Python script (python3 run_scraping.py) or with an app with an input file (impactx input_scraping.in).
Each can also be prefixed with an MPI executor, such as mpiexec -n 4 ... or srun -n 4 ..., depending on the system.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-
from impactx import ImpactX, Map6x6, distribution, elements
sim = ImpactX()
# set numerical parameters and IO control
sim.max_level = 1
sim.n_cell = [16, 16, 20]
sim.blocking_factor_x = [16]
sim.blocking_factor_y = [16]
sim.blocking_factor_z = [4]
sim.space_charge = False
sim.poisson_solver = "fft"
sim.dynamic_size = True
sim.prob_relative = [1.2, 1.1]
# beam diagnostics
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
# domain decomposition & space charge mesh
sim.init_grids()
# load a 2 GeV electron beam with an initial
# unnormalized rms emittance of 2 nm
kin_energy_MeV = 250 # reference energy
bunch_charge_C = 1.0e-9 # used with space charge
npart = 100000
# reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(1.0).set_mass_MeV(938.27208816).set_kin_energy_MeV(kin_energy_MeV)
# problem parameters
beam_radius = 2.0e-3
aperture_radius = 3.5e-3
correlation_k = 0.5
drift_distance = 6.0
# particle bunch
distr = distribution.Kurth4D(
lambdaX=beam_radius / 2.0,
lambdaY=beam_radius / 2.0,
lambdaT=beam_radius / 2.0,
lambdaPx=0.0,
lambdaPy=0.0,
lambdaPt=0.0,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# initialize the linear map
Iden = Map6x6.identity()
Rmat = Iden
Rmat[2, 1] = correlation_k
Rmat[4, 3] = correlation_k
# elements
drift1 = elements.Drift(
name="d1",
ds=drift_distance,
aperture_x=aperture_radius,
aperture_y=aperture_radius,
nslice=40,
)
map1 = elements.LinearMap(R=Rmat)
# design the accelerator lattice
sim.lattice.extend([monitor, map1, drift1, monitor])
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 100000
beam.units = static
beam.kin_energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = kurth4d
beam.lambdaX = 1.0e-3
beam.lambdaY = beam.lambdaX
beam.lambdaT = 1.0e-3
beam.lambdaPx = 0.0
beam.lambdaPy = 0.0
beam.lambdaPt = 0.0
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor map1 drift1 monitor
lattice.nslice = 40
map1.type = linear_map
map1.R21 = 0.5
map1.R43 = map1.R21
drift1.type = drift
drift1.ds = 6.0
drift1.aperture_x = 3.5e-3
drift1.aperture_y = 3.5e-3
monitor.type = beam_monitor
monitor.backend = h5
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
Analyze
We run the following script to analyze correctness:
Script analysis_scraping.py
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-
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, sigpx, sigpy, sigpt, 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,
sigpx,
sigpy,
sigpt,
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_beam = series.iterations[1].particles["beam"]
final_beam = series.iterations[last_step].particles["beam"]
initial = initial_beam.to_df()
final = final_beam.to_df()
# compare number of particles
num_particles = 100000
assert num_particles == len(initial)
# problem parameters
beam_radius = 2.0e-3
aperture_radius = 3.5e-3
correlation_k = 0.5
drift_distance = 6.0
print("Initial Beam:")
sigx, sigy, sigt, sigpx, sigpy, sigpt, emittance_x, emittance_y, emittance_t = (
get_moments(initial)
)
print(f" sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(f" sigpx={sigpx:e} sigpy={sigpy:e} sigpt={sigpt:e}")
print(
f" emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)
atol = 0.0 # ignored
rtol = 2.0 * num_particles**-0.5 # from random sampling of a smooth distribution
print(f" rtol={rtol} (ignored: atol~={atol})")
assert np.allclose(
[sigx, sigy, sigt],
[
beam_radius / 2.0,
beam_radius / 2.0,
beam_radius / 2.0,
],
rtol=rtol,
atol=atol,
)
atol = 1.0e-11 # ignored
print(f" atol={atol}")
assert np.allclose(
[sigpx, sigpy, sigpt, emittance_x, emittance_y, emittance_t],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
atol=atol,
)
# calculation of predicted final beam parameters
beam_radius_no_aperture = beam_radius * (1.0 + correlation_k * drift_distance)
beam_radius_with_aperture = min(beam_radius_no_aperture, aperture_radius)
fractional_loss = 1.0 - min(1.0, (aperture_radius / beam_radius_no_aperture) ** 2)
sigma_x_final = beam_radius_with_aperture / 2.0
sigma_px_final = correlation_k / (1.0 + correlation_k * drift_distance) * sigma_x_final
print("")
print("Predicted Final Beam:")
print(f" sigx={sigma_x_final:e} sigy={sigma_x_final:e} sigt={beam_radius / 2.0:e}")
print(f" sigpx={sigma_px_final:e} sigpy={sigma_px_final:e} sigpt=0.0")
print(f" fractional_loss={fractional_loss:e}")
print("")
print("Final Beam:")
sigx, sigy, sigt, sigpx, sigpy, sigpt, emittance_x, emittance_y, emittance_t = (
get_moments(final)
)
print(f" sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(f" sigpx={sigpx:e} sigpy={sigpy:e} sigpt={sigpt:e}")
print(
f" emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)
atol = 0.0 # ignored
rtol = 2.0 * num_particles**-0.5 # from random sampling of a smooth distribution
print(f" rtol={rtol} (ignored: atol~={atol})")
assert np.allclose(
[sigx, sigy, sigt, sigpx, sigpy],
[
sigma_x_final,
sigma_x_final,
beam_radius / 2.0,
sigma_px_final,
sigma_px_final,
],
rtol=rtol,
atol=atol,
)
charge_i = initial_beam.get_attribute("charge_C")
charge_f = final_beam.get_attribute("charge_C")
loss_pct = 100.0 * (charge_i - charge_f) / charge_i
print(f" fractional loss (%) = {loss_pct}")
atol = 0.2 # tolerance 0.2%
print(f" atol={atol}")
assert np.allclose(
[loss_pct],
[100 * fractional_loss],
atol=atol,
)