Final Report for NASA CAN NCCS5-153


A Multipurpose Code for 3-D Relativistic Astrophysics and Gravitational Wave Astronomy: Application to Coalescing Neutron Star Binaries



Overview of the project

The construction and validation of a multi-purpose 3D code for computational general relativistic astrophysics that makes efficient use of massively parallel computers are major challenges in computational physics, and are important to high energy astronomy and gravitational wave astronomy. We took on this challenge starting fall of 1996 with the support of the NASA Round-2 Grand Challenge Program. We used the simulation of coalescing neutron stars as a driver for the development of a massively parallel code that is capable of solving the Einstein equations coupled to the general relativistic hydrodynamic equations. We have now successfully concluded the project. The major achievements in this project are:


The application of the code to the headon and inspiral collisions with different equations of state, and the validation of
such applications are now in progress.  Neutrino transport and magneto-hydrodynamics, which are important for many
neutron star processes, are NOT included in this code.  They are items for future grand challenge projects.

Back to the top of page

Significant Milestones of the Project

                                          Achieved February 15, 1997                                           Achieved February 15, 1997                                           Achieved May 10, 1998                                           Achieved October 15, 2000

In the following we describe the project in more details.

Back to the top of page

Background

General Relativistic Astrophysics--astrophysics that involves strong and dynamical gravitational fields requiring the Einstein theory of relativity for its understanding--is becoming an exciting area of research, due to the large amount of data in high energy astronomy, and to the great promise of gravitational wave astronomy. In particular, Computational General Relativistic Astrophysics, with the recent stunning increases in computer power, may hold the key to the understanding of many observations in high energy astronomy and gravitational wave astronomy, two major frontiers of astronomy in the new century.

With the advent of the large amount of observational data in high energy astronomy and gravitational astronomy, there is an urgent need for the capability to carry out large scale simulations in relativistic astrophysics.  High energy (X- and Gamma-ray) astronomy is driven by X- and Gamma-ray satellites from NASA and other space agencies, presently and/or soon to be in orbit  [1]. Gravitational wave astronomy is driven by the gravitational wave observatories scheduled to be on line within the next few years [2]. The former probes the most energetic events of the universe, while the latter opens a completely new window into our universe. Both kinds of observation concern astrophysical phenomena in strong and dynamical gravitational fields. If we are to understand the observational data generated by non-linear and dynamical gravitational fields, detailed modeling taking dynamic general relativity into account must be carried out.

The numerical treatment of strong field astrophysical phenomena, such as those involving neutron stars and black holes, requires the integration of the Einstein theory with other tools of astrophysics, in particular relativistic hydrodynamics and nuclear astrophysics. The equations of the Einstein theory describing the dynamics of the spacetime itself form probably the most complicated set of equations in physics, with rich mathematical structure. It is highly nonlinear with literally thousands of terms coupling together tensor fields, vector fields, and scalar fields. Its integration with relativistic hydrodynamics and nuclear astrophysics for simulations of realistic astrophysical phenomena in full 3+1 dimension (3D) calls for state-of-the-art parallel computational, visualization and software engineering technologies. We took on this challenge starting fall of 1996 with the support of the NASA Round-2 Grand Challenge Program, using the simulation of the coalescence of neutron stars as the driver application of the code development effort.

The project is a collaboration between six institutes: the Argonne Lab, Livermore Lab, Max Planck Institute, SUNY at Stony Brook, University of Illinois and Washington University. The computer science groups at UI and Livermore led by Saylor and Ashby are responsible for developing parallel elliptic equation algorithms. The group at Argonne Lab led by Foster (together with Parashar now at Rutgers) is responsible for developing meta-computing capability and adaptive mesh refinement treatments. The nuclear physics group at Stony Brook led by Lattimer is responsible for developing realistic equation of state treatment. The astrophysics group at the National Center for Supercomputing Applications (NCSA) of UI led by Norman and Swesty is responsible for developing an artificial viscosity based hydro code. The Max Planck Institute group led by Seidel and Schutz, and the Washington University group led by Suen and Will are responsible for developing the spacetime code. The Washington University group is also responsible for the development of high resolution shock capturing hydrodynamic codes (Newtonian hydro and general relativistic hydro), general relativistic initial data preparations and neutron star simulations.

Back to the top of page

Relativity, Astrophysics, and Computer Sciences Aspects of the Code Development Effort

One major difficulty in the code development effort is that the multidisciplinary nature of the project requires inclusion of vastly different expertise, from computer and computational sciences to the different physics disciplines involved, including relativity, astrophysics, nuclear physics and computational hydrodynamics. This calls for building the simulation code on top of a collaborative infrastructure that enables the simultaneous development of the code by a large and geometrically distributed team of developers. Moreover, as we are aiming at a multi-purpose code that will enable a large variety of astrophysical applications, it is also important that the code is flexible and easy to use by astrophysicists with varying degrees of computational expertise. The Max Planck Institute/Wash U team in the NASA Grand Challenge collaboration, together with others in the community, have constructed a computational software infrastructure  called the "Cactus Computational Toolkit'', that enables collaborative code development. We built our code for solving the Einstein equations on this infrastructure.  It has a modular code structure that allows related code modules to be arranged into ``arrangements'', and be plugged in as ``thorns'' to the core part of the code (the ``flesh'').  This enables the interoperability of the code modules. The ``flesh'' contains the parallel domain decomposition layer based on message passing interface (MPI) for solving partial differential equations on massively parallel machines, including multigrid and other solvers for elliptic systems, parallel I/O, checkpointing, visualization tools and fixed mesh refinement. We have also constructed a consistency test suite library (to make sure that new thorns will not conflict with other parts of the code), and various code development tools that provide a complete environment for code development and testing. For a review of the collaborative infrastructure of the code, see [4].

Our spacetime evolution code contains the ADM [5], conformal [6,7,8,9] and hyperbolic formulations [10] of the Einstein equations. Many different explicit finite differencing schemes (MacCormack, Iterative Crank-Nicholson, and staggered and non-staggered Leapfrog) are included as options for the numerical integration of the Einstein evolution equations. It has the capability of solving the coupled set of the constraint equations for initial data. Different boundary conditions (static, flat, radiative and Robin) can be used. Various time slicings (maximal slicing and various algebraic slicings) and spatial coordinate conditions (minimal distortion and related shifts [11]) are included as options for evolving different spacetimes.

We have also constructed a code for evolving the general relativistic hydrodynamic (GRHydro) equations coupled to the Einstein equations (with perfect fluid source). Our GRhydro treatment makes use of the high resolution shock capturing (HRSC) schemes (for review, see [12]). The curved spacetime HRSC treatment is based on a spectral decomposition (eigen field structure)  we obtained [3] for a general metric, building on earlier work [13]. The choices of HRSC schemes in the code presently include the 2nd order flux limited flux split, Roe, Marquina, ENO and HLLE schemes. The GRHydro evolution is coupled to the spacetime evolution in a second order accurate manner in both space and time. In both the spacetime and GRHydro treatments we have included different treatments as options. Considering the complexity and generality of the code, along with the fact that the solution space of the set of differential equations is largely unexplored, it is imperative that we have the capability to compare results based on different mathematical formulations of the equations, and different numerical schemes for solving them. For full details, see papers published with the support of the grant cited in the last section of the report.

A largely unexplored area of research in computational general relativistic astrophysics is the gauge degrees of freedom in the simulations of dynamic spacetimes, especially those driven by hydrodynamic flow. In our studies, we find that suitable coordinate conditions are crucial for the long term stability of simulations involving neutron stars, black holes and gravitational waves. We have investigated various time slicings of the spacetimes, including various slicings in the class of "generalized harmonic slicing" and the maximal slicing condition in the presence of matter sources. Investigations of various spatial coordinate conditions, e.g., the minimal shear and minimal distortion conditions have also been carried out. Many of the coordinate conditions are enforced by solving elliptic equations, which call for efficient parallel elliptic equation treatments. In GR3D_EOS, two elliptic equations modules are included: "BAM" is specially developed for our spacetime code and has recently been extended to handle matter source terms and vector elliptic equations. "PETSc" is a general purpose elliptic solver developed by the Argonne PETSc group, and linked to our spacetime code.

