Summary of ACCSIM and ORBIT benchmarking simulations

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH CERN ⎯ BEAMS DEPARTMENT BE-Note-2009-011 Summary of ACCSIM and ORBIT benchmarking simulations Masamitsu AIBA Abstract We have performed a benchmarking study of ORBIT and ACCSIM which are accelerator tracking codes having routines to evaluate space charge effects. The study is motivated by the need of predicting/understanding beam behaviour in the CERN Proton Synchrotron Booster (PSB) in which direct space charge is expected to be the dominant performance limitation. Historically at CERN, ACCSIM has been employed for space charge simulation studies. A benchmark study using ORBIT has been started to confirm the results from ACCSIM and to profit from the advantages of ORBIT such as the capability of parallel processing. We observed a fair agreement in emittance evolution in the horizontal plane but not in the vertical one. This may be partly due to the fact that the algorithm to compute the space charge field is different between the two codes. Geneva, Switzerland February, 2009 1. Introduction We have performed a benchmarking study of ORBIT and ACCSIM which are accelerator tracking codes having routines to evaluate space charge effects. The study is motivated by the need of predicting/understanding beam behaviour in the CERN Proton Synchrotron Booster (PSB) in which direct space charge is expected to be the dominant performance limitation, and to evaluate the impact of the increase of the injection energy from 50 MeV to 160 MeV as part of the Linac4 project. A lot of efforts both from experimental side and analytical/simulation side, found in Ref [1, 2] and many other papers, have been made in order to ensure the best use of injectors toward the LHC. Historically at CERN, ACCSIM has been employed for space charge simulation studies. A benchmark study using ORBIT has been started to confirm the results from ACCSIM and to profit from the advantages of ORBIT such as the capability of parallel processing. In order to make numerical noise small enough, a large number of macro-particles must be tracked in Particle-In-Cell (PIC) method, and generally simulation takes long time. Parallel processing is a possible solution to reduce the simulation time. Although this kind of benchmark has been already performed for several cases [3], it is still worth to benchmark the two codes using the specific machine parameters of the PSB in which the tune spread is unusually large (up to ~0.5 for high intensity beams and more than 0.3 for LHC type beams) and overlap with the low order resonances. 2. Differences in codes The differences in physics modeling should be mentioned when we compare the results from two codes. Both codes have a so-called 2.5 dimensional (2.5D) space charge model. A 3 dimensional (3D) space charge model is also available in ORBIT but the results in this study are from the 2.5D model. In the 2.5D model, all macro-particles in a bunch are projected to a transverse plane. This means that the transverse charge distribution is assumed to have the same shape in any longitudinal slice. The transverse kicks are finally computed from the charge distribution and weighted by the longitudinal density. In ACCSIM, the weighting can be enabled or disabled. It is, however, pointed out that the weighting has to be turned on in PSB simulations to ensure a correct evaluation of space charge forces for a bunched beam [4]. In both codes several methods are implemented to find the space charge field. There are two methods in ACCSIM†: 1) Multiple-Fourier-Transform (MFT) and 2) Hybrid-Fast-Multipole (HFM), and three methods in ORBIT: 1) Pair-wise sum, 2) BruteForce PIC and 3) FFT-PIC. In terms of both computation time and accuracy of the computed space charge field, HFM is the better algorithm in ACCSIM [7] and has been used for this study. The FFT-PIC in ORBIT is the fastest solver of the three and is used for the simulations presented here. In conclusion, the algorithm to find the space charge field is different between ACCSIM and ORBIT. Another notable difference is the binning for PIC. The grid size is automatically adjusted to the current beam size in ORBIT whereas it is fixed in ACCSIM. This kind of flexible binning aims at following the beam envelope which changes according to the beta function of the ring. It should be however noted that halo particles having very large † The first space charge implementation in ACCSIM was based on the computation of amplitude dependent tune shift [5]. The PIC type space charge modeling (MFT and HFM) was developed later [6]. amplitude transverse oscillations could cause a considerable computation error since they make the grid size too large compared to the size of beam core. This problem can be avoided by introducing a physical aperture, removing particles with large oscillation amplitude. The tracking models to transport particles through accelerator elements have also some differences. Basically transfer matrices are adopted in both codes. One of the differences is the implementation of chromaticity. Since the tracking in ORBIT is based on Teapot style (see also Appendix), every transfer matrix takes into account the momentum deviation dP/P. On the other hand, the chromatic tune shift is approximated in ACCSIM, on top of linear tracking, by introducing an additional transfer matrix which shifts the betatron phase once per turn (or several times if specified) according to the momentum deviation and the ring chromaticity. 3. Benchmarking 3-1. Beam parameters and PSB lattice A beam similar to the LHC type beam, which will be provided with the future Linac4, is examined in the benchmark study. The beam parameters are summarized in Table 1. Gaussian and elliptic distributions matched to linear optics, are generated and used as the initial distribution. Table 1: Beam parameters for benchmark Parameter Kinetic energy Intensity Initial transverse distribution Transverse emittance Initial longitudinal profile Longitudinal emittance for bunched beam (100%) Value 160 MeV 3.25×1012 proton/ring Gaussian / Elliptic ~2.5 πmm-mrad. r.m.s. normalized Parabolic 0.9 eV-s For the benchmark simulations presented here, the beam is captured with an h=1 rf system with 8 kV only, whereas in the real machine a second harmonic RF system is also used for bunch flattening. The lattice of PSB is modeled as 16 periodic cells excluding injection bumpers so that the benchmarking is as simple as possible. Bare tunes are then (4.28, 5.45). Boundary condition and physical aperture are also not modeled, that is, the beam is in a free space. The number of transverse space charge nodes is 224 in ACCSIM and about 300 in ORBIT. 3-2. Benchmarking with Gaussian distribution A Gaussian distribution is approximated by applying Joho’s formalism [8] and setting the binomial factor m=100. The incoherent space charge tune shifts are computed for the Gaussian like distribution with ORBIT and shown in Fig. 1. 5.5 Qy 5.0 3.5 4.0 4.5 Qx Figure 1: Space charge tune shifts for Gaussian distribution. Norm. RMS emittance (πmm-mrad) 5.5 5.0 4.5 4.0 3.5 3.0 2.5 0 ORBIT, Hor. (Nx,y=32) ORBIT, Ver. (Nx,y=32) ORBIT, Hor. (Nx,y=64) ORBIT, Ver. (Nx,y=64) ACCSIM, Hor. ACCSIM, Ver. 2500 5000 7500 10000 12500 15000 Turn Figure 2: Emittance evolutions for Gaussian initial distribution. Np=99999. ACCSIM results are taken from the study of Ref. [2]. Emittance evolutions are computed for the number of PIC grids of Nx,y=32 and 64 as shown in Fig. 2. The number of macro-particles Np is 99999. The same initial distribution is used in both codes to have the same initial condition. For the Gaussian distribution, the emittance evolution does not depend significantly on the number of PIC grids, and we see rather good agreement between the codes. This means that the number of PIC grids of 32 is already good enough to find the space charge field for the Gaussian like distribution, in which the tail of distribution approaches to zero smoothly. Actually, the number of particles of 105 is the limit in ACCSIM. Though it is possible to increase it with small modification of the source code, it is also a practical limit in single processing with the present computation capability. Simulations with larger numbers of macro-particles have been carried out with ORBIT, and the results are shown in Fig. 3. Norm. RMS emittance ( πmm-mrad) 3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 0 2000 4000 6000 8000 Horizontal Np=99999 Np=199999 Np=499999 Np=999999 Np=1999999 Vertical Np=99999 Np=199999 Np=499999 Np=999999 Np=1999999 Turns Figure 3: Emittance evolutions for Gaussian initial distribution for various number of particles (Nx,y=64). Whereas the emittance evolution does not depend on the number of PIC grids, it depends on the number of macro-particles rather strongly. The number of particle of 105 seems not enough in this case, and the emittance evolutions show signs of convergence with more than 5×105 particles. 3-3. Benchmarking with elliptic distribution Elliptic distributions applying Joho’s formalism (m=1.5) are also used to generate initial distributions. Figure 4 shows the incoherent tune shifts for the elliptic distribution obtained with ORBIT. As expected from analytical formulae, the maximum tune shift (linear tune shift) is smaller than that of a Gaussian distribution due to the smaller density at the beam center. 5.5 Qy 5.0 3.5 4.0 4.5 Qx Figure 4: Space charge tune shifts for elliptic distribution The emittance evolutions for the elliptic distribution are shown in Fig. 5. The same initial distribution is again used in both codes to have the same initial distribution, and the number of macro-particles is 99,999. Norm. RMS emittance (πmm-mrad) 2.8 2.7 ORBIT, Hor. (Nx,y=32) ORBIT, Ver. (Nx,y=32) ORBIT, Hor. (Nx,y=64) ORBIT, Ver. (Nx,y=64) ACCSIM, Hor. ACCSIM, Ver. 2.6 2.5 2.4 0 1000 2000 3000 4000 Turn Figure 5: Emittance evolutions for elliptic initial distribution. Np=99999. ACCSIM results are taken from the study of Ref[2]. It is observed that the horizontal emittance evolution shows rather good agreement whereas the vertical one is quite different. For the elliptic distribution, the vertical emittance growth is particularly sensitive to the number of PIC grids. Generation of halo particles is observed together with a sudden blow-up in ORBIT. This generation of halo is also observed in ACCSIM but does not contribute significantly to the r.m.s. emittance. The number of PIC grids of 32 would not be good enough to describe the elliptic distribution, in which the tail of distribution approaches zero quickly, and the motion of halo particles. The emittance evolutions for lager number of macro-particles are simulated also for elliptic distribution with ORBIT as shown in Fig. 6. Again, the number of particles of 105 seems not enough, and the emittance evolutions show signs of convergence with more than 5×105 particles. Norm. RMS emittance (πmm-mrad) 2.60 2.55 2.50 Horizontal Np=99,999 Np=199,999 Np=499,999 Np=999,999 Vertical Np=99,999 Np=199,999 Np=499,999 Np=999,999 2.45 0 1000 2000 3000 4000 5000 Turn Figure 6: Emittance evolutions for elliptic initial distribution for various number of particles (Nx,y=64). In order to understand the vertical blow-up from possible observables, vertical phase space plots around the blow-up and the time evolution of second order moment of particle distribution are shown in Fig. 7 and Fig. 8. 6 4 2 turn 800 6 4 2 turn 1000 6 4 2 turn 1200 y' (mrad) y' (mrad) 0 -2 -4 -6 -20 0 -2 -4 -6 -20 y' (mrad) -15 -10 -5 0 5 10 15 20 0 -2 -4 -6 -20 -15 -10 -5 0 5 10 15 20 -15 -10 -5 0 5 10 15 20 y (mm) y (mm) y (mm) 6 4 2 turn 1400 6 4 2 turn 1600 y' (mrad) 0 -2 -4 -6 -20 y' (mrad) -15 -10 -5 0 5 10 15 20 0 -2 -4 -6 -20 -10 0 10 20 y (mm) y (mm) Figure 7: Evolution of vertical phase space during sudden blow-up in ORBIT. For Np=99,999, Nx,y=32. Particles of 10,000 are plotted turn 400 turn 800 turn 1200 (Blow-up) turn 1600 norm. 2nd order vertical moment (μm) 4.6 4.8 Norm. 2nd vertical moment (μm) 4.6 4.4 4.4 Amplitude of quadrupole oscillation 4.2 4.2 4.0 4.0 0 500 1000 Turn 1500 2000 0 20 40 60 80 100 120 140 s (m) (a) Turn by turn (b) Snap shots for every 400 turns Figure 8: Second vertical moment during sudden blow-up in ORBIT. As seen in Fig. 7, halo particles are appearing during only a few hundred turns. Although no significant change of the beam core is observed, the second moment shows the excitation of quadrupole oscillation plotted in Fig. 8 (a). Figure 8(b) shows the second moment over one turns observed every 400 turns. During the sudden blow-up (turn 1,200), there seems the harmonics of five overlapping the lattice periodicity of 16, thus implying that the resonance of 2Qy=10 might be one of sources of emittance blow-up. 4. Summary and outlook A benchmarking study of ACCSIM and ORBIT codes has been performed using a beam similar to the LHC beam, which will be provided by the PSB with Linac4. We observed a fair agreement in emittance evolution in the horizontal plane but not in the vertical one. Note that the algorithm to compute the space charge field is different between the two codes. The simulation parameters, like the number of PIC grids and the number of macro-particles, should be properly chosen. In ORBIT, the optimum number of PIC grids would depend on the particle distribution, and is relatively large when the incoherent motion of halo particles is concerned. From ORBIT simulations, we observed that a larger number of macro-particles results in smaller emittance blow-up. The number of macro-particles of 105, which is the practical limit in ACCSIM, seems not sufficient for the PIC-FFT algorithm used for ORBIT runs. It is, however, not clear whether it is enough or not for the HFM algorithm used for ACCSIM runs. Therefore, it would be interesting to install HFM into ORBIT and to examine a larger number of macro-particles with parallel processing. Further studies such as a test of 3D space charge model, space charge under boundary condition would be of some importance. Acknowledgements I would like to thank Drs. G. Arduini, C. Carli, M. Chanel, S. Cousineau (SNS), A. Franchi, R. Garoby, F. Gerigk, M. Martini, M. Medahhi, A. Shishlo (SNS) and F. Zimmermann for valuable discussions and a lot of supports. References [1] M. Chanel, “SPACE CHARGE MEASUREMENTS AT THE PSB”, Proc. of BEAM’07, pp106-110 (2008) [2] M. Martini, “STATUS OF THE SPACE CHARGE SIMULATION STUDIES OF THE PS BOOSTER PERFORMANCE WITH LINAC4 INJECTION USING ACCSIM”, CERN-AB-2008-005 ABP (2008) [3] S. Cousineau, “SIMULATION TOOLS FOR HIGH INTENSITY RINGS”, Proc. of PAC’03, pp259-263 (2003) [4] S. Cousineau, “Summary of Initial ORBIT Simulations of the CERN PSB”, SPL meeting #82, https://twiki.cern.ch/twiki/bin/view/SPL/SPLminutes13December2006 [5] H. Schonauer, “ADDITION OF TRANSVERSE SPACE CHARGE TO ACCSIM CODE”, TRIUMP design note, TRI-DN-89-K50 (1989) [6] F. W. Jones and H. O. Schonauer, “New Space-Charge Methods in Accsim and Their Application to Injection in the CERN PS Boster”, Proc. of PAC’99, pp2933-2935 (1999) [7] F. W. Jones, “A Hybrid Fast-Multiple Technique for Space-Charge Tracking with Halos”, AIP Conf. Proc. 448, pp359-370 (1998) [8] W. Joho, “Representation of Beam Ellipses for Transport Calculations”, SINREPORT TM-11-14 (1980) [9] L. Schachinger and R. Talman, “Teapot: A Thin-Element Accelerator Program for Optics and Tracking”, Particle Accel. 22, 35 (1987) Appendix: Tracking modes in ORBIT ORBIT provides two algorithms to track particles through accelerator elements: - tracking using transfer matrices generated with MAD8 (or DIMAD) - Teapot [9] style tracking They are called linear tracking and nonlinear tracking, respectively. In the first method, the second order transfer matrices are also taken into account by implicitly specifying the command to use them in the input script to run the code. It is, however, not recommended to use them, especially for a long-term tracking, because the tracking algorithm does not satisfy the simplectic condition. The second method treats multipole elements properly. The main difference between two methods for PSB simulation, in which the lattice does not have any multipole element, is the chromatic tune shift. In Teapot style tracking, a kick applied to particle varies in each accelerator element depending on the momentum deviation dP/P whereas it is not taken into account in the first method. Figure A shows the chromatic tune shift in PSB computed with Teapot style tracking. 5.50 Linear tracking Nonlinear tracking Qy 5.45 5.40 4.25 4.30 Qx Figure A: Chromatic tune shift For a beam with dP/P=0.34%. The chromaticity of the ring is (-3.7, -10.3).