3.3 State-of-the-art three-dimensional codes3 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. 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.

3.2.1 Smoothed particle hydrodynamics 

The Lagrangian particle method SPH, derived independently by Lucy [158] and Gingold and Monaghan [104], 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 multi-dimensional applications, when, e.g., modelling regions containing large voids. Experience in Newtonian hydrodynamics shows that SPH produces very accurate results with a small number of particles (tex2html_wrap_inline5070 or even less). However, if more than tex2html_wrap_inline5072 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, may exceed the capabilities of current supercomputers. 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.

Reviews of the classical SPH equations are abundant in the literature (see, e.g., [187, 191Jump To The Next Citation Point In The Article] and references therein). The reader is addressed to [191Jump 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., [57] and references therein). However, SPH has been scarcely applied to simulate relativistic flows in curved spacetimes. Relevant references include [137Jump To The Next Citation Point In The Article, 146Jump To The Next Citation Point In The Article, 147, 271Jump To The Next Citation Point In The Article].

Following [146Jump 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_inline5074, its mean smoothed value tex2html_wrap_inline5076 can be obtained from


where W is the smoothing kernel, h the smoothing length, and tex2html_wrap_inline5082 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 [105]. The kernel is required to satisfy a normalization condition,


which is assured by choosing tex2html_wrap_inline5084, with tex2html_wrap_inline5086, tex2html_wrap_inline5088 being a normalization function, and tex2html_wrap_inline5090 a standard spherical kernel.

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


and the approximation of the divergence of a vector reads


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


with tex2html_wrap_inline5098 . These approximations can then be used to derive the SPH version of the general relativistic hydrodynamic equations. Explicit formulae are reported in [146Jump To The Next Citation Point In The Article]. The time evolution of the final system of ODEs is performed in [146Jump 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 [160Jump To The Next Citation Point In The Article].

Recently, Siegler and Riffert [271Jump To The Next Citation Point In The Article] 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 [271] 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., [160, 137, 146]). The conservative character of their scheme has allowed the modelling 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 [110Jump To The Next Citation Point In The Article, 52Jump 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 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 have been obtained in one-dimensional supernova collapse simulations, both in the framework of a tensor-scalar theory of gravitation [209, 211Jump To The Next Citation Point In The Article] and in general relativity [210].

Following [41Jump 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_inline5102 a constant parameter. In the linear case (tex2html_wrap_inline5104), and assuming the function u to be periodic, spectral methods expand the function into a Fourier series:


From the numerical point of view, the series is truncated for a suitable value of k . Hence, Equation (59Popup Equation), with tex2html_wrap_inline5104, 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 tex2html_wrap_inline5112 and tex2html_wrap_inline5114 are given by the Fourier expansion of u (x, 0).

In the nonlinear case, tex2html_wrap_inline5118, 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 tex2html_wrap_inline5122 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 (60Popup 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 [110, 52] 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, as well as the associated violent phenomena of gravitational collapse and supernova explosion. They have developed a fully object-oriented library (based on the C++ computer language) called LORENE [157] to implement (multi-domain) spectral methods within spherical coordinates. A comprehensive summary of applications in general relativistic astrophysics is presented in [41]. The most recent ones deal with the computation of quasi-equilibrium configurations of either synchronized or irrotational binary neutron stars in general relativity [112, 284, 283]. Such initial data are currently being used by fully relativistic, finite difference time-dependent codes (see Section  3.3.2) to simulate the coalescence of binary neutron stars.

3.2.3 Flow field-dependent variation method 

Richardson and Chung [240Jump To The Next Citation Point In The Article] have recently proposed 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 general relativistic hydrodynamic equations are expanded in a special form of the Taylor series:


with tex2html_wrap_inline5126 and tex2html_wrap_inline5128 denoting the first-order and second-order variation parameters. Using the above expressions, the time update then reads:


Combining the conservation form of the equations and the rearranged Taylor series expansion, allows us to rewrite the time update into standard matrix (residual) form, which can then be discretized using either standard finite difference or finite element methods [240Jump To The Next Citation Point In The Article].

The physical interpretation of the coefficients tex2html_wrap_inline5126 and tex2html_wrap_inline5128 is the foundation of the FDV method. The first-order parameter tex2html_wrap_inline5126 is proportional to tex2html_wrap_inline5136, where tex2html_wrap_inline5138 is the convection Jacobian representing the change of convective motion. If the Lorentz factor remains constant in space and time, then tex2html_wrap_inline5140 . However, if the Lorentz factor between adjacent zones is large, tex2html_wrap_inline5142 . Similar assessments in terms of the Reynolds number can be provided for the diffusion and diffusion gradients, while the Froude number is used for the source term Jacobian tex2html_wrap_inline5144 . Correspondingly, the second-order FDV parameters tex2html_wrap_inline5128 are chosen to be exponentially proportional to the first-order ones.

Obviously, the main drawback of the FDV method is the dependence of the solution procedure on a large number of problem-dependent parameters, a drawback shared to some extent by artificial viscosity schemes. Richardson and Chung [240Jump To The Next Citation Point In The Article] have implemented the FDV method in a C++ code called GRAFSS (General Relativistic Astrophysical Flow and Shock Solver). The test problems they report are the special relativistic shock tube (problem 1 in the notation of [166Jump To The Next Citation Point In The Article]) and the Bondi accretion onto 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. Applications to non-ideal hydrodynamics and relativistic MHD are in progress [240].

3.2.4 Going further 

The finite element method is the preferred choice to solve multi-dimensional structural engineering problems since the late 1960s [318]. First steps in the development of the finite element method for modeling astrophysical flows in general relativity are discussed by Meier [176Jump To The Next Citation Point In The Article]. 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 [176Jump To The Next Citation Point In The Article] 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 to accurately compute the equilibrium stellar structure of Newtonian polytropes, either spherical or rotating. The main limitation of the finite element method, as Meier points out [176], is that it is generally fully implicit, which results in severe computer time and memory limitations for three- and four-dimensional problems.

3.3 State-of-the-art three-dimensional codes3 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