Bending Dipole with Zero Field Strength
This example tests the effect of a bending dipole with zero field strength, as represented using the elements Sbend and ExactSbend.
In the limiting case of zero field strength, the effect of the dipole should be identical to that of a drift of the same length.
To test this, the map for a dipole is applied, and this is followed by applying the inverse map for a drift of the corresponding length.
We use a 2 GeV electron beam, identical to the distribution used in the test “examples-rotation”.
The second moments of x, y, and t should 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_zero_bend_inverse.pyorImpactX executable using an input file:
impactx input_zero_bend_inverse.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: Ryan Sandberg, 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.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 nm
kin_energy_MeV = 2.0e3 # reference energy
bunch_charge_C = 1.0e-9 # used without 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=4.0e-3,
lambdaY=4.0e-3,
lambdaT=1.0e-3,
lambdaPx=3.0e-4,
lambdaPy=3.0e-4,
lambdaPt=2.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice
ns = 24
length1 = 10.0
length2 = 5.0
linear_segment = [
monitor,
elements.Sbend(name="bend1", ds=length1, rc=0.0, nslice=ns),
elements.Drift(name="inv_drift1", ds=-length1, nslice=ns),
]
nonlinear_segment = [
elements.ExactSbend(name="bend2", ds=length2, phi=0.0, nslice=ns),
elements.ExactDrift(name="inv_drift2", ds=-length2, nslice=ns),
monitor,
]
# assign a lattice segment
sim.lattice.extend(linear_segment)
sim.lattice.extend(nonlinear_segment)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 2.0e3
beam.charge = 1.0e-9
beam.particle = electron
beam.distribution = waterbag
beam.lambdaX = 4.0e-3
beam.lambdaY = beam.lambdaX
beam.lambdaT = 1.0e-3
beam.lambdaPx = 3.0e-4
beam.lambdaPy = beam.lambdaPx
beam.lambdaPt = 2.0e-3
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor bend1 inv_drift1 bend2 inv_drift2 monitor
lattice.nslice = 24
monitor.type = beam_monitor
monitor.backend = h5
bend1.type = sbend # a zero-field "sbend" is equivalent to a "drift"
bend1.ds = 10.0
bend1.rc = 0.0
bend2.type = sbend_exact # a zero-field "sbend_exact" is equivalent to a "drift_exact"
bend2.ds = 5.0
bend2.phi = 0.0
inv_drift1.type = drift
inv_drift1.ds = -bend1.ds
inv_drift2.type = drift_exact
inv_drift2.ds = -bend2.ds
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
Analyze
We run the following script to analyze correctness:
Script analysis_zero_bend_inverse.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],
[
4.0e-03,
4.0e-03,
1.0e-03,
1.2e-06,
1.2e-06,
2.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-12 # 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,
)