4 Hydrodynamical Simulations in Relativistic 3 Numerical Schemes3.2 Other techniques

3.3 State-of-the-art three-dimensional codes 

The most advanced time-dependent, finite difference, three-dimensional Cartesian codes to solve the system of Einstein and hydrodynamics equations are those developed by Shibata [258Jump To The Next Citation Point In The Article] and by the Washington University/NCSA/AEI-Golm Numerical Relativity collaboration [96Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article]. The main difference between both codes lies in the numerical methods used to solve the relativistic hydrodynamic equations: artificial viscosity in Shibata's code [258Jump To The Next Citation Point In The Article], and upwind HRSC schemes in the code of [96Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article]. We note, however, that very recently Shibata has upgraded his code to incorporate HRSC schemes (in particular, a Roe-type approximate Riemann solver and piecewise parabolic cell-reconstruction procedures) [260Jump To The Next Citation Point In The Article]. Further 3D codes are currently being developed by a group in the U.S. (Duez et al. [73Jump To The Next Citation Point In The Article]) and by a E.U. Research Training Network collaboration [119].

3.3.1 Shibata 

In Shibata's code [258Jump To The Next Citation Point In The Article], the hydrodynamics equations are formulated following Wilson's approach (Section  2.1.2) for a conformal-traceless reformulation of the spacetime variables (see below). An important difference with respect to the original system, Equations (20Popup Equation, 21Popup Equation, 22Popup Equation), is that an equation for the entropy is solved instead of the energy equation. The hydrodynamic equations are integrated using van Leer's [290Jump To The Next Citation Point In The Article] second order finite difference scheme with artificial viscosity, following the approach of a precursor code developed by [214Jump To The Next Citation Point In The Article].

The ADM Einstein equations are reformulated into a conformal traceless system, an idea originally introduced by Shibata and Nakamura [264Jump To The Next Citation Point In The Article] (see also [197Jump To The Next Citation Point In The Article]), and further developed by Baumgarte and Shapiro [25Jump To The Next Citation Point In The Article]. This ``BSSN'' formulation, which shows enhanced long-term stability compared to the original ADM system, makes use of a conformal decomposition of the 3-metric, tex2html_wrap_inline5148 and the trace-free part of the extrinsic curvature, tex2html_wrap_inline5150, with the conformal factor tex2html_wrap_inline5152 chosen to satisfy tex2html_wrap_inline5154 . In this formulation, as shown by [25Jump To The Next Citation Point In The Article], in addition to the evolution equations for the conformal 3-metric tex2html_wrap_inline5156 and the conformal-traceless extrinsic curvature variables tex2html_wrap_inline5158, there are evolution equations for the conformal factor tex2html_wrap_inline5152, the trace of the extrinsic curvature K and the ``conformal connection functions'' tex2html_wrap_inline5164 . Further details can be found in [25Jump To The Next Citation Point In The Article, 258Jump To The Next Citation Point In The Article].

The code uses an approximate maximal slicing condition, which amounts to solving a parabolic equation for tex2html_wrap_inline5166, and a minimal distortion gauge condition for the shift vector. It also admits tex2html_wrap_inline5168 -rotation symmetry around the z -axis, as well as plane symmetry with respect to the z =0 plane, allowing computations in a quadrant region. In addition, Shibata has also implemented in the code the ``cartoon'' method proposed by the AEI Numerical Relativity group [6Jump To The Next Citation Point In The Article], which makes possible axisymmetric computations with a Cartesian grid. ``Approximate'' outgoing boundary conditions are used at the outer boundaries; these do not completely eliminate numerical errors due to spurious back reflection of gravitational waves Popup Footnote . A staggered leapfrog method is used for the time evolution of the metric quantities. Correspondingly, the hydrodynamic equations are updated using a second-order two-step Runge-Kutta scheme. In each time step, the staggered metric quantities needed for the hydrodynamics update are properly extrapolated to intermediate time levels to reach the desired order of accuracy.

The code developed by Shibata [258Jump To The Next Citation Point In The Article, 260Jump To The Next Citation Point In The Article] has been tested in a variety of problems, including spherical collapse of dust to a black hole, signalled by the appearance of the apparent horizon (comparing 1D and 3D simulations), stability of spherical stars and computation of the radial oscillation period, quadrupole oscillations of perturbed spherical stars and computation of the associated gravitational radiation, preservation of the rotational profile of (approximate) rapidly rotating stars, and the preservation of a co-rotating binary neutron star in a quasi-equilibrium state (assuming a conformally flat 3-metric) for more than one orbital period. Improvements of the hydrodynamical schemes have been considered very recently [260Jump To The Next Citation Point In The Article], in order to tackle problems in which shocks play an important role, e.g., stellar core collapse to a neutron star. Shibata's code has allowed important breakthroughs in the study of the morphology and dynamics of various general relativistic astrophysical problems, such as black hole formation resulting from both the gravitational collapse of rotating neutron stars and rotating supermassive stars, and, perhaps most importantly, the coalescence of binary neutron stars, a long-standing problem in numerical relativistic hydrodynamics. These applications are discussed in Section  4 . The most recent simulations of binary neutron star coalescence [267Jump To The Next Citation Point In The Article] have been performed on a FACOM VPP5000 computer with typical grid sizes of (505, 505, 253) in (x, y, z).

