Swifter

Author: David E. Kaufmann

Date : 01 September 2005

INTRODUCTION

The Swifter subroutine package provided here is designed to integrate a set of mutually gravitationally interacting massive bodies together with a group of massless test particles which feel the gravitational influence of the massive bodies but do not affect each other or the massive bodies.
(NOTE: the SyMBA integrator supports a second class of massive bodies whose masses are less than some user-specified value, MTINY. These bodies gravitationally affect and are affected by the more massive bodies, but do not interact with themselves.)
Seven integration techniques are included thus far:

1. Wisdom-Holman Mapping (WHM). This is the N-body mapping method of Wisdom & Holman (1991; AJ, 102, 1528).
2. Regularized Mixed Variable Symplectic (RMVS) method. This is an extension of WHM that handles close approaches between test particles and massive bodies. See Levison & Duncan (1994; Icarus, 108, 18).
3. Democratic Heliocentric (DH, or HELIO) method. This is a basic symplectic integrator (i.e., no close encounters) that uses democratic heliocentric coordinates. See Duncan, Levison, & Lee (1998; AJ, 116, 2067).
4. Symplectic Massive Body Algorithm (SyMBA). This is an extension of HELIO that handles close approaches between massive bodies (with mass greater than or equal to MTINY) and any other type of object (massive body of any mass or massless test particle). This algorithm is also described in Duncan, Levison, & Lee (1998). See also Levison & Duncan (2000; AJ, 120, 2117).
5. A fourth-order T+U Symplectic (TU4) method. This is the T+U method of Candy & Rozmus (1991; J. Comp. Phys., 92, 230). Also see Gladman, Duncan, & Candy (1991; CeMDA, 52, 221).
6. A nonsymplectic fifteenth-order integrator that uses Gauss-Radau spacings (RA15). This algorithm is described by Everhart (1985; ASSL Vol. 115: IAU Colloq. 83: Dynamics of Comets: Their Origin and Evolution, 185).
7. A Bulirsch-Stoer (BS) method. This is a Bulirsch-Stoer from Press, Teukolsky, Vetterling, & Flannery (1992; Numerical Recipes in FORTRAN).

In order to get to this file you have presumably obtained, uncompressed and untar-ed the swifter.tar.Z file.
You will now find various files and subdirectories.
Most of the directories contain the subroutines that make up Swifter.
However, fully functional examples of working drivers are found in the subdirectory main.
The subroutines used are grouped by function in subdirectories with what we hope to be good internal documentation.

COMPILING Swifter

There is one file that must be edited in order to compile the Swifter library and drivers. In the top Swifter directory, you will find a file called Makefile.Defines. This file contains variable definitions that the make program uses to control the building of the Swifter library and drivers (via the supplied Makefile) as well as the FXDR library (via the supplied Makefile.fxdr). (The FXDR library and XDR, the eXternal Data Representation, are described further at the end of this file.) Edit Makefile.Defines. You may have to change the values of:

• SWIFTER_HOME : the top-level Swifter directory (this directory should be specified by its full pathname)

• USER_MODULES : space-separated list of any additional modules that you write and place in directory $(SWIFTER_HOME)/module for inclusion in the Swifter library

• FORTRAN : the command to run the Fortran 90/95 compiler on your machine

• FFLAGS : the flags you use for your Fortran (NOTE: this flag list should NOT include the switch, typically -c, used to cause the compiler to generate an object file only)

• CC : the command to run the C compiler on your machine

• CFLAGS : the flags you use for your C (NOTE: this flag list should NOT include the switch, typically -c, used to cause the compiler to generate an object file only)

• SHELL : the shell make uses to execute commands in a rule (the default is the Bourne shell /bin/sh and should NOT need to be changed in general)

• AR : the library archive command on your machine (the default is ar and should NOT need to be changed in general)

• RANLIB : the command to generate an index to a library archive (the default is 'ranlib', which should be available on most machines; an alternative would be 'ar -s')

• INSTALL : the command to install a file in a filesystem (the default is 'install' and should NOT need to be changed in general)