We have also made significant achievements in the Computer Science aspects of our project, including code performance (see below ) and meta-computing. We integrated the code with the Globus toolkit to assist us in executing the code on distributed resources. Globus consists of a set of tools for tasks such as communication, resource location, scheduling, access to remote data, and authentication (see "Globus metacomputing toolkit" ). We use several of these low-level services to execute our code using a Graphical Resource Co-allocator (GRC). The GRC allows us to specify in abstract terms the resources required for the application (e.g., architecture type, number of processors, amount of memory) and also the applications themselves (executable, arguments). For example, in the figure below, we show the specifications used to achieve remote job submission and rendering. The subpanel labeled "Compute1" specifies that the astrophysics code is to run on 64 nodes of any computer running the IRIX operating system, while the "Compute2" subpanel indicates that the visualization component should run on a specific display device (in this case, an ImmersaDesk virtual reality system located at the High Performance Distributed Computing Conference in Chicago). The subpanel to the right indicates that the two computers must be connected by a network with at least 10 Megabits per sec of bandwidth between the two computers and a latency of at most 20 milliseconds.

Before a user starts to interact with Globus, they authenticate themselves to the Globus environment. The user can then execute applications on any computer running Globus to which the user has access, without re-authentication. After authentication, the user starts the GRC, enters a description of their application, and asks the tool to find candidate sets of resources for the application. The tool queries the Globus information service to find sets of candidate resources and displays information about the candidates. The user selects one candidate and asks the high-level tool to execute the application on those resources. The GRC uses the common Globus interface to local scheduling systems to request that the local scheduling systems controlling the selected resources start the applications. The executable and data files are then staged to the remote systems, the applications are run, and any result files are staged to where the user wants them.

Together with the general relativistic code, we have also developed two Newtonian gravity codes for comparison
to the relativity code.  One uses the artificial viscosity scheme (by the NCSA team) while the other uses the same high resolution shock capturing scheme as in the relativistic version (by the Wash U team).

Back to the top of page

Code Validation

Considering the complexity of the code development effort, it is imperative to carry out many tests to insure the fidelity of the discretization to the original differential equations. This is especially so for our GR3D_EOS spacetime-GRHydro code (or any other code of this nature) because:

1. the solution space of the couple Einstein-GRHydro equations is to a large extent an un-explored territory,

2. the complexity of the code: the set of differential equations we are attempting to solve consists of very complicated, coupled partial differential equations involving thousands of terms,

3. the generality of the code: with the many different options in choosing the formulation, the solver, the time and space coordinate choices, the numerical methods and the boundary conditions,

4. the involvement of so much computational science infrastructure for massively parallel computation (one tiny memory allocation error can change the final result!).

Code testing and validation have been the major part of our work in the last few years. Our code has been subjected to the following classes of tests/checks:

I. Comparison with analytic solutions and model testbeds.

One obvious first step in validating the code is to carry out simulations of known solutions. We have carried out studies of the following systems [3,14]:

1. Riemann shock tubes, including ultra relativistic shocks at a random angle with respect to the Cartesian grid of the code.

2. Friedmann-Robertson-Walker cosmology tests.

3. Oppenheimer-Synder collapse.

4. Equilibrium self-gravitating compact (TOV) stars. It turns out that evolving a TOV star in its equilibrium configuration is quite non-trivial, especially for stars with a stiff equation of state. In such cases, the density profile near the surface of the star changes abruptly (the first derivative of the density is discontinuous at the surface), and the second derivatives of the metric has a kink there. The second derivatives of the metric enter explicitly in the evolution equations. The sharp changes demand high resolution and induce instability in numerical evolutions. For our treatment of this problem, see [3].

5. Evolution of rotating neutron stars in equilibrium, using shift vectors.

6. Relativistically boosted TOV stars. This is a particularly valuable testbed. The evolution of the TOV solution boosted at an angle with respect to the Cartesian coordinates (moving across the grid) involves ALL terms in the spacetime-GRHydro evolution equations [3]. We urge all researchers interested in constructing a spacetime-GRHydro code to go through this test!

