Testing Charge and Field Sign Consistency
ImpactX uses a normalized field (focusing) strength internally for tracking. For example, in a Quad a focusing strength k > 0 indicates horizontal focusing, independent of the charge of the particle (just as in MAD-X).
For elements that support an option for user-provided field strength in SI units (e.g. ChrQuad and ExactSbend), varying the sign of the charge should provide tracking consistent with the sign of the specified field.
This example tests the transport of a 1 GeV electron bunch through multiple field sign combinations. The transport occurs through pairs of elements. For example, the element
quad1 specifies a positive normalized field strength, and should provide horizontal focusing. The element quad2inv specifies a quadrupole field gradient in units of T/m, which
should provide an equivalent horizontal focusing. (For electrons, the required field gradient is negative.)
To test this, the map for quad1 is applied (ds > 0), and this is followed by applying the inverse map for quad2inv (ds < 0), which should result in the identity map. Similar
considerations apply to the other element pairs.
As a result, the second moments of x, y, and t and the associated emittances of the bunch should all be exactly unchanged.
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.
Run
This example can be run either as:
Python script:
python3 run_charge_sign.pyorImpactX executable using an input file:
impactx input_charge_sign.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 -*-
import numpy as np
from impactx import ImpactX, distribution, elements
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 initial beam
kin_energy_MeV = 1.0e3 # 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)
# 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-5,
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")
# problem parameters
rigidity = 3.337345025729098
g_gradient = 100.0
k_gradient = g_gradient / rigidity
length = 0.25
angle_deg = 45.0
angle_rad = np.pi * angle_deg / 180.0
field = rigidity * angle_rad / length
# design the accelerator lattice)
ns = 20 # number of slices per ds in the element
quad1 = elements.ChrQuad(name="quad1", ds=length, k=k_gradient, unit=0, nslice=ns)
quad2inv = elements.ChrQuad(
name="quad2inv", ds=-length, k=-g_gradient, unit=1, nslice=ns
)
quad3 = elements.ChrQuad(name="quad3", ds=length, k=-k_gradient, unit=0, nslice=ns)
quad4inv = elements.ChrQuad(
name="quad4inv", ds=-length, k=g_gradient, unit=1, nslice=ns
)
bend1 = elements.ExactSbend(name="bend1", ds=length, phi=angle_deg, B=field, nslice=ns)
bend2inv = elements.ExactSbend(name="bend2inv", ds=-length, phi=angle_deg, nslice=ns)
bend3 = elements.ExactSbend(name="bend3", ds=length, phi=angle_deg, B=-field, nslice=ns)
bend4inv = elements.ExactSbend(name="bend4inv", ds=-length, phi=-angle_deg, nslice=ns)
line = [
monitor,
quad1,
quad2inv,
quad3,
quad4inv,
bend1,
bend2inv,
bend3,
bend4inv,
monitor,
]
sim.lattice.extend(line)
# number of periods through the lattice
sim.periods = 1
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 1.0e3
beam.charge = 1.0e-9
beam.particle = electron
beam.distribution = waterbag
beam.lambdaX = 1.0e-3
beam.lambdaY = beam.lambdaX
beam.lambdaT = 0.3
beam.lambdaPx = 2.0e-4
beam.lambdaPy = beam.lambdaPx
beam.lambdaPt = 2.0e-5
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor quad1 quad2inv quad3 quad4inv bend1 bend2inv bend3 bend4inv monitor
lattice.nslice = 20
monitor.type = beam_monitor
monitor.backend = h5
quad1.type = quad_chromatic
quad1.ds = 0.25
quad1.k = 29.963938169129919
quad1.units = 0
quad2inv.type = quad_chromatic
quad2inv.ds = -0.25
quad2inv.k = -100.0
quad2inv.units = 1
quad3.type = quad_chromatic
quad3.ds = 0.25
quad3.k = -29.963938169129919
quad3.units = 0
quad4inv.type = quad_chromatic
quad4inv.ds = -0.25
quad4inv.k = 100.0
quad4inv.units = 1
bend1.type = sbend_exact
bend1.ds = 0.25
bend1.phi = 45.0
bend1.B = 10.484578615324974
bend2inv.type = sbend_exact
bend2inv.ds = -0.25
bend2inv.phi = 45.0
bend3.type = sbend_exact
bend3.ds = 0.25
bend3.phi = 45.0
bend3.B = -10.484578615324974
bend4inv.type = sbend_exact
bend4inv.ds = -0.25
bend4inv.phi = -45.0
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
Analyze
We run the following script to analyze correctness:
Script analysis_charge_sign.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],
[
1.0e-03,
1.0e-03,
0.3,
2.0e-07,
2.0e-07,
6.0e-06,
],
rtol=rtol,
atol=atol,
)
print("")
print("Final Beam:")
sigxf, sigyf, sigtf, emittance_xf, emittance_yf, emittance_tf = get_moments(final)
print(f" sigx={sigxf:e} sigy={sigyf:e} sigt={sigtf:e}")
print(
f" emittance_x={emittance_xf:e} emittance_y={emittance_yf:e} emittance_t={emittance_tf:e}"
)
atol = 0.0 # ignored
rtol = 1.0e-11 # exact to within roundoff tolerance
print(f" rtol={rtol} (ignored: atol~={atol})")
assert np.allclose(
[sigxf, sigyf, sigtf, emittance_xf, emittance_yf, emittance_tf],
[sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
rtol=rtol,
atol=atol,
)