Assessment of the MUSTA approach for numerical relativistic hydrodynamics
^{1} Laboratory for Scientific Computing, Department of Physics, Cavendish Laboratory, J J Thomson Avenue, Cambridge, CB3 0HE, UK
email: pmb39@cam.ac.uk; nn10005@cam.ac.uk
^{2} Department of Mathematical Sciences, Rensselaer Polytechnic Institute, NY 12180 Troy, USA
email: henshw@rpi.edu
Received: 20 October 2014
Accepted: 25 December 2014
Aims. We evaluate some approximations for solving the equations of special relativistic hydrodynamics within complex geometries. In particular, we assess the following schemes: the Generalized FORCE (GFORCE) and MUlti STAge (MUSTA) approaches which are used as the basis for a secondorderaccurate SlopeLImitedCentred (SLIC) method. These do not require detailed knowledge of the characteristic structure of the system, but have the potential to be nearly as accurate as more expensive schemes which do require this knowledge.
Methods. In order to treat complex geometries, we use multiple overlapping grids which allow the capturing of complex geometries while retaining the efficiencies associated with structured grids.
Results. The schemes are evaluated using a suite of one dimensional problems some of which have known exact solutions, and it is shown that the schemes can be used at CFL numbers close to the theoretical stability limit. We compare the effects of the MUSTA approach when applied to two different schemes. The scheme is further validated on a number of problems involving complex geometries with overlapping grids.
Key words: methods: numerical / hydrodynamics / relativistic processes / shock waves
© ESO, 2015
1. Introduction
The numerical solution of the relativistic hydrodynamical equations is of importance to the simulation of astrophysical phenomena such as gammaray bursts, supernova corecollapse, and relativistic wind accretion. Although approximations can be obtained by the use of a Newtonian fluid alone, relativistic effects must be taken into account in order to make accurate predictions that can be tested against observations.
These problems often require calculations to be performed in three dimensions on nonCartesian grids (for reasons to be explained later) and are thus computationally expensive. An ideal numerical scheme for evolving such problems would therefore have comparatively low computational cost, and be easily incorporated into existing codes capable of using threedimensional, nonCartesian grids.
There are various approaches to evolving the relativistic hydrodynamical equations. The most popular are finitevolume based schemes, but other methods such as smoothed particle hydrodynamics (Rosswog 2010), spectral methods (van Putten 1993) and the relativistic beam scheme (Sanders & Prendergast 1974) have been employed in the past. A full review of the field can be found in Martí & Müller (2003). Highorder and highaccuracy methods are often based on either an exact or an approximate Riemann solver, coupled with a highorder reconstruction procedure. However, the full solution of the relativistic hydrodynamical Riemann problem requires knowledge of the wave structure and the characteristics, and these are algebraically complicated, expensive to calculate, and their form depends on the equation of state. Therefore it could be advantageous to have a numerical scheme that does not require knowledge of the characteristic structure.
It has been demonstrated by Eulderink & Mellema (1995) that Riemann solver based methods are capable of providing accurate evolutions of special relativistic problems. However, we show here that the computational expense associated with a Riemann solver based method may not, in fact, be necessary, as a similar accuracy can be achieved using a method that requires less information about the characteristic structure. We further show that we can choose the timestep based on a CFL number of around 0.95, as opposed to values such as 0.5 as used by Del Zanna & Bucciantini (2002).
There do exist numerical schemes which do not require a Riemann solver or any other information about the equations, save for the conserved variables and fluxes given in the conservation law form, and the speed of the fastest moving wave, in order to maintain stability. These schemes are typically centred in that they do not take into account the direction of travel of the waves. In contrast, methods based on approximate or exact Riemann solutions tend to upwind the appropriate characteristic variables.
Our objective is to use numerical methods which are both accurate and robust and can be used for practical largescale calculations. Following the arguments above, an ideal method would have the advantages of centred schemes in terms of computational expense, and ease of implementation, and would also approach the accuracy and robustness of a traditional Riemann solver based method. Such a method should then be suitable for the solution of more complex systems of equations, such as those of special relativistic hydrodynamics (SRHD) or general relativistic hydrodynamics (GRHD), and could be readily incorporated into existing multidimensional frameworks.
To this end, the MUltiSTAged (MUSTA) approach proposed by Toro & Titarev (2006) has been considered as a methodology which meets the above objectives. This takes a basic centred scheme and attempts to improve the solution by repeatedly solving the local Riemann problem in a predictorcorrector fashion. The same paper proposed the Generalized FORCE (GFORCE) scheme, which relies on estimates of the local wavespeeds, and which can provide better results than the MUSTA approach when applied to the FORCE scheme. In particular, when implemented for SRHD flows, the GFORCE scheme appears to combine the accuracy of upwind schemes with the simplicity of centred methods. The MUSTA approach also displays these properties, but we give an explicit demonstration here as to why care should be taken when using the MUSTA approach.
The idea behind the MUSTA approach is to solve the local Riemann problem numerically rather than analytically, by means of a first order centred method which is simple and computationally inexpensive as compared to approximate or exact analytic solutions, but which makes full use of local wavepropagation speed information. The MUSTA approach seeks to improve upon this by evolving the local Riemann problem multiple times in a predictorcorrector approach, with the expectation that this will lead to a more accurate approximation to the flux, depending on the number of corrector steps used. The number of steps taken by the MUSTA approach can be varied, and the effect of using different numbers of steps upon the accuracy and robustness are explored in this paper.
We use the special relativistic system as a steppingstone towards the solution of the full GRHD system, to be presented in another communication. We use the SRHD system to evaluate these numerical methods since exact solutions are known for this system, and there are existing results showing its solution with various other numerical schemes. Therefore, our numerical solutions can be compared to both exact solutions and to existing numerical schemes.
The intention is that the code being developed here for SRHD will be extended to the evolution of a fluid on an arbitrary spacetime metric, specifically, in the vicinity of one or more black holes. The presence of a singularity in the metric leads to a severe reduction in the maximum stable timestep permitted for the explicit methods being considered. One approach to avoid this timestep restriction is to excise the singularity from the domain, ensuring that the excision boundary is within the apparent horizon of the black hole, so that causality will not allow any boundary effects to propagate into the region outside the horizon. Further, we must specify outer boundary conditions at some large distance from the black hole since it is hard to determine suitable boundary conditions in regions where the spacetime is not approximately flat. One way to address both these issues is to use curvilinear overlapping grids. These grids will allow the excision boundary to be smooth, thus reducing the likelihood of generating spurious numerical effects there. The approach also allows a spherical grid to be used at large distances from the black hole, which generally requires fewer grid points to reach a suitably large distance for the outer boundary than would a Cartesian grid. We have therefore based our code on Overture as described by Brown et al. (1997), which is a framework for developing numerical PDE solvers on curvilinear overlapping grids capable of efficiently treating complex, possibly moving, geometries.
In this paper we demonstrate the use of the MUSTA approach as developed by Toro & Titarev (2006) to solve the SRHD equations. Although the methodology is simpler, than traditional Riemann Problem based ones, we assess whether it is capable of producing results which are at least as accurate as previously used schemes. Further, we show that it can be implemented within a curvilinear overlapping grid framework, and give some sample results.
The rest of this paper is organised as follows: in Sect. 2, a discrete form of the SRHD equations in conservation law form is formulated. We outline the finitevolume approach in Sect. 3 and the FORCE and SlopeLImitedCentred (SLIC) schemes in Sect. 3.1. The MUSTA and GFORCE approaches are then described in Sects. 3.2 and 3.3, and some notation is defined. This includes some analysis on the MUSTA approach to explore its limitations when applied to the FORCE scheme, and we suggest that it may be more robust for GFORCE. In Sect. 4 we describe the Overture software and its applicability to astrophysical fluid dynamics. In Sect. 5 we discuss how to adapt the SLIC scheme to evolution on curvilinear grids.
In Sect. 6 we validate the multistaged FORCE and GFORCE schemes in 1D against exact solutions and comparisons to a well known approximate Riemann solver. We include a demonstration of an evolution on a nonuniform 1D grid. In Sect. 7.1 we demonstrate the evolution of a 1D problem across an overlapping curvilinear grid system.
We then proceed to some truly twodimensional problems, and in Sect. 7.2 evaluate our code against the fourquadrant problem proposed by Del Zanna & Bucciantini (2002). Finally, in Sects. 7.3–7.5, we demonstrate the application of our code to some examples where overlapping curvilinear grids are required for their accurate solution.
We then present our conclusions in Sect. 8 and evaluate the potential of this method for the solution of more complex systems of equations.
2. Equations of SRHD
Throughout this paper, we use units where the speed of light, c = 1. Greek indices run over space and time: μ,ν,... = 0,1,2,3, and Roman indices run over space only: i,j,... = 1,2,3. We use Minkowski space with metric η_{μν} = diag( − 1, + 1, + 1, + 1).
The evolution of a relativistic fluid in flat space is governed by the conservation law (1)where (2)In the case of an ideal gas, the conserved fluid variables U are defined in terms of the primitive fluid variables as (3)where ρ is the fluid density, p is the pressure, v^{i} is the velocity vector, and W is the Lorentz factor. The specific enthalpy h = h(ρ,p) characterises the heat content of a fluid per unit volume, and it defines the equation of state for the fluid when written as a function of the pressure and density. The conserved variables are the relativistic density D, the threemomentum S_{i}, and the relativistic energy E, although we evolve τ, being the difference of the total relativistic energy and the relativistic density, in place of E to allow for more accurate evolutions as otherwise variations in p can be swamped by the relatively large relativistic density D, following Martí & Müller (2003).
Given an equation of state defined through the specific enthalpy h, the energy momentum tensor for a perfect fluid is given by (4)where the fourvelocity u^{μ} is defined in terms of the threevelocity v^{i} via (5)Here, W = (1 − v^{i}v_{i})^{− 1/2}, calculated for Minkowski spacetime. The specific enthalpy of an ideal fluid is (6)where Γ (assumed to be constant) is the adiabatic index of the fluid (often written as γ, but this symbol is reserved for the determinant of the threemetric). The speed of sound in the fluid, c_{s}, is then given by (7)In the case of nonrelativistic idealgas, the recovery of the primitive variables from the conserved variables is an algebraic operation. However, for relativistic fluids the presence of the Lorentz factor complicates the issue, since the threemomentum components are no longer independent of each other, so we use a NewtonRaphson method to recover the primitive variables. Other techniques are available, including the solution of a quartic, but tend to be computationally involved and expensive. The NewtonRaphson approach we use follows the one given in Appendix A of Eulderink & Mellema (1995). We chose this over the approach of Martí & Müller (2003) as in our experience the latter appeared susceptible to numerical roundoff errors when evaluating W, and appeared to lead to problems if the conserved variables had accumulated errors that took them outside the physically valid range. Eulderink & Mellema’s method, on the other hand, requires no calculation of W and, although it may be susceptible to issues of numerical precision, or lead to very small densities and pressure, we believe that we have rectified these in a robust fashion, as follows:
Given the conserved variables given above, we solve the following for ξ in order to recover the primitive variables:
We note that, in order to retain numerical accuracy, the expression for calculating u^{0} = W in Eq. (A.12) of Eulderink & Mellema (1995) should be rewritten as (9)taking C_{0}ξ outside the brackets, since for large W the bracketed expression is of order unity and C_{0}ξ is large. We then use the following to recover the remaining primitive variables: (10)where we have chosen ρ_{min} = 10^{8}. Given that ρ may have been altered from its correct value of D/W, we recalculate the appropriate value of p from h, enforce p ≥ p_{min} if necessary, and then recalculate h before calculating v_{i} as (11)According to Eulderink & Mellema (1995), “with some effort it may be shown that the derivative of [the quartic (8)] is always at least 2/Γ in its physical root, thus guaranteeing quadratic convergence for the Newton[Raphson] method”. This gives us confidence that the approach is robust, and will remain so even when applied to more extreme flows, in cases where numerical errors are more likely, such as the evolution of a fluid on a strongly curved spacetime metric.
3. HRSC schemes and MUSTA
Finite volume schemes (FV) are often used for the solution of equations in conservationlaw form. The solution is approximated using averages over grid cells, and is advanced in time by the addition and subtraction of fluxes of the conserved quantities across the grid cell edges: (12)where we have a gridspacing Δx, a timestep Δt, is an approximation to the solution expressed in conserved variables in gridcell i at timestep n, and is an approximation to the flux of conserved quantities between cells i − 1 and i at timestep n.
Since the finite volume approach is based on the integral form of the conservation law, it is suited to problems where discontinuities are expected. The sine qua non of finitevolume schemes is therefore the accurate approximation of the fluxes of the conserved quantities through the cell edges. The way in which the solution is updated using Eq. (12)ensures conservation, but care must be taken to ensure that no spurious oscillations arise in the vicinity of shocks.
A subset of finite volume schemes is the set of highresolution shockcapturing (HRSC) schemes. These are capable of capturing shocks with very little smearing, even at low resolution, and are of at least secondorder accuracy in smooth regions of the flow.
3.1. FORCE and SLIC schemes
One of the most basic FV schemes is the First ORder CEntred (FORCE) scheme. It can be applied to any set of PDEs in conservation form (omitting source terms) with little effort. The flux approximation is calculated as follows: (13)
This is firstorder accurate, and is a centred scheme, meaning that it does not have a (solution dependent) bias to information from either the left or the right. As will become relevant later, FORCE is the simple average of two other numerical schemes: the LaxWendroff scheme, being the term , and the LaxFriedrichs scheme, being the term . The LaxWendroff scheme is secondorder but oscillatory, and the LaxFriedrichs scheme is firstorder, but very diffusive. It turns out that the FORCE scheme is the least diffusive nonoscillatory threepoint firstorder linear scheme that is not dependent on local solution information. The FORCE scheme is based on the piecewise constant reconstruction of the solution, i.e. the left and right states are taken directly from the solution vectors and .
One improvement over the constant piecewise reconstruction, known as the SLIC scheme (SlopeLImited Centred), is to use cells i − 1 and i + 2 to make a piecewise linear reconstruction of the solution at time n + 1/2 to provide initial states for the Riemann problem. However, we must take care that, in reconstructing the solution, we do not introduce any new extrema into the solution. We first reconstruct the slope of the solution, Δ_{i}, using data from either side of the current cell, taking the average of the forward and backward differences: (14)and then determine a limited slope by multiplying each component by a factor dependent on the local solution behaviour: (15)where (16)The preceding two equations are applied to each individual component of Δ_{i}, , and to form the limited slope . The SLIC flux is then calculated as
where we use and as the left and right states for the FORCE flux.
There are many choices for the slopelimiter ξ, even given the constraint that the overall numerical scheme should be Total Variation Diminishing (TVD) as defined by Toro (1999). If we restrict the limiter to be a function of the ratio of the forward and backward differences: (21)then the constraints on ξ(r) can be shown to be (Toro 1999, p. 505, setting ω = 0): where (22)with C_{CFL} being the local Courant number.
We see that the ratio r is a measure of the change in slope around cell i. It is negative if i lies at an extremum of the solution, as one of and will be positive and the other negative, and this is the reason that we require ξ(r) = 0 for r< 0 as any nonconstant reconstruction would result in an overshoot or undershoot of the extremum.
The slopelimiter that we will use in this paper is the vanLeer limiter given by: (23)This limiter is a good compromise between compressive (superbee) and diffusive (minmod) limiters, and we have found it suitable for our purposes, not leading to excessive smearing of shocks, nor any spurious oscillations.
3.2. MUltiSTAging approach (MUSTA)
When considering the accuracy of a numerical scheme, one must take into account both the order of accuracy (for example secondorder in both space and time) and the “errorconstant”. In other words, if the solution error is locally of the form (24)then the accuracy depends both on p and q and also the constant c_{i}. Thus the accuracy of the scheme can be improved by decreasing the error constant without decreasing the order of accuracy.
In order to improve on the accuracy of SLIC, we could consider Riemann solvers (as an alternative to FORCE to calculate the firstorder solution) which use information about the wave types making up the solution to calculate the flux. These tend to be computationally intensive, often requiring iterative or matrix methods to determine the flux at the cell face. In order to increase the order of accuracy, we could use a higherorder reconstruction than the linear one that SLIC uses. Again, this would require more computational expense. However, Toro (2003) has suggested a simple way of improving the accuracy of any numerical scheme, which can produce improved accuracy for less computational effort than would be required if the accuracy were to be attained by increasing resolution. This is called MUSTA, which stands for MUltiSTAging and can be understood as follows.
Fig. 1 How to multistage a simple numerical scheme (e.g. FORCE). The Riemann problem is solved for several time steps using the simple scheme on the twocell domain, with transmissive boundary conditions as shown by the empty arrows. The solid lines show the wave structure of the solution. As the number of stages increases, the waves move away from the centre, so that the flux at the centre tends to the correct limit, which is represented by the filled arrow at the top. This figure corresponds to calculating the FORCE^{3} scheme, if FORCE is used as the underlying simple scheme. 