Unfortunately, for the spacetime-GRHydro equations, only in cases with high degrees of symmetry one can have analytic solutions to compare with, and systems with symmetries can be less effective in revealing problems in the code. (Typically symmetric problems involve less terms in the equations. On the other hand, the advantage is that one can easily pick out "bugs'' that explicitly break the symmetry, see below.) This lack of analytic solutions for comparison highlight our ignorance of the solution space of the equations, and the need for performing many different classes of tests.

II. Comparison with lower dimension/linear perturbation/Newtonian results.

As there are not many known analytic solutions of the coupled spacetime-GRHydro equations representing dynamical spacetimes with matter, we want to compare our results to other numerical results which are established/reliable/accurate. We numerically evolve spherical and axisymmetric systems with our 3D Cartesian code, and compare to results of a 1D code in spherical coordinates, a 2D code in axisymmetry, and a 2D axisymmetric implementation of the 3D code. The lower dimensional simulations can be carried out in much higher resolution, and/or specially chosen coordinate system adapted to the symmetry to reduce/avoid difficulties coming from the gauge degrees of freedom in the full 3D code.

We have a spherical coordinate code that solves the spacetime equations in standard ADM formulation in polar slicing and radial gauge, coupled to the GRHydro equations solved with the same high-resolution shock capturing treatments as in our 3D Cartesian code (only the flux split and Roe treatments are presently included in the 1D code; for details of the spacetime numerical treatment in the 1D code, see [4]; for details of the GRHydro treatments see [3]). Comparisons of spherical collapses, static TOV stars and radially perturbed TOV stars have been carried out, and more tests are in progress.

For comparison to axisymmetric systems, we have two approaches. One is to construct axisymmetric codes in cylindrical coordinates $(\rho, z, \phi)$, and impose an explicit $\phi$ independence. We have such a vacuum spacetime evolution code which we used extensively in black hole and gravitational wave studies (see, e.g., [15,16]). GRHydro treatment has recent been added to this code [17]. We have carried out various comparisons of this code to the 3D Cartesian code. Testbeds carried out/being carried out include: (i) Radial pulsation of stable stars, (ii) axisymmetric perturbation of rotating stars, (iii) collapse of unstable neutron stars to black holes; and (iv) Migration of unstable neutron stars to the stable branch. These comparisons are reported in [14].

However, one disadvantage of using the cylindrical coordinates is that in longer time evolutions it tends to induce numerical instability along the $z$ axis, with singular terms containing factors of $1/ \rho$ in the dynamical evolution equations. Due to this axis instability problem, and the advantage of building our existing computational infrastructure for easiness in making direct comparison, we are presently constructing an axisymmetric treatment on top of our 3D Cartesian code [18]: For axisymmetric systems, one can evolve just one 2D slice in the 3D Cartesian grid of the 3D code, with the boundary conditions on the two neighboring slices set by symmetry. In this way, the resolution can be increased by more than an order of magnitude compared to the 3D code, while making use of the same parallel computational infrastructure. Since the calculation is based on the Cartesian coordinates in this treatment, there is no numerical axis instability. This treatment has been used in vacuum spacetime studies [18] and is being extended to spacetime-GRHydro studies. There are added complication in treating shocks in GRHydro in this approach. We have added this axisymmetric treatment to our Newtonian gravity-hydrodynamic 3D code, and are testing the propagation of shock waves in this treatment without the complication of general relativity.

We have built a Newtonian gravity-hydrodynamic 3D code  on the same Cactus computational infrastructure as the general relativity code to enable direct comparisons. It uses the same high resolution shock capturing (HRSC) schemes as our GRHydro code. Studies of shock propagation, static and boosted neutron stars, head-on collisions and inspiral coalescences of binary neutron stars have been carried out [19].  There is also an independent artificial viscosity based Newtonian gravity-hydrodynamic code in the Grand Challenge collaboration using a different parallel computational layer that we can compare our results to  [20]. Not only do the Newtonian codes provide direct comparisons to the general relativity code in the weak gravity limit, it also provides a way to identify the effects due to relativity and the point that such effects become important.

Another direction in which we are carrying out tests to validate our general relativity code is to compare it to results of perturbation calculations. Such comparisons have been successful in black hole collision studies [21,22]. We studied radial and nonradial perturbed TOV stars, extracted the gravitational wave output, and compared the frequencies of the stellar oscillations and gravitational waves to the frequencies obtained in linear perturbation theory [14].

III. Comparison with other treatments in the code.

We put much effort in building a flexible code with many options as discussed above. The many finite difference forms of the equations, choices of numerical evolution schemes, time slicing and coordinate choices, and boundary treatments enable many different combinations and cross checks. Not only does this flexibility enables treating different systems with methods most suitable for each of them (making our code a multi-purpose code), it also provides comparison tests for any particular choice. Any new treatment that we add to the code, e.g., a new hydrodynamic scheme, could be conveniently tested by comparing to existing treatments. We have set up ``standard test-suites'' that we insist new schemes and all upgrades of the code to pass through. There are over 25 test-suites on evolutions of gravitational waves, black holes and neutron stars that we use regularly.

IV. Short term convergence tests.

Convergence test is an effective tool in flushing out ``bugs'' in the construction of the finite differencing evolution code.

The point of short term convergence test is to check, with decreasing grid separations, whether the finite difference solutions obtained converge to a definite solution, with the expected convergence rate of the numerical scheme. Typically the differential operators in the evolution equations are represented by finite difference expressions in a second order fashion. (For example., $df/dx$ at the $j^{th}$ grid point is represented by $(f(j+1)-f(j-1))/(2 \Delta x)$, where $\Delta x$ is the grid separation, and f(j) represents$f$ at the $(j)^{th}$ grid point etc. Second order accuracy is often achieved through symmetrization, i.e., the symmetric use the two neighboring grid points about the $j^{th}$ point in this case.) In our 3D Cartesian code, except at special points (of measure zero, see below), the finite difference evolution equations are second order accurate in both space and time. This means that in the limit of the spatial grid separations $\Delta x , \Delta y,\Delta z$ going to zero (with a corresponding decrease in $\Delta t$ to satisfy the Courant stability requirement (see, e.g., [23])), the finite difference solutions of the evolution equations evaluated at a fixed time should converge to the limiting solution in a $\Delta x ^ 2$ manner, for small enough $\Delta x$. (That is, for an expression that should converge to zero, its value should decrease by a factor of four for each factor of two decrease in $\Delta x$.)

While in principle in convergence testing a solution one should check for convergence of all variables involved, with the large amount of variables involved in a spacetime-GRHydro simulation, we can only check a subset of the variables or certain combinations of them for convergence. Making good choices helps. The hamiltonian and momentum constraints (the time-time and space-time components of the Einstein equations) are among the four quantities that we check in all cases in our free evolution code. Not only are the constraints convenient to check with their known limiting values (namely, they converge to zero), they are also sensitive indicators of the accuracy of the simulation, as they involve the highest derivatives used in the code. The trace of the extrinsic curvature $K$, giving some summary information of the extrinsic curvature is also often a good choice (except in those simulations that we explicitly force $K$ to zero on each time step of the numerical evolution, e.g., in our so-called ``AF2'' and ``AFK'' treatments [24] in the conformal-trace free (CT) formulations [6,7,8,9]). For a strongly gravitating isolated object, the radial-radial component of the metric, formed using the evolved Cartesian metric components, is also often a good choice. Other quantities often used in our tests include the curvature invariants, and the Newman-Penrose quantities in gravitational wave systems.

We carried out extensive convergence tests in our code construction effort. In Ref. [3], twelve combinations of spacetime-GRHydro treatments were convergence tested in many model problems.

A few finer points to note concerning convergence tests: (i.) While our evolution code uses finite differencing schemes which are second order accurate in space and in time at a generic point, some parts of our hydrodynamic finite differencing schemes, namely the flux limiters, are first order at points where some of the hydro variables attend extremium values, e.g., at the center of a star. We have confirmed the convergence rates with and without the flux limiters turned on. (ii.) Some boundary conditions we used are only first order accurate. The effects of lower convergence rate at the boundary can often be identified and considered separately, especially for the short term convergence tests considered here (in contrast to the long term convergence test to be discussed below). A caution is that in cases where elliptic operators are involved, e.g., simulations using maximal slicing with the trace of $K$ assumed in the evolution equations, lower orderness of the boundary conditions can affect the whole evolution in a way that cannot be easily separated. (iii.) In checking for convergence rate, not only one has to make sure that a fixed spacetime point is used in the different resolution runs, but also that some spacetime points at a certain set of resolutions may not be suitable for such studies. For example, in the case of a ``loused'' instability  [25] existing in the numerical solution, one has to pick time steps that allow the ``loused'' solutions to be compared, e.g., at the maximum value points of the various resolutions used. For more details, see [3].

V. Consistency check.

In contrast to convergence, which is a relation of the finite difference solutions at different resolutions to their continuum limit, consistency is about the relation between a finite difference equation to a differential equation. Consistency is to require that the finite differencing operator converges in the continuum limit to the correct differential operator. This is to be true for the finite difference operator acting on arbitrary smooth functions, which are not necessarily solutions to the equations.

The importance of consistency check in a numerical code is highlighted by the Lax theorem [27] (and its various generalizations [28]), which basically states that for a well-posed initial value problem involving a linear hyperbolic system (and its generalization), as long as the finite difference operator is consistent with the differential operator, convergence of the numerical solution implies the stability of it, and vise versa. Here, stability means that all fourier components of the finite difference solution are bounded during evolution. While there is no generalization of the theorem to cover directly the spacetime-GRHydro equations, which are quasi-linear and non-strictly hyperbolic, we expect nevertheless the three properties (consistency, convergence and stability) to be closely linked.

We carried out consistency checks for all finite difference operators in our spacetime-GRHydro code. To check a finite difference operator, we use analytic test functions that activate all terms of the operator. The results of the finite difference operator acting on the test functions are compared to the analytic evaluation of the differential operator acting on the same test functions, and the difference between the two should converge at a rate determined by the scheme used. Indeed consistency checks are the first test we carry out in the construction of any finite difference operator in the code. As the test functions are not necessarily solutions to the equations, consistency checks are relatively easy but powerful checks that can be used in all cases. This is in contrast to long term convergence and stability studies (see below) which are computationally expensive to carry out. As the three (consistency, stability and convergence) are closely related (see below), consistency tests are particularly worthwhile. We advocate their usage especially in spacetime-GRHydro codes, in which very complicated differential operators are involved.

VI. Residual evaluators.

The use of ``residual evaluators'' is, in a sense, an extension of the convergence tests. For key quantities obtained through complicated procedures, we often construct a residual evaluator for each of them. For example, for the initial data we constructed through the post-newtonian formulation, we also construct an independent finite difference version of the post-newtonian operator (independent of that used in the solving of the post-newtonian equations) which in the continuum limit should annihilate the initial metric and extrinsic curvature. With finite resolution, the annihilation will not be exact. We demand the residual

(1) to be smaller than what is required by physics (in the post-newtonian consideration we want the residual to be less than the next newtonian order neglected, this sets a requirement on the resolution used in the simulation), and

(2) to converge to zero in increasing resolution with the expected convergence rate. This guarantees the solution satisfies the differential equations it is supposed to solve.

We have constructed residual evaluators also for
(i) the quasi-equilibrium initial data generated using the ``conformally flat'' approach by our Meudon collaborators [29],
(ii) minimal distortion shift [11], and
(iii) lapse, extrinsic curvature and its time derivative in maximal slicing.
The use of residual evaluators have been effective in identifying ``bugs'' in various parts of the code.; A residual evaluator directly for the dynamical variables in evolution will also be added to the code.

VII. Long term convergence test.

While short term convergence tests as discussed above have now been quite widely used in numerical relativity, we would like to particularly advocate long term convergent tests. By long term convergence tests we mean the convergence of the solution throughout the numerical evolution until the point that the physical results are extracted. In particular, the physics results themselves must pass the convergence tests. This is in contrast to short term convergence tests which test the convergence of the finite difference solution at the first few, or first few tens, time steps. No numerical results should be reported without these results themselves having been convergent tested.

Long term convergence tests for EACH specific problem is essential for the spacetime-GRHydro system, which is complicated and highly non-linear. Unlike linear theory, the convergence of a non-linear problem depends on which solution we are looking at, and at what spacetime point. Short term convergence tests of a different model problem, and even for the same problem of interest, are not enough to provide confidence in the numerical results obtained at a later time of the simulation.

The major difficulty in carrying out long term convergence tests is that the convergence rate would be at the desired order only for a small enough $\Delta x$. To get to the convergence regime in a 3D simulation one often needs a lot of grid points, and that in turns requires a computer with a large memory size. For some problems this may not be possible even with the largest parallel computer available to us at NASA or the National SuperComputing Centers. Unlike short term convergence tests using known analytic solutions in which one can test the evolution code by restricting the computational domain to a small region to get a high enough resolution, in long term convergence tests typically one cannot move the outer boundary too close.

Even in cases where we have large enough memory size to get to the convergent regime, for our explicit finite differencing code a small $\Delta x$ implies a small time step and hence a large number of iterations are needed to evolve the system to the desired time. This could prove prohibitively expensive computationally. (We have carried out tests of our code on a 1024 node T3E-1200 with half a terabyte of memory; but it is not available for long time scale studies.). The computational restrictions mean that in many cases we cannot afford to have enough grid points to make $\Delta x$ small enough to put the solutions in the second order convergence regime in long time scale simulations. Finite difference errors of order $\Delta x ^3$ or higher are then not negligible. One may see that the solutions are only ``converging to converge'', when $\Delta x$ is reduced. In such cases, we can only confirm the convergence of the code by matching to a converging rate of $ a \Delta x ^ 2 + b \Delta x ^3 $ for some constants $a$ and $b$ (or even higher order terms). The confirmation of this requires carrying out the convergence study at many different resolutions.

Not only does one need to confirm that the physics results pass the convergence tests, one must also insist that the numerical solutions converge from the initial time all the way to the final time at which the physics results are extracted. For non-linear evolution equations, even smooth initial data could develop high frequency components and bring the numerical solution out of the convergence regime (error in higher order of $\Delta x$ not negligible). While this often leads to instabilities crashing the simulation, this is not always the case depending on the finite differencing scheme used. It could return to the convergence regime at some later time and yield fictitious physics results.

While the usefulness and importance of long term convergence tests cannot be over-emphasized, we should keep in mind that converge tests are not leak-proof. Firstly the Lax theorem (and its presently known extensions) is not strictly applicable to the spacetime-GRHydro system. But more importantly, in practice one can only carry out convergence tests with a limited range of $\Delta x$. This is especially true for long term convergence test. Hence in practice one cannot claim convergence with full confidence. Two complications related to this limitation of finite $\Delta x$ are: (i.) There can be a term dominating the truncation error which converges at the correct order, covering up a "bug'' that can only be revealed by convergence tests with a smaller $\Delta x$. (ii.) For physical systems whose short wave length components can be generated during evolution, convergence tests may not be useful if $\Delta x$ cannot be small enough to cover the shortest wave length generated. This is especially so if the short wave length features are not limited to specific spacetime points.

VIII. Stability with respect to numerical choices.

The numerical solution should not depend sensitively on the numerical parameters used in the simulation. The numerical choices we pay special attention to include:

(i) Position of the outer boundary. The boundary position of the computational domain could affect physics results through reflection of waves, inaccuracy in the fall off rate, and through solving of elliptic equations. While the choice of the conditions to impose at the boundary is determined by the physics (e.g., outgoing wave boundary condition or Robin condition at an appropriate fall off rate) and is not a numerical choice, the position at which we implement the condition is. One should make sure that the outer boundary is far enough out that the physics results do not sensitively depend on its position. For example, in studies involving the existence/location of an apparent horizon as in the head-on collision study described below, we made sure that position of the apparent horizons determined in the simulations converge with respect to our moving the boundary further out. (Again we see the connection between stability and convergence just like in the Lax theorem.  Only in this case it is not with respect to $\Delta x$, but with respect to the total number of grid points used with $\Delta x$ fixed.)

(ii) Use of different numerical schemes: while one scheme may be more suitable than another for a particular system, they should yield converging results with increasing resolutions. Here different numerical schemes could mean leapfrog verses iterative Crank-Nickolson in the spacetime part of the evolution, or Roe verses Marquina in the GRHydro part of the evolution. While one may not be able to afford checking many different combinations in the high resolution production runs, it is useful to check that at least the qualitative results are the same in lower resolution runs with a few combinations at a few different resolutions.

(iii) Tolerance used in solving elliptic equations: Elliptic equations are used in the construction of initial data, and often in choosing the lapse and the shift. Too stringent a tolerance requirement would make the simulation computationally too expensive, while too loose a requirement can affect convergence tests, but not always. For example a loose tolerance elliptic solve for initial data may yield extra ``radiation content'' in the simulation that cannot be picked up by a convergence test.

There are many other numerical choices like the courant factor (ratio between $\Delta x$ and $\Delta t$ in the evolution), explicit enforcement of one or the other constraints (not only the hamiltonian and momentum constraints, but also constraints coming from gauge requirements, and the hyperbolic and conformal requirements), and even on what machine to run on (the main platforms we used presently are the Origin 2000, T3E and Linux clusters, problems like ``memory leak'' could happen on one system but not the others). Fortunately many problems in making inappropriate choices can be picked up by the convergence tests discussed above.

IX. Stability with respect to physical parameters.

For physical results that should not sensitively depend on some physical parameters, we have to carry out tests to make sure that this is the case. This is important especially in cases that we want to model a realistic physical system: In many cases, we may not be able to set initial data that represents precisely the physical system. In the standard York [30] approach in setting initial data, part of initial metric and extrinsic curvature is determined by the constraint equations and cannot be specified a priori. One needs to make sure that the physics results coming out from the simulation are not sensitive to the uncertainty in the initial setup. For example, in the simulation of the head-on collision of 1.4$M_{\odot}$ neutron stars freely falling in from infinity (to be discussed below), one needs to start with the two stars at finite separation with the correct infalling velocity. However, we do not have a way to specify this initial infall velocity field precisely. (Approximately, one can use using a post-newtonian treatment, or a full GR treatment with the point mass approximation.) In our investigation of whether the collision would lead to a prompt formation of a black hole, we need to make sure that the prompt collapse behavior we observed is not sensitive to this inability to specify the initial infall velocity. We need to carry out simulations covering the range of uncertainty in the velocity. (For discussion of this point, see [31].) Note that this stability/sensitivity requirement could strongly depend on what physics questions one is asking: If the aim is to get the gravitational waveform from the head-on collision process, the requirement would be much more severe than to determine whether a prompt collapse occurred.

Other physical parameters that one need to consider in neutron star processes include the mass of the objects, spin states, equation of state of the matter etc., but the basic idea of carrying out simulations to cover the range of unknown is the same.

A different class of stability tests is on the effects of using different slicings and spatial coordinates. The final physics results, in suitable measure, should be independent of the coordinate choices. We have built various coordinate choices in our code (as discussed above), and have constructed various coordinate invariant measures (including curvature invariants, geodesic integrators and a proper time integrator) for such tests (and other purposes). While we routinely carry out runs with different slicings and spatial coordinates to investigate what would give us most stable and accurate evolutions, we do not carry out systematic studies of the coordinate independent nature of the final results often enough. We plan to build tools to make such study more convenient in the future.

X. Special tests: Symmetry, Scaling, Principle of Equivalence and Others.

In many problems there are often special tests that one can carry out. While these tests may not be applicable in general, they are valuable nevertheless. Here we discuss various considerations we have in this direction.

One obvious test is symmetry. Since our 3D spacetime-GRHydro code does not assume any symmetry, for any symmetric system the preservation of the symmetry provides a useful test. In this regard, one needs to be careful that the numerical methods used respect the symmetry. For example one of the two implementations of the MacCormack finite differencing scheme in our code breaks parity.   We note that many finite differencing schemes break reflection and rotation symmetry just due to the ordering of terms. When the numerical schemes used respect the symmetry of the system, the symmetry should be exactly preserved even in low resolution simulation.  Symmetry considerations have enabled us to pick out various code problems in early stage of our code development effort.

Scaling is a particular kind of symmetry. For neutron star studies using a polytropic equation of state (EOS) $P = K \rho ^ \Gamma$ (which is what most of the neutron star numerical studies used), it turns out that the spacetime-GRHydro equations admit a scaling of $t -> t/c, x-> x/c$, $\rho -> c^2 \rho$ and $K -> c^{2(1-\Gamma)} K$. The scaling means that a run involving a 1.4 $M_{\odot}$ neutron star with a polytropic EOS of $K=80$ should be identical to a run using a 0.7 $M_{\odot}$ star with $K=20$ at the correspondingly scaled spacetime points. For the evolution equations themselves, the scaling test is just a special case of a consistency test. However the scaling test include more aspects of the code, and can potentially pick out problems in initial data and boundary conditions. We have carried out such scaling tests in our neutron star studies.

One basic problem in validating a spacetime-GRHydro code is our lack of knowledge of curved space phenomena. One remedy is that some aspects of the curved spacetime evolution are linked to aspects in flat space through the equivalence principle, namely, local processes in curved space should be the same as those in flat space in terms of properly measured quantities. For example, in the 1.4 $M_{\odot}$ neutron star head-on collision simulation, the validation of the shock propagation speed is important for the confirmation of the prompt collapse to a black hole [32]. As there are very few curved spacetime shock propagation testbeds one can carry out, testings of the curve space shock treatment are particularly desirable. We note that shock propagation is a local process that one can investigate through the principle of equivalence, and our shock treatment has been well established in flat space with many flat space shock tube tests. We hence extracted the locally measured fluid properties on the two sides of the shock front in the neutron star headon collision curve space simulations, re-created a flat space shock with the same fluid properties in the proper frame of the pre-shock fluid, and confirmed the shock propagation speed. Similar considerations can be useful in other cases.

Each of the above classes of tests (I to X) serve particular purposes. In practice, we do not carry out these tests in the order listed above, but rather in the following manner. In constructing a new piece of code, e.g., adding a new hydrodynamic numerical treatment, we will first confirm the consistency of the finite differencing equations with the differential equations (V. above) before carrying any numerical evolution. In the first tests of numerical evolutions, we confirm the symmetry properties (X. above), and carry out the short term convergent tests (IV. above), comparing to the analytic solutions and the lower dimensional/perturbation results (I. and II. above). Then we make sure that the new treatment is converging to the known results with the order expected, even in a long time evolution (VII. above). The new treatment will then be compared to results obtained by other numerical options already established in the code (III. above); hopefully the new treatment would yield better results for the class of problems it is designed for, e.g., for capturing some features of shock wave. When there are puzzling features observed, residual evaluators will be constructed to zero in on the problem (VI. Above). After establishing the code with the above, we will then apply it to new physical problems. Due to the non-linearity of the equations, there is no guarantee that we still have the same convergence properties for a new physical problem (convergence and stability is solution dependent). Especially for the physics results obtained after a certain time of evolution (e.g., the popping into existence of the apparent horizon in a gravitational collapse), we demand that the solution, and in particular the physics results, are convergent at this particular time (VII. above). The physics results will then be verified to be stable with respect to changes of the computational parameters, numerical methods, choices of boundary settings (VIII. above), and to change of physical parameters covering the range of uncertainty (IX. above). At this point we will also devise further checks for various testable features (X. above).

The validation of a physics result obtained with the spacetime-GRHydro numerical simulations is an extremely time consuming and computational expensive process. However, we see no short cut to this process especially in the present early stages of the development. We further note that no amount of tests can 100% guarantee the absolute correctness of a code, and the ``correctness'' could depend on the which physics system one is studying. Extra-ordinary results need an extra-ordinary amount of testings. Of course, at the end of all the tests, one still have to appeal to one's physical intuition and understanding as the final judge.

In the above we used often use as an example our neutron star head-on collision study, one of our first physics application of the GR-Hydro code .   More of this application appear in the next section.

Back to the top of page

Astrophysics Applications

While the main goal of the Grand Challenge Project are the construction and validation of the multi-purpose Spacetime-GRHydro code, we believe it is important to push for some astrophysical applications in the process of code development. One major application we are using as driver for the code development effort is the head-on collision of two neutron stars. In the following we first give an introduction to the background of the physics involved in this problem and present the results obtained using the code.

There is an intriguing conjecture recently put forth by Shapiro [33] stating that for neutron stars (NSs) colliding head-on in a free-fall from infinity, the merged object will be supported by the thermal pressure, and it will not collapse into a black hole before significant radiative cooling.; We called this a delayed collapse, in contrast to a prompt collapse that happens on the dynamical time scale. Delayed and prompt collapses are drastically different: radiative cooling time scale is of order seconds, much longer than the dynamical timescale which is of order $10^{-3} $ seconds. A delayed collapse will be spherical and generate much less gravitational wave. While for low mass neutron star collisions (in which the merged object is not much heavier than the critical TOV mass) one may not be surprised to see delayed collapse, the intriguing point is: The conjecture [33] states that it will always be a delayed collapse independent of the masses of the initial neutron stars with any polytropic equation of state $P = K(s) \rho ^ \Gamma$, where the polytropic coefficient $K$ is a function of entropy $s$, and polytropic index $\Gamma$ is a space-time constant.

The essence of the argument of [33] is that, with the polytropic EOS, there always exists a stable equilibrium TOV configuration with the same total mass and total energy as the two initially isolated TOV stars. As mass and energy are conserved to a good extent in the head-on collision process, it is concluded [33] that the merged object will be in this equilibrium state until energy is radiated away at a much longer timescale.

The argument based on conservation is appealing. However we feel that there is a potential danger which could invalidate the argument, namely, quasi-equilibrium consideration may not be applicable to this problem. If the evolution is not quasi-equilibrium, the existence of an equilibrium state with the same mass and energy does not imply that this state will be attained in the dynamical evolution process.

The comparison of the various time scales in the collision process is crucial. The equilibrium state can be attained if the time scale for thermo and dynamical equilibrium is short enough compare to the other time scales. The relevant time scales in the collision process are 1. The time scale of infall: $t_i = R/ V_i$, R=the radius of the NS, $V_i$= (infall velocity at the point of contact). 2. The time scale associated with the local sound velocity: $t_s = R/ V_s$ , $V_s$= (sound velocity). 3. The time scale associated with the velocity of the shock (velocity in the rest frame of fluid): $t_{sh} = R/V_{sh}$, $V_{sh}$=(shock velocity). 4. The time scale for the merged object to thermalize globally, in the sense of being describable by a single EOS with the same polytropic index $K$ everywhere: $t_e$. 5. the neutrino cooling time scale: $t_n$, 6. The gravitational radiation cooling time scale $t_g$. 7. The time scale for local thermal equilibrium $t_l$. 8. The gravitational collapse time scale $t_c$.

It is clear that $t_n$ and $t_g$ are orders of magnitudes longer than the other time scales in the head-on collision and are irrelevant to the question of whether equilibrium can be achieved. On the other hand, $t_l$ is shorter than all other time scales, and one can always assume a local equation of state, e.g., one in the form of $P = K(s) \rho ^ \Gamma$ as used in the conjecture. However, the time for GLOBAL equilibrium$t_e$ is not necessarily short, and $K(s)$ may not be a spatial constant throughout the merged object in a time scale short compare to the dynamical timescales. $t_e$ can at most be as short as $t_{sh}$, assuming that the shock is effective in converting the bulk kinetic energy into thermal energy. But the question is if it can be as short as the dynamical time scale.

In Ref. [32] we estimate these time scales using the case of two 1.4 $M_{\odot}$ NSs. We model the stars with a polytropic EOS $P = K \rho ^ \Gamma$ with a polytropic index of $\Gamma =2$. The initial polytropic coefficient $K$ is taken to be $1.16 \times 10^5 \; \frac {{cm}^5}{g \; s^2}$, a typical value used in NS collision studies (see e.g., [34] and references therein). As assumed in Shapiro's conjecture, the polytropic index $\Gamma$ is taken to be a constant in space and in time, whereas the polytropic coefficient $K$ is a function of entropy and hence would change in space and time. We note that the argument in [33] is applicable to all polytropic models.

It is straightforward to carry out an order of magnitude estimate of the various time scales for the head-on collision of the 1.4 $M_{\odot}$ NSs freely falling in from infinity as discussed in Ref. [32]. The bottom line is that there is no reason to expect that the shock wave can convert the bulk kinetic energy to thermal energy in a short enough time scale to justify the quasi-equilibrium assumption needed in the Shapiro conjecture. And apparently there is no other mechanism to achieve equilibrium in a shorter time scale. As quasi-equilibrium may not be maintained in the process, the argument in [33] based on the conservation laws is not applicable to the collision of the 1.4 $M_{\odot}$ NSs: the prompt formation of a black hole is also consistent with the conservation laws. The existence of an equilibrium TOV solution is just a necessary condition, but not a sufficient condition, for a delayed collapse.

As we believe that Shapiro's argument cannot give a conclusive indication of whether one has a delayed or prompt collapse in the case of 1.4 $M_{\odot}$ collision, what actually happens needs to be determined through a numerical simulation. The collapse time scale $t_c$ is governed by the non-linear process involving the spacetime-hydro coupling.; Its determination needs the full Einstein theory of gravity coupled to GR-hydrodynamics (but does not require radiation transport, which is absence in our treatment).;; This serves as a good physics problem to drive the code development effort.

Some of the results of our numerical study are reported in [32,31] for stars with polytropic EOS (as assumed in the Shapiro conjecture) . The main result is that for the 1.4 $M_{\odot}$ case we find an apparent horizon in the infalling timescale $t_i$. The numerical study yields the following picture.

When the two 1.4$M_{\odot}$ NSs touch after infalling from infinity, they approach one another at a fraction of the light speed. A strong shock is generated, converting about $90\%$ of the bulk kinetic energy to thermal energy in the post-shock material. However, the rest of the stars (pre-shock) cramp in so rapidly that an apparent horizon quickly forms, trapping everything, including the shock front, inside. This is in sharp contrast to the ``equilibrium picture'' of [33,35]. In the ``equilibrium picture'' one envisions that the shock wave bounces a couple of times across the whole coalesced object, heating it up approximately uniformly [so that the polytropic constant$K(s)$ becomes a spatial constant], and the coalesced object can be described by a TOV configuration in equilibrium.

The dynamical nature of the collision can be illustrated by the world line of the shock front and its relation to the AH. In Fig. 1 below, the solid line is the world line of the shock front in the $z$ direction (direction of infall), plotted in coordinate distance and coordinate time. The ``*'' represents the location of the apparent horizon found in the simulation. The dotted line represents one leg of the light cone extending backward from the apparent horizon. The backward light cone meets the shock front at $t=0.195 ms$, which is $0.048 ms$ after the two stars first touched. At this point the shock has just barely propagated outward and the material heated up by the shock wave (the post-shock material) made up of only a small fraction of the total amount of matter on this time slice. In terms of the number of baryons, it is less than 5%. (One could consider different time slices, but the basic picture is the same.) We do not expect that such an insignificant fraction of heated material could have much effect on the overall dynamics of the infalling matter, let alone providing the thermal pressure to prevent the collapse as envisioned in the Shapiro argument. At later times the shock heated material does make up a bigger fraction of the total mass. However, from the point of intersection of the dashed and solid line at $t=0.195 ms$ onward, the shock heated material is causally disconnected with the apparent horizon shown and could not affect its formation. This highlights how far away from thermal and dynamical equilibrium the collision process is: in the ``equilibrium picture'' the shock is supposed to heat up the whole star and brings it to a stable TOV configuration.


Figure: The causal relation between the apparent horizon and the shock front is displayed. The point marked by ``*'' represents the intersection of the apparent horizon with the $+z$ axis (the stars are falling towards one another along the $z$ axis with the center of collision at $z=0$). A null geodesic traveling along the $z$ axis is shown as a dotted line. The solid line represents the world line of the intersection of the shock front on the $+z$ axis. World lines are plotted in coordinate distance vs. coordinate time.
This ``dynamical picture'' of the 1.4 $M_{\odot}$ NS collision has been subjected to ALL tests mentioned in the previous section (test I. to X). For some of these tests, in particular, the dependence of this picture on the initial setup, see [31].

We have recently carried out various pilot studies of the head-on collisions of lower mass NSs, aiming at determining the transition point between prompt and delayed collapse.  We are also carrying out headon collision study with the final version of the Spacetime-GRHydro code GR3D_EOS, using the Lattimer-Swesty equation of state.  Results will be reported elsewhere.

Back to the top of page

Visualization

 In the following we show a few more test cases performed using GR3D. The runs are obtained with the ADM evolution scheme for the spacetime evolution, and the hyperbolic hydrodynamics scheme, MAHC, for the hydrodynamics evolution. Parameter files are available for each evolution, allowing the user to replicate these results (see code documentation for a full description "here" ). For better visualization, choose the "Quicktime" option.
[Movies/AHEpsRho.mpg]  Click here for the MOV version of the movie.

This is a movie of the head on collision of two neutron stars using the general relativistic code. The blue/green part is the logarithm of energy density in geometric unit with solar mass equals 1. The red/yellow part is the logarithm of the internal energy. The apparent horizon is the dotted surface that appears at around T=63.

[Movies/RealNewtonII.mpg]  Click here for the MOV version of the movie.

This is a simulation of two inspiralling neutron stars using the Newtonian high resolution shock capturing code.  The green/blue portion is the density and the red/orange portion is the volumetric isolevels of the internal energy.

[Movies/RogRho.fast.mpg]  Click here for the MOV version of the movie.

This is the logarithm of the density for the same simulation as the previous one.

[Movies/RhoAlpha.fast.mpg]  Click here for the MOV version of the movie.

This is a simulation of two inspiralling neutron stars using the general relativity code GR3D.  The two spheres at top are two isosurfaces of the density of the two netron stars, and the curved surface is the lapse of the slice that cuts through the center of the stars.

[Movies/RhoBeta.mpg]  Click here for the MOV version of the movie.

The brown surface is an isosurface of the density and the arros represent the shift vector for the same simulation as the previous one.

[Movies/ra_01.mpg]   Click here for the MOV version fo the movie.

The two spheres at top are two isosurfaces of the density of the two netron stars in a grazing collision, simulated using the general relativity code GR3D.  The curved surface is the lapse of the slice that cuts through the center of the stars.

[Movies/r_v_h.mpg]    Click here for the version of the movie.

This is the logarithm of the density of the grazing collision in geometric units with solar mass equals 1, for the same simulation as the previous one.

[Movies/Rho_01.mpg] Click here for the MOV version of the movie.

These are three isosurfaces of the density of the grazing collision with the axis in the background, for the same simulation as the previous one.

Back to the top of page

GR3D Code Performance (Last Performance Milestone Code)

Machine Specification: Date tested: Code tested:
Testbed  32 bit  64 bit 
Grid Size per Processor  84x84x84 66x66x66 
Processor topology  8 x8 x16  8 x8 x16 
Total Grid Size  644 x 644 x 1284  500 x 500 x 996
Single Proc MFlop/sec  144.35  118.33 
Aggregate GFlop/sec  142.2  115.8 
Scaling efficiency  96.2%  95.6% 

Back to the top of page

Credit for Code Development

The spacetime evolution code in GR3D_EOS is a version of the Cactus code specially adapted for this NASA neutron star grand challenge project. The Cactus code is jointly developed by past and present members of the NCSA/Potsdam/Wash U collaboration, including Miguel Alcubierre, Gabrielle Allen, Werner Benger, Bernd Brugmann, ChangXue Deng, Edwin Evans, Tom Goodale, Achameevda Gopakuma, Philip Gressman, Carsten Gundlach, Ken Hsieh, Sai Iyer, Gerd Lanfermann, Lap Ming Lin, Joan Masso, Mark Miller, Thomas Radke, John Shalf, Hisaaki Shinkai, Malcolm Tobias, and Paul Walker, together with many others in the numerical relativity community. In particular, for the various important subroutines included in this GR3D_EOS version, Maximal (for constructing maximal slicing) is written by Miguel Alcubierre, Joan Masso, Malcolm Tobias and Paul Walker, BAM_Elliptic (for solving elliptic equations with multi-grid methods) is by Bernd Brugmann, SOR and Elliptic (for solving elliptic equations with SOR and other methods) by Bernd Brugmann, Joan Masso and Paul Walker, BonaMasso_Evolution (for evolution of Einstein equations in the BM formulation) by Tom Goodale, Joan Masso, Mark Miller and Paul Walker, BSSN and Conf-ADM (for evolution of Einstein equations in the conformal-tracefree formulations) by Miguel Alcubierre and Mark Miller, LinearWaves (for constructing gravitational wave spacetimes) by Joan Masso and Malcolm Tobias, Util, StandardEinstein, BMUtil (various utilities needed for the code) by Joan Masso and Paul Walker, ADM_Evolution (for evolution of Einstein equations in the ADM formulation) by Joan Masso, Mark Miller, and Paul Walker, TSIVP (for setting up time symmetric initial data) by Mark Miller and Malcolm Tobias, PNinit (for constructing initial data based on the post-newtonian formulation) by Hisaaki Shinkai, Newton (for evolution in Newtonian gravity) by Philip Gressman, IEEEIO (for data input and output) by John Shalf, and parallelism software (PUGH) by Joan Masso, Mark Miller, Paul Walker, and others.

The hyperbolic hydrodynamic module in GR3D_EOS, MAHC, is written by Mark Miller, and developed by Jose A. Font, Philip Gressman, Mark Miller and Malcolm Tobias.

The coupling of the code to tabular equation of state, and the import of the Lattimer-Swesty EOS,  is made possible by
Teepanis Chachiyo, Edwin Evans, Tom Goodale, Achameevda Gopakuma, Philip Gressman and Mark Miller.  The Lattimer-Swesty equation of state is by Jim Lattimer and Doug Swesty.
 

We thank Tom Clune at Cray and Wolfgang Mertz at SGI for help in optimization and benchmarking.

Back to the top of page

Publication & Presented Papers

  1. ``Three Dimensional Numerical General Relativistic Hydrodynamics: Formulations, Methods, and Code Tests", J. A. font, M. Miller, W.-M. Suen and M. Tobias, Phys. Rev. D, Phys.Rev. D61 (2000) 041501.
  2. On the Numerical Stability of the Einstein Equations, Mark Miller, Phys. Rev. D, (submitted, gr-qc/0008017)
  3. ``Towards a Stable Numerical Evolution of Strongly Gravitating System: The Conformal Treatments", M. Alcubierre, G. Allen, B. Brugmann, T. Dramlitsch. J. A. Font, P. Papadopoulos, E. Seidel, N. Stergioulas, R. Takahashi, M. Tobias, Phys.Rev. D62 (2000) 044034
  4. ``Collision of 1.4 solar mass Neutron Stars: Dynamical or Quasi-Equilibrium?", Mark Miller, Wai-Mo Suen, Malcolm Tobias, gr-qc/9910022
  5. The Shapiro Conjecture: Prompt or Delayed Collapse in the Head-On Collision of Neutron Stars?, M. Miller, W.-M. Suen and M. Tobias, Rapid Comm. Phys. Rev. D (submitted), gr-qc/9904041)
  6. ``Numerical Relativity As A Tool For Computational Astrophysics", E. Seidel and W.-M. Suen, Journal of Computational and Applied Mathematics, invited review, (in press)
  7. ``Three Dimensional Numerical Relativity with a Hyperbolic Formulation", Carles Bona, Joan Masso, Edward Seidel, Paul Walker, Phys. Rev. D, (in press), gr-qc/9804052
  8. ``Dynamical Evolution of Boson Stars II: Excited States and Self-Interacting Fields", Jayashree Balakrishna, Edward Seidel, Wai-Mo Suen, Phys.Rev. D58 (1998) 104004, gr-qc/9712064
  9. ``Numerical Methods for Modeling Binary Neutron Star Systems", A. Calder, F. D. Swesty, and E. Wang, to appear in the Proceedings of the Second Oak Ridge Symposium on Atomic and Nuclear Astrophysics (gr-qc 9806019)
  10. ``Numerical Models of Newtonian and Post-Newtonian Binary Neutron Star Mergers", E. Wang, F. D. Swesty, and A. Calder , to appear in the Proceedings of the Second Oak Ridge Symposium on Atomic and Nuclear Astrophysics (gr-qc 9806022)
  11. ``The Development of a 3-D General Relativistic Self-Gravitating Hydrodynamics Code", F. D. Swesty and E. Wang, to appear in the proceedings of the NASA HPCCP/CAS Workshop"
  12. ``Numerical Relativity: Towards Simulations of 3D Black Hole Coalescence", Edward Seidel, plenary talk given at GR15, Poona, India, to appear in the proceedings.
  13. ``The Synergy between Numerical and Perturbative Approaches to Black Holes", Edward Seidel, to appear in "On the Black Hole Trail" eds. Bala Iyer and Biplab Bhawal (Kluwer)
  14. ``Perturbative Approach to Quasi-Normal Modes of Dirty Black Holes", P. T. Leung, Y. T. Liu, W.-M. Suen, C. Y. Tam and K. Young, Phys. Rev. D (in press)
  15. ``Numerical Evolution of Matter in Dynamical Axisymmetric Black Hole Spacetimes. I: Methods and Tests", S. Brandt, J.A. Font, J.M. Ibanez, J. Masso, and E. Seidel, Comput.Phys.Commun. 124 (2000) 169, gr-qc/9807017
  16. ``Numerical Relativity in a Distributed Environment", Werner Benger, Ian Foster, Edward Seidel, John Shalf, Warren Smith, Paul Walker, 1999 SIAM Conference on Parallel Processing (submitted).
  17. ``Test-Beds and Applications for Apparent Horizon Finders in Numerical Relativity", M. Alcubierre, S. Brandt, B. Bruegmann, C. Gundlach, J. Massó, E. Seidel, P. Walker, Potsdam Preprint, Phys. Rev D (submitted).
  18. ``Numerical Evolution of Matter in Dynamical Axisymmetric Black Hole Spacetimes. I. Methods and Tests", S. Brandt, J.A. Font, J.M. Ibanez, J. Massó, E. Seidel, Phys Rev D (submitted)
  19. ``Black Hole Spectroscopy: Determining Waveforms from 3D Excited Black Holes", Gabrielle Allen, Karen Camarda, Edward Seidel, Phys. Rev. Letters (submitted)
  20. ``Evolution of Distorted Black Holes: A Perturbative Approach", Gabrielle Allen, Karen Camarda, Edward Seidel, Phys. Rev. D (submitted)
  21. ``Three-Dimensional Simulations of Distorted Black Holes. I. Comparison with Axisymmetric Results", Karen Camarda, Edward Seidel, Phys.Rev. D59 (1999) 064019.
  22. ``Adaptive Computation of Gravitational Waves from Black Hole Interactions", Philippos Papadopoulos, Edward Seidel, Lee Wild, Phys.Rev. D58 (1998) 084002
  23. ``Towards an Understanding of the Stability Properties of the 3+1 Evolution Equations in General Relativity", M. Alcubierre, G. Allen, B. Brügmann, E. Seidel and W.-M. Suen, Phys. Rev. D, (in press) (gr-qc/9908079)
  24. ``Numerical Relativity as a Tool for Computational Astrophysics", E. Seidel and W.-M. Suen, J. of Computational and Applied Mathematics (1999) (in press). ``Computational General Relativistic Astrophysics: the Neutron Star Grand Challenge Project", W.-M. Suen, Yukawa International Conference on Gravitational Waves and Relativistic Astrophysics (in press)
  25. ``Gravitational Collapse of Gravitational Waves in 3D Numerical Relativity", Miguel Alcubierre, Gabrielle Allen, Bernd Bruegmann, Gerd Lanfermann, Edward Seidel, Wai-Mo Suen, Malcolm Tobias, Phys. Rev. D 61 041501 (2000)
  26. ``Numerical Relativity and Gamma Ray Burst", American Physical Society Meeting (invited talk), April 2000
  27. ``Neutron Star Coalescences and Gamma Ray Bursts", Wai-Mo Suen, (invited talk), Proceedings of the International Conference on Role of Physics in the New Millennium: Research, Education & Society, July, 2000, Hong Kong
  28. ``A Conformal Hyperbolic Formulation of the Einstein Equations", Miguel Alcubierre, Bernd Brugmann, Mark Miller and Wai-Mo Suen, Phys. Rev. D 60, 64017 (1999)
  29. ``Event Horizons in Numerical Relativity II: Analyzing the Horizon", Joan Massó, Edward Seidel, Wai-Mo Suen, Paul Walker, Phys. Rev. D 59 64022 (1999)
  30. ``Symmetric Hyperbolic System in The Ashtekar Formulation", Gen Yoneda and Hisa-aki Shinkai, Phys. Rev. Lett. 82, 263-266 (1999)
  31. ``Newtonian and Post-Newtonian Binary Neutron Star Mergers", Hisaaki Shinkai, Wai-Mo Suen, F. Douglas Swesty, Malcolm Tobias, Edward Y. M. Wang and Clifford M. Will, Proceedings of the 8th Marcel Grossmann Meeting, (World Scientific, 1999).
  32. ``Numerical Evolution of Dynamic 3D Black Holes: Extracting Waves", Karen Camarda and Edward Seidel, Phys. Rev. D 57, R3204 (1998)
  33. "Finding Apparent Horizons in Dynamic 3D Numerical Spacetimes", Peter Anninos, Karen Camarda, Joseph Libson, Joan Massó, Edward Seidel, and Wai-Mo Suen, Phys. Rev. D 58, 24003 (1998)
  34. "Fully Relativistic Approach to Neutron Star Collisions", E. Seidel, Bad Honnef Meeting on Relativistic Astrophysics, Ed. by Harald Riffert (Vieweg, Braunschweig/Wiesbaden, 1998), page 229-236.
  35. ``Waves in Open Systems Via Bi-Orthogonal Basis", P. T. Leung, W.-M. Suen, C. P. Sun and K. Young, Phys. Rev. E 57, 6101 (1998)
  36. ``Logarithmic Perturbation Theory for Quasinormal Modes", P. T. Leung, Y. T. Liu, W.-M. Suen, C. Y. Tam and K. Young, J. Phys. A31, 3271 (1998)
  37. ``Generation of Scalar-Tensor Gravity Effects in Equilibrium State Boson Stars", G. L. Comer and H. Shinkai, Class. Quantum Grav. 15, 669, (1998) [gr-qc/9708071]
  38. ``Boosted Three-Dimensional Black-Hole Evolutions with Singularity Excision With the Binary Black Hole Grand Challenge Alliance", Phys. Rev. Lett. 80, 2512 (1998)
  39. ``Dynamical Evolution of Boson Stars in Brans-Dicke Theory", Jayashree Balakrishna and Hisa-aki Shinkai, Phys. Rev. D 58, 044016 (1998)
  40. ``The Thermal Response of a Pulsar Glitch: The Non-spherical Symmetric Case", K. S. Cheng, Y. Li and W.-M. Suen, ApJ. Lett. 449, 45 (1998)
  41. ``Finding Apparent Horizons in Dynamic 3D Numerical Spacetimes", P. Anninos, Karen

    42. Camarda, Joseph Libson, J. Masso, E. Seidel, and W.-M. Suen, Phys. Rev. D 58, 24003 (1998)
  42. ``Dynamics of Boson Star II: Excited States and Self-Interacting Fields", J. Balakrishna, E. Seidel and W.-M. Suen, Phys. Rev. D 58, 104004 (1998)
  43. Quasinormal-mode Expansion for Waves in Open Systems, E. S. C. Ching, P. T. Leung, W.-M. Suen, S. S. Tong and K. Young, Rev. Mod. Phys. 70, 1545, (1998).
  44. ``NASA HPCC/ESS Grand Challenge Applications Project NCCS5-153 1998 Annual Report", with P. Saylor et al., (http://wugrav.wustl.edu/Relativ/annual98.rep.)
  45. ``Dynamical evolution of boson stars", Jayashree Balakrishna, G. L. Comer, Ed Seidel, Hisa-aki Shinkai and Wai-Mo Suen, Proceedings of Numerical Astrophysics 98 (March 98, Tokyo), (Kluwer Academic, 1998)
  46. ``Impact of a Multi-TeraFlop Machine to Gravitational Physics", Wai-Mo Suen, Report to the National Science Foundation, June, 1998 ``Scalable, Hydrodynamici and Radiation-Hydrodynamic studies of Neutron Star Mergers and Supernovae Explosion", P.Saylor, D. Smolarski, F. D. Swesty, and E. Wang, in the ``Supercomputing 97"
  47. ``Three Dimensional Distorted Black Holes: Initial Data and Evolution", S. R. Brandt, K. Camarda, E. Seidel, in the Proceedings of the Marcel Grossmann 8 (Jerusalem, 1997)
  48. ``Dynamics of Gravitational Waves in 3D: Formulations, Methods, and Tests", P. Anninos, J. Masso, E. Seidel, W.-M. Suen and M. Tobias, Phys. Rev. D 56, 842 (1997)
  49. ``Quasi-Normal Modes of Dirty Black Holes", P. T. Leung, Y. T. Liu, W.-M. Suen, C. Y. Tam and K. Young, Phys. Rev. Lett. 78, 2894 (1997)
  50. ``Trick for Passing Degenerate Points in the Ashtekar Formulation", G. Yoneda, H. Shinkai and A. Nakamichi, Phys. Rev. D56, 2086, (1997)
  51. ``Gravitational Waves in Brans-Dicke Theory: Analysis by Test Particles around Kerr Black Hole", M. Saijo, H. Shinkai and K. Maeda, Phys. Rev. D 56, 785 (1997)
  52. ``Can We Detect Brans-Dicke Scalar Gravitational Waves in Gravitational Collapse?", H. Shinkai, M. Saijo and K. Maeda, The 18th Texas Symposium on Relativistic Astrophysics, (1996 December, Chicago) Proceedings, (World Scientific, 1997)
  53. ``The Numerical Evolution of Gravitational Waves", Malcolm Tobias, Ph.D thesis, Washington University, 1997 (unpublished)
  54. ``Methods in Numerical Relativity", W.-M. Suen, Computational Physics Workshop, Hong Kong, 1997
  55. ``Lorentzian Dynamics in Ashtekar Gravity", H. Shinkai and G. Yoneda, The 8th Marcel Grossman Meeting (1997 June, Jerusalem) Proceedings, (World Scientific, 1997)
  56. ``Waves in Open Systems", E. S. C. Ching, P. T. Leung, A. Maassen Van den Brink, S. S Tong, W.-M. Suen, and K. Young, 2nd Joint Meeting of Chinese Physical Societies, Taipei, 1997
  57. ``Numerical Relativity and General Relativistic Astrophysics", Wai-Mo Suen, Paper Presented at the Second Joint Meeting of Chinese Physical Societies, Taipei, 1997
  58. ``Open Wave Systems", E. S. C. Ching, P. T. Leung, A. Maassen Van den Brink, S. S. Tong, W.-M. Suen, and K. Young, Invited Paper, 7th Asia Pacific Physics Conference, Beijing, 1997.
  59. ``Computational General Relativistic Astrophysics", E. Seidel and W.-M. Suen, Abstract, The 15th International Conference on General Relativity and Gravitation, India, 1997
  60. ``Event Horizons in Numerical Relativity: Understanding the Horizon", J. Massó, E. Seidel, W.-M. Suen and P. Walker, Abstract, The 15th International Conference on General Relativity and Gravitation, India, 1997
  61. ``NASA HPCC/ESS Grand Challenge Applications Project NCCS5-153 1997 Annual Report", with P. Saylor et al., (http://wugrav.wustl.edu/Relativ/annual.rep.)
  62. ``Computational General Relativistic Astrophysics and Coalescences of Neutron Star Binaries", Wai-Mo Suen, Paper Presented at the Second Joint Meeting of Chinese Physical Socities, Taipei, 1997
  63. ``When Black Holes Collide", P. Anninos, D. Hobill, E. Seidel and W.-M. Suen, Proceedings of the Sixth Canadian Conference on General Relativity and Relativistic Astrophysics, Fields Institute Communications, 54, 151, (1997)

    64.  

     
     
     


