Iteration of a User-Defined Linear Map
This example illustrates the application of a user-defined linear map via a matrix.
Here, the linear map represents an abstract symplectic transformation of the beam in 6D phase space. If desired, the user may interpret the matrix as the one-turn map of a storage ring or circular collider.
The (fractional) tunes (Qx, Qy, Qt) of the map are given by (0.139, 0.219, 0.0250). We use a 45.6 GeV electron beam that is invariant under the action of the linear map (matched). The horizontal and vertical unnormalized emittances are 0.27 nm and 1.0 pm, respectively.
These parameters are based on the single-beam parameters of FCC-ee (Z-mode). (backup).
The second moments of the phase space variables should be unchanged under application of the map.
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 tunes associated with a single particle orbit are extracted, and must agree with the values given above.
Run
This example can be run either as:
Python script:
python3 run_map.pyorImpactX executable using an input file:
impactx input_map.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 elements import LinearTransport
import numpy as np
from impactx import ImpactX, Map6x6, distribution, elements, twiss
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 2 GeV electron beam with an initial
# unnormalized rms emittance of 2 nm
kin_energy_MeV = 45.6e3 # 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(0.510998950).set_kin_energy_MeV(kin_energy_MeV)
# target beta functions (m)
beta_star_x = 0.15
beta_star_y = 0.8e-3
beta_star_t = 9.210526315789473
# particle bunch
distr = distribution.Waterbag(
**twiss(
beta_x=beta_star_x,
beta_y=beta_star_y,
beta_t=beta_star_t,
emitt_x=0.27e-09,
emitt_y=1.0e-12,
emitt_t=1.33e-06,
alpha_x=0.0,
alpha_y=0.0,
alpha_t=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
# desired tunes
Qx = 0.139
Qy = 0.219
Qt = 0.0250
# desired phase advance
phi_x = 2.0 * np.pi * Qx
phi_y = 2.0 * np.pi * Qy
phi_t = 2.0 * np.pi * Qt
# matrix elements for the horizontal plane
Rmat[1, 1] = np.cos(phi_x)
Rmat[1, 2] = beta_star_x * np.sin(phi_x)
Rmat[2, 1] = -np.sin(phi_x) / beta_star_x
Rmat[2, 2] = np.cos(phi_x)
# matrix elements for the vertical plane
Rmat[3, 3] = np.cos(phi_y)
Rmat[3, 4] = beta_star_y * np.sin(phi_y)
Rmat[4, 3] = -np.sin(phi_y) / beta_star_y
Rmat[4, 4] = np.cos(phi_y)
# matrix elements for the longitudinal plane
Rmat[5, 5] = np.cos(phi_t)
Rmat[5, 6] = beta_star_t * np.sin(phi_t)
Rmat[6, 5] = -np.sin(phi_t) / beta_star_t
Rmat[6, 6] = np.cos(phi_t)
# design the accelerator lattice
map = [
monitor,
elements.LinearMap(R=Rmat),
]
sim.lattice.extend(map)
# number of periods through the lattice
sim.periods = 4
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 45.6e3
beam.charge = 2.72370027e-8 #population 1.7e11
beam.particle = electron
beam.distribution = waterbag_from_twiss
beam.alphaX = 0.0
beam.alphaY = 0.0
beam.alphaT = 0.0
beam.betaX = 0.15
beam.betaY = 0.8e-3
beam.betaT = 9.210526315789473
beam.emittX = 0.27e-09
beam.emittY = 1e-12
beam.emittT = 1.33e-6
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.periods = 5
lattice.elements = monitor map1
monitor.type = beam_monitor
monitor.backend = h5
map1.type = linear_map
# horizontal plane
map1.R11 = 0.642252653176584
map1.R12 = 0.114973951021402
map1.R21 = -5.109953378728999
map1.R22 = 0.642252653176584
# vertical plane
map1.R33 = 0.193549468050860
map1.R34 = 0.0007848724139547
map1.R43 = -1226.363146804167548
map1.R44 = 0.193549468050860
# longitudinal plane
map1.R55 = 0.987688340595138
map1.R56 = 1.440843756949495
map1.R65 = -0.016984313347225
map1.R66 = 0.987688340595138
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = false
Analyze
We run the following script to analyze correctness:
Script analysis_map.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 = 2.2 * 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],
[
6.363961030678928e-6,
28.284271247461902e-9,
0.0035,
0.27e-9,
1.0e-12,
1.33e-6,
],
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 = 2.2 * 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],
[
6.363961030678928e-6,
28.284271247461902e-9,
0.0035,
0.27e-9,
1.0e-12,
1.33e-6,
],
rtol=rtol,
atol=atol,
)
# Specify time series for particle j
j = 5
print(f"output for particle index = {j}")
# Create array of TBT data values
x = []
px = []
y = []
py = []
t = []
pt = []
n = 0
for k_i, i in series.iterations.items():
beam = i.particles["beam"]
turn = beam.to_df()
x.append(turn["position_x"][j])
px.append(turn["momentum_x"][j])
y.append(turn["position_y"][j])
py.append(turn["momentum_y"][j])
t.append(turn["position_t"][j])
pt.append(turn["momentum_t"][j])
n = n + 1
# Output number of periods in data series
nturns = len(x)
print(f"number of periods = {nturns}")
print()
# Approximate the tune and closed orbit using the 4-turn formula:
# from x data only
argument = (x[0] - x[1] + x[2] - x[3]) / (2.0 * (x[1] - x[2]))
tunex = np.arccos(argument) / (2.0 * np.pi)
print(f"tune output from 4-turn formula, using x data = {tunex}")
# from y data only
argument = (y[0] - y[1] + y[2] - y[3]) / (2.0 * (y[1] - y[2]))
tuney = np.arccos(argument) / (2.0 * np.pi)
print(f"tune output from 4-turn formula, using y data = {tuney}")
# from t data only
argument = (t[0] - t[1] + t[2] - t[3]) / (2.0 * (t[1] - t[2]))
tunet = np.arccos(argument) / (2.0 * np.pi)
print(f"tune output from 4-turn formula, using t data = {tunet}")
rtol = 1.0e-3
print(f" rtol={rtol}")
assert np.allclose(
[tunex, tuney, tunet],
[
0.139,
0.219,
0.0250,
],
rtol=rtol,
)