The Lagrangian particle method SPH, derived independently by Lucy  and Gingold and Monaghan , has shown successful performance in modelling fluid flows in astrophysics and cosmology. Most studies to date consider Newtonian flows and gravity, enhanced with the inclusion of the fluid self-gravity.
In the SPH method a finite set of extended Lagrangian particles replaces the continuum of hydrodynamic variables, the finite extent of the particles being determined by a smoothing function (the kernel) containing a characteristic length scale . SPH-discretized representations of the equations of motion can be obtained after the introduction of mean-smoothed values for functions and for the gradient and divergence operators, as is briefly outlined below. Reviews of the classical SPH equations are abundant in the literature (see, e.g., [267, 271] and references therein).
The main advantage of the SPH method is that it does not require a computational grid, avoiding mesh tangling and distortion. Hence, compared to grid-based finite-volume methods, SPH avoids wasting computational power in multidimensional applications, when, e.g., modelling regions containing large voids. Among the limitations of SPH we note the difficulties in modelling systems with extremely different characteristic lengths and the fact that boundary conditions usually require a more involved treatment than in finite-volume schemes. Comparisons between SPH and HRSC methods have been reported by various authors (see, e.g., [271, 4] and references therein). In particular, the recent comparison study of state-of-the-art hydrodynamic codes carried out by Agertz et al.  has revealed fundamental differences between SPH and grid methods in modelling interacting multiphase fluids. More precisely, dynamic instabilities such as Kelvin–Helmholtz or Rayleigh–Taylor are much better resolved by Eulerian grid-based methods than by SPH techniques, the reason being the appearance of spurious pressure forces on particles in regions where shocks or steep gradients appear.
In the past few years, implementations of SPH to handle (special) relativistic (and even ultrarelativistic) flows have been developed (see, e.g.,  and references therein). However, SPH has been so far scarcely applied in simulating relativistic flows in curved spacetimes. Relevant references include [189, 212, 213, 381, 300, 122, 299].
Following , let us describe the implementation of an SPH scheme in general relativity. Given a function , its mean smoothed value can be obtained from. The kernel is required to satisfy a normalization condition
The smooth approximation of gradients of scalar functions can be written as
Discrete representations of these procedures are obtained after introducing the number density distribution of particles , with the collection of -particles where the functions are known. The previous representations then read:. The time evolution of the final system of ODEs is performed in  using a second-order Runge–Kutta time integrator with adaptive timestep. As in non-Riemann solver based finite-volume schemes, in SPH simulations involving the presence of shock waves, artificial viscosity terms must be introduced as a viscous pressure term .
Recently, Siegler and Riffert  have developed a Lagrangian conservation form of the general-relativistic hydrodynamic equations for perfect fluids (with artificial viscosity) in arbitrary background spacetimes. Within that formulation these authors  have built a general-relativistic SPH code using the standard SPH formalism as known from Newtonian fluid dynamics (in contrast to previous approaches, e.g., [232, 189, 212]). The conservative character of their scheme has allowed the modelling of ultrarelativistic flows including shocks with Lorentz factors as large as 1000.
The basic principle underlying spectral methods consists of transforming the partial differential equations into a system of ordinary differential equations by means of expanding the solution in a series on a complete basis. The mathematical theory of these schemes is presented in [156, 71] (see also  for a recent Living Review article on spectral methods for numerical relativity). Spectral methods are particularly well suited to the solution of elliptic and parabolic equations. Good results can also be obtained for hyperbolic equations as long as no discontinuities appear in the solution. When a discontinuity is present, some amount of artificial viscosity must be added to avoid spurious oscillations. In the specific case of relativistic problems, where coupled systems of elliptic equations (i.e., the Einstein constraint equations) and hyperbolic equations (i.e., hydrodynamics) must be solved, an interesting strategy is to use spectral methods to solve the elliptic equations and HRSC schemes for the hyperbolic ones. Using such combined methods, first results were obtained in one-dimensional supernova collapse simulations, both in the framework of a tensor-scalar theory of gravitation [292, 294] and in general relativity . Multidimensional approaches for core collapse and studies of neutron star dynamics are available in [97, 99].
Following  we illustrate the main ideas of spectral methods considering the quasi-linear one-dimensional scalar equation:
From the numerical point of view, the series is truncated for a suitable value of . Hence, Equation (92), with , can be rewritten as
Finding a solution of the original equation is then equivalent to solving an “infinite” system of ordinary differential equations, where the initial values of the coefficients and are given by the Fourier expansion of .
In the nonlinear case, , spectral methods proceed in a more convoluted way: first, the derivative of is computed in the Fourier space. Then, a step back to the configuration space is taken through an inverse Fourier transform. Finally, after multiplying by in the configuration space, the scheme returns again to the Fourier space.
The particular set of trigonometric functions used for the expansion in Equation (93) is chosen because it automatically fulfills the boundary conditions, and because a fast transform algorithm is available (the latter is no longer an issue for today’s computers). If the initial or boundary conditions are not periodic, Fourier expansion is no longer useful because of the presence of a Gibbs phenomenon at the boundaries of the interval. Legendre or Chebyshev polynomials are, in this case, the most common base of functions used in the expansions (see [156, 71] for a discussion on the different expansions).
Extensive numerical applications using (pseudo) spectral methods have been undertaken by the LUTH Relativity Group at the Observatoire de Paris in Meudon. The group has focused on the study of compact objects, providing stationary initial data for single and binary compact stars and black holes, as well as performing dynamic studies of gravitational collapse to neutron stars and black holes. The Meudon group has developed a fully object-oriented library (based on the C++ computer language) called LORENE  where (multiple domain) spectral methods have been implemented within spherical coordinates. This library is currently being used by a growing number of numerical relativity groups worldwide, as can be inferred from the comprehensive summary of applications in general-relativistic astrophysics presented in  (see also  for an earlier review). Additional numerical relativity codes based on spectral methods have been developed by Ansorg and collaborators, who have built remarkably accurate stationary data for relativistic stars and black holes [21, 22], by the Cornell/Caltech group, who are focused on the black-hole–binary problem, and by Faber et al. , who have investigated the black-hole–neutron-star binary system for a conformally-flat metric employing LORENE and an SPH code to solve for the hydrodynamics (see below).
The finite-element method is the preferred choice for solving multidimensional structural-engineering problems since the late 1960s . First steps in the development of the finite-element method for modeling astrophysical flows in general relativity are discussed by Meier . The method, designed in a fully covariant manner, is valid not only for the hydrodynamic equations but also for the Einstein and electromagnetic field equations. The most unique aspect of the approach presented in  is that the grid can be arbitrarily extended into the time dimension. Therefore, the standard finite-element method generalizes to a four-dimensional covariant technique on a spacetime grid, with the engineer’s “isoparametric transformation” becoming the generalized Lorentz transform. At present, the method has shown its suitability for accurately computing the equilibrium stellar structure of Newtonian polytropes, either spherical or rotating. The main limitation of the finite-element method, as Meier points out , is that it is generally fully implicit, which results in severe computer time and memory limitations for three and four-dimensional problems.
Richardson and Chung  proposed recently the flow field-dependent variation (FDV) method as a new approach for general relativistic (non ideal) hydrodynamics computations. In the FDV method, parameters characteristic of the flow field are computed in order to guide the numerical scheme toward a solution. The basic idea is to expand the conservation flow variables into a Taylor series in terms of the FDV parameters, which are related to variations of physical parameters such as the Lorentz factor, the relativistic Reynolds number and the relativistic Froude number. The main drawback of the FDV method is the dependence of the solution procedure on a large number of problem-dependent parameters. Richardson and Chung  implemented the FDV method in a C++ code called GRAFSS (General-Relativistic Astrophysical Flow and Shock Solver). The test problems reported are the special relativistic shock tube (problem 1 in the notation of ) and the Bondi accretion on to a Schwarzschild black hole. While their method yields the correct wave propagation, the numerical solution near discontinuities has considerably more diffusion than with upwind HRSC schemes. Nevertheless, the generality of the approach is worth considering.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.