Open with DEXTER 
The FORCE flux provides a simple way of estimating the flux at a cell face. We can define an approximate Riemann solver by using the FORCE (or any other) flux to evaluate the initial flux and then seeking to correct it using multiple evaluations of the flux on the resulting solution, and using the resulting flux on the face to define an improved numerical flux for the original scheme. The MUSTA approach uses this idea to numerically evolve the Riemann problem on a two cell grid, using some appropriate scheme, for a few timesteps, or stages, and to take the central flux at the last stage as the estimate for the final flux. This is represented in Fig. 1. When using the FORCE scheme as the basis for MUSTA, the MUSTA algorithm can be written as the iterative scheme where and give the initial conditions based on the original Riemann problem. This scheme uses the transmissive boundary conditions where the solution has zero derivative at the boundary.
The best improvements in accuracy for least computational expense are likely to be found when multistaging is applied to as simple and appropriate scheme as possible. Therefore, in order to improve the accuracy of SLIC, we do not multistage SLIC itself, but rather multistage the FORCE scheme on which SLIC is based. We write the multistaged FORCE scheme as FORCE^{k}, where k is the number of stages, so that FORCE^{1} ≡ FORCE. This notation is as in Toro (2003) but differs from that in Toro & Titarev (2006).
Although the preceding explanation seems plausible enough, we can examine the theoretical effect of MUSTA for the linear advection equation. As mentioned in Sect. 3.1, the FORCE scheme can be expressed as (28)with .
The local CFL number c_{L} is derived from the local wave speed and global time step Δt as (29)For linear advection, many numerical schemes with a stencil width of three can be written in the form (28)for some ω, which may depend on the solution. Therefore we can plot the dependence of ω on c_{L} for various schemes, as in Fig. 2. The Godunov scheme corresponds to the choice , the FORCE scheme to , LaxWendroff to ω = 1, and LaxFriedrichs to ω = 0. The Godunov scheme, based on the exact Riemann solution, falls on the curved line shown. This splits the diagram into two parts. Above the curve lie schemes that are oscillatory, and below it lie schemes that are nonoscillatory, but are more diffusive than the Godunov scheme. The firstorder Godunov scheme is the least diffusive nonoscillatory linear firstorder method for the linear advection equation. The FORCE scheme can therefore be seen to be the least diffusive nonoscillatory firstorder method that does not depend on the local wave speed information given by c_{L}.
Fig. 2 Demonstration of how different firstorder schemes for solving the linearadvection equation can be represented in the form 28 where ω depends on c_{L} the local CFL condition. 

