ORBIT Code Development Activities at FNAL Jean Francois Ostiguy by mikesanye

VIEWS: 8 PAGES: 30

									ORBIT Code Development
       at FNAL:
   A Progress Report

              Jean-Francois Ostiguy
       Acceleration Integration Department
                      FNAL
                             ostiguy@fnal.gov




 J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                Collaborators



Jeff Holmes                SNS/ORNL

Leo Michelotti BD/BP FNAL
Weiren Chou BD/BP FNAL




  J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                        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.
           J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
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.




 J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                  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)



           J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
    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
          J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
        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.
         J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
           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 ...



         J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                  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 ...




J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                       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 automatically
concatenated 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



                  J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
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.




       J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
      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.
        J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                                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.


              J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                       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




          J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                               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.

               J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
            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.


            J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                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.

              J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
         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 !



               J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                 Porting Strategy


●   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.

        J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                                     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)
            ;
}




                    J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                                  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)
           ;
}


                 J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
           ORBIT script vs PyORBIT script
       Booster.sc                                                     Booster.py
//…                                              # ...
// Add an RF Cavity                              ## Add an RF Cavity
Integer nRFHarms = 1;                            nRFHarms = 1
RealVector volts(nRFHarms);                      global volts
RealVector harmNum(nRFHarms);                    volts   = RealVector(nRFHarms)
RealVector RFPhase(nRFHarms);                    harmNum = RealVector(nRFHarms)
                                                 RFPhase = RealVector(nRFHarms)
harmNum(1) = 1;
RFPhase(1) = 0.;                                 harmNum.set(1,1.0)
Void PDIIVolts()                                 RFPhase.set(1,0.0)
  {                                              def PDIIVolts():
    tFactor = time/63.213;
    if(tFactor > 1.) tFactor = 1.;                     tFactor = Ring.time/63.213
    volts(1) = 800. + 0.0 * tFactor;                   if(tFactor > 1.):
  }                                                        tFactor = 1.0
addRampedRFCavity("RF 1", 4, nRFHarms, volts,          volts.set(1, (800.*tFactor))
                 harmNum, RFPhase, PDIIVolts);
                                                 Accelerate.addRampedRFCavity("RF 1", 4, nRFHarms, volts,
// Add a Longitudinal Space Charge Node                                      harmNum, RFPhase, PDIIVolts)

ComplexVector ZImped(32);                        ## Add a Longitudinal Space Charge Node
ZImped = Complex(0.,0.);
                                                 ZImped = ComplexVector(32)
// …                                             ZImped.setall(0+0j)

                                                 # …




                      J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                      Tune Computation
●   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.




             J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
     Phase Accumulation
                             Contour distorted by non-linearities
                                             Weakly nonlinear one turn map

X2,Y2
     2
                     X1,Y1              (X2,Y2) = M (X1, Y1)
                  1

                                       In presence of nonlinearities, averaging
                                       over multiple turns gives better results.




J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                            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
...)


             J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
                            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.


          J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
        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.


          J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
       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.


          J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003
            Conclusions and Status
●   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




           J.-F. Ostiguy – Midwest Accelerator Physics Collaboration Meeting – Sep 30 2003

								
To top