Collaborators

1. ORBIT Code Developmentat FNAL:A Progress ReportJean-Francois OstiguyAcceleration Integration DepartmentFNALostiguy@fnal.gov

2. Collaborators Jeff Holmes SNS/ORNL Leo Michelotti BD/BP FNAL Weiren Chou BD/BP FNAL

3. Why ORBIT ? • We considered a variety of existing codes • ORBIT was selected as our primary tool because of: • Source code and documentation publicly available • Well-developed diagnostics (tune footprints, moment evolution etc .. ) • Some efforts to validate the code at existing machines (PSR) • Support for both decoupled transverse and longitudinal (21/2D) and 3D space charge • Support for parallel execution (MPI) • Note: (Synergia) - essentially a derivative of IMPACT (developed to study halo in high intensity linacs) is also under development at FNAL. The focus is “full 3D” simulation in synchrotrons with high statitics. Cross checks are useful and important for validation of ORBIT and Synergia and other codes.

4. Beam Physics Dept Parallel Linux Cluster 32 2-CPU Nodes (1.4 Ghz AMD Athlon) 1 Gbyte RAM / Node Gigabit Ethernet Total Cost: 65 K$ (Note:adequate cooling: +25 K$ !) Adequate for 21/2 D simulations. 100-1000 turns with O(105) macro particles. We expect adding 16 more nodes in FY2004.

5. Code Development Recent development efforts at Fermilab have been focused on: • MAD parser • Internally generated high-order maps • A high level Python shell • Better support for Acceleration (e.g. energy-dependent multipoles, transition etc…) • Validation and bug fixes (e.g.: tune footprint Computation in presence of RF cavities)

6. A Lex/Yacc based MAD parser • A Lex/Yacc-based MAD parser was developed a few years ago for internal needs (BEAMLINE and MARS). • Designed as a generic system usable either for off-line translation to another human-readable description language or for dynamic definition of objects. • Few minor restrictions (no abbreviations, no action commands, no use of undefined variables) • Successfully validated on very large lattice files (e.g. complete Tevatron lattice) • BTW: This parser has been borrowed by Synergia

8. Internal Map Generation • Until recently, ORBIT had been relying on MAD or DIMAD to produce maps and lattice functions. The lattice information is read from MAD (ascii) output. • As a consequence, ORBIT could not internally recompute maps and/or relevant quantities such as lattice functions or chromaticities. • Map generation and particle propagation (map evaluation) using BEAMLINE/Mxyzptlk is now fully functional and validated.

9. Mxyzptlk and Beamline(L. Michelotti) • Mxyzptlk: a C++ class library to perform automatic differentiation to a user-specified order n. Provides overloaded operators, trig functions etc … • Beamline: a C++ class library build on top of mxyzptlk. Provides facilities to create, edit and manipulate beamlines hierarchies, compute lattice functions, chromaticities, maps (to order n), concatenate maps etc ...

10. Mxyzptlk ? Superman's foe from the 5th dimension. He will return to his own dimension only if he is made to spell his name backwards ...

11. Drifts, Slots etc … Beamline uses local coordinates attached to the magnet. No assumption is made about the closed orbit. Slots can be considered as drifts with (usually) one face that not perpendicular to the closed orbit. Slots are automaticallyconcatenated with CF_rbends to form “conventional” (MAD-style) maps. VL01:right 140 s 0.0375 2.07827 drift 19 QL01: 150 s 0.075 2.15327 quadrupole 20 SL01: 160 s 0.075 2.22827 quadrupole 21 L01H1: 170 s 0.125 2.35327 drift 22 SPACE0: 180 s 3.00961 5.36288 Slot CF_rbend Slot 23 SPACE1: 190 s 0.38 5.74288 drift 24 SPACE0: 200 s 3.00961 8.75249 Slot CF_rbend Slot 25 DRSEXTS1: 210 s 0.11 8.8625 drift 26 SEXTS: 220 s 0.3 9.16249 sextupole 27 DRBPML: 230 s 0.035 9.19749 drift

12. Why Use the BEAMLINE Class ? • Written in C++, just like ORBIT • The same code automatically supports 1st , 2nd and or order n maps if desired • Support for arbitrary misaligments(tilt, yaw, offset etc ..) • No “large ring” approximation(s). High intensity boosters are likely to be “small” rings. • Propagation physics and computation completely under user control. • Maps and thin kicks (a la TPOT) can freely be intermixed. • BTW: Synergia uses BEAMLINE to generate (first order) map coefficients; however propagation is handled by IMPACT.

13. Beamline Class Library Computational Performance • Initial attempts at using the BEAMLINE library resulted in discouraging performance (1-2 orders of magnitude !) • In MXYZPTLK, polynomials are implemented as doubly linked lists. Each list node contains a non-zero monomial coefficient as well as an integer which is uniquely mapped to actual monomial exponents. • evaluating a polynomial (map) implies traversal of a linked list (indirections) and recovery of the actual monomial exponents from unique integer encoding • A compact storage scheme is desirable, especially at high orders because of a rapid explosion in the no of polynomila terms. At low orders, it was found necessary to store map coefficients and invidual monomial exponents explicitly in linear arrays in order to get satisfactory map evaluation performance.

