4 Simulations of Astrophysical Phenomena3 Numerical Schemes3.1 Finite difference schemes

3.2 Other techniques 

Two of the most frequently employed alternatives to finite difference schemes in numerical hydrodynamics are Smoothed Particle Hydrodynamics (SPH) and, to a lesser extent, (pseudo-) Spectral Methods.

3.2.1 Smoothed particle hydrodynamics 

The Lagrangian particle method SPH, derived independently by Lucy [120] and Gingold and Monaghan [80], has shown successful performance to model 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 hydrodynamical variables, the finite extent of the particles being determined by a smoothing function (the kernel) containing a characteristic length scale h . The main advantage of this 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., modeling regions containing large voids. Experience in Newtonian hydrodynamics shows that SPH produces very accurate results with a small number of particles (tex2html_wrap_inline3881 or even less). However, if more than tex2html_wrap_inline3883 particles have to be used for certain problems, and self-gravity has to be included, the computational power, which grows as the square of the number of particles, would exceed the capabilities of current supercomputers. Among the limitations of SPH we note the difficulties in modeling systems with extremely different characteristic lengths and the fact that boundary conditions usually require a more involved treatment than in finite volume schemes.

Reviews of the classical SPH equations are abundant in the literature (see, e.g., [144, 146Jump To The Next Citation Point In The Article] and references therein). The reader is addressed to [146Jump To The Next Citation Point In The Article] for a summary of comparisons between SPH and HRSC methods.

Recently, implementations of SPH to handle (special) relativistic (and even ultrarelativistic) flows have been developed (see, e.g., [51] and references therein). However, SPH has been scarcely applied to simulate relativistic flows in curved spacetimes. Relevant references include [106Jump To The Next Citation Point In The Article, 112Jump To The Next Citation Point In The Article, 113, 207Jump To The Next Citation Point In The Article].

Following [112Jump To The Next Citation Point In The Article] let us describe the implementation of an SPH scheme in general relativity. Given a function tex2html_wrap_inline3885, its mean smoothed value tex2html_wrap_inline3887 can be obtained from


where W is the smoothing kernel, h the smoothing length and tex2html_wrap_inline3893 the volume element. The kernel must be differentiable at least once, and the derivative should be continuous to avoid large fluctuations in the force felt by a particle. Additional considerations for an appropriate election of the smoothing kernel can be found in [81]. The kernel is required to satisfy a normalization condition,


which is assured by choosing tex2html_wrap_inline3895, with tex2html_wrap_inline3897, tex2html_wrap_inline3899 a normalization function, and tex2html_wrap_inline3901 a standard spherical kernel.

The smooth approximation of gradients of scalar functions can be written as


and that corresponding to the divergence of a vector reads


Discrete representations of these procedures are obtained after introducing the number density distribution of particles tex2html_wrap_inline3903, with tex2html_wrap_inline3905 being the collection of N -particles where the functions are known. The previous representations then read:


with tex2html_wrap_inline3909 . These approximations can be used to derive the SPH version of the general relativistic hydrodynamic equations. Explicit formulae are reported in [112Jump To The Next Citation Point In The Article]. The time evolution of the final system of ODEs is performed in [112Jump To The Next Citation Point In The Article] using a second-order Runge-Kutta time integrator with adaptive time step. 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 [121Jump To The Next Citation Point In The Article].

Recently, Siegler and Riffert [207] 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 they have built a general relativistic SPH code using the standard SPH formalism as known from Newtonian fluid dynamics (in contrast to previous approaches [121, 106, 112]). The conservative character of their scheme allows the modeling of ultrarelativistic flows including shocks with Lorentz factors as large as 1000.

3.2.2 Spectral methods 

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 [86Jump To The Next Citation Point In The Article, 46Jump To The Next Citation Point In The Article]. 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 obtain a smooth solution. 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 solution strategy is to use spectral methods for the elliptic equations and HRSC schemes for the hyperbolic ones. Using such combined methods the first results have been obtained in one-dimensional supernova collapse in the framework of a tensor-scalar theory of gravitation [160, 161].

Following [36Jump To The Next Citation Point In The Article] we illustrate the main ideas of spectral methods considering the quasi-linear one-dimensional scalar equation:


with u = u (t, x), and tex2html_wrap_inline3913 a constant parameter. In the linear case (tex2html_wrap_inline3915), and assuming the function u to be periodic, spectral methods expand the function in a Fourier series:


From the numerical point of view, the series is truncated for a suitable value of k . Hence, Eq. (66Popup Equation), with tex2html_wrap_inline3915, 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 coefficients tex2html_wrap_inline3923 and tex2html_wrap_inline3811 are given by the Fourier expansion of u (x,0).

In the non-linear case, tex2html_wrap_inline3929, spectral methods, or, in this case, more correctly pseudo-spectral methods, proceed in a more convoluted way: First, the derivative of u is computed in the Fourier space, then they come back to the configuration space by an inverse Fourier transform, multiply tex2html_wrap_inline3933 by u in the configuration space and, finally, come back again to the Fourier space.

The particular set of trigonometric functions used for the expansion is chosen because it automatically fulfills the boundary conditions and because a fast transform algorithm is available. 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 [86, 46] for a discussion on the different expansions).

A comprehensive summary of applications of (pseudo-) spectral methods in general relativistic astrophysics is presented in [36]. Among those involving time-dependent hydrodynamical simulations we mention the spherically symmetric collapse computations of Gourgoulhon [87Jump To The Next Citation Point In The Article] as well as the three-dimensional Newtonian simulations of Bonazzola and Marck et al. [37Jump To The Next Citation Point In The Article, 122].

4 Simulations of Astrophysical Phenomena3 Numerical Schemes3.1 Finite difference schemes

image Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de