Soft-edge solenoid
Proton beam propagating through a 6 m region containing a soft-edge solenoid.
The solenoid model used is the default thin-shell model described in: P. Granum et al, “Efficient calculations of magnetic fields of solenoids for simulations,” NIMA 1034, 166706 (2022) DOI:10.1016/j.nima.2022.166706
The solenoid is a cylindrical current sheet with a length of 1 m and a radius of 0.1667 m, corresponding to an aspect ratio diameter/length = 1/3. The peak magnetic field on-axis is 3 T.
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.
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_solenoid_softedge.pyorImpactX executable using an input file:
impactx input_solenoid_softedge.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: Chad Mitchell, Axel Huebl
# 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 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)
# design the accelerator lattice
sol = elements.SoftSolenoid(
ds=6.0,
bscale=1.233482899483985,
cos_coefficients=[
0.350807812299706,
0.323554693720069,
0.260320578919415,
0.182848575294969,
0.106921016050403,
4.409581845710694e-002,
-9.416427163897508e-006,
-2.459452716865687e-002,
-3.272762575737291e-002,
-2.936414401076162e-002,
-1.995780078926890e-002,
-9.102893342953847e-003,
-2.456410658713271e-006,
5.788233017324325e-003,
8.040408292420691e-003,
7.480064552867431e-003,
5.230254569468851e-003,
2.447614547094685e-003,
-1.095525090532255e-006,
-1.614586867387170e-003,
-2.281365457438345e-003,
-2.148709081338292e-003,
-1.522541739363011e-003,
-7.185505862719508e-004,
-6.171194824600157e-007,
4.842109305036943e-004,
6.874508102002901e-004,
6.535550288205728e-004,
4.648795813759210e-004,
2.216564722797528e-004,
-4.100982995210341e-007,
-1.499332112463395e-004,
-2.151538438342482e-004,
-2.044590946652016e-004,
-1.468242784844341e-004,
],
sin_coefficients=[
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
],
mapsteps=200,
nslice=4,
)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
sim.lattice.extend(
[
monitor,
sol,
monitor,
]
)
# run simulation
sim.evolve()
# 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 sol1 monitor
lattice.nslice = 1
monitor.type = beam_monitor
monitor.backend = h5
sol1.type = solenoid_softedge
sol1.ds = 6.0
sol1.bscale = 1.233482899483985
sol1.mapsteps = 800
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = false
Analyze
We run the following script to analyze correctness:
Script analysis_solenoid_softedge.py
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Chad Mitchell, Axel Huebl
# 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 = 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,
)
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],
[
2.425578906459e-3,
2.654015302646e-3,
9.985897906860e-3,
1.365357890e-6,
1.634641555e-6,
1.000000000e-6,
],
rtol=rtol,
atol=atol,
)