Go to previous page Go up Go to next page

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

As mentioned in the introduction of the article, there have been a number of important key advances in computational relativistic astrophysics and numerical relativity in the very recent past. Not only black-hole–binary simulations have become finally possible but also nonvacuum spacetimes have started to be investigated numerically with increasing levels of complexity in the input physics, namely through the incorporation of both magnetic fields and microphysical, nonzero temperature equations of state. In addition, the steady increase in computational power is turning the realm of three-dimensional simulations routine. All such advances rely on state-of-the-art (in terms of accuracy and efficiency) numerical codes, the most salient of which are briefly summarized in this section. In the spirit of the topic of the article we focus the discussion on those codes, which deal with solving the equations of relativistic hydrodynamics and MHD for either background or dynamic spacetime metrics. Therefore, we do not consider any of those various codes and approaches, which have made it possible to solve numerically the black-hole–binary problem, the recent major breakthrough of numerical relativity (unless they can also handle matter sources). We instead direct the interested reader to [323Jump To The Next Citation Point] and references therein for an up-to-date overview of such approaches. We also consider only multidimensional (3D) codes.

4.4.1 Hydrodynamics

We start with those purely hydrodynamic codes which also consider dynamic spacetimes. Hydrodynamic codes, which only solve for the matter dynamics on fixed background metrics (quite more numerous), are listed in Table 1, where we give their main basic features and provide the relevant references for further reading. Correspondingly, Table 2 lists test-fluid GRMHD codes, which could also, in principle, deal with unmagnetized flows by simply setting to zero the magnetic field contribution. In addition, the reader is directed to the Living Review article of [240] where similar tables can be found for the particular case of special relativistic hydrodynamics.

CACTUS/WHISKY: The whisky code was the result of a collaboration among several European institutions [24627Jump To The Next Citation Point]. This code solves the general relativistic hydrodynamics equations formulated as in Section 2.1.3 on a 3D numerical grid with Cartesian coordinates. The code has been constructed within the framework of the Cactus Computational Toolkit (see [244Jump To The Next Citation Point] for details), developed at the Albert Einstein Institute (Potsdam-Golm, Germany) and at Louisiana State University (Baton Rouge, Louisiana). This public-domain code provides high-level facilities such as parallelization, input/output, portability on different platforms and several evolution schemes to solve general systems of partial differential equations. Special attention is dedicated to the solution of the Einstein equations, whose matter terms in nonvacuum spacetimes are handled by the whisky code. The initial development of whisky benefitted in part from the release of a public version of the general-relativistic hydrodynamics code gr_astro described in [137Jump To The Next Citation Point131Jump To The Next Citation Point], and developed mostly by the group at Washington University (St. Louis, Missouri). The main features of the whisky code are (a) cell-reconstruction procedures with minmod and monotonized central (MC) slope limiters, along with PPM and ENO methods; (b) full-wave decomposition and incomplete Riemann solvers to compute the numerical fluxes such as HLLE, Roe-like, and Marquina flux formula; (c) analytic expression for the right and left-eigenvectors of flux-vector Jacobian matrices [176] and compact flux formulae [10] for the Roe-Riemann solver and Marquina’s flux formula; (d) a “method of lines” (MoL) approach for the implementation of conservative, high-order time evolution schemes; (e) the possibility to couple the GRHD equations with a conformally-decomposed three-metric in the BSSN formulation[363Jump To The Next Citation Point39Jump To The Next Citation Point]. It is worth noting that whisky also incorporates excision methods adapted to HRSC schemes for studying black-hole formation, as described in [166]. In particular, a modified PPM reconstruction scheme was developed to handle the boundary of the excised region inside apparent horizons. The accuracy of the simulations performed with the whisky code can also benefit from the mesh refinement driver implemented in the cactus code, called carpet, as recently demonstrated in the gravitational collapse studies of [28Jump To The Next Citation Point305Jump To The Next Citation Point].

