silburt / swifter Goto Github PK

*******************************************************************************
                               Swifter README
                   (A Preliminary User's Guide to 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.

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 2*PI.  A third useful set of
units has lengths specified in AU and time in years, yielding a Solar mass of
4*PI^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).

A.S. note: tool_follow is an .exe in bin/ (like the integrators/drivers). You need to copy/move it to your example folder, run it, feed it param.in, and then use frequency of 1 (lowest integer value). But, need to do that for each body, so that could take a long time. Probably need to create a script that will do it all in one go for me.  

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.