An Introduction to the ST3D Numerical Relativity Code
Joan Massó, Mark Miller, Malcolm Tobias
March 21, 1997
Formulation
- Einstein Equations in 3+1 Dimensions
- Hyperbolic System
- First Order Flux Conservative
Numerical Techniques
- Parallelized FORTRAN using MPI
- Finite Difference Discretization
- Operator Splitting for First Order System
- MacCormack Method
Applications
- Gravitational wave evolution
Click here for a postscript version of the Documentation
The ST3D code can be downloaded here.
These instructions are for running on the Cray T3D at the Goddard
Space Flight Center. To run the code, carry out the following steps:
- Set the environmental variable NPES to the number of processors
you plan on using (i.e. type 'setenv NPES 1' for one processor).
- Type 'make'.
- Edit the inputxxx.par file, where xxx is the number of
processors you plan to use.
- Edit the batch script (titled ``batch'') and make sure that it
points to the directory in which your executable lives, and also that it
sets the variable mpp_p to be greater than or equal to the number of
processors you plan on using, and that the environmental variable NPES
agrees with that in step 1.
- Submit the job to the queuing system by typing ``qsub
batch''.
Introduction
ST3D is a Numerical Relativity code that
solves the full Einstein Equations in 3D using a
finite differencing conservative scheme.
ST3D is the Cray T3D version of H3expresso, a
``concentrated'' (``espresso'') and fast (``express'') version of
a general 3D hyperbolic code that solves the general
relativistic Einstein's Equations. The hyperbolic formulation of
the Einstein equations, by using a special class of
gauge choice, casts
the Einstein equations into a First Order Flux Conservative
Hyperbolic system and thus allows the use of advanced numerical techniques
[1]. This code is the first step in the numerical study
based on these hyperbolic Einstein equations.
The Harmonic Formulation of the Enistein Equations
ST3D is based on the work of Bona and
Massó[2] that casts the Einstein equations in an
explicitly first order, flux conservative, hyperbolic form.
In the standard 3+1 formulation of the Einstein theory, the spacetime
metric is given by
The goal of the numerical evolution is to obtain the 6 dynamical
variables 
as functions of the spatial coordinates
and time t, given a choice of the ``lapse
function'' 
and the shift functions 
, for a
given set of initial data, namely, the metric functions and
their first time derivatives at time t=0.
It was
shown in Ref. [2] that if one chooses the lapse
to satisfy the
harmonic slicing condition
one can write the Einstein
evolution equations as a hyperbolic first order system of balance laws.
In vacuum and in the case of zero shift, the equations take the form
The 
quantities are proportional to the extrinsic
curvature.
The connection coefficients 
, defined as
are constructed from the
first derivatives of the metric
These derivatives are evolved using Eq. (
5 ).
Eq. (
8
) is only used in the initial slice. Similarly, the
derivatives of the lapse are used on the initial slice to construct
and to derive the initial values of the momentum
constraint related variables:
These variables are evolved using Eq. (
6
)
while Eq. (
9
) is used to compute the 
during the evolution.
Basic Numerical Methods
Since in the hyperbolic formulation the Einstein equations take
the form of first order, flux
conservative, hyperbolic equations, a wide variety of modern numerical
techniques can be employed. Operator splitting allows
for the principal part of the system to be treated as a flux
conservative first order system. In this code, a flux conserving
MacCormack method
is used for the principal part of the evolution system. More advanced
differencing schemes and operator splitting treatments are allowed in the
formulation and can be added on to the code without difficulty.
Initial Data
The initial data in this code represents low amplitude quadrupole
gravitational waves [3, 4].
These waves are solutions to the linearized Einstein constraint equations.
For this problem, there are
2 user-specified the parameters, A the amplitude of the wave
(which should be much smaller than one) and m the azimuthal mode
number for the wave. For the formula relating the metric
to the parameters A and m and a more general
discussion of numerical evolutions of these waves, see
[5] and the references cited therein.
Files
Fortran Files
A list of all fortran subroutines can be found in SUBROUTINES IN THE ST3D CODE subsection.
Makefile
Makefile to build the code.
inputxxx.par
A parameter file that defines the parameters to us when running on xxx
number of processors. In this file, nx, ny, and nz stand for the number of gridpoints in each direction
on each processor, while px, py, and pz stand for the number of
processors being used in each dimension. Therefore the total
dimension of the gird is nx px X ny py X nz pz.
batch
A sample batch script for submitting jobs to the Cray T3D.
Variable list
All variables are 3-D arrays.
-
x, y, z, r: the x, y, z, and radius coordinate for
each grid point
-
psi: conformal factor 
-
alp: lapse function
-
cux, cuy, cuz: momentum constraint related
variables
-
uxx, uxy, uxz, uyy, uyz, uzz: Contravariant
(``upper'') metric 
-
gxx, gxy, gxz, gyy, gyz, gzz: Covariant metric
(note: these are not evolved, but computed from the contravariant
metric)
-
rg: square root of the determinant of the
covariant metric
-
qxx, qxy, qxz, qyy, qyz, qzz: Extrinsic curvature
related variables
-
dxuxx, dxuxy, dxuxz, dxuyy, dxuyz, dxuzz: x-derivatives
of the contravariant metric
-
dyuxx, dyuxy, dyuxz, dyuyy, dyuyz, dyuzz: y-derivatives
of the contravariant metric
-
dzuxx, dzuxy, dzuxz, dzuyy, dzuyz, dzuzz: z-derivatives
of the contravariant metric
Subroutines
- Analywave: Initialize the metric tensor with a linearized
Eppley style packet of incoming/outgoing waves. The metric
is computed in spherical coordinates and then transformed to
cartesian.
- Boundaries: Provide boundary conditions to all the metric
variables that
have been evolved in the interior. ST3D simplifies the
process by enforcing the same boundary condition on all
quantities thru a call to Onebound for each variable.
- Centralderiv: Compute the derivatives of the metric tensor
with a central finite differencing scheme that
assumes a regular grid spacing. The finite
differencing is implemented using triplet notation.
- Curelation: Compute the vector
based on the metric and
derivatives and the algebraic relation that holds at
the initial time.
- dummyentry: a break point for the compilers
- Fluxes: Compute the fluxes of the equations
- Goodbye: A goodbye.
- H3: Main code: it does the evolution by calling the routines
to set the initial data at iteration 0 (ST3D restriction)
and then calling the method (ST3D only implements Macormack)
and the boundaries at every time step. ST3D does not
have I/O so no output routine is called. Instead, some minimal
information is given at every iteration.
- Infoheader: Display a header for the Information routine
output
- Information: Provide some info on the evolution as it runs.
ST3D only shows the current iteration, the time and
the minimum and maximum values of the variable sent as an
argument.
- InitH3: Wrapper for the main code: i) reads the parameters
thru readinput, ii) calls the main code H3, iii) gives timing
info
- Initial: Initialize the metric structure with the initial data.
See the comments inside the routine for specific info, as some
special routines have to be called to supply initial data.
- Invert: Invert the metric tensor to compute covariant from
contravariant and vice versa. The square root of
the determinant is also computed.
- Method: Evolve the metric quantities one time step
using a simple finite differencing CFD-like method: MacCormack.
The finite differences are implemented using triplet notation
and assuming a regularly spaced grid.
- Onebound: Simple boundary conditions on one
variable.
- Readinput: Read from the text file
``input.par'' the control parameters.
Stops if the file is not found or the format is the right
one.
- Sources: Compute the sources of the equations.
- timeshow: Show some statistics based on the cpu
time and the spacetime grid.
- timestart: Start the timer..
- timestop: Stop the timer.
- Welcome: A welcome.
References
- 1
- See e.g., C. Bona, J. Massó, E. Seidel, and
J. Stela, Phys. Rev. Lett., 75, 600, (1995); A. Abrahams, A. Anderson,
Y. Choquet-Bruhat and J. W. York, Phys. Rev. Lett. 75, 3377
(1995), and references cited therein.
- 2
-
C. Bona and J. Massó, Phys. Rev. Lett. 68, 1097 (1992).
- 3
-
K. Eppley, in Sources of Gravitational Radiation, edited by L. Smarr
(Cambridge University Press, Cambridge, England, 1979), p. 275.
- 4
-
S. Teukolsky, Phys. Rev. D 26, 745 (1982).
- 5
-
P. Anninos, J. Massó, E. Seidel, W.-M. Suen and M. Tobias,
submitted to Phys. Rev. D.
This distribution of the code runs on the Cray T3D.
Malcolm Tobias
Physics Department
Washington University, St. Louis, MO 63130
mtobias@wugrav.wustl.edu
(314)-935-6249
The H3expresso code was originally developed by Joan Masso, with
the initial data routine by Malcolm Tobias, MPI
support and extensions by Rob Gjersten and Marc
Nardulli, and automatic cache optimization by Paul Walker.
The T3D version ST3D has been further developed by Tom Clune, Ian Foster,
and Mark Miller.