14. Propagator • In the beamline library, the Propagator is a functor ( i.e. a function object) that determines how a particle is propagated through a beamline element. • Each element is assigned a default propagator, which can be overridden by a user-supplied alternate • When the element is a Map, the propagator operator() trivially evaluates a polynomial in 6N variables for each phase space dimension . • Because the Propagator is an object, an alternate (private) polynomial representation can be instantiated by the propagator constructor. • Using this technique, tracking using the facilities provided by the beamline library has been verified to be as efficient as more conventional code.

15. 2nd Order Maps Validation Simple test: tune spread associated with 2 different momentum distributions (dp/p= 2.0 E-3, 4.0 E-3). Lattice-based tunes:6.935 H, 6.662 V Lattice-based chromaticities:-9.86H, -6.89V Tracking-based chromaticities: -9.3 H, -6.0 V

16. Python Shell Problem: ORBIT is build on top of SuperCode, an intepreted C++ like scripting language that is orphaned and minimally documented. Solution: use Python ! • Python is a mature scripting language • Its object model is semantically highly compatible with that of C++ • Operator overloading • Exceptions • Good tools are available for automatic interface code generation. • A wealth of high quality Python code is freely available for reuse.

17. ORBIT Code Structure • ORBIT is structured as “modules” controlled by a high level interpreter, SuperCode • SuperCode was designed to have a C++ -like syntax • interface code generated from special interface definition files by a program : MGen • In principle, modules could be written in f77, C or C++ • Exported interface from modules is the common denominator between all these languages: static functions and variables • BTW: In ORBIT, the shell is an integral part of the code. Input syntax checks and runtime diagnostics are generated by the shell.

18. Why Preserve the Interpreter/Modules Structure ? • A high level interpreter allows for rapid implementation of custom features. While the efficiency of a low-level language is required to propagate large no of particles, diagnostics and a posteriori analysis can benefit from high level language implementations because they are often problem-specific. • If efficiency becomes an issue, functionality implemented at the interpreter level can be reimplemented into a compiled module without affecting existing scripts. • The interpreter/module structure promotes well-defined interfaces. This makes it easy to contribute new functionality without deep knowledge of the entire code.

19. Python/C++ Interface Code Generation Currently are three systems available: • SWIG (www.swig.org) by David Beazley, U of Chicago. Comprehensive system, support for most scripting languages. While support for C is excellent, support for C++ constructs has serious limitations. • SIP (www.riverbankcomputing.co.uk) by Phil Thompson. Similar in philosophy to SWIG, but python/C++ specific. Not well documented, requires special interface files • Boost.python (www.boost.org) by David Abrahams. Python/C++ specific. Implemented as a generic C++ library (mostly header files). Uses template metaprogramming techniques to generate interface code; no special program needed beyond an ANSI standard C++ compiler. Python classes may inherit from C++ classes and override virtual functions !

20. Emulate existing SuperCode data types (e.g. Vector, matrix, 3D array, String) Use std C++ functionality to emulate existing (custom) RTTI and exception functionality As much as possible, preserve existing syntax to facilitate porting. Porting Strategy

21. Boost.python Example : SuperCode Data Type Emulation #include "sc-types.h" #include "sc-string.h" #include "Python.h" #include <boost/python/python.hpp> using std::complex; using namespace boost::python; void wrap_supercode() { // ** ComplexMatrix ** python::class_<Matrix<complex<Real> > >("ComplexMatrix", init<int,int>()) .def(init<const Matrix<complex<Real> > >()) .def("get", &Matrix<complex<Real> >::get) .def("set", &Matrix<complex<Real> >::set) .def("__repr__", &Matrix<complex<Real> >::print) .def("clear", &Matrix<complex<Real> >::clear) .def("resize", &Matrix<complex<Real> >::resize) .def(python::self + python::self) .def(python::self - python::self) .def(python::self * python::self) .def(python::self ^ python::self) ; }

22. Boost.python Example: Exporting an ORBIT class void wrap_lspacecharge() { class_<LSpaceCharge>("LSpaceCharge", no_init) .def("addFFTLSpaceCharge", &LSpaceCharge::addFFTLSpaceCharge, "Routine to add an FFT Longitudinal space charge node") .staticmethod("addFFTLSpaceCharge") .def("showLSpaceCharge", &showLSpaceCharge, "Routine to show the Longitudinal space charge info") .staticmethod("showLSpaceCharge") .def("dumpKick", &LSpaceCharge::dumpKick, "Routine to dump kick of lsc node n to at an angle") .staticmethod("dumpKick") .def("dumpPhiCount", &LSpaceCharge::dumpPhiCount, "Routine to dump phi Bin information to vector p") .staticmethod("dumpPhiCount") .def("dumpFFTStuff", &LSpaceCharge::dumpFFTStuff, "Routine to dump phi Bin information to vector p") .staticmethod("dumpFFTStuff") .def("dumpZImpedStuff", &LSpaceCharge::dumpZImpedStuff, "Routine to dump Z and chi information to vectors p1, p2") .staticmethod("dumpZImpedStuff") .def_readwrite("nLongSCs", &LSpaceCharge::nLongSCs) .def_readwrite("nFFTLongs", &LSpaceCharge::nFFTLongs) .def_readwrite("nLongBins", &LSpaceCharge::nLongBins) ; }