Open with DEXTER 
Fig. 3 Same as Fig. 2, with additional multistaged versions of FORCE ^{k} where k = 1...5 and k increases up the page (labels for k = 3,4,5 are omitted for reasons of clarity). 

Open with DEXTER 
Although the MUSTA approach does not depend explicitly on c_{L}, it does depend on it implicitly, through its dependence on the local data for the predictor step (25). This dependence can be calculated in a straightforward, but algebraically intensive manner, which we have done using Maple, and we plot the results for FORCE ^{n}, n = 1,..., 5 in Fig. 3. It is immediately clear that MUSTA may not be a suitable evolution scheme at high CFL numbers, even for linear advection, as all of the schemes now protrude into the oscillatory region. Using very high numbers of stages shows that FORCE^{k} tends to the LaxWendroff scheme as k → ∞ (except for c_{L} = 1 where for all k, so that the convergence is not uniform).
3.3. Generalised FORCE flux
The impetus for the development of MUSTA was the fact that even for the simplest of all conservation laws, the linear advection equation, the FORCE method was more diffusive than the Godunov scheme using the exact Riemann solver for CFL less than unity, and the approach was to improve the accuracy of a scheme, without requiring local wave speed information. However, we have seen that the MUSTA approach is not ideal as it can lead to oscillatory methods, even for linear advection.
Toro & Titarev (2006) presented a solution to this in the same paper as their MUSTA approach. In that paper, they calculate the generalised FORCE flux (GFORCE) as a weighted average of the LaxFriedrichs and LaxWendroff fluxes, now dependent on local wave speed information: (30)where c_{L} is the local CFL number, derived from the local wave speed and global time step Δt as (31)With these definitions, GFORCE replicates the Godunov upwind scheme for the linearadvection case. This suggests that the GFORCE flux could well provide improved accuracy for more complex systems of equations.
We note that the SR wave speeds in one dimension are given by (32)As the GFORCE scheme already coincides with the Godunov scheme for the linear advection equation, applying the MUSTA approach makes no difference in this case, simply maintaining the Godunov scheme. However, this is not the case for more complicated systems of equations, which have multiple wavespeeds, and we shall see that multistaging GFORCE provides some increase in accuracy.
3.4. HLL flux
The GFORCE flux, being equivalent to the flux from the firstorder Godunov solver in this case, is the most accurate nonoscillatory firstorder scheme for the linear advection equation. However, more complex systems of equations have more complex solutions to their Riemann problems. It is possible, in many cases, to derive a linearised or exact Riemann solver for a system of equations. However, this will usually involve either calculating a matrix inverse or using an iterative method to find the flux. This tends to be computationally expensive, and is therefore often avoided unless the payoff of computational time versus accuracy is worthwhile.
However, a simple Riemann solver, applicable to any system of equations, was suggested by Harten et al. (1983), and is known as the HLL Riemann solver. The rationale behind this approach is that if all the waves in the solution are moving in one direction only, then the flux is taken wholly from one side or the other, and an appropriately weighted average is taken otherwise: This is only dependent on approximations to the left and right local wave speeds a_{L} and a_{R}. These can be calculated from the fluid variables as suggested by Schneider et al. (1993) (33)where and are the arithmetic means of the initial velocities and sound speeds in cells i and i + 1.
In the case of a curvilinear grid, we note that we have to transform the variables into the computational grid’s local coordinate system in order to find the wave speed parallel to the grid coordinate direction. We elaborate on this in Sect. 5.
3.5. Recommendations for the MUSTA approach
In summary, therefore, the multistaging approach at first sight appears to produce a straightforward and cheap way of improving the accuracy of a simple flux approximation scheme. However, some investigation shows that applying it to the FORCE scheme alone may not be a viable approach as FORCE for k> 1 are oscillatory schemes for linear advection. An improvement on this appears to be to apply it to the GFORCE scheme, where GFORCE is not oscillatory for linear advection. Likewise, the HLL schemes are not oscillatory.
However, as the analysis has not been extended to more complex systems of equations, our recommendations are as follows: use of FORCE may be advantageous in that it captures peaks and discontinuities with less diffusion as k increases, but this method should be used with care since for k> 1 the scheme no longer satisfies the TVD property for the linear advection equation as illustrated in Fig. 3. Without further analysis the GFORCE and HLL schemes should also be used with caution although they are likely more robust than the FORCE scheme. The use of GFORCE^{k} or HLL^{k} does not, as yet, have any evidence of causing the loss of the TVD property although, given the analysis for FORCE, we suggest caution as to whether this result can be applied to more complex systems of equations.
We attach a caveat emptor here, therefore, to any use of multistaging, but suggest that it may be useful if used with care, and with only a few extra stages, perhaps around three or so. We present some numerical results and comparisons in Sect. 6. In particular, the examples shown in Fig. 7 show the effect of increasing the multistaging parameter.
In terms of a comparison to using an exact or approximate Riemann solver, the relative accuracy will depend on the system being solved. We have not made any direct comparisons with a Riemannsolverbased method, but in general centred schemes tend to smear contact discontinuities substantially more than Riemann solvers designed to capture these exactly. However, other regions of the flow should be captured with similar accuracy if sufficient MUSTA stages are used.
The MUSTA approach permits improved accuracy on top of that provided by FORCE by the simple addition of an extra containing loop within a single flux calculation. No further algebraic derivations are required as for a Riemann solver, and neither are any higherorder reconstructions, such as those required by CENO (Del Zanna & Bucciantini 2002). As demonstrated by Eulderink & Mellema (1995), constructing a Riemann solver takes substantial work.
Therefore, the advantage of our aproach is that it only requires the system fluxes and maximum wavespeeds to be specified, and affords an easy way of improving accuracy with minimal time penalty, without having to derive, implement, and test a full Riemann solver.
In terms of computational cost, we have not made a comparison with a relativistic Riemann solver. However, a comparison between FORCE and an exact Riemann solver for an ideal Eulerian gas in one dimension suggested that FORCE took about a third of the time, FORCE about twothirds of the time, and FORCE about the same time. We would expect that for SRHD, the Riemann solver would take comparatively substantially longer than this, and so the performance gain would be far greater.
4. Overture
Overture^{1} is an infrastructure that supports the development of PDE solvers for overlapping curvilinear grids, as described by Brown et al. (1997). Overture includes support for grid generation, for individual component grids, and for generating the overlapping grid, and associated connectivity information required for interpolation. Users of Overture may utilise existing PDE solvers or implement their own approximations.
Our main reason for using Overture is that we wish to capture a spherical geometry, with a smooth inner excision boundary, so that we can include general relativistic effects and ultimately simulate fluid flow into a black hole. Relevant results for the GRHD equations will be presented in a future paper.
There are two main issues with regard to dealing with overlapping grids. The first is that of constructing the full grid structure. This is one of the tasks that Overture has been designed to perform even for very complicated grids. Given a set of logically rectangular grids, along with mappings to the physical grid coordinates, as well as boundary conditions, Overture can remove parts of grids that are obscured by other grids and calculate which cells on which grids should be used to interpolate data onto grid cells that are now near to cutout portions of the grid.
In order to interpolate solution variables between grids, we have used polynomial interpolation. Overture uses Lagrange polynomials for interpolation, and takes into account the grid geometries involved. One potential drawback of using Lagrange polynomials is that, except for the linear case, it is not guaranteed that the interpolated values will be contained in the same range as the values from which they are interpolated. This could lead to overshoots developing near interpolation boundaries. However, although monotonicity preserving interpolation schemes can be derived, it was not found that the extra expense was worthwhile.
We used a fifthorder accurate interpolation, based on a fivepoint stencil, which required two lines of interpolation points at the edge of each grid. We required two lines of interpolation points as the slopelimiting methods we use require two ghost cells at the grid boundaries.
5. Curvilinear grids
When solving a PDE on overlapping grids, each with its own curvilinear coordinate system, it is usual to transform the governing equations to the new grid coordinate systems but retain the Cartesian components of the dependent variables. This simplifies interpolation between different component grids.
The curvilinear grid is defined by a smooth transformation, x^{i} = φ(r^{i}), from Cartesian coordinates x^{i} to the unit square parameter space coordinates r^{i}. This means that we transform the governing equations as follows (34)where J is the Jacobian of the transformation: (35)The third term in Eq. (34)is added so that at the discrete level we maintain the geometric conservation law (this is also known as freestream preservation). When the flux is not spatially dependent, the second term in Eq. (34)is zero analytically, but may not be zero when evaluated numerically. The third term is added to ensure that, when the flux is constant spatially, the flow remains constant. The term is composed of an evaluation of f on a cellcentred value of u, multiplied by a geometric term. This geometric term must be calculated in a manner consistent with the second term, so that, for a constant state u, the geometric terms will cancel exactly (at least up to roundoff error).
The process of evaluating the third term is fairly straightforward in practice. For any given HRSC scheme, we can do this by starting with a fluid with constant density, velocity, and pressure, and compute the fluxes analytically, ensuring that the flow is preserved, regardless of the geometry or other parameters, such as slopelimiter. In the case of SLIC, the result is that the flux for each cell face is composed of the product of a flux, derived from the given constant state, and a geometry term of the form . Since the final flux for a given coordinate direction r^{i} is then simply a difference of these, the term that cancels these is (36)where the values of the Jacobian and coordinate derivatives have to be calculated in exactly the same way as in the rest of the HRSC scheme being used.
5.1. Curvilinear SLIC scheme
The SLIC scheme requires the reconstruction of states at cellfaces. We have changed the SLIC algorithm slightly to take account of the spatial of the Jacobian on curvilinear grids.
(37)where J_{i} is the Jacobian of cell i, and . The derivatives ∂r^{i}/∂x^{j}, however, are averaged differently since, in one dimension, we know that J^{1} = ∂r/∂x, and so we want to ensure that (38)where there is no sum over j.
One way to ensure that this is always true is to average the derivatives of the inverse mapping using the harmonic mean (39)where the matrix D_{i} is defined as (40)Even in the case of multidimensions, we never need to calculate the fluxes across more than one cellface at a time, so the above averaging is only ever done in one computational grid direction.
5.2. HLL flux
In order to calculate the HLL flux in curvilinear coordinates, we first transform the Cartesian components of the solution variables into the local grid coordinate basis. The wave speeds are then calculated for the required coordinate direction from the transformed variables. In the same way as SLIC, we calculate the fluxes in the local grid coordinate direction, multiplied by the local Jacobian. In order to achieve second order convergence, we adjust the HLL flux to read: (41)where the fluxes are calculated parallel to the local grid coordinate i. We note that the local CFL number defined in Eq. (31)required for the calculation of GFORCE is also derived from the wave speeds calculated along the local grid coordinate directions.
5.3. Determination of the timestep
On Cartesian grids, the maximum timestep that retains stability for the SLIC scheme is given by (42)where C_{CFL} ≤ 1 is the CFL number, Δx is the grid spacing, and S_{max} is the maximum wavespeed, obtained from (32). On curvilinear grids, we determine the global timestep as the minimum of local timesteps determined on a cellbycell basis, where Δx is chosen appropriately from the local cell size. The alternative, taking the minimum cell size and the maximum wave speed over the whole grid, can lead to time steps that are more restrictive than necessary. Further, when we wish to calculate the wavespeeds, we in fact transform the numerical solution into the local grid basis so that the correct wave speeds along the coordinate directions can be found, as these may not be aligned with Cartesian directions, which could lead to an underestimate of the local wavespeed.
However, a drawback of using curvilinear grids is that the cells are not of similar size everywhere. Even with the above method of calculating the timestep, this still leads to an effectively lower CFL number for evolution on larger grid cells, as the timestep has been restricted by the stability condition on smaller grid cells. Therefore waves may become somewhat smeared on a coarser area of the grid only because of finer resolution elsewhere.
We can reduce the effect of this problem by using as accurate a scheme as possible, namely GFORCE or HLL, which provide accurate solutions even at lower CFL numbers.
6. Onedimensional test problems
In order to test the stability and accuracy of the GFORCE scheme and the MUSTA approach, we have applied them to some standard onedimensional testcases which have been used previously in the literature. The tests we consider are carried out using a onedimensional uniform resolution code using the same numerical schemes as the Overture code, but without the overhead imposed by Overture.
An exact solution to the Riemann problem for the onedimensional special relativistic hydrodynamical equations exists, and so we can compare our numerical results to this. In order to generate the exact solutions, we have used the code riemann.f from Martí & Müller (2003).
The testcases we use are taken from Del Zanna & Bucciantini (2002), so that a comparison can be made with their numerical results, which were obtained using a Convex ENO reconstruction and the HLL Riemann solver. The three testcases below all use an ideal fluid with adiabatic index Γ = 5/3 and are performed on a domain of x ∈ [0,1].
6.1. Mildly relativistic shock
The initial conditions for this test are (43)It is fairly easy to obtain good results for this test, as it contains no extreme velocities or densities. The solution is evolved to a time t = 0.4. It has been noticed by other researchers (e.g. Del Zanna & Bucciantini 2002) that, for certain schemes, the initial zero pressure has to be set to some small nonzero value. We use a minimum pressure of p_{min} = 10^{8}, assuming reference pressures of order 1.
Fig. 4 Results for a mildly relativistic shock wave, (ρ,u,p)^{L} = (10, 0, 13.3) and (ρ,u,p)^{R} = (1, 0, 0), using 400 cells, a CFL number of 0.99, and the van Leer limiter, with the SLIC scheme based on the GFORCE^{3} flux, simulated to time t = 0.4. The exact solution is given by the solid line and is composed of a rarefaction, a contact discontinuity, and a shock. 

