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 testfluid 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 [246, 27]. 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 [244] for details), developed at the Albert Einstein Institute (PotsdamGolm, Germany) and at Louisiana State University (Baton Rouge, Louisiana). This publicdomain code provides highlevel 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 generalrelativistic hydrodynamics code gr_astro described in [137, 131], and developed mostly by the group at Washington University (St. Louis, Missouri). The main features of the whisky code are (a) cellreconstruction procedures with minmod and monotonized central (MC) slope limiters, along with PPM and ENO methods; (b) fullwave decomposition and incomplete Riemann solvers to compute the numerical fluxes such as HLLE, Roelike, and Marquina flux formula; (c) analytic expression for the right and lefteigenvectors of fluxvector Jacobian matrices [176] and compact flux formulae [10] for the RoeRiemann solver and Marquina’s flux formula; (d) a “method of lines” (MoL) approach for the implementation of conservative, highorder time evolution schemes; (e) the possibility to couple the GRHD equations with a conformallydecomposed threemetric in the BSSN formulation[363, 39]. It is worth noting that whisky also incorporates excision methods adapted to HRSC schemes for studying blackhole 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 [28, 305].
CoCoNuT: This code, described in detail in [97], was designed with the aim of studying generalrelativistic 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 cellreconstruction 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 [97] 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 (Poissonlike) 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 [98, 305]. Extensions of the CF approximation also implemented in the CoCoNuT code are considered by [74].
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 [353]) and in a conservative way similar to the approach of Section 2.1.3 (see [355]). In the former case, an important difference with respect to the original system given by Equations (20) – (22), is that an equation for the entropy is solved instead of the energy equation. The hydrodynamic equations are integrated using van Leer’s [403] secondorder 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 highorder central schemes are implemented in the code, along with thirdorder parabolic cellreconstruction 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 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 [355] for details particular to the formulation of the spacetime variables and their solution). In addition, [355] uses the transport velocity throughout (instead of ) and in the energy equation the continuity equation is not subtracted (contrary to the approach of [36]). The choice of coordinates in Shibata’s code depends on the dimensionality of the simulation under study. In the full 3D case, Cartesian coordinates 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 , with and . However, the Einstein’s equations are solved in the plane using Cartesian coordinates. To be able to compute derivatives in this plane the code implements the “cartoon” method [6], which makes possible axisymmetric computations with a Cartesian grid. Next, the hydrodynamics equations are rewritten in Cartesian coordinates using appropriate relations at the plane (see [355] for details). Comprehensive tests of the code are described in [353, 355, 361]. Shibata’s code has allowed important breakthroughs in the study of the morphology and dynamics of various generalrelativistic astrophysical problems, such as blackhole 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 neutronstar binaries, a longstanding problem in numerical relativistic hydrodynamics. These applications are discussed in Section 5.
Duez et al. [109]: The code of [109] shares many features with that of Shibata just discussed. As in the case of [355] the Einstein equations are formulated in the BSSN approach and solved using an iterative Crank–Nicholson scheme for the time update and secondorder centered differencing for the spatial derivatives. Correspondingly, the hydrodynamics equations are formulated in the same conservative approach followed by [355], 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 lowdensity 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 [109]. Extensions of this code to account for viscous fluids through the solution of the relativistic Navier–Stokes equations are reported in [107].

Oechslin et al. [299]: As in the CoCoNuT code described above, in the code of [299] the relativistic hydrodynamics equations are solved together with the Einstein equations in the conformallyflat 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 neutronstar matter through two finitetemperature 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 [299] 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 [378, 379], 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, highorder cellreconstruction 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 [46], using the same conservative formulation of Section 2.2.2. However, instead of relying on Riemannsolverbased methods for their solution, the hydrodynamics equations are solved in the Pitt code with a dissipative central scheme designed by [83]. Using this code, short time evolutions of a selfgravitating 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 lightcone approach developed by Winicour et al. [419] or the conformal formalism of Friedrich [142].
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 testfluid 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, characteristicinformationfree 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 fullwavedecomposition approximate Riemann solvers (such as those used in the codes of [201, 24]) 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 [366]: As described in Section 3.1.2 the GRMHD equations are formulated in the code of [366] in a conservative way. For their solution the semidiscrete highresolution central scheme of Kurganov–Tadmor [211] is used (see also [361]), together with a piecewise parabolic interpolation of the primitive variables for the cellreconstruction 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 [366] follow the convenient prescription proposed by [149] and approximate the fourthorder equation by a quadratic equation for which there exists an analytic solution. The magnetic field divergencefree 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 [355] for further details). The code is axisymmetric and, as in the hydrodynamic case, the “cartoon” method is employed [6] 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 gammaray burst mechanisms, as we discuss in Section 5.
Duez et al. [108]: 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 [105, 106]. 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 [109] the Einstein equations are formulated in the BSSN approach and solved using a threestep iterative Crank–Nicholson scheme for the time update, which yields secondorder 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 shockcapturing capabilities. More precisely, the GRMHD evolution equations are solved with the HLL approximate Riemann solver, together with different possibilities for the cellreconstruction step (either MC, CENO, or PPM). It is worth noting that, contrary to the noatmosphere approach used in the purely hydrodynamic code [109] to handle vacuum regions outside stars, the GRMHD code of [108] is not suitable for such an approach and a very small positive density must be maintained outside the stars.

WHISKYMHD: 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 threedimensional 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 whiskyMHD 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 [363, 39]. The GRMHD equations are implemented using the conservative formulation of [24] discussed in Section 3.1.2 and integrated in time employing the method of lines. The code uses an HLLE approximate Riemann solver, with secondorder TVD slope limiters for the cellreconstruction procedure, and the magneticfield divergencefree constraint is enforced using the fluxCT scheme briefly discussed in Section 4.3.
Anderson et al. [12]: The most distinctive feature of this code, which is mainly described in [12] (see also [286, 14] for further details and [13] 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 neutronstar binary mergers and the computation of the gravitational radiation in the wave zone far from the merger [13, 12]. 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 divergencefree constraint. Correspondingly, the Einstein equations, cast in firstorder symmetric hyperbolic form, are solved in the generalized harmonic decomposition.
CerdáDurán et al. [75]: The GRMHD code developed by [75] has been designed to study magnetorotational, relativistic, stellar core collapse [76]. It is an extension of the axisymmetric hydrodynamics code developed by [95] (whose 3D extension constitutes the CoCoNuT code described above), in which magnetic fields are included following the approach laid out in [24]. Einstein’s equations are formulated using the conformallyflat condition, which has proved very accurate for studying rotational core collapse [96, 98] 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 fluxCT scheme for the magneticfield divergencefree constraint. As a first step towards magnetorotational core collapse simulations, an early version of the code assumed a passive (test) magnetic field, a justifiable assumption since weaklymagnetized fluids are present in such astrophysical scenarios. An extension of the code to relax this assumption has recently been finished.
http://www.livingreviews.org/lrr20087 
This work is licensed under a Creative Commons AttributionNoncommercialNo Derivative Works 2.0 Germany License. Problems/comments to 