Kurth Distribution in a Periodic Focusing Channel

Matched Kurth distribution in a periodic focusing channel (without space charge).

The distribution is radially symmetric in (x,y,t) space, and matched to a radially symmetric periodic linear focusing lattice with a phase advance of 121 degrees.

We use a 2 GeV proton beam with initial unnormalized rms emittance of 1 um in all three phase planes.

The particle distribution should remain unchanged, to within the level expected due to numerical particle noise. This is tested using the second moments of the distribution.

In this test, the initial and final values of \(\lambda_x\), \(\lambda_y\), \(\lambda_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_kurth_periodic.py or

  • ImpactX executable using an input file: impactx input_kurth_periodic.in

For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Listing 20 You can copy this file from examples/kurth/run_kurth_periodic.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Ryan Sandberg, 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
# unnormalized rms emittance of 1 um in each
# coordinate plane
kin_energy_MeV = 2.0e3  # reference energy
bunch_charge_C = 1.0e-8  # 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.Kurth6D(
    lambdaX=1.11e-3,
    lambdaY=1.11e-3,
    lambdaT=3.74036839224568e-4,
    lambdaPx=9.00900900901e-4,
    lambdaPy=9.00900900901e-4,
    lambdaPt=2.6735334467940146e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)

# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")

# design the accelerator lattice
constf1 = elements.ConstF(ds=2.0, kx=0.7, ky=0.7, kt=0.7)
drift1 = elements.Drift(ds=1.0)
sim.lattice.extend([monitor, drift1, constf1, drift1, monitor])

# run simulation
sim.evolve()

# clean shutdown
sim.finalize()
Listing 21 You can copy this file from examples/kurth/input_kurth_periodic.in.
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 2.0e3
beam.charge = 1.0e-8
beam.particle = proton
beam.distribution = kurth6d
beam.lambdaX = 1.11e-3
beam.lambdaY = 1.11e-3
beam.lambdaT = 3.74036839224568e-4
beam.lambdaPx = 9.00900900901e-4
beam.lambdaPy = 9.00900900901e-4
beam.lambdaPt = 2.6735334467940146e-3


###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor drift1 constf1 drift1 monitor

monitor.type = beam_monitor
monitor.backend = h5

drift1.type = drift
drift1.ds = 1.0

constf1.type = constf
constf1.ds = 2.0
constf1.kx = 0.7
constf1.ky = 0.7
constf1.kt = 0.7


###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false

amr.n_cell = 40 40 32

Analyze

We run the following script to analyze correctness:

Script analysis_kurth_periodic.py
Listing 22 You can copy this file from examples/kurth/analysis_kurth_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()

# 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.7 * 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.11e-03,
        1.11e-03,
        3.74036839224568e-04,
        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.7 * 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.11e-03,
        1.11e-03,
        3.74036839224568e-04,
        1.000000000e-006,
        1.000000000e-006,
        1.000000000e-006,
    ],
    rtol=rtol,
    atol=atol,
)

Kurth Distribution in a Periodic Focusing Channel with Space Charge

Matched Kurth distribution in a periodic focusing channel with space charge.

The distribution is radially symmetric in (x,y,t) space, and matched to a radially symmetric constant linear focusing.

We use a 2 GeV proton beam with initial unnormalized rms emittance of 1 um in all three phase planes. The bunch charge is set to 10 nC, depressing the phase advance from 121 degrees to 74 degrees.

The particle distribution should remain unchanged, to within the level expected due to numerical particle noise. This is tested using the second moments of the distribution.

In this test, the initial and final values of \(\lambda_x\), \(\lambda_y\), \(\lambda_t\), \(\epsilon_x\), \(\epsilon_y\), and :math:`

Run

This example can be run as a Python script (python3 run_kurth_10nC_periodic.py) or as an app with an input file (impactx input_kurth_10nC_periodic.in). Each can also be prefixed with an MPI executor, such as mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Listing 23 You can copy this file from examples/kurth/run_kurth_10nC_periodic.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Ryan Sandberg, 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]  # use [72, 72, 72] for increased precision
sim.particle_shape = 2  # B-spline order
sim.space_charge = True
# 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
# unnormalized rms emittance of 1 um in each
# coordinate plane
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.Kurth6D(
    lambdaX=1.46e-3,
    lambdaY=1.46e-3,
    lambdaT=4.9197638312420749e-4,
    lambdaPx=6.84931506849e-4,
    lambdaPy=6.84931506849e-4,
    lambdaPt=2.0326178944803812e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)

# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")

# design the accelerator lattice
nslice = 20  # use 30 for increased precision
constf1 = elements.ConstF(ds=2.0, kx=0.7, ky=0.7, kt=0.7, nslice=nslice)
drift1 = elements.Drift(ds=1.0, nslice=nslice)
sim.lattice.extend([monitor, drift1, constf1, drift1, monitor])

# run simulation
sim.evolve()

# clean shutdown
sim.finalize()
Listing 24 You can copy this file from examples/kurth/input_kurth_10nC_periodic.in.
###############################################################################
# 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 = kurth6d
beam.lambdaX = 1.46e-3
beam.lambdaY = 1.46e-3
beam.lambdaT = 4.9197638312420749e-4
beam.lambdaPx = 6.84931506849e-4
beam.lambdaPy = 6.84931506849e-4
beam.lambdaPt = 2.0326178944803812e-3


###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor drift1 constf1 drift1 monitor
lattice.nslice = 20
#lattice.nslice = 30 #optional for increased precision

monitor.type = beam_monitor
monitor.backend = h5

drift1.type = drift
drift1.ds = 1.0

constf1.type = constf
constf1.ds = 2.0
constf1.kx = 0.7
constf1.ky = 0.7
constf1.kt = 0.7


###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = true

amr.n_cell = 48 48 40
#amr.n_cell = 72 72 72  # optional for increased precision
geometry.prob_relative = 3.0

Analyze

We run the following script to analyze correctness:

Script analysis_kurth_10nC_periodic.py
Listing 25 You can copy this file from examples/kurth/analysis_kurth_10nC_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()

# 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.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.46e-3,
        1.46e-3,
        4.9197638312420749e-4,
        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 = 2.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.46e-3,
        1.46e-3,
        4.9197638312420749e-4,
        1.000000000e-006,
        1.000000000e-006,
        1.000000000e-006,
    ],
    rtol=rtol,
    atol=atol,
)