• INSTALL_PROGRAM : the command to install an executable file in a filesystem (the default is '$(INSTALL) -m 755' and should NOT need to be changed in general)

• INSTALL_DATA : the command to install a non-executable file in a filesystem (the default is '$(INSTALL) -m 644' and should NOT need to be changed in general)

Also present in Makefile.Defines are definitions of the variables F77CMD, F77OPTS, CCCMD and CCOPTS.
These are needed for compatibility with the FXDR library makefile and SHOULD NOT BE CHANGED.

After editing Makefile.Defines, type make or make all. This will cause the following actions to be taken:

1. The modules in directory $(SWIFTER_HOME)/module will be compiled to object code and added to the Swifter library (libswifter.a) located in $(SWIFTER_HOME)/lib. Also, all the .mod files generated will be installed to $(SWIFTER_HOME)/lib as well. These files are needed for subsequent compilation of other source files that reference the modules. (NOTE: DO NOT CHANGE the ordering of the modules in the definition of $(SWIFTER_MODULES) in the Makefile, as the successful compilation of some of the modules depends on this ordering.)
2. The remaining Swifter library source code located in the various subdirectories of $(SWIFTER_HOME) will be compiled to object code and added to libswifter.a.
3. The FXDR library (libfxdr.a) will be regenerated from source, tested, then installed in $(SWIFTER_HOME)/lib.
4. Any Swifter drivers located in $(SWIFTER_HOME)/main will be compiled and linked with libswifter.a and libfxdr.a, and the resulting executables installed in $(SWIFTER_HOME)/bin. The name of the resulting executable will by default be the name of the original source file minus the .f90 suffix (e.g., swifter_whm.f90 -> swifter_whm)
5. Any tools located in $(SWIFTER_HOME)/tool will be compiled, linked, installed and named in exactly the same manner as were the drivers in Step (4) above.

One by-product of the initial call to make or make all is the creation of soft links to the files Makefile and Makefile.Defines from each Swifter subdirectory containing source code. This facilitates subsequent compilation directly from these subdirectories, although the make processes described below should work correctly, unless otherwise noted, from any directory having access to these two files.

1. make mod : recompiles the module source code located in $(SWIFTER_HOME)/module, replaces these objects within libswifter.a, located in $(SWIFTER_HOME)/lib, and installs the .mod files to $(SWIFTER_HOME)/lib
2. make lib : rebuilds the entire Swifter library libswifter.a from source, replacing all the old archive members
3. make libdir : rebuilds only the libswifter.a source code LOCATED IN THE CURRENT WORKING DIRECTORY, replacing the old archive members
4. make fxdr : rebuilds the FXDR library from source, tests it, and installs it to $(SWIFTER_HOME)/lib
5. make drivers : recompiles all the driver source code located in $(SWIFTER_HOME)/main, links to the Swifter and FXDR libraries, and installs the resulting executables in $(SWIFTER_HOME)/bin
6. make tools : performs exactly the same steps as for the driver source code above, but on the tools source code located in $(SWIFTER_HOME)/tools
7. make bin : recompiles all the source code LOCATED IN THE CURRENT WORKING DIRECTORY, links to the Swifter and FXDR libraries, and installs the resulting executables in $(SWIFTER_HOME)/bin
8. make clean : removes any soft links to the files Makefile and Makefile.Defines from subdirectories of $(SWIFTER_HOME), removes all executables from $(SWIFTER_HOME)/bin, removes all libraries (lib*.a) and compiled modules (*.mod) from $(SWIFTER_HOME)/lib and the FXDR include file fxdr.inc from $(SWIFTER_HOME)/include

THE DRIVERS

All of the drivers are designed to take essentially the same input files and to produce the same type of output. The basic step for all of them begins with heliocentric positions and velocities and advances them a timestep dt.

