Thin Dipole
This test involves tracking a 5 MeV proton beam through a 90 degree sector bend, using two different methods:
#. Using a sequence of drifts and thin kicks as described by: G. Ripken and F. Schmidt, “A Symplectic Six-Dimensional Thin-Lens Formalism for Tracking,” CERN/SL/95-12 (1995).
#. Using the exact nonlinear transfer map for a sector bend as described by: D. Bruhwiler et al, “Symplectic Propagation of the Map, Tangent Map and Tangent Map Derivative through Quadrupole and Combined-Function Dipole Magnets without Truncation,” THP41C, EPAC98, pp. 1171-1173 (1998).
To compare the two models, the beam is tracked using 1), followed by the inverse of 2). The beam moments and emittances should coincide with those of the initial distribution.
Run
This example can be run either as:
Python script:
python3 run_thin_dipole.pyorImpactX executable using an input file:
impactx input_thin_dipole.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: 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 = True
# domain decomposition & space charge mesh
sim.init_grids()
# load initial beam
kin_energy_MeV = 5.0 # reference energy
bunch_charge_C = 1.0e-9 # used with space charge
npart = 10000 # number of macro particles
# 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)
# particle bunch
distr = distribution.Waterbag(
lambdaX=1.0e-3,
lambdaY=1.0e-3,
lambdaT=0.3,
lambdaPx=2.0e-4,
lambdaPy=2.0e-4,
lambdaPt=2.0e-4,
muxpx=-0.0,
muypy=0.0,
mutpt=0.0,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice)
ns = 1 # number of slices per ds in the element
segment = [
elements.Drift(ds=0.003926990816987, nslice=ns),
elements.ThinDipole(theta=0.45, rc=1.0),
elements.Drift(ds=0.003926990816987, nslice=ns),
]
bend = 200 * segment
inverse_bend = elements.ExactSbend(ds=-1.570796326794897, phi=-90.0)
sim.lattice.append(monitor)
sim.lattice.extend(bend)
sim.lattice.append(inverse_bend)
sim.lattice.append(monitor)
# run simulation
sim.evolve()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 5.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.0e-3
beam.lambdaY = 1.0e-3
beam.lambdaT = 0.3
beam.lambdaPx = 2.0e-4
beam.lambdaPy = 2.0e-4
beam.lambdaPt = 2.0e-4
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor bend inverse_bend monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
#90 degree sbend using drift-kick-drift:
bend.type = line
bend.elements = dr kick dr
bend.repeat = 200
dr.type = drift
dr.ds = 0.003926990816987
kick.type = thin_dipole
kick.theta = 0.45
kick.rc = 1.0
#inverse of a 90 degree sbend using exact nonlinear map:
inverse_bend.type = sbend_exact
inverse_bend.ds = -1.570796326794897
inverse_bend.phi = -90.0
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
Analyze
We run the following script to analyze correctness:
Script analysis_thin_dipole.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, 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.8 * 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],
[
1.0e-3,
1.0e-3,
3.0e-1,
2.0e-7,
2.0e-7,
6.0e-5,
],
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.8 * 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],
[
1.0e-3,
1.0e-3,
3.0e-1,
2.0e-7,
2.0e-7,
6.0e-5,
],
rtol=rtol,
atol=atol,
)
Visualize
Note
TODO :)