CoCoNuT: This code, described in detail in [97Jump To The Next Citation Point], was designed with the aim of studying general-relativistic astrophysical scenarios such as rotational core collapse to neutron stars and black holes, as well as pulsations and instabilities of the formed compact objects. Contrary to most 3D codes, which are based on Cartesian coordinates, the CoCoNuT code employs spherical coordinates. Upwind HRSC methods are used for the hydrodynamic equations in the formulation of Section 2.1.3 (Marquina’s flux formula and PPM cell-reconstruction procedures) and spectral methods (through the LORENE library) for solving the metric equations which are formulated in the conformal flatness (CF) approximation (an approach, which [97Jump To The Next Citation Point] coined Mariage des Maillages, i.e., French for grid wedding). In this approximation, and under the maximal slicing condition, the 3+1 ADM equations reduce to a set of five coupled elliptic (Poisson-like) nonlinear equations for the metric components (lapse function, shift vector and conformal factor), which are optimally suited to be solved using spectral methods (see Section 4.2.2). The CF approximation has proved to agree remarkably well with BSSN in simulations performed with the whisky code for rotating neutron stars and stellar core collapse [98Jump To The Next Citation Point305Jump To The Next Citation Point]. Extensions of the CF approximation also implemented in the CoCoNuT code are considered by [74Jump To The Next Citation Point].

Shibata’s GRHD code: In this code the hydrodynamics equations are formulated both in a nonconservative way (following Wilson’s approach of Section 2.1.2; see [353Jump To The Next Citation Point]) and in a conservative way similar to the approach of Section 2.1.3 (see [355Jump To The Next Citation Point]). In the former case, an important difference with respect to the original system given by Equations (20View Equation) – (22View Equation), is that an equation for the entropy is solved instead of the energy equation. The hydrodynamic equations are integrated using van Leer’s [403Jump To The Next Citation Point] second-order finite difference scheme with artificial viscosity, following the approach of a precursor code developed by [303]. For the conservative system, both upwind HRSC schemes and high-order central schemes are implemented in the code, along with third-order parabolic cell-reconstruction procedures. As mentioned before, the choice of conserved variables in the conservative formulation is, however, slightly different to that of Section 2.1.3, mainly due to the presence of a common factor e6φ in all quantities, φ being a spacetime conformal factor, motivated by the formulation used in the code for solving Einstein’s equations. Those are written using the BSSN formulation and solved with finite differences (see [355Jump To The Next Citation Point] for details particular to the formulation of the spacetime variables and their solution). In addition, [355Jump To The Next Citation Point] uses the transport velocity throughout (instead of &tidle;vi) and in the energy equation the continuity equation is not subtracted (contrary to the approach of [36Jump To The Next Citation Point]). The choice of coordinates in Shibata’s code depends on the dimensionality of the simulation under study. In the full 3D case, Cartesian coordinates (x,y, z) are employed. While this is also the case in axisymmetry there is an important subtle difference worth mentioning; the hydrodynamics equations are first written using cylindrical coordinates (ϖ, φ,z), with ϖ = √x2-+--y2- and φ = arctan (y∕x). However, the Einstein’s equations are solved in the y = 0 plane using Cartesian coordinates. To be able to compute y derivatives in this plane the code implements the “cartoon” method [6Jump To The Next Citation Point], which makes possible axisymmetric computations with a Cartesian grid. Next, the hydrodynamics equations are rewritten in Cartesian coordinates using appropriate relations at the y = 0 plane (see [355Jump To The Next Citation Point] for details). Comprehensive tests of the code are described in [353Jump To The Next Citation Point355Jump To The Next Citation Point361Jump To The Next Citation Point]. 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, dynamic instabilities of rotating neutron stars, and, perhaps most importantly, the coalescence of neutron-star binaries, a long-standing problem in numerical relativistic hydrodynamics. These applications are discussed in Section 5.

