Quadrupole with Alignment Errors

A 2 GeV proton beam propagates through a single quadrupole with 3 mm horizontal misalignment and 30 degree rotation error.

The first and second moments of the particle distribution before and after the quadrupole should coincide with analytical predictions, to within the level expected due to noise due to statistical sampling.

In this test, the initial and final values of \(\mu_x\), \(\mu_y\), \(\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_alignment.py or

• ImpactX executable using an input file: impactx input_alignment.in

For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Python: Script

Listing 12 You can copy this file from examples/alignment/run_alignment.py.

#!/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 a 2 GeV proton beam
kin_energy_MeV = 2.0e3  # reference energy
bunch_charge_C = 1.0e-9  # used with space charge
npart = 100000  # 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.16098260008648811e-3,
    lambdaY=1.16098260008648811e-3,
    lambdaT=1.0e-3,
    lambdaPx=0.580491300043e-3,
    lambdaPy=0.580491300043e-3,
    lambdaPt=2.0e-3,
    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
lattice = [
    monitor,
    elements.Quad(ds=1.0, k=0.25, dx=0.003, dy=0.0, rotation=30.0, nslice=ns),
    monitor,
]
sim.lattice.extend(lattice)

# run simulation
sim.evolve()

# clean shutdown
sim.finalize()

Executable: Input File

Listing 13 You can copy this file from examples/alignment/input_alignment.in.

###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 100000
beam.units = static
beam.kin_energy = 2.0e3
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.16098260008648811e-3
beam.lambdaY = 1.16098260008648811e-3
beam.lambdaT = 1.0e-3
beam.lambdaPx = 0.580491300043e-3
beam.lambdaPy = 0.580491300043e-3
beam.lambdaPt = 2.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0


###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor quad_err monitor
lattice.nslice = 1

monitor.type = beam_monitor
monitor.backend = h5

quad_err.type = quad
quad_err.ds = 1.0
quad_err.k = 0.25
quad_err.dx = 0.003
quad_err.rotation = 30.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_alignment.py

Listing 14 You can copy this file from examples/alignment/analysis_alignment.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, tmean


def get_moments(beam):
    """Calculate standard deviations of beam position & momenta
    and emittance values

    Returns
    -------
    meanx, meany, sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
    """
    meanx = tmean(beam["position_x"])
    meany = tmean(beam["position_y"])
    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 (meanx, meany, 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 = 100000
assert num_particles == len(initial)
assert num_particles == len(final)

print("Initial Beam:")
meanx, meany, sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(
    initial
)
print(f"  meanx={meanx:e} meany={meany:e}")
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 = 3.0 * num_particles ** (-0.5) * sigx
print(f" atol~={atol}")

assert np.allclose(
    [meanx, meany],
    [
        0.0,
        0.0,
    ],
    atol=atol,
)

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, emittance_x, emittance_y, emittance_t],
    [
        1.160982600086e-3,
        1.160982600086e-3,
        1.0e-3,
        6.73940299e-7,
        6.73940299e-7,
        2.0e-6,
    ],
    rtol=rtol,
    atol=atol,
)


print("")
print("Final Beam:")
meanx, meany, sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(
    final
)
print(f"  meanx={meanx:e} meany={meany:e}")
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 = 3.0 * num_particles ** (-0.5) * sigx
print(f" atol~={atol}")

assert np.allclose(
    [meanx, meany],
    [
        1.79719761842e-4,
        3.24815908981e-4,
    ],
    atol=atol,
)

atol = 0.0  # ignored
rtol = 3.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, emittance_x, emittance_y, emittance_t],
    [
        1.2372883901369e-3,
        1.3772750218080e-3,
        1.027364e-03,
        7.39388142e-7,
        7.39388142e-7,
        2.0e-6,
    ],
    rtol=rtol,
    atol=atol,
)