Back to the top of page

Acknowledgement

We thank our collaborators and especially those mentioned by name in the code credit section for their contribution to the project. We thank Tom Clune at Cray and Wolfgang Mertz at SGI for help in optimization and benchmarking. We thank Ed Evans, Achameevda Gopakuma, Ken Hsieh and Pranoat Priesmeyer at Washington University for their helps in preparing this document.

Back to the top of page

Bibliography

1
For a summary of the present and planned high energy missions, see e.g., http://heasarc.gsfc.nasa.gov/docs/heasarc/missions.html.
2
K. Thorne, in Proceedings of the Eight Nishinomiya-Yukawa Symposium on Relativistic Cosmology, edited by M. Sasaki (Universal Academy Press, Japan, 1994).
3
J. A. Font, M. Miller, W. M. Suen, and M. Tobias, (1998), Phys. Rev. D, in press.
4
E. Seidel and W.-M. Suen, J. Comp. Appl. Math. (1999), in press.
5
R. Arnowitt, S. Deser, and C. W. Misner, in Gravitation: An Introduction to Current Research, edited by L. Witten (John Wiley, New York, 1962), pp. 227-265.
6
T. Nakamura and K. Oohara, in Frontiers in Numerical Relativity, edited by C. Evans, L. Finn, and D. Hobill (Cambridge University Press, Cambridge, England, 1989), pp. 254-280.
7
M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995).
8
T. W. Baumgarte and S. L. Shapiro, Physical Review D 59, 024007 (1999).
9
M. Alcubierre et al., (1999), gr-qc/9908079.
10
C. Bona, J. Massó, E. Seidel, and J. Stela, Phys. Rev. D 56, 3405 (1997).
11
L. Smarr and J. York, Phys. Rev. D 17, 2529 (1978).
12
J. Ibanez, J. Font, J. Marti, and J. Miralles, in Proceedings of The 18th Texas Symposium on Relativistic Astrophysics, edited by J. F. A. Olinto and D. Schramm (World Scientific, Singapore, 1998).
13
F. Banyuls et al., ApJ 476, 221 (1997).
14
J. A. Font, et al., (1999), in preparation.
15
P. Anninos et al., Phys. Rev. D 52, 2044 (1995).
16
P. Anninos et al., in Computational Astrophysics: Gas Dynamics and Particle Methods, edited by W. Benz, J. Barnes, E. Muller, and M. Norman (Springer-Verlag, New York, 1997), in press.
17
S. Brandt et al., Phys. Rev. D (1998), submitted, gr-qc/9807017.
18
M. Alcubierre et al., (1999), gr-qc/9908012.
19
P. e. a. Gressman, (1999), in preparation.
20
E. Wang and F. D. Swesty, in Proceedings of The 18th Texas Symposium on Relativistic Astrophysics, edited by J. F. A. Olinto and D. Schramm (World Scientific, Singapore, 1998).
21
P. Anninos et al., Phys. Rev. D 52, 4462 (1995).
22
J. Baker et al., Phys. Rev. D 55, 829 (1997).
23
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (Cambridge University Press, Cambridge, England, 1986).
24
M. Alcubierre et al., in preparation (1999).
25
K. C. B. New, K. Watt, C. W. Misner, and J. M. Centrella, Phys. Rev. D 58, 064022 (1998), gr-qc/9801110.
26
M. Miller, Ph.D. thesis, University of Syracuse, 1996.
27
P. Lax, Comm. Purr Appl. Math 9, 135 (1956).
28
R. D. Richtmyer and K. Morton, Difference Methods for Initial Value Problems (Interscience Publishers, ADDRESS, 1967).
29
S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys.Rev.Lett. 82, 892 (1999), gr-qc/9810072.
30
J. York, in Sources of Gravitational Radiation, edited by L. Smarr (Cambridge University Press, Cambridge, England, 1979).
31
M. Miller, W.-M. Suen, and M. Tobias, , gr-qc/9904041.
32
M. Miller, W.-M. Suen, and M. Tobias, , submitted to Rapid Commun. Phys. Rev. D; gr-qc/9904041.
33
S. Shapiro, Phys. Rev. D 58, (1998).
34
M. Ruffert and H.-T. Janka, Astron. Astrophys. 338, 535 (1998).
35
S. Shapiro, (1999), gr-qc/9909059.


Back to the top of page