Duez et al. [109Jump To The Next Citation Point]: The code of [109Jump To The Next Citation Point] shares many features with that of Shibata just discussed. As in the case of [355Jump To The Next Citation Point] the Einstein equations are formulated in the BSSN approach and solved using an iterative Crank–Nicholson scheme for the time update and second-order centered differencing for the spatial derivatives. Correspondingly, the hydrodynamics equations are formulated in the same conservative approach followed by [355Jump To The Next Citation Point], yet they are not solved using HRSC schemes, but through an artificial viscosity (either quadratic or linear), which allows it to handle shock waves. A noteworthy feature of this code is the incorporation of an algorithm, which makes unnecessary the addition of a low-density atmosphere required in most Eulerian hydrodynamic codes, either Newtonian or relativist, when studying the dynamics of matter sources of compact support (as, e.g., the pulsations of stars). Technical details of boundary conditions and gauge choices, along with a comprehensive list of tests passed by the code, are available in [109Jump To The Next Citation Point]. Extensions of this code to account for viscous fluids through the solution of the relativistic Navier–Stokes equations are reported in [107].


Table 1: State-of-the-art GRHD multidimensional numerical codes for test-fluid evolutions.
Code name or authors Main features
   
Riemann Solver HRSC
   
EM95 [117Jump To The Next Citation Point] 2D. Conservative formulation. Linearized Riemann solver based
on Roe averaging. Second-order accurate.
   
ToniK [138Jump To The Next Citation Point] 2D. Conservative formulation. Roe-like, HLLE, and Marquina solvers.
Minmod, ENO and PPM cell-reconstruction procedures.
   
genesis [9Jump To The Next Citation Point, 7Jump To The Next Citation Point] 3D. Conservative formulation. Roe-like, HLLE, and Marquina solvers.
PPM cell-reconstruction procedures.
   
pizaa [188Jump To The Next Citation Point] 3D. Conservative formulation. Modified HRSC scheme optimized
for quasistationary neutron-star spacetimes.
   
Symmetric HRSC
   
cosmos [19Jump To The Next Citation Point] 3D. Conservative energy formulation. MUSCL-type piecewise linear
cell-reconstruction. Second-order accuracy.
   
wham [396Jump To The Next Citation Point] 3D. Conservative total energy formulation well suited for high Mach
number flows. Fifth-order weighted ENO scheme for cell-reconstruction.
Runge–Kutta time stepping.
   
Artificial Viscosity
   
HSW84 [170Jump To The Next Citation Point, 171Jump To The Next Citation Point]; CW [73] 2D. Nonconservative formulation. Second-order accurate.
   
cosmos [19] 3D. Internal energy formulation. Both staggered and nonstaggered
AV methods. Dimensional splitting. Second-order spatial finite differencing.
Single-step time integration.
   
SPH
   
LMZ93 [212] SPH. Specific internal energy equation. AV extra terms in Euler
and energy equations. Second-order time stepping with Runge–Kutta
scheme.
   
SR00 [381] SPH. Conservative formulation. Consistent AV formulation.

Oechslin et al. [299Jump To The Next Citation Point]: As in the CoCoNuT code described above, in the code of [299Jump To The Next Citation Point] the relativistic hydrodynamics equations are solved together with the Einstein equations in the conformally-flat approximation. The code has evolved from an early Newtonian version, which was designed with the definite aim of studying the merger of neutron star binaries. It has gradually been improved to provide ever more accurate descriptions of such mergers. In its latest version the code incorporates microphysical treatment of neutron-star matter through two finite-temperature nuclear EOSs along with a modern treatment to extract gravitational waveform information. Compared to the other codes mentioned before, the most noticeable feature of the code of [299Jump To The Next Citation Point] is the use of SPH techniques to solve for the relativistic hydrodynamics equations, which are formulated in conservation form. The code implements the new SPH artificial viscosity formalism of [79] in which the artificial viscosity interaction is determined by approximately solving a Riemann problem between particle neighbors. Binary merger simulations performed with this code are discussed in Section 5.