Open with DEXTER 
An example result is shown in Fig. 4 using the SLIC scheme with the van Leer slope limiter, based upon the GFORCE ^{3} scheme, and using a CFL number of 0.99. The rightmost discontinuity, corresponding to a shock wave, is captured with only four points. In this, and in the capturing of other discontinuities, we see that our results are as good as those of Del Zanna & Bucciantini (2002). However, the capturing of the density is seen to be somewhat better in LucasSerrano et al. (2004). We also note that the kink we find in the velocity at x ≈ 0.6 (near where the solution has a discontinuity in the slope) is present to a slightly lesser extent in LucasSerrano et al. (2004), and not at all in Del Zanna & Bucciantini (2002).
Fig. 5 Results for the density of a mildly relativistic shock wave, (ρ,u,p)^{L} = (10, 0, 13.3) and (ρ,u,p)^{R} = (1, 0, 0), using 400 cells, a CFL number of 0.99, and the van Leer limiter, with the SLIC scheme based on the fluxes: FORCE (circles), GFORCE (squares), GFORCE^{3} (crosses), and HLL (+ signs). The exact solution is given by the straight line. We restrict the region plotted so that the differences between the methods can be seen. 

Open with DEXTER 
Figure 5 shows the result near the leading shock using four different numerical schemes for the same problem. All the methods use the SLIC scheme, based on four different fluxes: FORCE, GFORCE, GFORCE^{3}, and HLL. It can be seen that there is very little difference indeed between the methods. This is due to the mild nature of the problem. However, we can make the differences between the methods more apparent by evaluating the ℒ_{1} error of the density for each of several methods, using a range of CFL numbers. The results are shown in Table 1. The error norm is taken over the whole domain and so includes the discontinuities. We see that increasing the number of stages used for a scheme reduces the error. In all cases, the GFORCE scheme is better than the FORCE scheme, although the FORCE^{2} scheme is always better than the GFORCE^{2} scheme, except for the case of C_{CFL} = 0.4, where the numerical evolution blows up. We note this data as evidence for the caveat emptor in Sect. 3.5. The reason for the improvement of FORCE ^{k} over GFORCE is that the FORCE ^{k} scheme is usually above the GFORCE scheme as shown in Fig. 3, as a result of its lessening diffusion at the expense of potentially introducing oscillations into the solution. The GFORCE scheme is of similar accuracy to the HLL scheme, although it is slightly worse at CFL 0.4. There is a slight improvement in accuracy when increasing the number of stages used for the GFORCE scheme. We also note that the error decreases for higher CFL numbers, both as a result of the reduced number of timesteps required to reach the final time and the fact that schemes are often more accurate for CFL numbers close to 1.
However, in general, the FORCE and GFORCE schemes are more accurate than the HLL scheme, at least for k ≥ 2. This is roughly in line with what we would expect from the analysis of the various schemes on the linear advection equation in Sect. 3.
Error norms of the density using ∥ ·∥ _{1} for various base schemes for the SLIC method applied to the mildly relativistic shock test.
6.2. Strongly relativistic BlastWave
The initial data for this problem is given by (44)This is a more difficult testcase to evolve, as the solution contains a narrow peak of high density, which is hard to capture at low resolution, and also contains velocities which approach c = 1. This problem is also evolved to a time t = 0.4. The exact solution for all variables is shown as part of Fig. 6. This shows a narrow region of high density moving close to the speed of light. The maximum density is very difficult to capture at a moderate resolution such as we are using due to the narrowness of the peak. A good test of the accuracy of a scheme, therefore, is the fraction of the peak value that is attained.
An example result, obtained using the SLIC scheme based on GFORCE, using the van Leer limiter, with a CFL of 0.99, and using 400 cells, is shown in Fig. 6. The solution is already quite accurate, and the density peak is at 65.6% of its correct value. Del Zanna & Bucciantini (2002) achieved peaks of 70.1% and 62.4% for two different versions of their method, and LucasSerrano et al. (2004) attained a density peak of 76% of the exact value (the comparison is not quite valid as Del Zanna & Bucciantini (2002) only evolved to time t = 0.35, rather than t = 0.4 as used here and in LucasSerrano et al. 2004).
On comparing different underlying firstorder methods for SLIC, we find that once again there is little visual difference between the results, and we therefore give the density peaks in Table 2. These results are in line with what we have already seen, that FORCE ^{k} tends to be more accurate than GFORCE ^{k} for k ≥ 2, and that HLL is more accurate than GFORCE when no multistaging is used.
Fig. 6 Results for the blastwave Riemann problem given by (44)using 400 cells, a CFL number of 0.99, and the van Leer limiter, with the SLIC scheme based on the GFORCE flux, to time t = 0.4. The exact solution is given by the solid lines, and is composed of a rarefaction, a contact discontinuity, and a shock. The maximum velocity is v ≈ 0.960, corresponding to a Lorentz factor of 3.59. Numerical results are given for pressure (squares), velocity (circles), and density (crosses), all scaled by appropriate factors to fit them all onto one plot. 

