Error scaling test for nonlinear dipole fringe field map
This benchmark tests the use of the nonlinear DipEdge model for integrating through a dipole fringe field.
Six distinct initial conditions are tested for a nominal proton beam with 800 MeV kinetic energy. The values of the field integrals K0-K6 are set to the default values, as used in:
Hwang and S. Y. Lee, “Dipole fringe field map for compact synchrotrons,” Phys. Rev. Accel. Beams 18, 122401 (2015)
The initial conditions are chosen with increasing distance from the origin in phase space. The value of the Lie generator is a dynamical invariant of the ideal fringe field map. In reality, there is an error in the final variables that scales with (x,px,y,py,t,pt,g) as degree 3. As a result, the Lie generator is not an exact invariant of the numerically-computed fringe field map.
In this test, the change in the Lie generator for each initial condition should coincide with its (small) nominal value.
Run
This example can be run either as:
Python script:
python3 run_dipedge.py
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: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-
import math
import pandas as pd
import amrex.space3d as amr
from impactx import Config, ImpactX, elements
sim = ImpactX()
# set numerical parameters and IO control
sim.space_charge = False
sim.slice_step_diagnostics = True
# domain decomposition & space charge mesh
sim.init_grids()
# load initial beam parameters
kin_energy_MeV = 0.8e3 # 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)
qm_eev = 1.0 / 938.27208816 / 1e6 # electron charge/mass in e / eV
pc = sim.particle_container()
df_initial = pd.read_csv("./initial_coords.csv", sep=" ")
dx = df_initial["x"].to_numpy()
dpx = df_initial["px"].to_numpy()
dy = df_initial["y"].to_numpy()
dpy = df_initial["py"].to_numpy()
dt = df_initial["t"].to_numpy()
dpt = df_initial["pt"].to_numpy()
# dg = df_initial["gap"].to_numpy()
if not Config.have_gpu: # initialize using cpu-based PODVectors
dx_podv = amr.PODVector_real_std()
dy_podv = amr.PODVector_real_std()
dt_podv = amr.PODVector_real_std()
dpx_podv = amr.PODVector_real_std()
dpy_podv = amr.PODVector_real_std()
dpt_podv = amr.PODVector_real_std()
else: # initialize on device using arena/gpu-based PODVectors
dx_podv = amr.PODVector_real_arena()
dy_podv = amr.PODVector_real_arena()
dt_podv = amr.PODVector_real_arena()
dpx_podv = amr.PODVector_real_arena()
dpy_podv = amr.PODVector_real_arena()
dpt_podv = amr.PODVector_real_arena()
for p_dx in dx:
dx_podv.push_back(p_dx)
for p_dy in dy:
dy_podv.push_back(p_dy)
for p_dt in dt:
dt_podv.push_back(p_dt)
for p_dpx in dpx:
dpx_podv.push_back(p_dpx)
for p_dpy in dpy:
dpy_podv.push_back(p_dpy)
for p_dpt in dpt:
dpt_podv.push_back(p_dpt)
# print(dt)
print(dpt)
# print(dg)
pc.add_n_particles(
dx_podv, dy_podv, dt_podv, dpx_podv, dpy_podv, dpt_podv, qm_eev, bunch_charge_C
)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice)
ns = 1 # number of slices per ds in the element
edge_angle = math.pi / 8.0
dipedge0 = elements.DipEdge(
name="dipedge0", psi=edge_angle, rc=10.0, g=1.0e-3, model="linear", location="entry"
)
dipedge1 = elements.DipEdge(
name="dipedge1",
psi=edge_angle,
rc=10.0,
g=1.0e-3,
model="nonlinear",
location="entry",
)
line = [
monitor,
# dipedge0,
dipedge1,
monitor,
]
# assign a segment
sim.lattice.extend(line)
# run simulation
sim.track_particles()
# clean shutdown
sim.finalize()
Analyze
We run the following script to analyze correctness:
Script analysis_dipedge.py
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
import math
import numpy as np
import openpmd_api as io
# 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()
beam_final = series.iterations[last_step].particles["beam"]
final = beam_final.to_df()
# Basic input parameters
g = 1.0e-3
phi = math.pi / 8.0
rc = 10.0
R = 1.0
K0 = math.pi**2 / 6.0
K3 = 1.0 / 6.0
Kar = [K0, 0, 0, K3, 0, 0, 0]
delta = 0.0
# Derived quantities
cs = math.cos(phi)
sn = math.sin(phi)
tn = sn / cs
sc = 1.0 / cs
# Lie generator coefficients
c1 = g * Kar[1] / (rc * cs)
c2 = sn * g**2 * Kar[0] / rc**2 * 1.0 / (2.0 * cs**3 * (1 + delta))
c3 = g**2 / rc * Kar[0] / (cs**2 * (1 + delta))
c4 = 1 / (1 + delta) * g / rc * Kar[1] * sn / cs**2
c5 = sn / cs * 1.0 / (2 * rc)
c6 = g * Kar[1] / rc * sn**2 / (4 * rc * (1 + delta) * cs**3)
c7 = (
1
/ (2 * cs**3 * (1 + delta))
* (g * Kar[1] / (2 * rc**2) + (1 + sn**2) * g / rc**2 * Kar[2])
)
c8 = 1 / 6 * tn**3 / (2 * rc**2 * (1 + delta))
c9 = 1 / 2 * (tn * sc**2 / (2 * rc**2 * (1 + delta)))
c10 = 1 / (2 * (1 + delta)) * tn**2 / rc
c11 = 1 / (2 * rc * (1 + delta))
c12 = 1 / 24 * (4 / cs - 8 / cs**3) * Kar[3] / (rc**2 * g * (1 + delta))
c13 = sn**2 / (2 * cs**3) * g**2 / (rc * R) * Kar[4]
c14 = 1 / 2 * sn / cs**3 * g / (rc * R) * Kar[5]
c15 = Kar[6] / (rc * R) * 1 / cs**3
xi = initial["position_x"]
pxi = initial["momentum_x"]
yi = initial["position_y"]
pyi = initial["momentum_y"]
ti = initial["position_t"]
pti = initial["momentum_t"]
Omega_initial = (
xi * c1
- xi * c2
+ pxi * c3
+ (xi * pxi - yi * pyi) * c4
+ (xi**2 - yi**2) * c5
- xi**2 * c6
+ yi**2 * c7
- xi**3 * c8
+ xi * yi**2 * c9
+ (xi**2 * pxi - yi**2 * pxi - 2 * xi * yi * pyi) * c10
- yi**2 * pxi * c11
+ yi**4 * c12
+ xi * c13
+ (yi**2 - xi**2) * c14
+ (xi * yi**2 / 2 - xi**3 / 6) * c15
)
xf = final["position_x"]
pxf = final["momentum_x"]
yf = final["position_y"]
pyf = final["momentum_y"]
tf = final["position_t"]
ptf = final["momentum_t"]
Omega_final = (
xf * c1
- xf * c2
+ pxf * c3
+ (xf * pxf - yf * pyf) * c4
+ (xf**2 - yf**2) * c5
- xf**2 * c6
+ yf**2 * c7
- xf**3 * c8
+ xf * yf**2 * c9
+ (xf**2 * pxf - yf**2 * pxf - 2 * xf * yf * pyf) * c10
- yf**2 * pxf * c11
+ yf**4 * c12
+ xf * c13
+ (yf**2 - xf**2) * c14
+ (xf * yf**2 / 2 - xf**3 / 6) * c15
)
Delta_Omega = (Omega_final - Omega_initial).abs()
dx = (xf - xi).abs().max()
dpx = (pxf - pxi).abs().max()
dy = (yf - yi).abs().max()
dpy = (pyf - pyi).abs().max()
dt = (tf - ti).abs().max()
dpt = (ptf - pti).abs().max()
print("Change in the coordinates and momenta:")
print("dx", dx)
print("dpx", dpx)
print("dy", dy)
print("dpy", dpy)
print("dt", dt)
print("dpt", dpt)
print("Change in the Lie generator, for each initial condition:")
print(Delta_Omega)
atol = 1.5e-10
print(f" atol={atol}")
assert np.allclose(
[Delta_Omega.max()],
[
0.0,
],
atol=atol,
)