Constant Focusing Channel
Stationary beam in a constant focusing channel (without space charge).
The matched Twiss parameters at entry are:
\(\beta_\mathrm{x} = 1.0\) m
\(\alpha_\mathrm{x} = 0.0\)
\(\beta_\mathrm{y} = 1.0\) m
\(\alpha_\mathrm{y} = 0.0\)
We use a 2 GeV proton beam with initial unnormalized rms emittance of 1 um. The longitudinal beam parameters are chosen so that the bunch has radial symmetry when viewed in the beam rest frame.
The particle distribution should remain unchanged, to within the level expected due to numerical particle noise. This fact is independent of the length of the channel. This is tested using the second moments of the distribution.
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_cfchannel.py
orImpactX executable using an input file:
impactx input_cfchannel.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: Marco Garten, 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 with an initial
# normalized transverse rms emittance of 1 um
kin_energy_MeV = 2.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(938.27208816).set_kin_energy_MeV(kin_energy_MeV)
# particle bunch
distr = distribution.Waterbag(
lambdaX=1.0e-3,
lambdaY=1.0e-3,
lambdaT=3.369701494258956e-4,
lambdaPx=1.0e-3,
lambdaPy=1.0e-3,
lambdaPt=2.9676219145931020e-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.ConstF(ds=2.0, kx=1.0, ky=1.0, kt=1.0),
monitor,
]
)
# run simulation
sim.evolve()
# 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 = proton
beam.distribution = waterbag
beam.lambdaX = 1.0e-3
beam.lambdaY = 1.0e-3
beam.lambdaT = 3.369701494258956e-4
beam.lambdaPx = 1.0e-3
beam.lambdaPy = 1.0e-3
beam.lambdaPt = 2.9676219145931020e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor constf1 monitor
monitor.type = beam_monitor
monitor.backend = h5
constf1.type = constf
constf1.ds = 2.0
constf1.kx = 1.0
constf1.ky = 1.0
constf1.kt = 1.0
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false
Analyze
We run the following script to analyze correctness:
Script analysis_cfchannel.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 # a big number
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-003,
1.0e-003,
3.369701494258956e-4,
1.0e-006,
1.0e-006,
1.0e-006,
],
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 # a big number
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-003,
1.0e-003,
3.369701494258956e-4,
1.0e-006,
1.0e-006,
1.0e-006,
],
rtol=rtol,
atol=atol,
)
Constant Focusing Channel with Space Charge
RMS-matched beam in a constant focusing channel with space charge.
The matched Twiss parameters at entry are:
\(\beta_\mathrm{x} = 1.477305\) m
\(\alpha_\mathrm{x} = 0.0\)
\(\beta_\mathrm{y} = 1.477305\) m
\(\alpha_\mathrm{y} = 0.0\)
We use a 2 GeV proton beam with initial unnormalized rms emittance of 1 um. The longitudinal beam parameters are chosen so that the bunch has radial symmetry when viewed in the beam rest frame. The bunch charge is set to 10 nC, resulting in a transverse tune depression ratio of 0.67. The initial distribution used is a 6D waterbag.
The beam second moments should remain nearly unchanged, except for some small emittance growth due to nonlinear space charge. This is tested using the second moments of the distribution.
In this test, the initial and final values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and :math:`
Run
This example can be run as a Python script (python3 run_cfchannel_10nC.py
) or as an app with an input file (impactx input_cfchannel_10nC.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: Marco Garten, 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.n_cell = [48, 48, 40] # [72, 72, 64] for increased precision
sim.particle_shape = 2 # B-spline order
sim.space_charge = True
sim.prob_relative = [3.0]
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = True
# domain decomposition & space charge mesh
sim.init_grids()
# load a 2 GeV proton beam with an initial
# normalized transverse rms emittance of 1 um
kin_energy_MeV = 2.0e3 # reference energy
bunch_charge_C = 1.0e-8 # used with space charge
npart = 10000 # number of macro particles; use 1e5 for increased precision
# 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.2154443728379865788e-3,
lambdaY=1.2154443728379865788e-3,
lambdaT=4.0956844276541331005e-4,
lambdaPx=8.2274435782286157175e-4,
lambdaPy=8.2274435782286157175e-4,
lambdaPt=2.4415943602685364584e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice
nslice = 50 # use 1e5 for increased precision
# design the accelerator lattice
sim.lattice.extend(
[
monitor,
elements.ConstF(ds=2.0, kx=1.0, ky=1.0, kt=1.0, nslice=nslice),
monitor,
]
)
# run simulation
sim.evolve()
# clean shutdown
sim.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
#beam.npart = 100000 # optional for increased precision
beam.units = static
beam.kin_energy = 2.0e3
beam.charge = 1.0e-8
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.2154443728379865788e-3
beam.lambdaY = 1.2154443728379865788e-3
beam.lambdaT = 4.0956844276541331005e-4
beam.lambdaPx = 8.2274435782286157175e-4
beam.lambdaPy = 8.2274435782286157175e-4
beam.lambdaPt = 2.4415943602685364584e-3
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor constf1 monitor
lattice.nslice = 50
#lattice.nslice = 60 # optional for increased precision
monitor.type = beam_monitor
monitor.backend = h5
constf1.type = constf
constf1.ds = 2.0
constf1.kx = 1.0
constf1.ky = 1.0
constf1.kt = 1.0
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = true
amr.n_cell = 48 48 40
#amr.n_cell = 72 72 64 # optional for increased precision
geometry.prob_relative = 3.0
Analyze
We run the following script to analyze correctness:
Script analysis_cfchannel_10nC.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.2154443728379865788e-003,
1.2154443728379865788e-003,
4.0956844276541331005317e-004,
1.000000000e-006,
1.000000000e-006,
1.000000000e-006,
],
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 = 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.2154443728379865788e-003,
1.2154443728379865788e-003,
4.0956844276541331005317e-004,
1.000000000e-006,
1.000000000e-006,
1.000000000e-006,
],
rtol=rtol,
atol=atol,
)