Open with DEXTER 
In Table 3 we give the ℒ_{1} errors in the density for various resolutions and methods. This shows near firstorder accuracy over the whole problem, which is the best that we can expect for a solution containing discontinuities.
Density peak for strongly relativistic shock as predicted by various methods at a resolution of 400 cells.
ℒ_{1} norm errors in density for the blastwave test, evaluated over the whole domain.
Fig. 7 Results for the densityperturbation initial conditions. The first plot shows a full solution obtained using the SLIC scheme based on FORCE^{1}. The other plots show closeups of the density when replacing FORCE by other schemes. For all results, a CFL of 0.99 and the van Leer limiter were used, on a grid of 400 cells. 

Open with DEXTER 
6.3. Perturbed density test
The initial data for this problem is given by (45)This problem provides a good indicator of how well a numerical method can capture extrema. There is no analytic solution to this, and instead we use a converged numerical solution, generated at highresolution. This test is evolved to t = 0.35, since the results of Del Zanna & Bucciantini (2002) are consistent with this, not t = 0.4 as stated in their main text. The converged solution, evolved using the SLIC ^{1} scheme with van Leer limiter, and 20 000 cells, is shown as the first part of Fig. 7. We note that this is visually indistinguishable from the result using 10 000 cells. Superimposed on this is an example evolution, using 400 cells and the SLIC scheme based on FORCE. Points to note are that the three sharp peaks are not particularly well captured, and that the final peak at the right of the figure is clipped. We note that Del Zanna & Bucciantini (2002) manage to capture the peaks somewhat better than we do whilst using half the number of cells.
The rest of Fig. 7 shows the results for various schemes, focusing on the region 0.6 ≤ x ≤ 1. Now we see that increasing the number of stages used improves the capturing of the three peaks and the knees at x ≈ 0.7 and x ≈ 0.85. The ordering of the plots has been chosen to show a progressively increase in accuracy. We see that multistaging GFORCE outperforms HLL, and that multistaged FORCE still produces better accuracy than multistaged GFORCE. However, in our experience, multistaged FORCE appears to be more robust, and we note that although increasing multistaging generally improves accuracy, the payoff after about three stages is small relative to the extra computational cost.
7. Validation in two spatial dimensions
In this section we give details of various validation tests we have carried out in two spatial dimensions, making use of the overlapping grid capabilities of our algorithm.
7.1. Onedimensional shock across overlapping grids
It is important to test the accuracy of the method on curvilinear, overlapping grids since use of these grids will be integral to our later work. To this end, we evolve a onedimensional Riemann problem with a known solution across the twodimensional domain, with outer vertices at (± 0.5, ± 0.5) shown in Fig. 8, where the Cartesian grid has resolution 400 × 400, and the annulus matches its resolution to that of the Cartesian grid at its inner circumference. This is designed to test the accuracy of the interpolation between a curvilinear grid and a regular Cartesian grid, as well as the accuracy of the evolution taking place on an annular geometry, rather than a Cartesian grid. From a physical point of view, this is just evolving the shocktest on a square domain, but the presence of the embedded curvilinear grid adds numerical and algorithmic complexity. As initial conditions we use the Riemann problem (46)so that the shock is propagating diagonally across the grid.
As this is only a 1D test evolved on a 2D grid, the exact solution can be found as in Sect. 6. We perform the evolution using the grid described in Fig. 8, and sample the solution at 400 points (corresponding to the Cartesian grid’s resolution) along the line x = y, perpendicular to the shock direction, at a time of t = 0.4. Figure 9 shows the results from using the GFORCE flux, with a CFL of 0.95 and the van Leer limiter. We can see that the solutions on the plain 2D Cartesian grid and on the overlapping grid are almost indistinguishable.
7.2. The twodimensional Riemann problem
We can also test the multidimensionalisation of our schemes in a simple twodimensional test, essentially a twodimensional Riemann problem. Due to the twodimensional nature of the setup, there is no known exact solution to this problem, but we can compare our results to those of other authors such as Del Zanna & Bucciantini (2002). The problem is defined on a region [0,1] × [0,1] and is evolved to a time t = 0.4, with initial conditions (47)This setup is symmetric about the line x = y, so that the evolved solution should maintain this symmetry.
Figure 10 shows the results for two different schemes, FORCE and GFORCE, both using the van Leer limiter. The resolution used is 400 cells in each direction, using a CFL of 0.95. Although there is no exact solution for this initial data, the results compare well with those of Del Zanna & Bucciantini (2002), in that the same features are captured, and with similar widths. We note that the long density peak in the upperright of the figure is more pronounced in our scheme than in Del Zanna & Bucciantini (2002), but that the feature just below and to the left of the centre is somewhat less welldefined in our plot than in Del Zanna & Bucciantini (2002). It is in fact the case that, within this region, the density drops to a very small value indeed, of the order of 10^{6}. It is also the case, although not entirely clear on the plots shown, that the discontinuities are slightly narrower for the GFORCE scheme than for the FORCE scheme.
Fig. 8 Grid for testing curvilinear and overlapping routines. An annulus is superimposed on a Cartesian grid and used to replace part of that grid. The Cartesian grid has its vertices at (± 0.5, ± 0.5), a resolution of 400 × 400, and the annulus resolution matches that of the Cartesian grid at the inner circumference. We only show every fourth gridline in this plot for clarity. 