23. ORBIT script vs PyORBIT script Booster.py Booster.sc # ... ## Add an RF Cavity nRFHarms = 1 global volts volts = RealVector(nRFHarms) harmNum = RealVector(nRFHarms) RFPhase = RealVector(nRFHarms) harmNum.set(1,1.0) RFPhase.set(1,0.0) def PDIIVolts(): tFactor = Ring.time/63.213 if(tFactor > 1.): tFactor = 1.0 volts.set(1, (800.*tFactor)) Accelerate.addRampedRFCavity("RF 1", 4, nRFHarms, volts, harmNum, RFPhase, PDIIVolts) ## Add a Longitudinal Space Charge Node ZImped = ComplexVector(32) ZImped.setall(0+0j) # … //… // Add an RF Cavity Integer nRFHarms = 1; RealVector volts(nRFHarms); RealVector harmNum(nRFHarms); RealVector RFPhase(nRFHarms); harmNum(1) = 1; RFPhase(1) = 0.; Void PDIIVolts() { tFactor = time/63.213; if(tFactor > 1.) tFactor = 1.; volts(1) = 800. + 0.0 * tFactor; } addRampedRFCavity("RF 1", 4, nRFHarms, volts, harmNum, RFPhase, PDIIVolts); // Add a Longitudinal Space Charge Node ComplexVector ZImped(32); ZImped = Complex(0.,0.); // …

24. To efficiently produce a (tranverse) tune footprint, ORBIT computes tunes by accumulating phase advance in normalized (Floquet) coordinates. This is an approximation since constant action contours are circular only in the linear approximation. The normalized coordinates do not include dispersion. (dp/p is a perturbation !). To compute the tune of an off-momentum particle, the dispersive contribution to the trajectory must first be subtracted off. All collective fields are “frozen” during a tune computation i.e. although they influence the macro particles, they are not updated. All RF cavities must effectively be turned OFF during a tune computation, otherwise the momentum upstream and downstream of a cavity is different, leading to erroneous dispersive corrections and erroneous tunes. This was not automatically enforced in ORBIT. Tune Computation

25. Phase Accumulation Contour distorted by non-linearities Weakly nonlinear one turn map X2,Y2 2 (X2,Y2) = M (X1, Y1) X1,Y1 1 In presence of nonlinearities, averaging over multiple turns gives better results.

26. Acceleration ORBIT currently provides some somewhat simple-minded support for acceleration. However, because SNS is an accumulation ring, investing time and efforts to further develop the existing functionality has not been the highest priority of the main developers (justifiably so !). The code reliance on externally computed maps, implicitly prevented changes during acceleration. While this may be an acceptable first approximation, a more realistic simulation should include: Saturation effects (energy dependent field defects) Remanent field effects (constant field defects) Transition crossing mechanics (phase jump, pulsed quadrupoles etc ...)

27. Saturation • In some machines, magnetostatic saturation effects are important. • Optical aberrations due to magnetic saturation increase rapidly with energy once the field reaches some threshold value. • Computer model: interpolation tables associated with each element. Map coefficients need to be recomputed as E varies. • Observation: Booster uses CF magnets which are mostly not pushed into saturation. Saturation effects in Booster are less important than in typical separated function synchrotron.

28. Remanent Magnetization • Remanent magnetization distorts the field significantly at injection. • Often sets the minimum injection energy and “dynamic range” of a synchrotron • Simple model: Remanent magnetization contributions to the total magnetic field are fixed and independent of energy. Therefore: optically, the aberrations due to remanent field are energy dependent. Their relative importance scales like ~1/E • Computer model: fixed contribution. Maps need to be recomputed as E varies. • In the FNAL Booster, remanent magnetization induced sextupole is known to be important.

29. Implementing Acceleration • Every FNALNode is a self-contained object which includes a pre-computed map as well as sufficient info to recompute a new map. • FNALNode inherits from ORBIT Node and shares its interface. • Maps are typically updated before acceleration begins (misalignements) or once/turn when necessary (e.g. saturation, gamma-t jump quads). • Typically, space charge field computations completely dominate execution time.

30. FNAL MAD parser integrated into ORBIT BEAMLINE based map generation functional and validated. 1st and 2nd order propagation optimized. Execution speed on par with the original code based on externally generated maps . Python shell work completed and ready for testing. Bug fixes: tune footprint computation in presence of RF voltage. incorrect scaling of longitudinal coordinates. Obscure memory management problems in mxyzpltk. Work on improved support for acceleration in progress FNAL Orbit now in sync with SNS ORBIT CVS as of Aug 20 2003 J. Holmes visiting FNAL Oct 2-3 Sources and more info available from unofficial FNAL Web:www-ap.fnal.gov/~ostiguy/ORBIT Conclusions and Status