3.3.2 CACTUS/GR_ASTRO 

A three-dimensional general relativistic hydrodynamics code developed by the Washington University/NCSA/AEI-Golm collaboration for the NASA Neutron Star Grand Challenge Project [296] is discussed in Refs. [96Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article]. The code is built upon the Cactus Computational Toolkit [171]. A version of this code that passed the milestone requirement of the NASA Grand Challenge project was released to the community. This code has been benchmarked at over 140 GFlop/sec on a 1024 node Cray T3E, with a scaling efficiency of over 95%, showing the potential for large scale 3D simulations of realistic astrophysical systems.

The hydrodynamics part of the code uses the conservative formulation discussed in Section  2.1.3 . A variety of state-of-the-art Riemann solvers are implemented, including a Roe-like solver [242] and Marquina's flux formula [70Jump To The Next Citation Point In The Article]. The code uses slope-limiter methods to construct second-order TVD schemes by means of monotonic piecewise linear reconstructions of the cell-centered quantities for the computation of the numerical fluxes. The standard ``minmod'' limiter and the ``monotonized central-difference'' limiter [289] are implemented. Both schemes provide the desired second-order accuracy for smooth solutions, while still satisfying the TVD property. In addition, third-order piecewise parabolic (PPM) reconstruction has been recently implemented and used in [279Jump To The Next Citation Point In The Article].

The Einstein equations are solved using the following different approaches: (i) the standard ADM formalism, (ii) a hyperbolic formulation developed by [40], and (iii) the BSSN formulation of [197Jump To The Next Citation Point In The Article, 264, 25]. The time evolution of both the ADM and the BSSN systems can be performed using several numerical schemes [96Jump To The Next Citation Point In The Article, 4Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article]. Currently, a leapfrog (non-staggered in time), and an iterative Crank-Nicholson scheme have been coupled to the hydrodynamics solver. The code is designed to handle arbitrary shift and lapse conditions, which can be chosen as appropriate for a given spacetime simulation. The AEI Numerical Relativity group [170] is currently developing robust gauge conditions for (vacuum) black hole spacetimes (see, e.g., [7] and references therein). Hence, all advances accomplished here can be incorporated in future versions of the code for non-vacuum spacetime evolutions. Similarly, since it is a general purpose code, a number of different boundary conditions can be imposed for either the spacetime variables or for the hydrodynamical variables. We refer the reader to [96Jump To The Next Citation Point In The Article, 4Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article] for additional details.

The code has been subjected to a series of convergence tests [96Jump To The Next Citation Point In The Article], with many different combinations of the spacetime and hydrodynamics finite differencing schemes, demonstrating the consistency of the discrete equations with the differential equations. The simulations performed in [96Jump To The Next Citation Point In The Article] include, among others, the evolution of equilibrium configurations of compact stars (solutions to the TOV equations) and the evolution of relativistically boosted TOV stars (v =0.87 c) transversing diagonally across the computational domain - a test for which an exact solution exists. In [4Jump To The Next Citation Point In The Article, 5] the improved stability properties of the BSSN formulation of the Einstein equations were compared to the ADM system. In particular, in [4] a number of strongly gravitating systems were simulated, including weak and strong gravitational waves, black holes, boson stars and relativistic stars. While the error growth-rate can be decreased by going to higher grid resolutions, the BSSN formulation requires grid resolutions higher than the ones needed in the ADM formulation to achieve the same accuracy. Furthermore, it was shown in [89Jump To The Next Citation Point In The Article] that the code successfully passed stringent long-term evolution tests, such as the evolution of both spherical and rapidly rotating, stationary stellar configurations, and the formation of an apparent horizon during the collapse of a relativistic star to a black hole. The high accuracy of the hydrodynamical schemes employed has allowed the detailed investigation of the frequencies of radial, quasi-radial and quadrupolar oscillations of relativistic stellar models, and use them as a strong assessment of the accuracy of the code. The frequencies obtained have been compared to the frequencies computed with perturbative methods and in axisymmetric nonlinear evolutions [88Jump To The Next Citation Point In The Article]. In all of the cases considered, the code reproduces these results with excellent accuracy and is able to extract the gravitational waveforms that might be produced during non-radial stellar pulsations.



4 Hydrodynamical Simulations in Relativistic 3 Numerical Schemes3.2 Other techniques

image Numerical Hydrodynamics in General Relativity
José A. Font
http://www.livingreviews.org/lrr-2003-4
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de