Positron Channel
Acceleration of a positron beam with large (10%) energy spread, from 10 GeV to 2.5 TeV, for parameters based on possible staging of a laser-wakefield accelerator.
The lattice consists of 250 periods, each consisting of a quadrupole triplet followed by 10 GeV energy gain in a uniform field.
We use a 190 pC positron beam with initial normalized rms emittance of 10 nm, rms beam size of 5 microns, and a triangular current pulse with end-to-end pulse length of 0.12 ps.
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_positron.py
orImpactX executable using an input file:
impactx input_positron.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 = False
# domain decomposition & space charge mesh
sim.init_grids()
# load a 10 GeV positron beam with an initial
# normalized rms emittance of 10 nm
kin_energy_MeV = 10.0e3 # reference energy
bunch_charge_C = 190.0e-12 # 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.Triangle(
lambdaX=5.054566450e-6,
lambdaY=5.054566450e-6,
lambdaT=8.43732950e-7,
lambdaPx=1.01091329e-7,
lambdaPy=1.01091329e-7,
lambdaPt=1.0e-2,
muxpx=0.0,
muypy=0.0,
mutpt=0.995037190209989,
)
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
period = [
monitor,
elements.ChrQuad(ds=0.1, k=-6.674941, units=1, nslice=ns),
elements.ChrDrift(ds=0.3, nslice=ns),
elements.ChrQuad(ds=0.2, k=6.674941, units=1, nslice=ns),
elements.ChrDrift(ds=0.3, nslice=ns),
elements.ChrQuad(ds=0.1, k=-6.674941, units=1, nslice=ns),
elements.ChrDrift(ds=0.1, nslice=ns),
elements.ChrAcc(ds=1.8, ez=10871.950994502130424, bz=1.0e-12, nslice=ns),
elements.ChrDrift(ds=0.1, nslice=ns),
monitor,
]
sim.lattice.extend(period)
# number of periods through the lattice
sim.periods = 250
# run simulation
sim.evolve()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 10.0e3
beam.charge = 190.0e-12
beam.particle = positron
beam.distribution = triangle
beam.lambdaX = 5.054566450e-6
beam.lambdaY = 5.054566450e-6
beam.lambdaT = 8.43732950e-7
beam.lambdaPx = 1.01091329e-7
beam.lambdaPy = 1.01091329e-7
beam.lambdaPt = 1.0e-2
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.995037190209989
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.periods = 250
lattice.elements = monitor quad1 drift1 quad2 drift1 quad1 drift2 unifacc drift2 monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
quad1.type = quad_chromatic
quad1.ds = 0.1
quad1.k = -6.674941
quad1.units = 1
drift1.type = drift_chromatic
drift1.ds = 0.3
quad2.type = quad_chromatic
quad2.ds = 0.2
quad2.k = 6.674941
quad2.units = 1
drift2.type = drift_chromatic
drift2.ds = 0.1
unifacc.type = uniform_acc_chromatic
unifacc.ds = 1.8
unifacc.ez = 10871.950994502130424
unifacc.bz = 1.0e-12
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false
#amr.n_cell = 56 56 48
#geometry.prob_relative = 1.0
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = false
Analyze
We run the following script to analyze correctness:
Script analysis_positron.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 = 2.1 * 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],
[
5.05456645029333603e-006,
5.05456645029333603e-006,
8.47941120001532497e-006,
5.1097284000861952e-13,
5.1097284000861952e-13,
8.47941120001532497e-008,
],
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.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],
[
5.26e-006,
5.26e-006,
8.47941120001532497e-006,
2.0439953822082447e-15,
2.0439953822082447e-15,
3.3919371010764148e-10,
],
rtol=rtol,
atol=atol,
)