Open with DEXTER 
Fig. 9 Crosssection of the evolution of the mildly relativistic shockwave across the grid shown in Fig. 8. The run was performed using the GFORCE^{1} scheme at a CFL of 0.95, and using the van Leer slopelimiter. The green + signs are the numerical solution on the overlapping grid, the red circles are from an equivalent simulation on a single 2D Cartesian grid, and the blue line is the exact solution for the problem. The two numerical solutions have been interpolated from along the line x = y. 

Open with DEXTER 
Fig. 10 Quadrants problem evolved at a resolution of 400 × 400, with a CFL of 0.95, using the FORCE^{1} and GFORCE^{1} schemes with the van Leer slopelimiter. We plot the log of the density ρ, using 30 equally spaced contours from ln0.1 = −2.303 to ln3.927 = 1.368. 

Open with DEXTER 
7.3. Frontfacing step
In order to demonstrate the overlapping grid approach, we apply our code to the problem of a fluid shock forming ahead of a solid frontfacing step, similar to the one presented in LucasSerrano et al. (2004). Since this benchmark has a gridaligned geometry, it can be performed with far simpler meshgeneration approaches, such as by a blockstructured methodology based on orthogonal cells. In order to demonstrate the versatility of Overture, we smooth the corner of the step, as shown in Fig. 11. The fluid has an initial density of 1.4, with an adiabatic index of 1.4. The fluid is initially steady, and the lefthand boundary has an inflow of fluid at velocity v_{x} = 0.8, at Mach 3, and the system is evolved to time t = 10.
The resulting flow pattern is shown in Fig. 12, where we plot isocontours of the log of the density. This compares qualitatively well to the results of LucasSerrano et al. (2004).
7.4. Reflective cylinder
Fig. 11 Grid used for the frontfacing step problem. The grid is made up of three separate grids: one for the incoming channel, one for the rest of the channel, and one for the step itself. The last has stretched grid cells that become narrower towards the edge of the step. For clarity we only show every other grid line. 