• swifter_whm : The basic WHM integrator when particles are to be removed at the time of close planetary encounters. Arbitrarily close solar encounters can occur.
• swifter_rmvs : The RMVS integrator that is like WHM, except that it can handle arbitrarily close encounters between massive bodies and massless test particles.
• swifter_helio: The basic HELIO integrator (uses democratic heliocentric coordinates) when particles are to be removed at the time of close planetary encounters. Close encounters with the sun are NOT allowed by the HELIO integrator.
• swifter_symba: The SyMBA integrator, based on HELIO, can follow arbitrarily close encounters between the massive bodies (i.e., planets with masses greater than MTINY) and any of the other bodies in the simulation (except the Sun). We are working to include the version that handles close encounters with the sun, but it is not available yet in this beta release. swifter_tu4 : Symplectic fourth-order T+V integrator when particles are to be removed at the time of close planetary or solar encounters. swifter_ra15 : Nonsymplectic fifteenth-order RADAU integrator. This integrator can follow arbitrarily close encounters between any of the objects, but is much less efficient than the included symplectic integrators. The accuracy is controlled by an error tolerance input by the user. swifter_bs : Nonsymplectic Bulirsch-Stoer integrator. This integrator can follow arbitrarily close encounters between any of the objects, but is much less efficient than the included symplectic integrators. The accuracy is controlled by an error tolerance input by the user.

INPUT/OUTPUT

Swifter takes input from three files: a parameter file, a planet file, and a test particle file. The user inputs the name of the parameter file at the prompting of the program. The names of the other two files are contained in the parameter file.

NEW The user may also enter the name of the input file as a command line argument, rather than at the prompt.

Parameter File: This file contains all the run-time parameters

The structure of the parameter file is a list of PARAMETER VALUE pairs. The list of recognized parameters is given below. The leftmost column indicates the type of parameter (R=real, I=integer, S=string, F=flag [flag values are the strings YES and NO. If the value is YES, the flag will be set internally to .TRUE., if NO, then .FALSE.]) The second column gives the parameter name, and the third column describes the parameter. Default values, if applicable, are enclosed in [brackets]. The parameter names and values are NOT case sensitive.

I     NPLMAX         maximum number of planets [-1] (not yet used)
I     NTPMAX         maximum number of test particles [-1] (not yet used)
R     T0             initial time (REQUIRED) [0.0]
R     TSTOP          time to stop the integration (REQUIRED) [0.0]
R     DT             time step (REQUIRED) [0.0]
S     PL_IN          planet data filename (REQUIRED) [""]
S     TP_IN          test particle data filename [""]
S     IN_TYPE        format of PL_IN, TP_IN ["ASCII"] | "XDR8"
I     ISTEP_OUT      number of time steps between outputs [-1]
S     BIN_OUT        binary output filename (REQUIRED if ISTEP_OUT > 0) [""]
S     OUT_TYPE       format of binary output file (REQUIRES BIN_OUT)
                                      "REAL4" | "REAL8" | ["XDR4"] | "XDR8"
S     OUT_FORM       data stored in binary output file (REQUIRES BIN_OUT)
                                      "EL" | ["XV"] | "FILT"
S     OUT_STAT       binary output file status (REQUIRES BIN_OUT)
                                      ["NEW"] | "UNKNOWN" | "APPEND"
I     ISTEP_DUMP     number of time steps between dumps [-1]
R     J2             spher. harm. term J_2*R^2, R = central body radius [0.0]
R     J4             spher. harm. term J_4*R^4, R = central body radius [0.0]
F     CHK_CLOSE      check for test particle/planet close encounters
                     "YES" | ["NO"]
R     CHK_RMIN       heliocentric distance at which a test particle is stopped
                     as being too close to the central body [-1.0]
R     CHK_RMAX       heliocentric distance at which a test particle is stopped
                     as being too distant from the central body [-1.0]
R     CHK_EJECT      heliocentric distance at which an energetically unbound
                     test particle is stopped as being too distant from the
                     central body [-1.0]
R     CHK_QMIN       pericenter distance at which a test particle is stopped
                     as being too close to the pericenter [-1.0] (the
                     pericenter is interpreted either as the center of the
                     central body or the system barycenter, depending on the
                     value of CHK_QMIN_COORD)
