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 (
or even less). However, if more than
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, 146] and references therein). The reader is addressed to [146] 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 [106, 112, 113, 207].

Following [112] let us describe the implementation of an SPH scheme in general
relativity. Given a function
, its mean
*smoothed*
value
can be obtained from

where
*W*
is the smoothing kernel,
*h*
the smoothing length and
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 , with , a normalization function, and 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
, with
being the collection of
*N*
-particles where the functions are known. The previous
representations then read:

with . These approximations can be used to derive the SPH version of the general relativistic hydrodynamic equations. Explicit formulae are reported in [112]. The time evolution of the final system of ODEs is performed in [112] 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 [121].

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.

Following [36] we illustrate the main ideas of spectral methods considering the quasi-linear one-dimensional scalar equation:

with
*u*
=
*u*
(*t*,
*x*), and
a constant parameter. In the linear case (), 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. (66), 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 coefficients
and
are given by the Fourier expansion of
*u*
(*x*,0).

In the non-linear case,
, 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
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 [87] as well as the three-dimensional Newtonian simulations of Bonazzola and Marck et al. [37, 122].

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