FODO Cell with RF
Stable FODO cell with short RF (buncher) cavities added for longitudinal focusing. The phase advance in all three phase planes is between 86-89 degrees.
The matched Twiss parameters at entry are:
\(\beta_\mathrm{x} = 9.80910407\) m
\(\alpha_\mathrm{x} = 0.0\)
\(\beta_\mathrm{y} = 1.31893788\) m
\(\alpha_\mathrm{y} = 0.0\)
\(\beta_\mathrm{t} = 4.6652668782\) m
\(\alpha_\mathrm{t} = 0.0\)
We use a 250 MeV proton beam with initial unnormalized rms emittance of 1 mm-mrad in all three phase planes.
The second moments of the particle distribution after the FODO cell should coincide with the second moments of the particle distribution before the FODO cell, 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_fodo_rf.py) or with an app with an input file (impactx input_fodo_rf.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, 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
# unnormalized rms emittance of 1 mm-mrad in all
# three phase planes
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=3.131948925200e-3,
sigmaY=1.148450209423e-3,
sigmaT=2.159922887089e-3,
sigmaPx=3.192900088357e-4,
sigmaPy=8.707386631090e-4,
sigmaPt=4.62979491526e-4,
)
sim.add_particles(bunch_charge_C, distr, npart)
# design the accelerator lattice
# Quad elements
quad1 = elements.Quad(ds=0.15, k=2.5)
quad2 = elements.Quad(ds=0.3, k=-2.5)
# Drift element
drift1 = elements.Drift(ds=1.0)
# Short RF cavity element
shortrf1 = elements.ShortRF(V=0.01, k=15.0)
lattice_no_drifts = [quad1, shortrf1, quad2, shortrf1, quad1]
# set first lattice element
sim.lattice.append(lattice_no_drifts[0])
# intersperse all remaining elements of the lattice with a drift element
for element in lattice_no_drifts[1:]:
sim.lattice.extend([drift1, element])
# 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 = 3.131948925200e-3
beam.sigmaY = 1.148450209423e-3
beam.sigmaT = 2.159922887089e-3
beam.sigmaPx = 3.192900088357e-4
beam.sigmaPy = 8.707386631090e-4
beam.sigmaPt = 4.62979491526e-4
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = quad1 drift1 shortrf1 drift1 quad2 drift1 shortrf1 drift1 quad1
quad1.type = quad
quad1.ds = 0.15
quad1.k = 2.5
drift1.type = drift
drift1.ds = 1.0
shortrf1.type = shortrf
shortrf1.V = 0.01
shortrf1.k = 15.0
quad2.type = quad
quad2.ds = 0.3
quad2.k = -2.5
###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false
Analyze
We run the following script to analyze correctness:
Script analysis_fodo_rf.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],
[
3.145694e-03,
1.153344e-03,
2.155082e-03,
9.979770e-07,
1.008751e-06,
1.000691e-06,
],
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],
[
3.112318e-03,
1.153322e-03,
2.166501e-03,
9.979770e-07,
1.008751e-06,
1.000691e-06,
],
rtol=rtol,
atol=atol,
)