Parallelized Vlasov-Fokker-Planck solver for desktop personal computers

The numerical solution of the Vlasov-Fokker-Planck equation is a well established method to simulate the dynamics, including the self-interaction with its own wake field, of an electron bunch in a storage ring. In this paper we present Inovesa, a modularly extensible program that uses opencl to massively parallelize the computation. It allows a standard desktop PC to work with appropriate accuracy and yield reliable results within minutes. We provide numerical stability-studies over a wide parameter range and compare our numerical findings to known results. Simulation results for the case of coherent synchrotron radiation will be compared to measurements that probe the effects of the microbunching instability occurring in the short bunch operation at ANKA. It will be shown that the impedance model based on the shielding effect of two parallel plates can not only describe the instability threshold, but also the presence of multiple regimes that show differences in the emission of coherent synchrotron radiation.

Figure 1

Unshielded (free space) CSR impedance (FS) and CSR impedance shielded by parallel plates (PP), calculated for the case of an accelerator with $f_{\text{rev}}=8.582\mathrm{GHz}$, and (for the shielded case) $g=32\mathrm{mm}$. For high frequencies both impedances converge, as short wavelength are not effected by the shielding. For $f\to 0\mathrm{Hz}$ both impedances approach $Z=0\mathrm{\Omega }$.

Figure 2

Computational time needed by an Intel Core i5 4258U for 1000 cubically interpolated rotation steps using different implementations. In this test case, the grid has 512 points per axis, no optimizations (besides SM) were used. The first bar represents the computational time a direct implementation of the rotation takes. The second and third bar show the time the implementation using the source map formalism take—running sequentially on the CPU or parallel on the integrated graphics processor.

Figure 3

Exponential damping of the energy spread ${\sigma }_E$ as a function of time for different data types. There is an initial jitter due to quantization noise in the test pattern that has been read in from a 16 bit PNG file and due to the fact that the initial distribution is not fully covered by the grid. Even with this problematic starting conditions, after $T=60{\mathrm{T}}_{\mathrm{s}}$, all data types have converged to ${\sigma }_p\approx 1$. Note that all simulation curves are almost perfectly overlapping, there is no noticeable difference coming from the used data type.

Figure 4

Results of simulation runs with different grid sizes and numbers of time steps to simulate one synchrotron period. Every run is represented by two tiles: The reconstructed damping time (${\tau }_{d,s}$) is shown in the upper row, the natural energy spread (${\sigma }_{E,0,s}$) can be found in the lower row. For the left two columns, the set damping time is ${\tau }_d=50T_s$ ($\beta =0.02$), for the others, it is ${\tau }_d=500T_s$ ($\beta =0.002$). The value corresponding to the analytic result is white, higher values are shown in red, lower ones in blue. It can be seen that slow damping is numerically more challenging—it requires bigger time steps (less steps per $T_s$) for good convergence. Independently from the numerical settings, Manhattan rotation shows to be more robust in reconstructing the set values.

Figure 5

Evolution of the energy spread over time (left) and bunch profiles of the simulation run using quadratic interpolation (right) at the points in time marked by the disks. For the case where cubic interpolation is used (left, solid black line), the energy spread stays at ${\sigma }_E\approx 1.01{\sigma }_{E,0}$. When using quadratic interpolation (left, dashed blue line), an increase in energy spread can be observed at $T=1T_s$. The right-hand side reveals that this increase is a numerical artifact. The bunch profiles computed during the initial increase of the energy spread show large ripples with a period length of two grid cells. The earlier and later profiles do not show such structures. This implies that the increase of energy spread is driven by a numerical instability.

Figure 6

Example for a simulated (left) and a measured (right) bursting spectrogram. For the simulated spectrogram the axis is scaled with a factor of $2\pi R/C$ to correct the mismatch due to the isomagnetic approximation. There are small differences, e.g. in the threshold current and in the frequencies. However, keeping in mind the very simple model, the general structure of the spectrograms matches quite well: There is an isolated finger pointing down (at $f\approx 3\mathrm{kHz}$); the fingertip ($I\approx 0.21\mathrm{mA}$) marks the instability threshold. For slightly higher currents, fluctuations in the lower frequency range ($f<10\mathrm{kHz}$) start and the finger broadens. For the highest currents displayed here, there is a regime showing parallel frequency lines that stay approximately constant with changing current.

Figure 7

Evolution of the RMS energy spread over a time of 500 synchrotron periods for $I=1.5\mathrm{mA}$. The width of the line is caused by high frequency oscillations. For the starting distributions the Haïssinski distribution (equilibrium at $I=675\mu \mathrm{A}$), a Gaussian distribution with ${\sigma }_p={\sigma }_q=2.25$, and the final distribution of the previous run (“adiabatic,” $I=1.54\mathrm{mA}$) were used. Note that the Haïssinski distribution is immediately blown up and becomes larger than the Gaussian distribution, and that the high frequency oscillation is systematically higher until $T\approx 300{\mathrm{T}}_{\mathrm{s}}$.