Swifter-OMP
Author: David A. Minton. Based on Swifter by David E. Kaufmann and Hal Levison
Date : 19 July 2010
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-OMP integrator supports a second class of massive bodies whose masses are less than some user-specified value, MTINY.
SyMBA-OMP is parallelized using OpenMP. Compile with OpenMP support and set your environment variable OMP_NUM_THREADS
to the number of threads you wish to run on prior to starting.
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:
- Wisdom-Holman Mapping (WHM). This is the N-body mapping method of Wisdom & Holman (1991; AJ, 102, 1528).
- 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).
- 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).
- 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).
- 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).
- 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).
- 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 shellmake
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 isar
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:
- 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.) - 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. - The FXDR library (
libfxdr.a
) will be regenerated from source, tested, then installed in$(SWIFTER_HOME)/lib
. - Any Swifter drivers located in
$(SWIFTER_HOME)/main
will be compiled and linked withlibswifter.a
andlibfxdr.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
) - 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.
make mod
: recompiles the module source code located in$(SWIFTER_HOME)/module
, replaces these objects withinlibswifter.a
, located in$(SWIFTER_HOME)/lib
, and installs the.mod
files to$(SWIFTER_HOME)/lib
make lib
: rebuilds the entire Swifter librarylibswifter.a
from source, replacing all the old archive membersmake libdir
: rebuilds only the libswifter.a source code LOCATED IN THE CURRENT WORKING DIRECTORY, replacing the old archive membersmake fxdr
: rebuilds the FXDR library from source, tests it, and installs it to$(SWIFTER_HOME)/lib
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
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
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
make clean
: removes any soft links to the filesMakefile
andMakefile.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 filefxdr.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.
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.
)
- (planet radius on line 2 only if
- 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.