Aperture Collimation
Proton beam undergoing collimation by a rectangular boundary aperture.
We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a vertical rms beam size of 2.21 mm.
After a short drift of 0.123, the beam is scraped by a 1 mm x 1.5 mm rectangular aperture.
In this test, the initial values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must agree with nominal values. The test fails if:
any of the final coordinates for the valid (not lost) particles lie outside the aperture boundary or
any of the lost particles are inside the aperture boundary or
if the sum of lost and kept particles is not equal to the initial particles or
if the recorded position \(s\) for the lost particles does not coincide with the drift distance.
Run
This example can be run as a Python script (python3 run_aperture.py) or with an app with an input file (impactx input_aperture.in).
Each can also be prefixed with an MPI executor, such as 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.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
sim.particle_lost_diagnostics_backend = "h5"
# domain decomposition & space charge mesh
sim.init_grids()
# load a 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
kin_energy_MeV = 250.0 # 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(938.27208816).set_kin_energy_MeV(kin_energy_MeV)
# particle bunch
distr = distribution.Waterbag(
lambdaX=1.559531175539e-3,
lambdaY=2.205510139392e-3,
lambdaT=1.0e-3,
lambdaPx=6.41218345413e-4,
lambdaPy=9.06819680526e-4,
lambdaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice
sim.lattice.extend(
[
monitor,
elements.Drift(name="drift", ds=0.123),
elements.Aperture(
name="collimator", aperture_x=1.0e-3, aperture_y=1.5e-3, shape="rectangular"
),
monitor,
]
)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.559531175539e-3
beam.lambdaY = 2.205510139392e-3
beam.lambdaT = 1.0e-3
beam.lambdaPx = 6.41218345413e-4
beam.lambdaPy = 9.06819680526e-4
beam.lambdaPt = 1.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor drift collimator monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
drift.type = drift
drift.ds = 0.123
collimator.type = aperture
collimator.shape = rectangular
collimator.aperture_x = 1.0e-3
collimator.aperture_y = 1.5e-3
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
diag.backend = h5
Analyze
We run the following script to analyze correctness:
Script analysis_aperture.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()
series_lost = io.Series("diags/openPMD/particles_lost.h5", io.Access.read_only)
particles_lost = series_lost.iterations[0].particles["beam"].to_df()
# compare number of particles
num_particles = 10000
assert num_particles == len(initial)
# we lost particles in apertures
assert num_particles > len(final)
assert num_particles == len(particles_lost) + 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.8 * 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.559531175539e-3,
2.205510139392e-3,
1.0e-3,
1.0e-6,
2.0e-6,
1.0e-6,
],
rtol=rtol,
atol=atol,
)
# particle-wise comparison against the rectangular aperture boundary
xmax = 1.0e-3
ymax = 1.5e-3
# kept particles
dx = abs(final["position_x"]) - xmax
dy = abs(final["position_y"]) - ymax
print()
print(f" x_max={final['position_x'].max()}")
print(f" x_min={final['position_x'].min()}")
assert np.less_equal(dx.max(), 0.0)
print(f" y_max={final['position_y'].max()}")
print(f" y_min={final['position_y'].min()}")
assert np.less_equal(dy.max(), 0.0)
# lost particles
dx = abs(particles_lost["position_x"]) - xmax
dy = abs(particles_lost["position_y"]) - ymax
print()
print(f" x_max={particles_lost['position_x'].max()}")
print(f" x_min={particles_lost['position_x'].min()}")
assert np.greater_equal(dx.max(), 0.0)
print(f" y_max={particles_lost['position_y'].max()}")
print(f" y_min={particles_lost['position_y'].min()}")
assert np.greater_equal(dy.max(), 0.0)
# check that s is set correctly
lost_at_s = particles_lost["s_lost"]
drift_s = np.ones_like(lost_at_s) * 0.123
assert np.allclose(lost_at_s, drift_s)
Aperture Collimation with Periodic Masking
Proton beam undergoing collimation by a periodic array of rectangular apertures, such as those used in a pepperpot emittance measurement.
We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a vertical rms beam size of 2.21 mm.
After a short drift of 0.123 m, the beam intercepts a periodic array of apertures of period 1 mm (in the horizontal and vertical), each of which is 0.15 mm x 0.1 mm in size. Every 2nd row of the aperture pattern (rows -3, -1, 1, 3, …) is shifted horizontally by half the horizontal period, creating a hexagonal / triangular pattern.
In this test, the initial values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must agree with nominal values. The test fails if:
any of the final coordinates for the valid (not lost) particles lie outside the aperture boundary or
any of the lost particles lie inside the aperture boundary or
if the sum of the numbers of lost and retained particles is not equal to the number of initial particles
Run
This example can be run as a Python script (python3 run_aperture_periodic.py) or with an app with an input file (impactx input_aperture_periodic.in).
Each can also be prefixed with an MPI executor, such as 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.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
sim.particle_lost_diagnostics_backend = "h5"
# domain decomposition & space charge mesh
sim.init_grids()
# load a 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
kin_energy_MeV = 250.0 # 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.559531175539e-3,
lambdaY=2.205510139392e-3,
lambdaT=1.0e-3,
lambdaPx=6.41218345413e-4,
lambdaPy=9.06819680526e-4,
lambdaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice
sim.lattice.extend(
[
monitor,
elements.Drift(name="drift", ds=0.123),
elements.Aperture(
name="pepperpot",
aperture_x=1.5e-4,
aperture_y=1.0e-4,
repeat_x=1.0e-3,
repeat_y=1.0e-3,
shift_odd_x=True,
shape="rectangular",
),
monitor,
]
)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 100000
beam.units = static
beam.kin_energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.559531175539e-3
beam.lambdaY = 2.205510139392e-3
beam.lambdaT = 1.0e-3
beam.lambdaPx = 6.41218345413e-4
beam.lambdaPy = 9.06819680526e-4
beam.lambdaPt = 1.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor drift pepperpot monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
drift.type = drift
drift.ds = 0.123
pepperpot.type = aperture
pepperpot.shape = rectangular
pepperpot.aperture_x = 1.5e-4
pepperpot.aperture_y = 1.0e-4
pepperpot.repeat_x = 1.0e-3
pepperpot.repeat_y = 1.0e-3
pepperpot.shift_odd_x = true
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
diag.backend = h5
Analyze
We run the following script to analyze correctness:
Script analysis_aperture_periodic.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()
series_lost = io.Series("diags/openPMD/particles_lost.h5", io.Access.read_only)
particles_lost = series_lost.iterations[0].particles["beam"].to_df()
# compare number of particles
num_particles = 100000
assert num_particles == len(initial)
# we lost particles in apertures
assert num_particles > len(final)
assert num_particles == len(particles_lost) + 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.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.559531175539e-3,
2.205510139392e-3,
1.0e-3,
1.0e-6,
2.0e-6,
1.0e-6,
],
rtol=rtol,
atol=atol,
)
# particle-wise comparison against the periodic rectangular aperture boundary
xmax = 1.5e-4
ymax = 1.0e-4
repeat_x = 1.0e-3
repeat_y = 1.0e-3
# kept particles, shifted to the fundamental domain
xshifted = abs(final["position_x"]) + xmax
yshifted = abs(final["position_y"]) + ymax
dx = repeat_x * 0.5
dy = repeat_y * 0.5
xshifted[np.fmod(np.floor((yshifted + dy) / repeat_y), 2) == 1] += dx
u = np.fmod(xshifted, repeat_x) - xmax
v = np.fmod(yshifted, repeat_y) - ymax
# difference from maximum aperture
dx = abs(u) - xmax
dy = abs(v) - ymax
print()
print(f" fundamental x_max={u.max()}")
print(f" fundamental x_min={u.min()}")
assert np.less_equal(dx.max(), 0.0)
print(f" fundamental y_max={v.max()}")
print(f" fundamental y_min={v.min()}")
assert np.less_equal(dy.max(), 0.0)
# lost particles, shifted to the fundamental domain
xshifted = abs(particles_lost["position_x"]) - xmax
yshifted = abs(particles_lost["position_y"]) - ymax
u = np.fmod(xshifted, repeat_x) - xmax
v = np.fmod(yshifted, repeat_y) - ymax
# difference from maximum aperture
dx = abs(u) - xmax
dy = abs(v) - ymax
print()
print(f" fundamental x_max={u.max()}")
print(f" fundamental x_min={u.min()}")
assert np.greater_equal(dx.max(), 0.0)
print(f" fundamental y_max={v.max()}")
print(f" fundamental y_min={v.min()}")
assert np.greater_equal(dy.max(), 0.0)
Collimation Using an Absorber
Proton beam undergoing collimation through partial absorption by a rectangular domain.
This test is the exact negative of the previous test, and illustrates the absorb option of the Aperture element.
We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a vertical rms beam size of 2.21 mm.
After a short drift of 0.123, the beam hits a 1 mm x 1.5 mm rectangular structure, resulting in particle loss.
In this test, the initial values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must agree with nomin> The test fails if:
any of the final coordinates for the valid (not lost) particles lie inside the absorber boundary or
any of the lost particles are outside the absorber boundary or
if the sum of lost and kept particles is not equal to the initial particles or
if the recorded position \(s\) for the lost particles does not coincide with the drift distance.
Run
This example can be run as a Python script (python3 run_absorber.py) or with an app with an input file (impactx input_absorber.in).
Each can also be prefixed with an MPI executor, such as 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.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
sim.particle_lost_diagnostics_backend = "h5"
# domain decomposition & space charge mesh
sim.init_grids()
# load a 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
kin_energy_MeV = 250.0 # 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(938.27208816).set_kin_energy_MeV(kin_energy_MeV)
# particle bunch
distr = distribution.Waterbag(
lambdaX=1.559531175539e-3,
lambdaY=2.205510139392e-3,
lambdaT=1.0e-3,
lambdaPx=6.41218345413e-4,
lambdaPy=9.06819680526e-4,
lambdaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice
sim.lattice.extend(
[
monitor,
elements.Drift(name="drift", ds=0.123),
elements.Aperture(
name="collimator",
aperture_x=1.0e-3,
aperture_y=1.5e-3,
shape="rectangular",
action="absorb",
),
monitor,
]
)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.559531175539e-3
beam.lambdaY = 2.205510139392e-3
beam.lambdaT = 1.0e-3
beam.lambdaPx = 6.41218345413e-4
beam.lambdaPy = 9.06819680526e-4
beam.lambdaPt = 1.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor drift collimator monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
drift.type = drift
drift.ds = 0.123
collimator.type = aperture
collimator.shape = rectangular
collimator.aperture_x = 1.0e-3
collimator.aperture_y = 1.5e-3
collimator.action = absorb
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
diag.backend = h5
Analyze
We run the following script to analyze correctness:
Script analysis_absorber.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()
series_lost = io.Series("diags/openPMD/particles_lost.h5", io.Access.read_only)
particles_lost = series_lost.iterations[0].particles["beam"].to_df()
# compare number of particles
num_particles = 10000
assert num_particles == len(initial)
# we lost particles in apertures
assert num_particles > len(final)
assert num_particles == len(particles_lost) + 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.8 * 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.559531175539e-3,
2.205510139392e-3,
1.0e-3,
1.0e-6,
2.0e-6,
1.0e-6,
],
rtol=rtol,
atol=atol,
)
# particle-wise comparison against the rectangular aperture boundary
xmax = 1.0e-3
ymax = 1.5e-3
# kept particles
dx = abs(final["position_x"]) - xmax
dy = abs(final["position_y"]) - ymax
print()
print(f" x_max={final['position_x'].max()}")
print(f" x_min={final['position_x'].min()}")
assert np.greater_equal(dx.max(), 0.0)
print(f" y_max={final['position_y'].max()}")
print(f" y_min={final['position_y'].min()}")
assert np.greater_equal(dy.max(), 0.0)
# lost particles
dx = abs(particles_lost["position_x"]) - xmax
dy = abs(particles_lost["position_y"]) - ymax
print()
print(f" x_max={particles_lost['position_x'].max()}")
print(f" x_min={particles_lost['position_x'].min()}")
assert np.less_equal(dx.max(), 0.0)
print(f" y_max={particles_lost['position_y'].max()}")
print(f" y_min={particles_lost['position_y'].min()}")
assert np.less_equal(dy.max(), 0.0)
# check that s is set correctly
lost_at_s = particles_lost["s_lost"]
drift_s = np.ones_like(lost_at_s) * 0.123
assert np.allclose(lost_at_s, drift_s)
Aperture Collimation for a Thick Element
Proton beam in a drift, undergoing collimation by a rectangular boundary aperture.
We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a vertical rms beam size of 2.21 mm.
The beam is scraped by a 1 mm x 1.5 mm rectangular aperture. For this test, the parameter nslice = 1, so application at the aperture boundary is equivalent to using a thin aperture element located at the exit of the drift.
In this test, the initial values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must> The test fails if:
any of the final coordinates for the valid (not lost) particles lie outside the aperture boundary or
any of the lost particles are inside the aperture boundary or
if the sum of lost and kept particles is not equal to the initial particles or
if the recorded position \(s\) for the lost particles does not coincide with the drift distance.
Run
This example can be run as a Python script (python3 run_aperture_thick.py) or with an app with an input file (impactx input_aperture_thick.in).
Each can also be prefixed with an MPI executor, such as 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.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
sim.particle_lost_diagnostics_backend = "h5"
# domain decomposition & space charge mesh
sim.init_grids()
# load a 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
kin_energy_MeV = 250.0 # 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(938.27208816).set_kin_energy_MeV(kin_energy_MeV)
# particle bunch
distr = distribution.Waterbag(
lambdaX=1.559531175539e-3,
lambdaY=2.205510139392e-3,
lambdaT=1.0e-3,
lambdaPx=6.41218345413e-4,
lambdaPy=9.06819680526e-4,
lambdaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice
sim.lattice.extend(
[
monitor,
elements.Drift(name="drift", ds=0.123, aperture_x=1.0e-3, aperture_y=1.5e-3),
monitor,
]
)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.559531175539e-3
beam.lambdaY = 2.205510139392e-3
beam.lambdaT = 1.0e-3
beam.lambdaPx = 6.41218345413e-4
beam.lambdaPy = 9.06819680526e-4
beam.lambdaPt = 1.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor drift monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
drift.type = drift
drift.ds = 0.123
drift.aperture_x = 1.0e-3
drift.aperture_y = 1.5e-3
###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
diag.backend = h5
Analyze
We run the following script to analyze correctness:
Script analysis_aperture_thick.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()
series_lost = io.Series("diags/openPMD/particles_lost.h5", io.Access.read_only)
particles_lost = series_lost.iterations[0].particles["beam"].to_df()
# compare number of particles
num_particles = 10000
assert num_particles == len(initial)
# we lost particles in apertures
assert num_particles > len(final)
assert num_particles == len(particles_lost) + 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.8 * 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.559531175539e-3,
2.205510139392e-3,
1.0e-3,
1.0e-6,
2.0e-6,
1.0e-6,
],
rtol=rtol,
atol=atol,
)
# particle-wise comparison against the rectangular aperture boundary
xmax = 1.0e-3
ymax = 1.5e-3
# kept particles
dx = abs(final["position_x"]) - xmax
dy = abs(final["position_y"]) - ymax
print()
print(f" x_max={final['position_x'].max()}")
print(f" x_min={final['position_x'].min()}")
assert np.less_equal(dx.max(), 0.0)
print(f" y_max={final['position_y'].max()}")
print(f" y_min={final['position_y'].min()}")
assert np.less_equal(dy.max(), 0.0)
# lost particles
dx = abs(particles_lost["position_x"]) - xmax
dy = abs(particles_lost["position_y"]) - ymax
print()
print(f" x_max={particles_lost['position_x'].max()}")
print(f" x_min={particles_lost['position_x'].min()}")
assert np.greater_equal(dx.max(), 0.0)
print(f" y_max={particles_lost['position_y'].max()}")
print(f" y_min={particles_lost['position_y'].min()}")
assert np.greater_equal(dy.max(), 0.0)
# check that s is set correctly
lost_at_s = particles_lost["s_lost"]
drift_s = np.ones_like(lost_at_s) * 0.123
assert np.allclose(lost_at_s, drift_s)