Go to previous page Go up Go to next page

4.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. This section contains a brief description of both approaches. In addition, we also mention the field-dependent variation method and the finite element method. We note, however, that both of these approaches have barely been used so far in relativistic hydrodynamics.

4.2.1 Smoothed particle hydrodynamics

The Lagrangian particle method SPH, derived independently by Lucy [231] and Gingold and Monaghan [152], 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 h. 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., [267271Jump To The Next Citation Point] 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., [271Jump To The Next Citation Point4Jump To The Next Citation Point] and references therein). In particular, the recent comparison study of state-of-the-art hydrodynamic codes carried out by Agertz et al. [4] 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., [79Jump To The Next Citation Point] and references therein). However, SPH has been so far scarcely applied in simulating relativistic flows in curved spacetimes. Relevant references include [189Jump To The Next Citation Point212Jump To The Next Citation Point213381Jump To The Next Citation Point300Jump To The Next Citation Point122Jump To The Next Citation Point299Jump To The Next Citation Point].

Following [212Jump To The Next Citation Point], let us describe the implementation of an SPH scheme in general relativity. Given a function f(x), its mean smoothed value ⟨f (x )⟩,(x = (x,y,z)) can be obtained from

∫ ∘ -- ⟨f(x )⟩ ≡ W (x,x ′;h )f (x′) g′d3x′, (85 )
where W is the smoothing kernel, h the smoothing length and √g-′d3x′ 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 [153]. The kernel is required to satisfy a normalization condition
∫ ∘ -- W (x,x ′;h) g′d3x′ = 1, (86 )
which is assured by choosing W (x, x′;h) = ξ(x)Ω (v ), with v = |x − x′|∕h, ξ(x) a normalization function and Ω (v), a standard spherical kernel.

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

⟨∇f (x)⟩ = ∇ ⟨f(x )⟩ − ⟨f(x )⟩∇ ln ξ(x), (87 )
and the approximation of the divergence of a vector reads
⟨∇ ⋅ A (x )⟩ = ∇ ⋅ ⟨A (x)⟩ − ⟨A (x )⟩ ⋅ lnξ(x ). (88 )

Discrete representations of these procedures are obtained after introducing the number density distribution of particles ∑N √ -- ⟨n(x)⟩ ≡ a=1 δ(x − xa)∕ g, with {xa}a=1,...,N the collection of N-particles where the functions are known. The previous representations then read:

∑N f(xb) ⟨f(xa)⟩ = ξ(xa) -------Ωab, (89 ) b=1 ⟨n(xb)⟩ ∑N f(x ) ⟨∇f (xa)⟩ = ξ(xa) ----b--∇xa Ωab, (90 ) b=1 ⟨n(xb)⟩ ∑N A (x ) ⟨∇ ⋅ A (xa)⟩ = ξ(xa) -----b-⋅ ∇xaΩab, (91 ) b=1 ⟨n(xb)⟩
with Ω ≡ Ω(x ,x ;h) ab a b. These approximations can then be used to derive the SPH version of the general-relativistic hydrodynamic equations. Explicit formulae are reported in [212Jump To The Next Citation Point]. The time evolution of the final system of ODEs is performed in [212Jump To The Next Citation Point] 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 [232Jump To The Next Citation Point].

Recently, Siegler and Riffert [381Jump To The Next Citation Point] 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 [381Jump To The Next Citation Point] 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., [232189212Jump To The Next Citation Point]). The conservative character of their scheme has allowed the modelling of ultrarelativistic flows including shocks with Lorentz factors as large as 1000.

4.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 [156Jump To The Next Citation Point71Jump To The Next Citation Point] (see also [159Jump To The Next Citation Point] 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 [292294Jump To The Next Citation Point] and in general relativity [293]. Multidimensional approaches for core collapse and studies of neutron star dynamics are available in [97Jump To The Next Citation Point99Jump To The Next Citation Point].

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

∂u ∂2u ∂u --- = ---2 + λu ---, t ≥ 0, x ∈ [0,1], (92 ) ∂t ∂x ∂x
with u = u(t,x) and λ a constant parameter. In the linear case (λ = 0), and assuming the function u to be periodic, spectral methods expand the function in a Fourier series:
∑∞ u(x, t) = [ak(t)cos(2πkx ) + bk(t)sin(2πkx )]. (93 ) k=0

From the numerical point of view, the series is truncated for a suitable value of k. Hence, Equation (92View Equation), with λ = 0, can be rewritten as

dak 2 dbk 2 ----= − k ak(t), ----= − k bk(t). (94 ) dt dt

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 ak and bk are given by the Fourier expansion of u(x, 0).

In the nonlinear case, λ ⁄= 0, spectral methods proceed in a more convoluted way: first, the derivative of u is computed in the Fourier space. Then, a step back to the configuration space is taken through an inverse Fourier transform. Finally, after multiplying ∂u ∕∂x by u in the configuration space, the scheme returns again to the Fourier space.

The particular set of trigonometric functions used for the expansion in Equation (93View Equation) 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 [15671] 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 [227Jump To The Next Citation Point] 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 [159] (see also [55] 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 [2122], by the Cornell/Caltech group, who are focused on the black-hole–binary problem, and by Faber et al. [122Jump To The Next Citation Point], 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).

4.2.3 Going further

The finite-element method is the preferred choice for solving multidimensional structural-engineering problems since the late 1960s [436]. First steps in the development of the finite-element method for modeling astrophysical flows in general relativity are discussed by Meier [253Jump To The Next Citation Point]. 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 [253Jump To The Next Citation Point] 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 [253], 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 [335Jump To The Next Citation Point] 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 [335] 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 [239Jump To The Next Citation Point]) 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.

  Go to previous page Go up Go to next page