Solenoid channel
Proton beam undergoing 90 deg X-Y rotation in an ideal solenoid channel.
The matched Twiss parameters at entry are:
\(\beta_\mathrm{x} = 2.4321374875\) m
\(\alpha_\mathrm{x} = 0.0\)
\(\beta_\mathrm{y} = 2.4321374875\) m
\(\alpha_\mathrm{y} = 0.0\)
We use a 250 MeV proton beam with initial unnormalized rms emittance of 1 micron in the horizontal plane, and 2 micron in the vertical plane.
The solenoid magnetic field corresponds to B = 2 T.
The second moments of the particle distribution after the solenoid channel are rotated by 90 degrees: the final horizontal moments should coincide with the initial vertical moments, and vice-versa, to within the level expected due to noise due to statistical sampling.
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 as a Python script (python3 run_solenoid.py) or with an app with an input file (impactx input_solenoid.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 -*-
import amrex
from impactx import ImpactX, RefPart, 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 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
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_energy_MeV(energy_MeV)
# particle bunch
distr = distribution.Waterbag(
sigmaX=1.559531175539e-3,
sigmaY=2.205510139392e-3,
sigmaT=1.0e-3,
sigmaPx=6.41218345413e-4,
sigmaPy=9.06819680526e-4,
sigmaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)
# design the accelerator lattice
sim.lattice.append(elements.Sol(ds=3.820395, ks=0.8223219329893234))
# run simulation
sim.evolve()
# clean shutdown
del sim
amrex.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.sigmaX = 1.559531175539e-3
beam.sigmaY = 2.205510139392e-3
beam.sigmaT = 1.0e-3
beam.sigmaPx = 6.41218345413e-4
beam.sigmaPy = 9.06819680526e-4
beam.sigmaPt = 1.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = sol1
lattice.nslice = 1
sol1.type = solenoid
sol1.ds = 3.820395
sol1.ks = 0.8223219329893234
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
Analyze
We run the following script to analyze correctness:
Script analysis_solenoid.py
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
import glob
import numpy as np
import pandas as pd
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["x"], moment=2) ** 0.5 # variance -> std dev.
sigpx = moment(beam["px"], moment=2) ** 0.5
sigy = moment(beam["y"], moment=2) ** 0.5
sigpy = moment(beam["py"], moment=2) ** 0.5
sigt = moment(beam["t"], moment=2) ** 0.5
sigpt = moment(beam["pt"], moment=2) ** 0.5
epstrms = beam.cov(ddof=0)
emittance_x = (sigx**2 * sigpx**2 - epstrms["x"]["px"] ** 2) ** 0.5
emittance_y = (sigy**2 * sigpy**2 - epstrms["y"]["py"] ** 2) ** 0.5
emittance_t = (sigt**2 * sigpt**2 - epstrms["t"]["pt"] ** 2) ** 0.5
return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
def read_all_files(file_pattern):
"""Read in all CSV files from each MPI rank (and potentially OpenMP
thread). Concatenate into one Pandas dataframe.
Returns
-------
pandas.DataFrame
"""
return pd.concat(
(
pd.read_csv(filename, delimiter=r"\s+")
for filename in glob.glob(file_pattern)
),
axis=0,
ignore_index=True,
).set_index("id")
# initial/final beam on rank zero
initial = read_all_files("diags/beam_000000.*")
final = read_all_files("diags/beam_final.*")
# 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 = 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,
)
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 = 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],
[
2.205510139392e-3,
1.559531175539e-3,
6.404930308742e-3,
2.0e-6,
1.0e-6,
1.0e-6,
],
rtol=rtol,
atol=atol,
)