Open with DEXTER 
Fig. 12 Solution at t = 10 for a shock hitting a frontfacing step. The fluid has adiabatic index Γ = 1.4, velocity v = 0.8 and is travelling at Mach 3. The scheme used was GFORCE^{1} with the van Leer slopelimiter, and with CFL 0.95. We plot 30 contours of the log of the density. 

Open with DEXTER 
Fig. 13 Solution for a reflective cylinder initially immersed in a constant velocity flow of Mach 3, velocity 0.8, adiabatic index Γ = 1.4, and density 1.4. The main channel has length 10 and width 5, composed of 100 × 200 cells. The annulus has inner radius 0.5 and outer radius 1.25, with 168 cells in the angular dimension, and 32 cells in the radial direction, where an exponential stretching function has been applied in the radial direction to restrict the radial size of cells nearer the inner boundary. For clarity, we only show every other grid line. In the top half of the plot, contours of ρ at time t = 10 are shown. The scheme used was GFORCE^{1}, with the van Leer slopelimiter, at a CFL of 0.95. There are 30 contours, equally spaced between 0 and 10. 

Open with DEXTER 
Since we intend to apply this code to problems in general relativity, we demonstrate a problem that is somewhat similar in its grid setup and initial conditions to that of a fluid impacting on a black hole. In this problem, we immerse a solid cylinder in a constant velocity flow and allow the numerical solution to evolve to steadystate. This general setup has been much studied in the nonrelativistic case, and we expect to see qualitatively similar results here. Our grid is constructed from a Cartesian grid for the channel of size 10 by 5, containing an annulus with inner radius 0.5, the inner circumference of which is a reflective boundary. The channel’s walls are given transmissive boundary conditions, with the exception of the left boundary, which is an inflow boundary. The grid is shown in the bottom half of Fig. 13. The state of the simulation at time t = 10 is shown in Fig. 13. The solution consists of a strong bow shock in front of the cylinder, and there is some evidence of an unstable wake beginning to form downstream, as would be expected. The solution has qualitatively similar characteristics to those exhibited by the nonrelativistic problem.
7.5. Multiple reflective cylinders
In Fig. 14 we show an extension of the problem from the previous section, to the case where there is more than one reflective cylinder in the domain. The same constant state initial data was used, this time incoming from the right. Since the numerical solver we developed can handle general overlapping grids, it was relatively simple to go from the single cylinder case to the multiple cylinder case. We show both the grid as used by Overture and the result which, while it cannot be validated, demonstrates qualitatively plausible behaviour.
Fig. 14 Simulation of flow past two solid cylinders. The first plot shows the grid and in the second we plot 30 contours of ρ from 0 to 10 for two solid cylinders initially immersed in a fluid of constant state, Mach 3 and velocity 0.8, with density 1.4 and adiabatic index 1.4. The solution is shown at time 10. The GFORCE^{1} scheme was used, with the van Leer slopelimiter, and a CFL of 0.95. 

Open with DEXTER 
8. Conclusions and further work
In this paper, we developed a numerical scheme to solve the SRHD equations in complex geometries. The approach was based on the secondorder SLIC method, with fluxes obtained through the MUSTA approach applied to the FORCE and GFORCE schemes. The properties of these new schemes were analysed and it was shown that the MUSTA approach can lead to more accurate results potentially at the cost of being less robust, although in practice the GFORCE scheme coupled with MUSTA presented no difficulties. The new schemes were then compared to each other and to the HLL scheme when applied to a number of onedimensional test problems. We found that the MUSTA approach applied to the GFORCE scheme generally gave more accurate results than the HLL scheme. The scheme was then applied to a number of two dimensional problems using overlapping grids, which demonstrated the flexibility of the scheme to handle complex geometries.
In a future paper, we shall therefore extend our scheme to evolve general relativistic hydrodynamical problems, and present the way that we have adapted the numerical schemes presented in this paper to allow for a nonflat metric. It is our intention to apply our code to problems involving a stationary blackhole metric, in particular, to BondiHoyleLyttleton accretion onto a Kerr black hole.
Freely available from http://www.overtureframework.org
Acknowledgments
This work was funded by an EPSRC Doctoral Training Grant. The authors would like to thank Eleuterio Toro for comments on the implementation of the MUSTA approach.
References
 Brown, D. L., Chesshire, G. S., Henshaw, W. D., & Quinlan, D. J. 1997, in Proc. of the Eighth SIAM Conf. on Parallel Processing for Scientific Computing [Google Scholar]
 Del Zanna, L., & Bucciantini, N. 2002, A&A, 390, 1177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eulderink, F., & Mellema, G. 1995, A&AS, 110, 587 [NASA ADS] [Google Scholar]
 Harten, A., Lax, P. D., & van Leer, B. 1983, SIAM Rev., 25, 35 [CrossRef] [MathSciNet] [Google Scholar]
 LucasSerrano, A., Font, J., Ibáñez, J., & Martí, J. 2004, A&A, 428, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martí, J. M., & Müller, E. 2003, Liv. Rev. Rel., 6, 7 [Google Scholar]
 Rosswog, S. 2010, Class. Quant. Grav., 27, 114108 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H., & Prendergast, K. H. 1974, ApJ, 188, 489 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, V., Katscher, V., Rischke, D., et al. 1993, J. Comput. Phys., 105, 92 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Toro, E. 1999, Riemann solvers and numerical methods for fluid dynamics, 2nd edn. (SpringerVerlag) [Google Scholar]
 Toro, E. 2003, Multistage predictorcorrector fluxes for hyperbolic equations http://www.newton.cam.ac.uk/preprints/NI03037.pdf [Google Scholar]
 Toro, E., & Titarev, V. 2006, J. Comput. Phys., 216, 403 [NASA ADS] [CrossRef] [Google Scholar]
 van Putten, M. H. P. M. 1993, ApJ, 408, L21 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Error norms of the density using ∥ ·∥ _{1} for various base schemes for the SLIC method applied to the mildly relativistic shock test.