Characteristic numerical relativity codes: Although most numerical relativity codes are based on the Cauchy problem, there also exist a couple of characteristic numerical relativity codes, which can handle hydrodynamics. On the one hand, there is the axisymmetric code of [378Jump To The Next Citation Point379Jump To The Next Citation Point], which has been used to perform dynamic evolutions of neutron stars, both pulsating and those formed after a core collapse. In this code the hydrodynamics equations are implemented following the conservative approach outlined in Section 2.2.2 and solved using upwind HRSC schemes, high-order cell-reconstruction procedures, and conservative Runge–Kutta schemes for the time update. Regarding full 3D, the Pittsburgh characteristic vacuum code has recently been upgraded to account for perfect fluid matter sources [46Jump To The Next Citation Point], using the same conservative formulation of Section 2.2.2. However, instead of relying on Riemann-solver-based methods for their solution, the hydrodynamics equations are solved in the Pitt code with a dissipative central scheme designed by [83Jump To The Next Citation Point]. Using this code, short time evolutions of a self-gravitating star in close orbit around a black hole are discussed in [46].

As an important note we point out that all codes based on the Cauchy approach place the outer boundary of the grid at a finite distance, which may potentially lead to numerical problems caused by unwanted reflections of outgoing waves back into the computational domain. Further work on the development of sophisticated boundary conditions is needed to solve (or at least ameliorate) this type of problem. Alternative solutions, which allow for compactified spacetimes that include future null infinity on the computational grid, are provided by the light-cone approach developed by Winicour et al. [419Jump To The Next Citation Point] or the conformal formalism of Friedrich [142].

4.4.2 Magnetohydrodynamics

As in the previous Section 4.4.1, we focus on magnetohydrodynamic codes, which also consider dynamic spacetimes. GRMHD codes which “only” solve for the matter dynamics on fixed background metrics are listed in Table 2, where we give their main basic features and provide references for further reading. As mentioned before, the development of such test-fluid GRMHD codes has witnessed spectacular progress in the last few years. This has been possible not only thanks to the various conservative formulations, which have been put forward but also, and perhaps most importantly, because of the use of accurate, characteristic-information-free HRSC schemes. Such schemes have rendered possible systematic investigations of strongly magnetized scenarios of relativistic astrophysics, which had thus far remained out of numerical reach. We note that in the characterization followed in Table 2 to classify the various existing codes we have made the explicit distinction between full-wave-decomposition approximate Riemann solvers (such as those used in the codes of [201Jump To The Next Citation Point24Jump To The Next Citation Point]) and incomplete Riemann solvers in which only a subset of wave speeds are used. In this regard, we place the widely used HLL (and HLLE) scheme in the latter group, together with symmetric schemes, despite the fact that HLL was originally developed as an approximate Riemann solver scheme [164] with a particularly simple numerical flux.

Shibata’s GRMHD code [366Jump To The Next Citation Point]: As described in Section 3.1.2 the GRMHD equations are formulated in the code of [366Jump To The Next Citation Point] in a conservative way. For their solution the semi-discrete high-resolution central scheme of Kurganov–Tadmor [211] is used (see also [361]), together with a piecewise parabolic interpolation of the primitive variables for the cell-reconstruction scheme. In order to obtain the maximum characteristic speeds needed for the central scheme, instead of solving the quartic equation in the characteristic polynomial, Shibata and Sekiguchi [366Jump To The Next Citation Point] follow the convenient prescription proposed by [149Jump To The Next Citation Point] and approximate the fourth-order equation by a quadratic equation for which there exists an analytic solution. The magnetic field divergence-free constraint is enforced with a constrained transport scheme. Correspondingly, and as in the purely hydrodynamic code, Einstein’s equations are also written in the GRMHD code using the BSSN formulation and solved with finite differences (see [355Jump To The Next Citation Point] for further details). The code is axisymmetric and, as in the hydrodynamic case, the “cartoon” method is employed [6Jump To The Next Citation Point] for evolving the BSSN equations, which makes possible the use of Cartesian coordinates in axisymmetric simulations, and a cylindrical grid is used for the GRMHD equations. It has recently been applied in the study of the collapse of magnetized hypermassive neutron stars to black holes in the context of gamma-ray burst mechanisms, as we discuss in Section 5.

