F. Douglas Swesty
dswesty@pulsar.astro.uiuc.edu
http://www.ncsa.uiuc.edu/~dswesty/
Paul Saylor
saylor@cs.uiuc.edu
http://www.uiuc.edu/ph/www/p-saylor
Dennis C. Smolarski, S.J.
DSMOLARSKI@SCUACC.SCU.EDU
http://www-acc.scu.edu/~dsmolarski/homepage.html
E. Y. M. Wang
eymw@deia.ncsa.uiuc.edu
http://nsnsgc.ncsa.uiuc.edu/nsnsgc/people.html
Scalable high performance computing has brought about a revolution in astrophysical modeling. One aspect of this revolution is the ability we now possess to conduct large scale 3-D fluid dynamics simulations for use in astrophysical modeling. Another aspect, that is only now emerging, is our ability to model the flow of radiation in multi-dimensions as well. In this paper we describe two ongoing astrophysical modeling projects that stress both of these aspects. One project involves the numerical modeling of the merger of two neutron stars in a binary system. This project has been selected as a NASA Earth and Space Sciences (ESS) High Performance Computing and Communication (HPCC) Grand Challenge. The other project involves the numerical modeling of convectively driven supernova explosions via multi-dimensional radiation hydrodynamics explosions. Both of these projects contain unique high performance computing challenges. In the subsequent sections we describe the scientific features of the problem. Then in succeeding sections we discuss high performance computing problems, results, and plans.
years the stars
will spiral rapidly together over a period of about 15 minutes. During
this brief period the gravitational waves will sweep up in frequency from
about 10Hz to about 1500Hz, when the final coalescence occurs. This
contains the range of the bandwidths of the LIGO and VIRGO observatories,
which should be operational around the year 2000. In the final coalescence
stage general relativistic effects, hydrodynamic instabilities, the
equation of state, and neutrino transport become extremely important
governors of the gravitational wave signal that is generated. It is during
this final phase that we intend to concentrate our efforts, where we can
bring a careful treatment of all these effects that has not been possible
previously. A number of independent calculations [2] show that
these final coalescence events will occur quite frequently in the universe,
the most conservative of which place the observable event rate at between
several and several dozen per year.
We briefly recount here some of the astrophysical context in which this problem is set and explain how our work fits in. Our calculations will provide important information about observational and theoretical issues in several major areas of astrophysics and relativity. Each of the items is discussed briefly in the following paragraphs.
First, gravitational waves from the merger of a neutron star with a neutron star (NSNS) and the merger of a neutron star with a black hole (NSBH) are likely to be detected by LIGO and VIRGO within 5--10 years, at a rate of dozens of events per year [2]. There is much interesting information contained in the gravitational wave signals. With the parameters of the coalescing objects, e.g., masses, radii, spin, etc. extracted from the signal through template matching, one can obtain information on the mass to radius relation, the distribution and range of neutron star masses etc. By tracking the system through the merger phase, the waveforms will tell us a great deal about the overall mass motions of the merging neutron star matter. An exciting possibility exists that the information contained in the waveform, in conjunction with detailed models for the dynamics of the merger, will allow us to constrain the equation of state (EOS) of dense matter. Such constraints would nicely complement neutron star thermal evolution constraints that can be placed on the EOS by observing data on galactic neutron stars from orbiting high energy astronomy satellites. Taken together, these combined constraints give important information about both neutron star structure and the physics of dense matter. All of these ideas call for detailed numerical models of the coalescence. Our intent is to carry out numerical simulations of the mergers to show in detail how variations in the EOS might produce discernible effects in the gravitational wave signal.
A second reason for studying NSNS mergers is that the coalescence of NSNS binary systems has been discussed as a potential site for the gamma--ray bursts [3] detected by the Burst and Transient Source Experiment (BATSE) on board the Compton GRO satellite. Despite the widespread discussion of NSNS mergers as a source of these bursts there have been few numerical calculations to assess the viability of this mechanism. These few calculations have left many unresolved questions about this model for bursts. Our proposed studies will direct attention to determining whether or not NS-NS mergers are compatible with the time scales and energetics of these bursts.
Still a third reason for studying NSNS mergers is that they have been cited as one of two possible sites of the astrophysical r--process, which is responsible for the production of many of the heavy neutron rich elements: The merger process of compact binaries may eject extremely neutron rich matter, which then undergoes the classical r--process nucleosynthesis. This could be a major source of r-process nuclei in the galaxy [4]. However, a proper treatment of this problem should include neutrinos and a realistic equation of state, as our calculations will.
The numerical modeling of neutron stars mergers presents some unique challenges. The problems requires not only the simulation of the motion of the fluid matter which constitutes the neutron stars, but also the gravitational fields which derives from the matter itself. Since the fluid bodies are self-gravitating, both the fluid and gravitational phenomena are intertwined in a complex synergy. Since the gravitational field of a neutron stars is sufficiently strong, it is more accurately described by Einstein's relativistic theory of gravity (general relativity, hereafter GR) as opposed to classical Newtonian gravitational theory. In GR, one must solve an extremely complex set of tensor field equations to describe the gravitational field instead of the relatively simple Poisson equation that arises from Newtonian gravity. Fortunately, the relativistic analog of the Euler equations of classical fluid dynamics mimic the structure of their Newtonian counterparts. The differences between the relativistic and non-relativistic fluid dynamics equations mainly arise from Lorentz contraction and space-time curvature effects. Because of the close similarity of the GR and Newtonian hydrodynamics equations we are able to adapt a variety of standard computational fluid dynamics methods to their solution. The work we describe in the remainder of this paper is based on finite-difference techniques for solving the fluid dynamics and gravitational field equations.
solar masses, and has central densities
in excess of 
g/cm
. At these densities the
neutrinos do not move freely through the matter, but instead diffuse
slowly in a random walk process. In contrast, at lower densities, i.e.
below 
g/cm
, the neutrinos move freely without
interacting with the matter. This wide range of behavior of the
neutrinos and the steep density gradients in the newly born neutron
star make the phenomenon hard to simulate. In addition to modeling
the flow of the matter, which requires solving the Euler equations of
hydrodynamics, one must also model the flow of radiation, in the form
of neutrinos, as well. The latter task is far more daunting then the
first.
Once the collapse of the iron-nickel core of the pre-supernova star has halted, the core ``bounces,'' i.e. it expands outward having overshot the equilibrium configuration during the collapse. The outer part of the core is still infalling supersonicly and is unaware that the inner portion of the core has reversed its motion. A shock wave forms near the sonic point in the flow and begins to propagate outward through the outer core. The shock wave drives a complex reactive flow involving the nuclear dissociation of heavy nuclei into free neutrons and protons. This process, which soaks up energy from the shock, eventually causes a shock to stall before reaching the edge of the iron-nickel core. The major challenge to supernova modelers is to explain the explosion of the star despite the fact that the shock wave stalls before escaping the core.
The weakening shock wave leaves behind another legacy: as the shock traverses the core it imprints a negative entropy gradient on the post-shocked core which is convectively unstable. Thus in the post bounce epoch of the supernova we also see convection occurring. The current thinking is that this convection revitalizes the shock and causes it to continue outward to eventually result in the explosion of the star. Thus the modeling of this phenomena requires that one simulate a radiating, reactive, convective flow in at least two (and preferably three) dimensions. Finally, the energy spectrum of the neutrinos must be discretized and this adds another energy dimension to the problem which becomes either three or four dimensional depending on how many spatial dimensions are modeled.
order
post-Newtonian corrections as described by Blanchet, Damour, &
Schäfer [10]. And the GR scheme is similar to that of
Hawley, Smarr, & Wilson [6].
Employing these methods we have been able to achieve high flop rates and reduced zone update times. The Newtonian code, V3D, currently achieves nearly 500 Mflops (497 Mflops to be precise) on a single Cray C90 processor. Using the F77 implementation of V3D we have been conducting preliminary simulations of neutron star mergers on the NCSA SGI Power Challenge. We have found that the F77 version of V3D can scale with 70% parallel efficiency across the entire 16-processor machine when we are employing polytropic equations of state. When we employ a more realistic equation the scaling improves beyond this. We are currently re-implementing this algorithm as a message passing code using the MPI message passing interface. We expect the scaling to improve still further when this is complete.
The bulk of the computation in V2D is spent in solving a linear system of equations arising from the implicit finite differencing of the radiation transport equations. The system of linear equations is block tridiagonal or pentadiagonal in nature and lends itself well to the use of Krylov subspace methods as we discuss later. Besides the cost of the linear systems solving the additional computational effort in this multigroup method additional computational cost of the multigroup approach is in the embarrassingly parallel operations of evaluating the neutrino-matter interaction terms. The V2D code currently exists as an f90 code and is being ported into an MPI implementation to take advantage of the finer grained parallelism on the Cray T3E. We have already implemented a linear system solver using f90 and MPI and we discuss the performance of this algorithm in a later section.
Recent algorithms proposed for approximate inverses, for example, [1], incorporate ingenious ideas that yield effective preconditioners, even compared to incomplete LU preconditioners on serial machines. We have taken a different, simplified approach however. We have selected block submatrices along the diagonal of our given matrix to come up with a nonsingular block diagonal matrix. We then compute the inverse of this and use it as our block inverse preconditioner. We have found this to be an effective method. (A separate report is under preparation.)
Using the aforementioned block inverse preconditioning method we have
implemented BiCG in MPI using a 1-D Cartesian process topology.
This particular topology is optimal for the solution of the transport
equation along radial rays through the supernova core. A typical
problem size would be approximately 
radial zones and
energy groups. This problem size will result in a system of
linear equations that must be solved. The matrix will
be block tridiagonal with the block size determined by the number of
energy groups, 
in this case. The domain
decomposition strategy results in each processor ``owning'' a number
of rows of blocks of the matrix. The vectors are also distributed
accordingly.
The bulk of the computational effort in this algorithm is due to
matrix-vector multiplies and the computation of the block inverses for
both the matrix and its transpose. The block inverse operation is
embarrassingly parallel while the matrix and transpose matrix
multiplies require communication of ghost rows from the neighboring
processes in the topology. The BiCG algorithm requires one matrix
multiply and one transpose multiply per iteration. Using asynchronous
sends and receives much of this communication can be hidden while work
is done on the interior rows of the process. The BiCG
algorithm also requires at least two dot products per iteration. A third
dot product is required if one wishes to employ the 
norm of the
residual as a stopping criterion. Each of these dot products requires
a costly global reduction. In fact, a performance
analysis reveals that this particular operation dominates the
communication costs of this algorithm. Unfortunately, MPI does not
provide an asynchronous global reduction operation so there exists
no means of hiding this particular portion of the message passing.
An example of the scaling we have achieved with this algorithm is
displayed in figure 1. We compare the scaling for a 
by
problem on the Cray T3E and the Silicon Graphics Origin 2000.
This particular problem took 13 iterations to converge. For problems
requiring fewer iterations the embarrassingly parallel computation of
the block inverse dominates the work and the scaling will increase.
Figure 1: Preliminary scaling and parallel efficiency
results for the MPI implementation of BiCG on the Cray T3E and
the Silicon Graphics Origin 2000. The solid lines indicate
the T3E results while the dashed lines indicate the Origin 2000
results. These results may not reflect the true capabilities of the
machines.
In order to address these needs we have developed a multi-grid Poisson solver that implements a full multi-grid W--cycle (FMW) algorithm. This code was implemented as a f90 code and is currently being used to carry out both Newtonian and post-Newtonian models of neutron star mergers on the NCSA SGI Power Challenge and Origin 2000. The current version of this code which relies on automatically parallelizing f90 compilers do not scale well. However, we are currently working on a port of this code into MPI so that we can implement it in a finer grained parallel environment such as the Cray T3E or the SGI Origin 2000.
The port of a multigrid algorithm to MPI is non-trivial. At grid levels where the grid size becomes comparable to the size of the processes topology some other algorithm must be adopted for distributing the mesh among the processors. We are currently investigating several alternate possibilities for dealing with this problem.
The motivation for developing a tabular EOS code was threefold. First, we wanted to provide a means of incorporating equations of state other than the LS EOS into our hydrodynamic codes without having to modify the interface. Also, this tabular EOS allows us to include extant EOS tables, as well as the results of numerically burdensome EOS calculations, directly into our code. The second motivation for the development of such a code was speed. A directly indexed table is usually much faster than most complex EOS codes. The third motivation for developing this code is that it allows us to easily match equations of state, which describe different ranges of density and temperature, in a smooth manner. Using this TCTEOS code we can easily incorporate almost any current dense matter EOS into our numerical simulations in a straightforward fashion.
Since the EOS evaluations usually dominate the flop count for our hydrodynamic simulations we have gone to great lengths to minimize the number of floating point operations needed to implement this tabular EOS scheme. We have also been able to construct the implementation in a fashion such that we are able to achieve nearly 50% of the processors theoretical peak speed on the MIPS R10000 processor which is used in the SGI Power Challenge and the SGI Origin 2000.
We have recently begun to carry out some preliminary studies of Newtonian neutron star mergers using the V3D code we described in the previous section. The reason for these simulations are fourfold: They give us estimates of the rate of gravitational energy radiated and the estimates of the waveforms (see figures 3 and 4), they give us a baseline to compare subsequent Post-Newtonian and relativistic simulations to, they allow us to begin to understand the physical environment that is present, and finally they allow us to conduct resolution studies. Using Newtonian studies we can also begin to evaluate issues like what the r-process yields from the merger are and how the EOS influences the dynamics of the coalescence.
We have found a number of intriguing results in the simulations that
we have conducted so far. The simulations have studied the Newtonian
coalescence, driven by tidal instabilities, of two 
neutron stars (modeled as n=1 polytropes) on a 
grid covering
a computational domain with a length of 
km on a side. Our
simulations yield gravitational wave signals that agree with the
results of other calculations [11,12] The radiated energy,
as predicted by the time variation of the quadrapole moment of the
mass, also yields radiated energies that are in good agreement with
these calculations.
Finally, it seems as though several tenths of solar mass of material is ejected from the outer layers of the stars into wide orbits. This is quite an exciting result from the standpoint of the r-process since such low density neutron rich material is exactly what is needed to make the r-process occur.
However, we have found that the onset of the tidal instability first
predicted by [7] seems to occur outside the 
separation. We now understand this result to be a consequence of the need
for better resolution in these calculations. By studying the resolution
needed to accurately represent static single neutron stars on a grid we
have discovered that the reason for our differences from the result of Lai
et al. is that higher resolution is required to adequately resolve the
stars advecting on a Cartesian mesh. The calculations of Lai et al. have
employed an SPH code, an inherently Lagrangian method, in a co-rotating
frame. However, since we ultimately would like to achieve fully
relativistic calculations on Cartesian grids, we needed to establish that
we could indeed achieve sufficient resolution needed to accurately
represent the stars on the grid. Also, since we are very much interested
in the r-process nucleosynthesis in the low density material that seems to
be ejected from the stars as they coalesce we want to use an Eulerian
method that can best track the flow of this low density material
accurately.
We have verified the resolution needed to represent stars in the
orbiting configuration by advecting single stars across the mesh to
make sure that they remain stable as they should. For a single star
we were able to employ a much smaller domain that yielded higher
resolution. Unfortunately, for the merger we cannot restrict the size
of our computational domain if we are to study how much mass is shed
during the merger, which is in turn important for calculating the
nucleosynthetic yields of r-process material. We will need to reach a
resolution of at least 
in order to clearly resolve the stars
and eliminate this problem. The 
resolution of the 
km
domain is less than adequate for resolving the neutron stars
on the 
grid we employed.
This agrees with the conclusions of
Ruffert et al.[11].
We have also begun testing of our V2D code with low angular resolution supernova models. These low resolution runs are effectively 1+1-D models since convection never develops under such conditions. However, they are important scientifically as they allow us to understand what the contributions of convection is to the whole phenomena. By comparing these low angular resolution model with high angular resolution models that are forthcoming we will be able to understand which aspects of the dynamics are due to convection and which are unrelated. The completion of our port of the V2D code to MPI will also enhance our ability to conduct these simulations at high resolution.
F. Douglas Swesty
has been working in the fields of Computational Astrophysics
and Nuclear Astrophysics for 10 years. He received his Ph.D. in Physics
from SUNY Stony Brook in 1993. He currently holds positions as a Visiting
Research Assistant Professor in the University of Illinois Department of
Astronomy and as a Research Scientist at the National Center for
Supercomputing Applications. His research interests are primarily focused
on the physics of compact objects (black holes, neutron stars, supernovae,
white dwarfs), radiation hydrodynamics, numerical relativity, and
computational astrophysics. Some of his currently active projects are:
research on the equation of state of hot, dense matter; 2-D and 3-D
modeling of core collapse supernovae; modeling the merger of binary neutron
star systems in 3-D; the development of adaptive mesh radiation
hydrodynamic algorithms; and the development of parallel methods for
radiation transport.
Paul Saylor has been a member of the Computer Science Department at
Illinois since 1967. Currently affiliated with Lawrence Livermore
Laboratory, and in the past with Lawrence Berkeley Lab, Los Alamos Lab
and Oak Ridge Lab. Principal Investigator for the NASA Earth and Space
Grand Challenge on Simulating the Merger of Binary Neutron Stars. Interests
are in high performance computing and the solution of large linear and
nonlinear systems with applications to astrophysics, computational
electromagnetics and groundwater flow.
Dennis C. Smolarski, S.J. is an associate professor in the Department of
Mathematics at Santa Clara University, specializing in Computer Science.
His research areas focus on iterative solutions of linear systems.
E. Y. M. Wang is a third year graduate student at Washington University in St. Louis in the department of physics. He is currently working on his Ph.D. thesis at the National Center for Supercomputing Applications as a research assistant. The thesis endeavours to examine compact object astrophysics, in particular the numerical modelling of neutron stars and their mergers in binary systems. This study will involve a knowledge and understanding of the structure of compact objects, the equation of state of dense matter and the dynamical behaviour Newtonian and general relativistic spacetimes. An understanding of numerical relativity and computational astrophysics will be necessary for this work. A tentative date for the completion of his Ph.D. work is July 1999.