Density peak for strongly relativistic shock as predicted by various methods at a resolution of 400 cells.
ℒ_{1} norm errors in density for the blastwave test, evaluated over the whole domain.
All Figures
Fig. 1 How to multistage a simple numerical scheme (e.g. FORCE). The Riemann problem is solved for several time steps using the simple scheme on the twocell domain, with transmissive boundary conditions as shown by the empty arrows. The solid lines show the wave structure of the solution. As the number of stages increases, the waves move away from the centre, so that the flux at the centre tends to the correct limit, which is represented by the filled arrow at the top. This figure corresponds to calculating the FORCE^{3} scheme, if FORCE is used as the underlying simple scheme. 

Open with DEXTER  
In the text 
Fig. 2 Demonstration of how different firstorder schemes for solving the linearadvection equation can be represented in the form 28 where ω depends on c_{L} the local CFL condition. 

Open with DEXTER  
In the text 
Fig. 3 Same as Fig. 2, with additional multistaged versions of FORCE ^{k} where k = 1...5 and k increases up the page (labels for k = 3,4,5 are omitted for reasons of clarity). 

Open with DEXTER  
In the text 
Fig. 4 Results for a mildly relativistic shock wave, (ρ,u,p)^{L} = (10, 0, 13.3) and (ρ,u,p)^{R} = (1, 0, 0), using 400 cells, a CFL number of 0.99, and the van Leer limiter, with the SLIC scheme based on the GFORCE^{3} flux, simulated to time t = 0.4. The exact solution is given by the solid line and is composed of a rarefaction, a contact discontinuity, and a shock. 

Open with DEXTER  
In the text 
Fig. 5 Results for the density of a mildly relativistic shock wave, (ρ,u,p)^{L} = (10, 0, 13.3) and (ρ,u,p)^{R} = (1, 0, 0), using 400 cells, a CFL number of 0.99, and the van Leer limiter, with the SLIC scheme based on the fluxes: FORCE (circles), GFORCE (squares), GFORCE^{3} (crosses), and HLL (+ signs). The exact solution is given by the straight line. We restrict the region plotted so that the differences between the methods can be seen. 

Open with DEXTER  
In the text 
Fig. 6 Results for the blastwave Riemann problem given by (44)using 400 cells, a CFL number of 0.99, and the van Leer limiter, with the SLIC scheme based on the GFORCE flux, to time t = 0.4. The exact solution is given by the solid lines, and is composed of a rarefaction, a contact discontinuity, and a shock. The maximum velocity is v ≈ 0.960, corresponding to a Lorentz factor of 3.59. Numerical results are given for pressure (squares), velocity (circles), and density (crosses), all scaled by appropriate factors to fit them all onto one plot. 

Open with DEXTER  
In the text 
Fig. 7 Results for the densityperturbation initial conditions. The first plot shows a full solution obtained using the SLIC scheme based on FORCE^{1}. The other plots show closeups of the density when replacing FORCE by other schemes. For all results, a CFL of 0.99 and the van Leer limiter were used, on a grid of 400 cells. 

Open with DEXTER  
In the text 
Fig. 8 Grid for testing curvilinear and overlapping routines. An annulus is superimposed on a Cartesian grid and used to replace part of that grid. The Cartesian grid has its vertices at (± 0.5, ± 0.5), a resolution of 400 × 400, and the annulus resolution matches that of the Cartesian grid at the inner circumference. We only show every fourth gridline in this plot for clarity. 

Open with DEXTER  
In the text 
Fig. 9 Crosssection of the evolution of the mildly relativistic shockwave across the grid shown in Fig. 8. The run was performed using the GFORCE^{1} scheme at a CFL of 0.95, and using the van Leer slopelimiter. The green + signs are the numerical solution on the overlapping grid, the red circles are from an equivalent simulation on a single 2D Cartesian grid, and the blue line is the exact solution for the problem. The two numerical solutions have been interpolated from along the line x = y. 

Open with DEXTER  
In the text 
Fig. 10 Quadrants problem evolved at a resolution of 400 × 400, with a CFL of 0.95, using the FORCE^{1} and GFORCE^{1} schemes with the van Leer slopelimiter. We plot the log of the density ρ, using 30 equally spaced contours from ln0.1 = −2.303 to ln3.927 = 1.368. 

Open with DEXTER  
In the text 
Fig. 11 Grid used for the frontfacing step problem. The grid is made up of three separate grids: one for the incoming channel, one for the rest of the channel, and one for the step itself. The last has stretched grid cells that become narrower towards the edge of the step. For clarity we only show every other grid line. 

Open with DEXTER  
In the text 
Fig. 12 Solution at t = 10 for a shock hitting a frontfacing step. The fluid has adiabatic index Γ = 1.4, velocity v = 0.8 and is travelling at Mach 3. The scheme used was GFORCE^{1} with the van Leer slopelimiter, and with CFL 0.95. We plot 30 contours of the log of the density. 

Open with DEXTER  
In the text 
Fig. 13 Solution for a reflective cylinder initially immersed in a constant velocity flow of Mach 3, velocity 0.8, adiabatic index Γ = 1.4, and density 1.4. The main channel has length 10 and width 5, composed of 100 × 200 cells. The annulus has inner radius 0.5 and outer radius 1.25, with 168 cells in the angular dimension, and 32 cells in the radial direction, where an exponential stretching function has been applied in the radial direction to restrict the radial size of cells nearer the inner boundary. For clarity, we only show every other grid line. In the top half of the plot, contours of ρ at time t = 10 are shown. The scheme used was GFORCE^{1}, with the van Leer slopelimiter, at a CFL of 0.95. There are 30 contours, equally spaced between 0 and 10. 

Open with DEXTER  
In the text 
Fig. 14 Simulation of flow past two solid cylinders. The first plot shows the grid and in the second we plot 30 contours of ρ from 0 to 10 for two solid cylinders initially immersed in a fluid of constant state, Mach 3 and velocity 0.8, with density 1.4 and adiabatic index 1.4. The solution is shown at time 10. The GFORCE^{1} scheme was used, with the van Leer slopelimiter, and a CFL of 0.95. 

Open with DEXTER  
In the text 