Duez et al. [108Jump To The Next Citation Point]: The development of this code has proceeded simultaneously to that of [366], to the point that they share some important features and both have been used to study similar astrophysical problems, whose results have been presented in joint publications [105Jump To The Next Citation Point106]. This code solves the GRMHD equations in both two dimensions (axisymmetry) and three dimensions, using a Cartesian grid for the latter case and the “cartoon” method in the former (although a cylindrical grid is used for the induction and MHD evolution equations). As in the purely hydrodynamic code of the same group [109Jump To The Next Citation Point] the Einstein equations are formulated in the BSSN approach and solved using a three-step iterative Crank–Nicholson scheme for the time update, which yields second-order accuracy in time. The main modification of the GRMHD code with respect to its predecessor is the complete reformulation of the hydrodynamic sector, whose accuracy has been improved through the use of shock-capturing capabilities. More precisely, the GRMHD evolution equations are solved with the HLL approximate Riemann solver, together with different possibilities for the cell-reconstruction step (either MC, CENO, or PPM). It is worth noting that, contrary to the no-atmosphere approach used in the purely hydrodynamic code [109Jump To The Next Citation Point] to handle vacuum regions outside stars, the GRMHD code of [108Jump To The Next Citation Point] is not suitable for such an approach and a very small positive density must be maintained outside the stars.


Table 2: State-of-the-art GRMHD multidimensional numerical codes for test fluid evolutions.
Code Main features
   
Full-wave decomposition Riemann Solver HRSC
   
K05 [199, 201Jump To The Next Citation Point] 2D. Special relativistic Riemann solver extended to GRMHD following the approach
of [320Jump To The Next Citation Point]. Right eigenvectors in primitive variables. No left eigenvectors. Jumps in
characteristic variables from physical conditions across discontinuities. Switch
between eigenvector sets for degenerate and nondegenerate states. Extra artificial
viscosity and resistivity. 2nd-order accuracy in space and time. CT scheme for
divergence-free constraint.
   
A06 [24Jump To The Next Citation Point] 2D. Right and left eigenvectors in conserved variables from covariant ones. Single set
of left/right eigenvectors for both degenerate and nondegenerate states. 2nd-order
accuracy in space and time. Extension to GRMHD follows the approach of [320Jump To The Next Citation Point].
  Flux-CT scheme for divergence-free constraint.
   
Symmetric HRSC / Incomplete Riemann solvers
   
K98 [196Jump To The Next Citation Point] 2D. Symmetric TVD scheme with nonlinear numerical dissipation [83Jump To The Next Citation Point]. 2nd-order accuracy
in space and time. No method to preserve the magnetic field divergence.
   
HARM [149Jump To The Next Citation Point] 2D. Kerr–Schild-type coordinates (Kerr black-hole spacetime). HLL approximate Riemann
solver. MC, minmod, and van Leer slope limiters. Flux-CT scheme for divergence-free
constraint. Approximation for maximum wave speeds (quadratic equation). 2nd-order
accuracy in space and time.
   
cosmos++ [20Jump To The Next Citation Point] 3D. Conservative energy formulation. Kurganov–Tadmor central scheme. Minmod, van
Leer, and superbee slope limiters. Unstructured grids and local AMR capabilities.
Divergence cleaning of magnetic field constraint (parabolic or hyperbolic). 2nd-order
accuracy.
   
raishin [265] 3D. Follows conservative formulation of [24Jump To The Next Citation Point]. HLL approximate Riemann solver. Flux-CT
scheme for divergence-free constraint. Cell reconstruction schemes: minmod, MC, CENO,
and PPM. 2nd and 3rd-order TVD Runge–Kutta schemes for time update.
   