S     CHK_QMIN_COORD coordinate frame to use for CHK_QMIN (REQUIRES CHK_QMIN)
                     ["HELIO"] | "BARY"
R R   CHK_QMIN_RANGE lower and upper boundaries of semimajor axis range to
                     perform the CHK_QMIN check (REQUIRES CHK_QMIN) [-1.0,-1.0]
S     ENC_OUT        encounter filename [""]
F     EXTRA_FORCE    use additional user-specified force routines
                     "YES" | ["NO"]
F     BIG_DISCARD    include data for all bodies > MTINY for each discard
                     record "YES" | ["NO"]
F     RHILL_PRESENT  Hill sphere radius is included for each planet in PL_IN
                     "YES" | ["NO"]

Sample Parameter File: contents between the = lines

===============================================================================
!
! start the run at time = 0
!
!NPLMAX         51
!NTPMAX         1001
T0             0.0E0
TSTOP          3.6525E7
DT             36.525E0            ! stepsize is 1/10 year
PL_IN          pl.in
TP_IN          tp.in
IN_TYPE        ASCII
ISTEP_OUT      10000
BIN_OUT        bin.dat
OUT_TYPE       XDR4
OUT_FORM       EL
OUT_STAT       NEW
ISTEP_DUMP     10000
!J2             0.01
!J4             0.001
CHK_CLOSE      yes
CHK_RMIN       1.0
CHK_RMAX       60.0
CHK_EJECT      -1.0
CHK_QMIN       1.0
CHK_QMIN_COORD HELIO
CHK_QMIN_RANGE 0.5 100.0
ENC_OUT        enc.dat
!EXTRA_FORCE    no
BIG_DISCARD    yes
!RHILL_PRESENT  no
===============================================================================

In the parameter file, portions of input lines following an exclamation point ('!') are treated as comments and ignored. Thus parameters NPLMAX, NTPMAX, J2, J4, EXTRA_FORCE and RHILL_PRESENT are not seen by the code and take their default values, -1, -1, 0.0, 0.0, .FALSE. and .FALSE., respectively.

Every ISTEP_OUT timesteps, the code outputs various quantities to the binary file specified by BIN_OUT. The exact data (and format thereof) written to this file are specified by OUT_TYPE and OUT_FORM.

If
     OUT_TYPE = "REAL4", data is written in 4-byte native Fortran binary format
              = "REAL8", data is written in 8-byte native Fortran binary format
              = "XDR4",  data is written in 4-byte XDR format
              = "XDR8",  data is written in 8-byte XDR format
     OUT_FORM = "EL", osculating orbital elements are written
              = "XV", heliocentric position and velocity components are written
              = "FILT", TBD filtered values are written (not yet implemented)

Every ISTEP_DUMP timesteps, the code dumps all of the information needed to resume the integration at that time in case of power failures or in case one wishes to resume an integration from its endpoint. The information is in 3 files called dump_pl1.bin, dump_tp1.bin, and dump_param1.dat (OR dump_pl2.bin, dump_tp2.bin, and dump_param2.dat). The first dump is written to the first set of files (set 1), and subsequent dumps alternate between the two. This is done so that at least one set of dump files will be preserved intact should the program die during the writing of the dump files. The format of dump_param#.dat is ASCII, that of dump_pl#.bin and dump_tp#.bin is XDR8, regardless of the format of the original input files. Note that TO in dump_param#.dat records the time of the dump and that OUT_STAT is changed to APPEND, so that these files can be used to restart a stopped integration. Depending on the situation, one may wish to increase TSTOP in order to extend an integration.

Planet file (PL_IN): This file contains all the initial planet data

The code requires units in which the gravitational constant G is unity. Any combination of lengths, masses, and times that keeps that true is OK. For example, one could use lengths specified in AU and time in days, thus forcing the Solar mass to be approximately 2.96E-4. Alternatively, one could use units in which lengths are in AU and the Solar Mass is unity, but then the orbital period of a test particle at r = 1 AU would be 2PI. A third useful set of units has lengths specified in AU and time in years, yielding a Solar mass of 4PI^2. The format is simple:

first the # of bodies on the first line (INCLUDING the Sun), then

3 lines for the Sun giving

• ID and mass on the first line
• heliocentric x,y,z on the next line and
• heliocentric vx,vy,vz on the third (NOTE: x,y,z and vx,vy,vz for the Sun MUST be 0!!)

3 (or 4) lines for each subsequent massive body giving

• ID and mass (and Hill sphere radius if RHILL_PRESENT = .TRUE.) on line 1
  • (planet radius on line 2 only if CHK_CLOSE = .TRUE.)
• heliocentric x,y,z on line 2 (or 3 if CHK_CLOSE = .TRUE.)
• heliocentric vx,vy,vz on line 3 (or 4 if CHK_CLOSE = .TRUE.)

If PL_IN is in XDR8 format, the ordering of the data in the file is identical to that described above even though line numbers lose their meaning. The ID for each body is an integer tag to help the code identify the body, thus no two planets can have the same ID. Furthermore, no planet can have the same ID as a test particle.

Test particle file (TP_IN): This file contains all the initial test particle data

In the same units as PL_IN, the first line is the number of test particles. The test particles are assumed to be massless, so for each particle there are 3 lines giving

• test particle ID on the first
• heliocentric x,y,z on the second and
• heliocentric vx,vy,vz on the third

No test particle can have the same ID as any planet or any other test particle.

Binary output file (BIN_OUT): This file, if defined in the input parameter file, contains snapshot frames of the system every ISTEP_OUT time steps, starting with the system prior to the first time step.

The data stored and its format are determined by the OUT_FORM and OUT_TYPE parameters defined above. Generally, each frame consists of a header record containing the time, the number of planets in the frame, the number of active test particles in the frame, and an integer identifier as to what type of data is being stored. Next, data for the planets are written. These are the ID, the mass, and the six quantities specifying its heliocentric orbit, either the osculating orbital elements or the heliocentric positions and velocities. Finally, data for the active test particles are written. These are identical to the planet data except that the mass, which is zero, is omitted. The tool tool_follow in the $(SWIFTER_HOME)/tool subdirectory shows how to access the data in this file. This particular tool will output to an ASCII file the data for a given body (specified by the ID value given by the user at the command line when the tool is run).

Discard file (discard.out): This ASCII formatted file, which has a fixed name in the current release, stores information on all discarded planets and test particles. First, the time, the number of bodies discarded this time step, and the value of BIG_DISCARD is output. Then, for each discarded body, it stores -1 if the body has been discarded, +1 if the body has been added (this only occurs for newly merged bodies in SyMBA), the body ID, an integer code (defined in module_parameters.f90) giving the reason for the discard, and the heliocentric position and velocity components of the body at the time of discard. If BIG_DISCARD is defined to be yes in the input parameter file, then similar data for the remaining active planets are also output after the discarded bodies.

Encounter file (ENC_OUT): This file, if defined in the input parameter file, contains output information on all close encounters that occurred during the run. (NOTE: this file is currently only used by the swifter_rmvs and swifter_symba integrators.

The tool tool_encounter_read in the $(SWIFTER_HOME)/tool subdirectory shows how to access the data in the encounter file.
For each encounter, the time of encounter, the ID's of the two bodies, their masses, heliocentric positions and velocities are given.
The format of the encounter file is determined by the OUT_TYPE parameter.
The tool_encounter_read tool will output all of the encounter data to an ASCII file.

XDR FILE FORMAT

XDR is a platform independent binary file format.
That is, it writes binary files that any machine can read.
So, the user does not have to worry about things like big endian vs. little endian.
We have adopted a package for reading and writing XDR called FXDR (for FORTRAN XDR).
FXDR was written by David W. Pierce of Scripps Institution of Oceanography.
Information about FXDR can be found at:

http://meteora.ucsd.edu/~pierce/fxdr_home_page.html

We are currently using version 2.1c.

EXAMPLE

An example set of input files and an associated README file can be found in the $(SWIFTER_HOME)/example subdirectory.