ECHO [91Jump To The Next Citation Point] 2D. HLL approximate Riemann solver. High-order reconstruction methods: minmod, MC,
and ENO-like schemes (up to 5th-order accurate for smooth flows). Upwind constrained
transport for divergence-free constraint. Runge–Kutta time stepping. Able to handle the
limiting case of magnetodynamics.
   
Artificial Viscosity
   
DVH03 [86Jump To The Next Citation Point] 3D. Nonconservative scheme. Boyer–Lindquist coordinates (Kerr black-hole spacetime).
Staggered grid, time-explicit, operator-split finite difference scheme. Method of
characteristics CT scheme for divergence-free constraint.
   
cosmos++ [20] 3D. Internal energy formulation and dual energy∕flux-conserving formulation. Latter
accurate for ultrarelativistic flows.
   

WHISKY-MHD: As the name indicates, this code is the extension of the hydrodynamic code whisky discussed above to solve for the full set of GRMHD equations in a dynamic spacetime. The code is described in [151]. It is a fully three-dimensional code in Cartesian coordinates, based on HRSC techniques on domains with adaptive mesh refinement, accomplished as in the whisky code through the AMR driver implemented in the cactus code, called carpet. The whisky-MHD code uses, hence, the cactus computational toolkit which provides the infrastructure for the parallelization and the I/O of the code, along with methods for the solution of the Einstein equations formulated in the BSSN approach [36339]. The GRMHD equations are implemented using the conservative formulation of [24Jump To The Next Citation Point] discussed in Section 3.1.2 and integrated in time employing the method of lines. The code uses an HLLE approximate Riemann solver, with second-order TVD slope limiters for the cell-reconstruction procedure, and the magnetic-field divergence-free constraint is enforced using the flux-CT scheme briefly discussed in Section 4.3.

Anderson et al. [12Jump To The Next Citation Point]: The most distinctive feature of this code, which is mainly described in [12Jump To The Next Citation Point] (see also [28614] for further details and [13Jump To The Next Citation Point] for the assessment of the purely hydrodynamic module of the code), is the fact that it uses a sophisticated computational infrastructure, which provides distributed AMR. This has allowed the investigation of both unmagnetized and magnetized neutron-star binary mergers and the computation of the gravitational radiation in the wave zone far from the merger [13Jump To The Next Citation Point12Jump To The Next Citation Point]. As mentioned in Section 3.1.2 the code uses a conservative formulation of the GRMHD equations. These are solved by means of the HLLE approximate Riemann solver with PPM cell reconstruction, in conjunction with divergence cleaning to control the divergence-free constraint. Correspondingly, the Einstein equations, cast in first-order symmetric hyperbolic form, are solved in the generalized harmonic decomposition.

Cerdá-Durán et al. [75Jump To The Next Citation Point]: The GRMHD code developed by [75Jump To The Next Citation Point] has been designed to study magneto-rotational, relativistic, stellar core collapse [76Jump To The Next Citation Point]. It is an extension of the axisymmetric hydrodynamics code developed by [95Jump To The Next Citation Point] (whose 3D extension constitutes the CoCoNuT code described above), in which magnetic fields are included following the approach laid out in [24Jump To The Next Citation Point]. Einstein’s equations are formulated using the conformally-flat condition, which has proved very accurate for studying rotational core collapse [96Jump To The Next Citation Point98Jump To The Next Citation Point] and are solved using spectral methods. In the code, both the HLLE approximate Riemann solver and the Kurganov–Tadmor central scheme are implemented to solve for the MHD evolution equations, along with the flux-CT scheme for the magnetic-field divergence-free constraint. As a first step towards magneto-rotational core collapse simulations, an early version of the code assumed a passive (test) magnetic field, a justifiable assumption since weakly-magnetized fluids are present in such astrophysical scenarios. An extension of the code to relax this assumption has recently been finished.


  Go to previous page Go up Go to next page