<figure><img src="image_1.png" /></figure>

HAL
open science

# XTROEM-FV: a new code for computational astrophysics based on very high order finite-volume methods - II. Relativistic hydro- and magnetohydrodynamics

Jonatan Núñez-de La Rosa, Claus-Dieter Munz

► To cite this version:

Jonatan Núñez-de La Rosa, Claus-Dieter Munz. XTROEM-FV: a new code for computational astrophysics based on very high order finite-volume methods - II. Relativistic hydro- and magnetohydrodynamics. 2016. [hal-01275178](https://hal.science/hal-01275178)

HAL Id: hal-01275178

https://hal.science/hal-01275178

Preprint submitted on 17 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# XTROEM-FV: a new code for computational astrophysics based on very high order finite-volume methods – II. Relativistic hydro- and magnetohydrodynamics

Jonatan Núñez-de la Rosa*, and Claus-Dieter Munz†

Institut für Aerodynamik und Gasdynamik, Universität Stuttgart. Pfaffenwaldring 21, D-70569 Stuttgart, Germany

Tuesday 16th February, 2016

## ABSTRACT

In this work we discuss the extension of the XTROEM-FV code to relativistic hydrodynamics and magnetohydrodynamics. XTROEM-FV is a simulation package for computational astrophysics based on very high order finite-volume methods on Cartesian coordinates. Arbitrary spatial high order of accuracy is achieved with a WENO reconstruction operator, and the time evolution is carried out with a strong-stability preserving Runge-Kutta scheme. In XTROEM-FV has been implemented a cheap, robust, and accurate shock capturing strategy for handling complex shock waves problems, typical in an astrophysical environment. The divergence constraint of the magnetic field is tackled with the generalized Lagrange multiplier divergence cleaning approach. Numerical computations of smooth flows for the relativistic hydrodynamics and magnetohydrodynamics equations are performed and confirm the high order accuracy of the main reconstruction algorithm for such kind of flows. XTROEM-FV has been subject to a comprehensive numerical benchmark, especially for complex flows configurations within an astrophysical context. Computations of problems with shocks with very high order reconstruction operators up to seventh order are reported. For instance, one-dimensional shock tubes problems for relativistic hydrodynamics and magnetohydrodynamics, as well as two-dimensional flows like the relativistic double Mach reflection problem, the interaction of a shock wave with a bubble, the relativistic Orszag-Tang vortex, the cylindrical blast wave problem, the rotor problem, the Kelvin-Helmholtz instability, and an astrophysical slab jet. XTROEM-FV represents a new attempt to simulate astrophysical flow phenomena with very high order numerical methods.

**Key words:** MHD – hydrodynamics – relativity – shock waves – methods: numerical.

## 1 Introduction

High energy astrophysical phenomena involve, in many cases, relativistic flows. Typical examples are superluminal motion of relativistic jets in extragalactic radio sources, accretion flows around massive compact objects, pulsar winds and gamma-ray bursts. When the magnetic fields are not considered in a first approximation, such flows are very well described with the equations of the relativistic hydrodynamics (SRHD). The special relativistic magnetohydrodynamics (SRMHD) equations provide a more accurate description of the underlying dynamics of many astrophysical plasma moving with speeds close to the speed of light. SRMHD plays a very important role when describing astrophysical jets emerging from super-massive black holes (Begelman et al. 1984). In fact, astronomical observations suggest that astrophysical jets emerging from compact objects involve significant magnetic fields and flows travelling

with relativistic speeds. General relativistic effects can be neglected when no strong gravitational fields are involved. The equations of the SRMHD describe the behaviour of relativistic, conducting fluids subject to electromagnetic fields. The main assumption is that the mean free-path of the electrons is much smaller than the characteristic length scale of the problem. This leads to a high-collisional frequency of the electrons. For numerically solving the SRHD and SRMHD equations, the so-called high-resolution shock-capturing (HRSC) schemes have provided the necessary tools in developing stable and robust relativistic fluid dynamical codes (for an excellent review, see Martí & Müller (2003, 2015)). Also see the references Schneider et al. (1993); Dolezal & Wong (1995); Martí et al. (1996); Aloy et al. (1999a,b); Mignone & Bodo (2005); Choi & Ryu (2005); Ryu et al. (2006); Dönmez & Kayali (2006); Dönmez (2006). The common feature of these schemes is that they are at most third-order accurate in space when they use the piecewise-parabolic method as reconstruction operator (Woodward & Colella 1984; Martí & Müller 1996). The SRHD has received a lot of attention over the last years from the high-order methods community. Some effort

<footnote>* E-mail: nunez@iag.uni-stuttgart.de</footnote>

<footnote>† E-mail: munz@iag.uni-stuttgart.de</footnote>
## 2 J. Núñez-de la Rosa, and C.-D. Munz

has been put into solving the equations with methods based on the WENO reconstruction procedure, as in, for example, the finite-volume code by Tchekhovskoy et al. (2007), the finite-difference code by Radice & Rezzolla (2012), or the discontinuous Galerkin by Dumbser et al. (2008) and by Radice & Rezzolla (2011). The numerical solution of the SRMHD equations received a boost with the development of finite-volume based codes by Komissarov (1999), and Balsara (2001). The authors used second order TVD schemes for solving the SRMHD equations written in conservation form. Further work for solving these equations has been done for several authors, which enormously improved the quality of the employed Riemann solvers (Mignone & Bodo 2006; Mignone & McKinney 2007; Mignone et al. 2009; Beckwith & Stone 2011), implemented high-order reconstruction operators in order to enhance the accuracy of the schemes (del Zanna et al. 2003; Anderson et al. 2006; del Zanna et al. 2007; Mignone et al. 2007; Dumbser et al. 2008), or developed robust adaptive mesh refinement (AMR) algorithms for SRMHD (Mignone et al. 2012; Keppens et al. 2012).

High-order numerical methods for conservation laws are numerical techniques with accuracy higher than second-order (Wang et al. 2013), and are used for solving hyperbolic systems of equations. Several schemes are available in computational fluid dynamics for structured and unstructured meshes. The conservative finite-difference method (Shu & Osher 1988, 1989), and the finite-volume method (Godunov 1959; van Leer 1979; Woodward & Colella 1984; Shu 2009) can achieve arbitrary spatial high-order accuracy provided the reconstruction operator is a high-order operator. This high-order accuracy can be only achieved in those regions with smooth flows. Among the available high order reconstruction operators, the WENO schemes are the most known and the most used in the computational fluid dynamics community. The WENO schemes are based on the ENO schemes from Harten et al. (1987), and were initially introduced by Liu et al. (1994). They have been extensively used in the last years for solving a wide spectrum of conservation laws (Shu 2009). They have been constructed for structured meshes (Liu et al. 1994; Balsara & Shu 2000) and also unstructured meshes (Dumbser & Käser 2007; Dumbser et al. 2007, 2013). XTROEM-FV makes use of the high-order WENO operators on structured meshes in order to achieve very high-order accuracy in the finite-volume scheme for solving the SRHD and SRMHD equations.

One of the main features of hyperbolic conservation laws is the presence of solution discontinuities. For high-order schemes, the appearance of discontinuities represents a major difficulty in order to achieve a higher-order accuracy. For that reason high-order methods are not widely used in computational fluid dynamics and computational astrophysics. The spurious oscillations originated around discontinuities produce unphysical states. Shock capturing strategies are necessary for handling such flow features if a high-order scheme is used. The slope/flux limiters used in the context of finite-difference and finite-volume schemes (Harten 1983; Harten et al. 1983) are mainly used for second-order methods. In this work, the higher-order methods use an adaptive strategy based on shock detection and order reduction of the reconstruction operator. This approach seems to be more robust and efficient, because it does not require the reconstruction of the characteristic variables, and the states in the region around the discontinuities are reconstructed with a very robust operator. Numerical computations of flow problems in SRHD and SRMHD confirm the reliability of the strategy (Mignone et al. 2007; Tchekhovskoy et al. 2007; Beckwith & Stone 2011; Radice & Rezzolla 2011, 2012).

The finite-volume scheme requires an additional layer for controlling the solenoidal constraint $\nabla \cdot \mathbf{B} = 0$ because this is not satisfied

from the numerical point of view. Actually, the spatial reconstruction and the time discretization generate errors in $\nabla \cdot \mathbf{B}$. These errors grow with time, and lead to unphysical states (Brackbill & Barnes 1980; Powell 1994). XTROEM-FV controls the solenoidal constraint with the generalized Lagrange multiplier (GLM) approach (Dedner et al. 2002; Mignone et al. 2010). This divergence cleaning scheme removes errors in $\nabla \cdot \mathbf{B}$ in a very cheap manner, in comparison with the Hodge projection (Chorin 1967; Brackbill & Barnes 1980), because this scheme has to solve a Poisson equation at each time step in order to maintain the divergence constraint up to machine accuracy. The GLM approach introduces a scalar field $\psi$ which couples the solenoidal constraint with the evolution equation of the magnetic field. The implementation of this approach is straightforward, and at the same time, conservation of all physical variables is maintained. We avoid the use of the constrained transport scheme (Evans & Hawley 1988) due to the complexity of the treatment of the staggered fields.

This paper is organized as follows: in Section 2 we present the SRHD and SRMHD equations. The eigenvalues of the Jacobian of the physical fluxes and the mapping from conservative to primitive variables are briefly discussed for both equation systems. In Section 3 an overview of the high order finite-volume method along with the WENO schemes for discretizing conservation laws is provided. In Section 4 are shown and discussed numerical computations for the SRHD and SRMHD equations in one- and two-dimensional spaces. The standard benchmark consists mainly of one-dimensional Riemann problems, and two-dimensional flows with shocks. Finally, in Section 5 the conclusions of this work are presented.

## 2 Governing equations

XTROEM-FV solves systems of conservation laws of special interest in astrophysics. In Núñez-de la Rosa & Munz (2016) we presented the XTROEM-FV code for the classical magnetohydrodynamics (MHD) equations. In this work we consider the extension to the SRHD and SRMHD equations. These equations describe the flow of fluid moving where the effects of special relativity cannot be neglected (Lichnerowicz 1967; Anile 1989; Martí & Müller 2003; Rezzolla & Zanotti 2013).

### 2.1 SRHD equations

The SRHD equations written as a system of conservation laws are given by (Font et al. 1994; Mignone et al. 2005; Ryu et al. 2006)

$$ \frac{\partial D}{\partial t} + \nabla \cdot (D\mathbf{v}) = 0, \qquad (1a) $$

$$ \frac{\partial S}{\partial t} + \nabla \cdot (S \otimes \mathbf{v} + p\mathbf{I}) = 0, \qquad (1b) $$

$$ \frac{\partial E}{\partial t} + \nabla \cdot \mathbf{S} = 0. \qquad (1c) $$

The five conserved quantities $D$, $S_x$, $S_y$, $S_z$ and $E$ are the mass density, the three components of the momentum density, and the total energy density, respectively. They are all measured in the laboratory frame, and are related to quantities in the local rest frame of the fluid (the primitive variables $\rho$, $v_x$, $v_y$, $v_z$, and $p$) through

$$ D = \rho\Gamma, \qquad (2a) $$

$$ S = \rho h \Gamma^2 \mathbf{v}, \qquad (2b) $$

$$ E = \rho h \Gamma^2 - p. \qquad (2c) $$
The three-velocity $v$ of the fluid is related with the four-velocity $u$ through

$$u = \Gamma(1, v), \qquad (3)$$

where $\Gamma$ is the Lorentz factor defined by

$$\Gamma = \frac{1}{\sqrt{1 - v^2}}, \qquad (4)$$

$\rho$ is the rest-mass density, $p$ is the isotropic gas pressure and $h$ is the specific enthalpy. The velocity $v = (v_x, v_y, v_z)$ of a fluid element is defined as the velocity of the rest frame of this fluid element with respect to the laboratory frame. The fluid velocity is a function of $(t, x, y, z)$, as are the thermodynamic quantities $\rho, p$.

The system (1) of partial differential equations is closed with an equation of state (EOS) $h = h(p, \rho)$ or $p = p(\rho, \varepsilon)$. The system (1) is hyperbolic for causal equations of state (Anile 1989; Rezzolla & Zanotti 2013), i.e., for those where the local sound speed satisfies $c_s < 1$, where $c_s$ is defined by

$$hc_s^2 = -\rho \frac{\partial h}{\partial \rho} \left( \rho \frac{\partial h}{\partial p} - 1 \right)^{-1}. \qquad (5)$$

In this work, we employ the most commonly used EOS, the ideal gas EOS, which is given by

$$p = (\gamma - 1)(\varepsilon - \rho), \quad \text{or} \quad h = 1 + \frac{\gamma}{\gamma - 1} \frac{p}{\rho}. \qquad (6)$$

Here $\gamma = c_p/c_v$ is the ratio of specific heats and $\varepsilon$ is the sum of the internal and rest-mass energy densities in the local frame and is related to the specific enthalpy as

$$h = \frac{\varepsilon + p}{\rho}. \qquad (7)$$

In this way, the speed of sound is given by

$$c_s^2 = \frac{\gamma p}{\rho h}. \qquad (8)$$

### 2.1.1 Eigenvalues of the SRHD system

The Jacobian matrices $A_i$ for the SRHD equations are defined by

$$A_i = \frac{\partial F_i(u)}{\partial u}, \quad (i = x, y, z), \qquad (9)$$

where $u$ is the state vector of conservative variables, and $F = [f, g, h]$ the tensor of physical fluxes. Assuming an ideal gas EOS, the eigenvalues of the matrix $A_x$ are given by

$$\lambda_1 = \frac{1}{1 - v^2 c_s^2} \left\{ \begin{aligned} & v_x (1 - c_s^2) \\ & - c_s \sqrt{(1 - v^2) [1 - v_x v_x - (v^2 - v_x v_x) c_s^2]} \end{aligned} \right\} \qquad (10)$$

$$\lambda_2 = v_x, \quad \lambda_3 = v_x, \quad \lambda_4 = v_x \qquad (11)$$

$$\lambda_5 = \frac{1}{1 - v^2 c_s^2} \left\{ \begin{aligned} & v_x (1 - c_s^2) \\ & + c_s \sqrt{(1 - v^2) [1 - v_x v_x - (v^2 - v_x v_x) c_s^2]} \end{aligned} \right\}. \qquad (12)$$

The cases for the matrices $A_y$ and $A_z$ easily follows from symmetry. The eigenvalues of $A_x$ represent the five characteristic speeds associated with two sound wave modes ($\lambda_1$ and $\lambda_5$) and three entropy modes ($\lambda_2, \lambda_3, \lambda_4$). The eigenvalues satisfy the following relation: $\lambda_1 < \lambda_2 = \lambda_3 = \lambda_4 < \lambda_5$.

### 2.1.2 Conservative to primitive variables conversion

Compared to the classical hydrodynamics equations, in the SRHD system, the relation between the primitive and conservative variables is not so simple. In this work, two strategies for expressing the primitive variables in terms of the conservative ones were followed: the first one, by solving iteratively an equation for the pressure and then determining the other variables from this (Mignone et al. 2005). The second one, which is more robust, consists in solving a quartic equation on the velocity $v$ (as it is outlined in Schneider et al. (1993); Duncan & Hughes (1994), and Ryu et al. (2006)). Though this last one is more computationally expensive, it provides more accurate results and does not have the problem of the iterative method with the initial guess of the pressure. In our algorithm, we combine both strategies in such way, that when the first one fails to iteratively find a root, we proceed to use the second algorithm. Only in very rare cases, for instance, when strong shocks are involved, the second algorithm is employed. This approach seems to be the most robust even for multidimensional flow simulations.

Solving numerically with Newton-Raphson. Here we follow the approach described in Mignone et al. (2005). The primitive variables $(\rho, v_x, v_y, v_z, p)$ are computed from the conservative quantities by solving with the one-dimensional Newton-Raphson algorithm

$$\frac{\partial f^n(p)}{\partial p} \delta p^{n+1} = f^n(p), \qquad (13)$$

the implicit equation

$$f(p) = D h \Gamma - E - p = 0, \qquad (14)$$

with $f'(p)$ given by

$$\frac{d f(p)}{d p} = \frac{\gamma}{\gamma - 1} \Gamma^2 - \frac{S^2 \Gamma^3}{(E + p)^3} \left( D + \frac{\gamma}{\gamma - 1} 2 p \Gamma \right) - 1, \qquad (15)$$

where the Lorentz factor is defined by

$$\Gamma_* = \frac{1}{\sqrt{1 - v_* \cdot v_*}}, \qquad (16)$$

and the components of the velocity are

$$v_* = \frac{S}{E + p}. \qquad (17)$$

Solving analytically a quartic polynomial on $|v|$. Here we follow the approach initially presented in Schneider et al. (1993) and further used in Choi & Ryu (2005), and Ryu et al. (2006). The quartic polynomial in $v$ we solve is given by

$$f(v) = \left[ \gamma v (E - S v) - S (1 - v^2) \right]^2 - (1 - v^2) v^2 (\gamma - 1)^2 D^2 = 0. \qquad (18)$$

Schneider et al. (1993) showed that the physically meaningful solution for $v$ is between the lower limit, $v_1$, and the upper limit, $v_2$,

$$v_1 = \frac{\gamma E - \sqrt{(\gamma E)^2 - 4 (\gamma - 1) S^2}}{2 (\gamma - 1) S}, \quad v_2 = \frac{S}{E}, \qquad (19)$$

and that the solution is unique. Once $v$ is known, the quantities $\rho, v_i$, and $\varepsilon$ can be straightforwardly calculated from the following relations

$$\rho = \frac{D}{\Gamma}, \quad v_x = \frac{S_x}{S} v, \quad v_y = \frac{S_y}{S} v, \quad v_z = \frac{S_z}{S} v, \qquad (20a)$$

$$p = (\gamma - 1) [(E - S_x v_x - S_y v_y - S_z v_z) - \rho]. \qquad (20b)$$
## 2.2 SRMHD equations

Considering the Minkowski spacetime with Cartesian coordinates $(t, x, y, z)$, the equations of the SRMHD form a system of conservation laws, which can be written as follows (Komissarov 1999; Balsara 2001)

$$ \frac{\partial D}{\partial t} + \nabla \cdot (D\mathbf{v}) = 0, \qquad (21a) $$

$$ \frac{\partial S}{\partial t} + \nabla \cdot (S \otimes \mathbf{v} + \mathbb{P}) = \mathbf{0}, \qquad (21b) $$

$$ \frac{\partial E}{\partial t} + \nabla \cdot \mathbf{S} = 0, \qquad (21c) $$

$$ \frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{B} \otimes \mathbf{v} - \mathbf{v} \otimes \mathbf{B}) = \mathbf{0}, \qquad (21d) $$

where the tensor $\mathbb{P}$ is defined as

$$ \mathbb{P} = \left( p + \frac{|\mathbf{B}|^2}{2\Gamma^2} + \frac{(\mathbf{v} \cdot \mathbf{B})^2}{2} \right) \mathbb{I} - \left( \frac{\mathbf{B}}{\Gamma^2} + (\mathbf{v} \cdot \mathbf{B})\mathbf{v} \right) \otimes \mathbf{B}. \qquad (22) $$

The quantity in the first bracket is the total pressure (which has contributions from the thermal and magnetic pressure). An EOS is used to close the system. We make use of the ideal gas EOS with adiabatic exponent $\gamma$

$$ p = (\gamma - 1)(\varepsilon - \rho) = (\gamma - 1) \left( E - \frac{1}{2} \rho |\mathbf{v}|^2 - \frac{1}{2} |\mathbf{B}|^2 \right). \qquad (23) $$

The system (21) must satisfy an additional constraint: the solenoidal property of the magnetic field,

$$ \nabla \cdot \mathbf{B} = 0, \qquad (24) $$

which is the same as in classical MHD. The eight conserved quantities $D, S_X, S_Y, S_Z, E, B_X, B_Y$, and $B_Z$ are the mass density, the three components of the momentum density, the total energy density, and the three components of the magnetic field, respectively. They are all measured in the laboratory frame, and are related to the quantities in the local rest frame of the fluid, the so-called primitive variables $\rho, v_X, v_Y, v_Z, p, B_X, B_Y, B_Z$, through

$$ D = \rho\Gamma, \qquad (25a) $$

$$ S = (\rho h \Gamma^2 + |\mathbf{B}|^2) \mathbf{v} - (\mathbf{v} \cdot \mathbf{B}) \mathbf{B}, \qquad (25b) $$

$$ E = \rho h \Gamma^2 - p + \frac{|\mathbf{B}|^2}{2} + \frac{|\mathbf{v}|^2 |\mathbf{B}|^2}{2} - \frac{(\mathbf{v} \cdot \mathbf{B})^2}{2}, \qquad (25c) $$

$$ \mathbf{B} = \mathbf{B}. \qquad (25d) $$

The rest-mass density is the quantity $\rho, p$ is the thermal pressure of the gas, and $h$ is the specific enthalpy. The system (21) of partial differential equations is closed with an EOS $h = h(\rho, \rho)$ or $p = p(\rho, \varepsilon)$. Like the SRHD equations, the system (21) is hyperbolic for a causal EOS (Anile 1989; Rezzolla & Zanotti 2013). In this work, we employ the ideal gas EOS, which is given by (6).

Others quantities of interest are the relativistic total (gas + magnetic) enthalpy

$$ h_{\text{tot}} = (\rho h + b^2), \qquad (26) $$

and the magnetic field in the fluid frame

$$ b^{\alpha} = \Gamma \left\{ \mathbf{v} \cdot \mathbf{B}, \frac{B^{i}}{\Gamma^2} + v^{i} (\mathbf{v} \cdot \mathbf{B}) \right\}. \qquad (27) $$

### 2.2.1 Eigenvalues of the SRMHD system

The characteristic structure of the SRMHD system was first studied by Anile & Pennisi (1987); more details can be found in Anile

(1989). These authors derived the eigenvalues and eigenvectors of the associated Jacobian $\mathbf{A}_X(\mathbf{u}) = \partial f / \partial \mathbf{u}$, with $f$ the physical flux in $x$-direction, by using the covariant notation.

XTROEM-FV does not require the left and right eigenvectors of the Jacobian matrices. Only the fastest and slowest characteristic speeds are needed. The eigenvectors have been reported in Anile & Pennisi (1987); Anile (1989); Antón et al. (2006), and Antón et al. (2010). Following del Zanna et al. (2003), the eigenvalues are: one entropy wave

$$ \lambda_0 = v_x, \qquad (28) $$

two Alfvén waves

$$ \lambda_A^{\pm} = \frac{u^x \pm \tilde{b}^x}{u^0 \pm \tilde{b}^0}, \qquad (29) $$

and four magneto-sonic waves (two fast and two slow waves). These are found by solving the non-linear quartic equation

$$ (1 - \varepsilon^2)(u^0 \lambda - u^x)^4 + (1 - \lambda^2)[(c_s^2(\tilde{b}^0 \lambda - \tilde{b}^x)^2 - \varepsilon^2(u^0 \lambda - u^x)^2] = 0, \qquad (30) $$

where $c_s^2 = \gamma p / \rho h$ is the sound speed squared, $\tilde{b}^\alpha = b^\alpha / \sqrt{h_{\text{tot}}}$ ($|\tilde{b}|^2 = \tilde{b}_\alpha \tilde{b}^\alpha = |b|^2 / h_{\text{tot}}$), and $\varepsilon^2 = c_s^2 + |\tilde{b}|^2 - c_s^2 |\tilde{b}|^2$.

### 2.2.2 Conservative to primitive variables conversion

Although the numerical scheme evolves the conservative variables $\mathbf{u}$, for the calculation of several quantities necessary for the execution of the code (for example, the physical fluxes, the fastest waves, the characteristic variables, etc), the primitive variables have to be calculated at least once every time step. The equations (25) express the relation between the conservative variables and the primitive variables. As it is clearly noticed from these equations, such relation is highly non-linear, and a procedure to find the primitive variables has to be developed. Following the algorithms discussed in Noble et al. (2006), and Mignone & McKinney (2007), we employ in our codes the strategies presented in those works. In this work we solve iteratively an equation for the pressure and then we determine the other variables from this (Mignone & McKinney 2007).

Solving numerically with Newton-Raphson. Here we will outline the approach presented by Mignone & McKinney (2007), considering only an ideal gas EOS and omitting the more general case when non-relativistic speeds are present (more details can be found in Mignone & McKinney (2007)). From equations (25), we take in consideration those for the energy and for the square modulus of the momentum. Defining the scalars $\hat{x} = Dh\Gamma$ and $\hat{y} = S \cdot \mathbf{B}$, we get

$$ E = \hat{x} - p + \frac{1 + |\mathbf{v}|^2}{2} |\mathbf{B}|^2 - \frac{\hat{y}^2}{2\hat{x}^2}, \qquad (31) $$

and

$$ |\mathbf{S}|^2 = (\hat{x} + |\mathbf{B}|^2)^2 |\mathbf{v}|^2 - \frac{\hat{y}^2}{\hat{x}^2} (2\hat{x} + |\mathbf{B}|^2). \qquad (32) $$

From the last equation, one can express $|\mathbf{v}|^2$ as a function of the unknown $\hat{x}$

$$ |\mathbf{v}|^2 = \frac{\hat{y}^2(2\hat{x} + |\mathbf{B}|^2) + |\mathbf{S}|^2 \hat{x}^2}{(\hat{x} + |\mathbf{B}|^2)^2 \hat{x}^2}. \qquad (33) $$

The next step is to insert this result into the expression for the energy (31)

$$ E = \hat{x} - p + \frac{|\mathbf{B}|^2}{2} + \frac{|\mathbf{B}|^2 |\mathbf{S}|^2 - \hat{y}^2}{2(|\mathbf{B}|^2 + \hat{x})^2}. \qquad (34) $$
This equation will be solved using a Newton-Raphson iterative scheme, where the $(k+1)$-th approximation to $\hat{x}$ is calculated as

$$ \hat{x}^{(k+1)} = \hat{x}^{(k)} - \left. \frac{f(\hat{x})}{df(\hat{x})/d\hat{x}} \right|_{\hat{x}=\hat{x}^{(k)}}. \quad (35) $$

The functional $f(\hat{x})$ is derived from the equation (34)

$$ f(\hat{x}) = \hat{x} - E - p + \frac{|\mathbf{B}|^2}{2} + \frac{|\mathbf{B}|^2|\mathbf{S}|^2 - \hat{y}^2}{2(|\mathbf{B}|^2 + \hat{x})^2}. \quad (36) $$

Computing the derivative $df(\hat{x})/d\hat{x} \equiv dE/d\hat{x}$, yields

$$ \frac{dE}{d\hat{x}} = 1 - \left. \frac{dp}{d\hat{x}} \right|_{\hat{x}} - \frac{|\mathbf{B}|^2|\mathbf{S}|^2 - \hat{y}^2}{(|\mathbf{B}|^2 + \hat{x})^3}. \quad (37) $$

Let us assume that $p = p(\rho, \hat{z})$, with

$$ \hat{z} \equiv \rho(h-1) = \frac{\rho\Gamma}{\Gamma}(h-1) = \frac{Dh\Gamma - D\Gamma}{\Gamma^2} = \frac{\hat{x} - D\Gamma}{\Gamma^2}. \quad (38) $$

Applying the chain rule, we are able to calculate the derivative $dp/d\hat{x}$,

$$ \frac{dp}{d\hat{x}} = \left. \frac{\partial p}{\partial \hat{z}} \right|_{\rho} \frac{d\hat{z}}{d\hat{x}} + \left. \frac{\partial p}{\partial \rho} \right|_{\hat{z}} \frac{d\rho}{d\hat{x}}. \quad (39) $$

Calculating the derivatives $d\hat{z}/d\hat{x}$, $d\rho/d\hat{x}$, yields

$$ \frac{d\hat{z}}{d\hat{x}} = \frac{1}{\Gamma^2} - \frac{\Gamma(D + 2\hat{z}\Gamma)}{2} \frac{d|\hat{v}|^2}{d\hat{x}}, \quad (40) $$

and

$$ \frac{d\rho}{d\hat{x}} = -\frac{D\Gamma}{2} \frac{d|\hat{v}|^2}{d\hat{x}}, \quad (41) $$

where

$$ \frac{d|\hat{v}|^2}{d\hat{x}} = -\frac{2}{\hat{x}^3} \frac{\hat{v}^2[3\hat{x}(\hat{x} + |\mathbf{B}|^2) + |\mathbf{B}|^4] + |\mathbf{S}|^2\hat{x}^3}{(\hat{x} + |\mathbf{B}|^2)^3}. \quad (42) $$

For the ideal gas EOS, we have

$$ p(\rho, \hat{z}) = \frac{\gamma - 1}{\gamma} \hat{z}. \quad (43) $$

From the expression $p(\hat{z}, \rho)$ we get easily

$$ \frac{\partial p}{\partial \hat{z}} = \frac{\gamma - 1}{\gamma}, \quad \frac{\partial p}{\partial \rho} = 0. \quad (44) $$

Once $\hat{x}$ has been determined to some accuracy with the Newton-Raphson procedure, we complete the inversion process by computing the velocities from an inversion of equation (21b) to obtain

$$ v_i = \frac{1}{\hat{x} + |\mathbf{B}|^2} \left( S_i + \frac{\hat{y}}{\hat{x}} B_i \right), \quad (45) $$

From the velocity, we can calculate the Lorentz factor $\Gamma$, and then we get directly $\hat{z}$ from equation (38),

$$ \hat{z} = \frac{\hat{x} - D\Gamma}{\Gamma^2}. \quad (46) $$

The pressure is obtained from the relation

$$ p(\rho, \hat{z}) = \frac{\gamma - 1}{\gamma} \hat{z}. \quad (47) $$

The rest-mass density is obtained from

$$ \rho = \frac{D}{\Gamma}, \quad (48) $$

and the magnetic field are the same as in the primitive variables.

## 2.3 Divergence cleaning with the GLM method

In the SRMHD equations, the solenoidal constraint of the magnetic field (24) is satisfied analytically but not numerically. Therefore a strategy for maintaining this constraint from the numerical point of view has to be used. The divergence cleaning proposed by Munz et al. (1999), and Dedner et al. (2002) is used in this work, specifically the mixed hyperbolic/parabolic approach. Within this scheme, the solenoidal constraint is coupled with the induction equation (21d) through the potential $\psi$. The induction equation and the solenoidal constraint are replaced by

$$ \frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{B} \otimes \mathbf{v} - \mathbf{v} \otimes \mathbf{B} + \psi \mathbf{I}) = \mathbf{0}, \quad (49a) $$

$$ \frac{\partial \psi}{\partial t} + \nabla \cdot (c_h^2 \mathbf{B}) = - \frac{c_h^2}{c_p^2} \psi. \quad (49b) $$

Equations (21a), (21b), (21c), (49a), and (49b), constitute the so-called GLM-SRMHD system.

# 3 Numerical methods

In this section we are going to briefly review the main ingredients of the numerical scheme: the finite-volume method, the WENO reconstruction schemes, the divergence cleaning, and the shock capturing approach followed in this work.

## 3.1 Finite-volume methods

We are interested in discretizing systems of conservation laws. We start by considering the system

$$ \frac{\partial \mathbf{u}}{\partial t} + \frac{\partial f(\mathbf{u})}{\partial x} + \frac{\partial g(\mathbf{u})}{\partial y} + \frac{\partial h(\mathbf{u})}{\partial z} = \mathbf{0}. \quad (50) $$

Finite-volume methods are numerical methods based on the weak solution of conservation laws in integral form. The semi-discrete scheme is obtained after integration of the Eq. (50) over the cell $\Omega_{ijk} = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}] \times [y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}] \times [z_{k-\frac{1}{2}}, z_{k+\frac{1}{2}}]$,

$$ \frac{d\mathbf{u}_{ijk}}{dt} = - \frac{\hat{f}_{i+\frac{1}{2},jk} - \hat{f}_{i-\frac{1}{2},jk}}{\Delta x} - \frac{\hat{g}_{i,j+\frac{1}{2},k} - \hat{g}_{i,j-\frac{1}{2},k}}{\Delta y} - \frac{\hat{h}_{i,j,k+\frac{1}{2}} - \hat{h}_{i,j,k-\frac{1}{2}}}{\Delta z}. \quad (51) $$

The quantity $\mathbf{u}_{ijk}$ is the spatial average of $\mathbf{u}$ in the cell $\Omega_{ijk}$ at the time $t$

$$ \mathbf{u}_{ijk} = \frac{1}{|\Omega_{ijk}|} \int_{\Omega_{ijk}} \mathbf{u}(x, y, z) \, dz \, dy \, dx, \quad (52) $$

with the volume of the cell $\Omega_{ijk}$ given by $|\Omega_{ijk}| = \Delta x \Delta y \Delta z$. The spatial averages of the physical fluxes $\hat{f}_{i \pm \frac{1}{2},jk}$, $\hat{g}_{i,j \pm \frac{1}{2},k}$, and $\hat{h}_{i,j,k \pm \frac{1}{2}}$, over the cell faces $x_{i \pm \frac{1}{2}}$, $y_{j \pm \frac{1}{2}}$, and $z_{k \pm \frac{1}{2}}$, respectively are defined by

$$ \begin{align} \hat{f}_{i \pm \frac{1}{2},jk} &= \frac{1}{|\sigma_{jk}|} \int_{\sigma_{jk}} f(\mathbf{u}(x_{i \pm \frac{1}{2}}, y, z)) \, dz \, dy, \nonumber \\ \hat{g}_{i,j \pm \frac{1}{2},k} &= \frac{1}{|\sigma_{ik}|} \int_{\sigma_{ik}} g(\mathbf{u}(x, y_{j \pm \frac{1}{2}}, z)) \, dz \, dx, \tag{53} \\ \hat{h}_{i,j,k \pm \frac{1}{2}} &= \frac{1}{|\sigma_{ij}|} \int_{\sigma_{ij}} h(\mathbf{u}(x, y, z_{k \pm \frac{1}{2}})) \, dy \, dx, \nonumber \end{align} $$
6 J. Núñez-de la Rosa, and C.-D. Munz

The surface elements in a Cartesian mesh are defined by $\sigma_{ij} = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}] \times [y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}]$, $\sigma_{jk} = [y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}] \times [z_{k-\frac{1}{2}}, z_{k+\frac{1}{2}}]$, and $\sigma_{ik} = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}] \times [z_{k-\frac{1}{2}}, z_{k+\frac{1}{2}}]$. The area of the faces are then $|\sigma_{ij}| = \Delta x \Delta y$, $|\sigma_{jk}| = \Delta y \Delta z$, and $|\sigma_{ik}| = \Delta x \Delta z$. Observe that the numerical fluxes defined as

$$
\begin{align}
\hat{f}_{i+\frac{1}{2},jk} &= \hat{f}(\mathbf{u}_{i-p,jk}, \dots, \mathbf{u}_{i+q,jk}), \nonumber \\
\hat{g}_{i,j+\frac{1}{2},k} &= \hat{g}(\mathbf{u}_{i,j-p,k}, \dots, \mathbf{u}_{i,j+q,k}), \tag{54} \\
\hat{h}_{ij,k+\frac{1}{2}} &= \hat{h}(\mathbf{u}_{ij,k-p}, \dots, \mathbf{u}_{ij,k+q}), \nonumber
\end{align}
$$

constitute an approximation of the physical fluxes and they are con-
sistent with them in the sense that $\hat{f}(\mathbf{u}, \dots, \mathbf{u}) = f(\mathbf{u})$, $\hat{g}(\mathbf{u}, \dots, \mathbf{u}) = g(\mathbf{u})$, and $\hat{h}(\mathbf{u}, \dots, \mathbf{u}) = h(\mathbf{u})$ (see LeVeque (1992); Bressan (2000); LeVeque (2002), and Toro (2009)).

The semi-discrete scheme (51) is written in conservative form. The Lax–Wendroff theorem guarantees that if the scheme (51) is convergent, then it converges to the weak solution (Lax & Wendroff 1960). Observe that the equation (51) is an exact relation as well as the averaged quantities (52) and (53).

In the finite-volume scheme the fluxes (53) have to be approxi-
mated. These integrals are discretized with a Gaussian quadrature.
In order to achieve a very high order of accuracy in space, a high-
order Gaussian quadrature is required. For Cartesian meshes, the
integrals are then given by

$$
\begin{align}
\hat{f}_{i \pm \frac{1}{2}, jk} &= \frac{1}{\Delta y} \frac{1}{\Delta z} \sum_{\alpha=1}^{N_{GP}} \sum_{\beta=1}^{N_{GP}} f(\mathbf{u}(x_{i \pm \frac{1}{2}}, y_\alpha, z_\beta)) \omega_\alpha \omega_\beta, \nonumber \\
\hat{g}_{i, j \pm \frac{1}{2}, k} &= \frac{1}{\Delta x} \frac{1}{\Delta z} \sum_{\alpha=1}^{N_{GP}} \sum_{\beta=1}^{N_{GP}} g(\mathbf{u}(x_\alpha, y_{j \pm \frac{1}{2}}, z_\beta)) \omega_\alpha \omega_\beta, \tag{55} \\
\hat{h}_{ij, k \pm \frac{1}{2}} &= \frac{1}{\Delta x} \frac{1}{\Delta y} \sum_{\alpha=1}^{N_{GP}} \sum_{\beta=1}^{N_{GP}} h(\mathbf{u}(x_\alpha, y_\beta, z_{k \pm \frac{1}{2}})) \omega_\alpha \omega_\beta. \nonumber
\end{align}
$$

The point-wise values of **u** at the Gaussian integration points at the faces are obtained through a high-order reconstruction. The reconstruction operator takes as input the cell averages **u**<sub>ijk</sub> at the cell barycenter, and produces the reconstructed values at the cell faces. Observe that at a given face, two sets of reconstructed values are available: those obtained through the use of the cell Ω<sub>ijk</sub> as the main cell in the reconstruction procedure and those obtained with the neighbouring cell. These are the so-called left and right values at the face interface: **u**<sub>L</sub> and **u**<sub>R</sub>. The fluxes are computed by replacing **f**. **g**, and **h** by a monotone flux (or an exact or approximated Riemann solver (Toro 2009)). The Riemann solvers used in XTROEM-FV for the SRHD and SRMHD equations are the very well known Rusanov numerical flux (Rusanov 1961), and the HLL Riemann solver (Toro 2009). The Rusanov numerical flux is given by

$$
\hat{f}(\mathbf{u}_L, \mathbf{u}_R) = \frac{1}{2} (f(\mathbf{u}_L) + f(\mathbf{u}_L) - |\lambda_{\max}|(\mathbf{u}_L - \mathbf{u}_R)), \quad (56)
$$

where $\lambda_{\max}$ is the largest local wave speed, which guarantees the
stability of the scheme.

**3.2 Discretization of the GLM method**

From the equations (49), we can observe that the quantities $B_x$ and $\psi$ are decoupled from the GLM-SRMHD system, yielding the linear system

$$
\frac{\partial}{\partial t} \begin{pmatrix} B_x \\ \psi \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ c_h^2 & 0 \end{pmatrix} \frac{\partial}{\partial x} \begin{pmatrix} B_x \\ \psi \end{pmatrix} = \begin{pmatrix} 0 \\ -\frac{c_h^2}{c_p^2} \psi \end{pmatrix}. \quad (57)
$$

The solution for the local Riemann problem with left and right states
$(B_{x,l}, \psi_l)^T$, $(B_{x,r}, \psi_r)^T$, respectively, is given by

$$
B_{x,m} = \frac{1}{2}(B_{x,r} + B_{x,l}) - \frac{1}{2c_h}(\psi_r - \psi_l), \quad (58a)
$$

$$
\psi_m = \frac{1}{2}(\psi_r + \psi_l) - \frac{c_h}{2}(B_{x,r} - B_{x,l}). \quad (58\text{b})
$$

Following Dedner et al. (2002), we use the exact solution (58) of the linear Riemann problem (57) as input for the Riemann solver used in the solution of the other conserved quantities. In order to deal with the source term in equation (49b), we follow the idea presented in Dedner et al. (2002): we first solve the homogeneous GLM–SRMHD system in a so-called hyperbolic step, and then we consider the source term in the source step. The scalar field $\psi$ is then

$$
\psi^{(\Delta t)} = \psi^{(0)} \exp \left( -\alpha_p \frac{c_h}{\Delta h / \Delta t} \right), \quad \text{with } \alpha_p = \Delta h \frac{c_h}{c_p^2}, \qquad (59)
$$

where $\psi^{(0)}$ has been computed in the hyperbolic step, and $\Delta h = \min(\Delta x, \Delta y, \Delta z)$ is the minimum mesh size. This approach is very simple to implement and is unconditionally stable (Dedner et al. 2002).

3.3 WENO schemes

In this section we are going to outline the principles of the WENO schemes. For a more detailed description, see Shu (2009). The WENO schemes are numerical methods that produce adaptive high-order reconstruction polynomials from cell average data in order to compute pointwise values at appropriate points, as for example at Gaussian integration points in the context of finite-volume methods. Provided a main stencil around a main cell, this is subdivided in to several small stencils, each containing the main cell. For every sub-stencil a reconstruction polynomial is built from the cell average of the cells in the substencl. These polynomials are then weighted in order to construct a higher-order polynomial. The final reconstruc-tion polynomial has the following properties: in smooth regions, the polynomial is the highest order polynomial obtained from the cell averages; in regions with discontinuities, it remains essentially non-oscillatory. The weights used in the reconstruction are got from each polynomial in each substencl, and they provide information about the smoothness of every polynomial.

In the case of structured meshes, the dimension-by-dimension WENO reconstruction procedure is the most efficient among all WENO approaches (Casper & Atkins 1993; Shu 2009), although the methodology introduced in Buchmüller & Helzel (2014) in the context of finite-volume methods makes the reconstruction procedure as simple as for finite difference methods.

In the three-dimensional case, the cell averages of the function
u = u(x, y, z) in the cell Ωijk is given by

$$
u_{ijk} = \frac{1}{|\Omega_{ijk}|} \int_{\Omega_{ijk}} u(x, y, z) \, dz \, dy \, dx. \quad (60)
$$

The WENO scheme reconstructs point values of *u* at the Gaussian integration points ($x_{i \pm \frac{1}{2}}$, $y_{j+\alpha}$, $z_{k+\beta}$), ($x_{i+\alpha}$, $y_{j \pm \frac{1}{2}}$, $z_{k+\beta}$), ($x_{i+\alpha}$, $y_{j+\beta}$, $z_{i \pm \frac{1}{2}}$) from the cell averages $u_{ijk}$. We employ a two-points Gaussian quadrature rule (Titarev & Toro 2004, 2005)

$$
\int_{-1}^{+1} f(x) \, dx \approx \sum_{i=1}^{2} w_i f(x_i) = f\left(-\frac{1}{\sqrt{3}}\right) + f\left(+\frac{1}{\sqrt{3}}\right), \quad (61)
$$

for calculating the surface integrals (55). These integration
points are given by $(x_{i\pm\frac{1}{2}}, y_{j+\alpha})$ and $(x_{i+\alpha}, y_{j\pm\frac{1}{2}})$ in 2D and $(x_{i\pm\frac{1}{2}}, y_{j+\alpha}, z_{k+\beta})$, $(x_{i+\alpha}, y_{j\pm\frac{1}{2}}, z_{k+\beta})$, $(x_{i+\alpha}, y_{j+\beta}, z_{i\pm\frac{1}{2}})$ in 3D, with $\alpha, \beta = \pm 1/2\sqrt{3}$. Titarev & Toro (2004, 2005) claim that the high-order accuracy of the scheme is provided by the high-order reconstruction of the function values at the Gaussian integration points. The two-points Gaussian quadrature employed by XTROEM-FV leads to a formal fourth order of accuracy. But several numerical experiments performed by Titarev & Toro (2004), and Titarev & Toro (2005) show that the accuracy is not affected with this choice.

For completeness reasons, we describe the dimension-by-dimension algorithm for the reconstruction at the Gaussian integration points $(x_{i\pm\frac{1}{2}}, y_{j+\alpha}, z_{k+\beta})$, also the points at faces $x_{i\pm\frac{1}{2}}$ of the cell $\Omega_{ijk}$. A similar procedure is used for reconstructing the point values of the function $u(x, y, z)$ at points $(x_{i+\alpha}, y_{j\pm\frac{1}{2}}, z_{k+\beta})$, and $(x_{i+\alpha}, y_{j+\beta}, z_{i\pm\frac{1}{2}})$, which correspond to the faces $y_{j\pm\frac{1}{2}}$ and $z_{k\pm\frac{1}{2}}$. The dimension-by-dimension algorithm employs three sweeps in the three-dimensional case and two sweeps in the two-dimensional case. Following Titarev & Toro (2004), the stencils used to reconstruct $u_{i\pm\frac{1}{2}, j+\alpha, k+\beta}$ with a WENO scheme of $(2N+1)$th order (where the polynomials used in every subStencil are of degree $N$), we require the stencil to be formed by the cells $\Omega_{i_x i_y i_z}$, with the indices $i_x, i_y$, and $i_z$ satisfying

$$ 
\begin{aligned}
i - N \le i_x \le i + N, \\
j - N \le i_y \le j + N, \\
k - N \le i_z \le k + N.
\end{aligned} 
\qquad (62) $$

The sweeps for the WENO scheme are the followings:

(1) First Sweep: From the cell averages $u_{ijk}$, a one-dimensional reconstruction in the x-direction is carried out for all values of the indices $i_y, i_z$ from the stencil. Two-dimensional averages at faces $x_{i\pm\frac{1}{2}}$ are obtained from this procedure

$$ \bar{u}_{i_y i_z} \Big|_{x_{i\pm\frac{1}{2}}} = \frac{1}{\Delta y} \frac{1}{\Delta z} \int_{y_{i_y-\frac{1}{2}}}^{y_{i_y+\frac{1}{2}}} \int_{z_{i_z-\frac{1}{2}}}^{z_{i_z+\frac{1}{2}}} u(x_{i\pm\frac{1}{2}}, y, z) \, dz \, dy. $$

(2) Second Sweep: From the obtained two-dimensional averages $\bar{u}_{i_y i_z}$, a one-dimensional reconstruction in the y-direction is carried out for all values of the index $i_z$ from the stencil. One-dimensional averages at lines $y_{j\pm\frac{1}{2\sqrt{3}}}$ on the faces $x_{i\pm\frac{1}{2}}$ are obtained from this procedure

$$ \bar{u}_{i_z} \Big|_{\substack{x_{i\pm\frac{1}{2}} \\ y_{j\pm\frac{1}{2\sqrt{3}}}}} = \frac{1}{\Delta z} \int_{z_{i_z-\frac{1}{2}}}^{z_{i_z+\frac{1}{2}}} u(x_{i\pm\frac{1}{2}}, y_{j\pm\frac{1}{2\sqrt{3}}}, z) \, dz. $$

(3) Third Sweep: In the last sweep, a one-dimensional reconstruction in the z-direction is carried out from all line averages $\bar{u}$ of every line $(x = x_{i\pm\frac{1}{2}}, y = y_{j\pm\frac{1}{2\sqrt{3}}})$. In this step, all point-wise values $u(x_{i\pm\frac{1}{2}}, y_{j\pm\frac{1}{2\sqrt{3}}}, z_{k\pm\frac{1}{2\sqrt{3}}})$ are reconstructed.

## 3.4 Shock capturing for high order finite-volume methods

In Núñez-de la Rosa & Munz (2016) we presented a shock capturing strategy for high order finite-volume schemes. Shock capturing algorithms are required in high-order numerical methods when shocks waves, or strong rarefactions are present in the simulation. If this kind of features are not properly tackled from the numerical point of view, it may typically generate unphysical states, like negative densities or pressures, or even superluminal speeds in relativistic

hydrodynamics. These unphysical states are found especially in the reconstruction step, and the conservative-to-primitive variables algorithm. Here we show again the shock capturing steps included in the base algorithm of the numerical scheme.

(1) The first step consists in checking quantities with physical restrictions and applying a correction mechanism. For example, we check the positivity of the density and the pressure, and check that there are no superluminal speeds

(2) In the second step we proceed with the detection of regions with strong shocks or discontinuities. For that purpose we use a shock indicator, and flag those troubled cells and their direct neighbours.

(3) In the third and last step, a special treatment in the troubled regions is applied. We employ robust Riemann solvers, and apply a robust second/third order reconstruction scheme.

We have to point out that a similar approach has been used by Mignone et al. (2007); Tchekhovskoy et al. (2007); Beckwith & Stone (2011); Radice & Rezzolla (2011), and Radice & Rezzolla (2012). This approach is usually refereed as the *fallback approach* because the reconstruction order is reduced in order to properly handle the shocks and discontinuities. After flagging the cells containing the shocks or discontinuities, we reduce the order of the reconstruction operator, and typically we employ the WENO3 algorithm and the MUSCL scheme. The shock indicator we use is the Jameson indicator in the pressure (Jameson et al. 1981),

$$ \eta_i = \frac{|p_{i+1} - 2p_i + p_{i-1}|}{|p_{i+1}| + 2|p_i| + |p_{i-1}|}. \qquad (63) $$

If the Jameson indicator is larger than a certain threshold value (in XTROEM-FV this value is $\eta = 5.0 \times 10^{-3}$ for all computations), then the cell is flagged, and the robust scheme is employed.

## 3.5 Time discretization

The time discretization is explicitly carried out with the family of strong-stability preserving Runge-Kutta (SSPRK) methods (Shu & Osher 1988; Shu 1988; Gottlieb & Shu 1998). The semi-discrete scheme (51) comprises a system of ordinary differential equations that is solved using the *method of lines* with a fourth-order SSPRK scheme. The SSPRK methods are used in hyperbolic problems where shocks and discontinuities arise in a natural way. As a time integrator, the SSPRK methods do not introduce additional oscillations because they preserve strong stability under the TVD norm (Ruuth & Spiteri 2002; Spiteri & Ruuth 2002, 2003; Hundsdorfer et al. 2003; Gottlieb 2005; Hesthaven & Warburton 2008; Gottlieb et al. 2009; Shu 2009).

If we write the system of ordinary differential equations (51) as

$$ \frac{du}{dt} = L(u, t), \quad u(t_0) = u_0, \quad t \in [t_0, t_f], \qquad (64) $$

where $-L(u, t)$ is the spatial discretization operator, then the fourth-order SSPRK with five stages as derived in Spiteri & Ruuth (2002)
is given by

$$
\begin{align}
\boldsymbol{u}^{(0)} &= \boldsymbol{u}^n \nonumber \\
\boldsymbol{u}^{(1)} &= \boldsymbol{u}^{(0)} + 0.39175222700392 \Delta t \boldsymbol{L}(\boldsymbol{u}^{(0)}, t^n) \nonumber \\
\boldsymbol{u}^{(2)} &= 0.44437049406734 \boldsymbol{u}^{(0)} + 0.55562950593266 \boldsymbol{u}^{(1)} \nonumber \\
&\quad + 0.36841059262959 \Delta t \boldsymbol{L}(\boldsymbol{u}^{(1)}, t^n + 0.39175222700392 \Delta t) \nonumber \\
\boldsymbol{u}^{(3)} &= 0.62010185138540 \boldsymbol{u}^{(0)} + 0.37989814861460 \boldsymbol{u}^{(2)} \nonumber \\
&\quad + 0.25189177424738 \Delta t \boldsymbol{L}(\boldsymbol{u}^{(2)}, t^n + 0.58607968896780 \Delta t) \nonumber \\
\boldsymbol{u}^{(4)} &= 0.17807995410773 \boldsymbol{u}^{(0)} + 0.82192004589227 \boldsymbol{u}^{(3)} \nonumber \\
&\quad + 0.54497475021237 \Delta t \boldsymbol{L}(\boldsymbol{u}^{(3)}, t^n + 0.47454236302687 \Delta t) \nonumber \\
\boldsymbol{u}^{(5)} &= 0.00683325884039 \boldsymbol{u}^{(0)} + 0.51723167208978 \boldsymbol{u}^{(2)} \nonumber \\
&\quad + 0.12759831133288 \boldsymbol{u}^{(3)} + 0.34833675773694 \boldsymbol{u}^{(4)} \nonumber \\
&\quad + 0.08460416338212 \Delta t \boldsymbol{L}(\boldsymbol{u}^{(3)}, t^n + 0.47454236302687 \Delta t) \nonumber \\
&\quad + 0.22600748319395 \Delta t \boldsymbol{L}(\boldsymbol{u}^{(4)}, t^n + 0.93501063100924 \Delta t) \nonumber \\
\boldsymbol{u}^{n+1} &= \boldsymbol{u}^{(5)}. \tag{65}
\end{align}
$$

This scheme has an SSP coefficient $c = 1.50818004975927$. The SSP coefficient satisfies $\Delta t \le c\Delta t_E$, where $\Delta t_E$ is the time step restricted by the Courant-Friedrichs-Levy (CFL) condition for the first-order forward Euler method. The CFL condition is a necessary condition for stability of the numerical scheme, and it is applied to constraint the time step (Courant et al. 1928; Blazek 2005). In the one-dimensional case this is given by

$$ \Delta t = C_{CFL} \min_k \left( \frac{\Delta x}{|\lambda_{k,x}|} \right), \quad (66) $$

with $\lambda_k$ the maximum eigenvalue of the physical $x$-flux Jacobian over all computational cells. The number $C_{CFL}$ depends on the spatial discretization and on the time-stepping scheme used for solving the hyperbolic equation, and besides it satisfies $C_{CFL} \le 1$. Only approximated values of the CFL condition are available for non-linear equations in the multidimensional case. Following Titarev & Toro (2005), the CFL condition is given by

$$ \Delta t = \frac{C_{CFL}}{2} \min_k \left( \frac{\Delta x}{|\lambda_{k,x}|}, \frac{\Delta y}{|\lambda_{k,y}|} \right), \quad (67) $$

for a two-dimensional space, and

$$ \Delta t = \frac{C_{CFL}}{3} \min_k \left( \frac{\Delta x}{|\lambda_{k,x}|}, \frac{\Delta y}{|\lambda_{k,y}|}, \frac{\Delta z}{|\lambda_{k,z}|} \right), \quad (68) $$

for a three-dimensional space. Observe that the number $C_{CFL}$ corresponds to the one-dimensional case.

# 4 Numerical Computations

## 4.1 Relativistic hydrodynamics

### 4.1.1 Convergence test

With this test we want to check the order of accuracy of the XTROEM-FV code for solving the SRHD equations when the RKFV method with WENO reconstruction is used. For this purpose, we solve a two-dimensional smooth problem where a wave is propagating in the physical domain $[0, 2/\sqrt{2}] \times [0, 2]$ at an angle $\theta = 30^\circ$

<table><caption>**Table 1.** Convergence rates for the SRHD equations with initial condition given by the smooth flow problem. The base numerical methods is the RKFV scheme. In these calculations we have employed the Rusanov Riemann solver. Results for $L_2$ norm error of the rest-mass density are provided, with reconstruction operator WENO3, WENO5 and WENO7. Simulation time was set to $t = 1.0$.</caption><thead><tr><th>Method</th><th>Cells</th><th><i>L</i><sub>2</sub> error</th><th><i>L</i><sub>2</sub> order</th></tr></thead><tbody><tr><td rowspan="4">WENO3</td><td>40 × 20</td><td>2.895 × 10<sup>-04</sup></td><td rowspan="2">2.86</td></tr><tr><td>80 × 40</td><td>3.984 × 10<sup>-05</sup></td></tr><tr><td>160 × 80</td><td>5.278 × 10<sup>-06</sup></td><td>2.92</td></tr><tr><td>320 × 160</td><td>6.383 × 10<sup>-07</sup></td><td>3.05</td></tr><tr><td rowspan="4">WENO5</td><td>40 × 20</td><td>3.542 × 10<sup>-06</sup></td><td rowspan="2">4.98</td></tr><tr><td>80 × 40</td><td>1.123 × 10<sup>-07</sup></td></tr><tr><td>160 × 80</td><td>3.543 × 10<sup>-09</sup></td><td>4.99</td></tr><tr><td>320 × 160</td><td>9.698 × 10<sup>-11</sup></td><td>5.19</td></tr><tr><td rowspan="4">WENO7</td><td>40 × 20</td><td>9.274 × 10<sup>-08</sup></td><td rowspan="2">7.08</td></tr><tr><td>80 × 40</td><td>6.845 × 10<sup>-10</sup></td></tr><tr><td>160 × 80</td><td>5.345 × 10<sup>-12</sup></td><td>7.00</td></tr><tr><td>320 × 160</td><td>3.856 × 10<sup>-14</sup></td><td>7.12</td></tr></tbody></table>

relative to the horizontal axis. Following He & Tang (2012a), the initial profile is given by

$$ \rho = 1 + A \sin(\phi(0)), \quad (69) $$

$$ v_x = v_0, \quad (70) $$

$$ v_y = 0, \quad (71) $$

$$ p = 1, \quad (72) $$

with $A = 0.2$, and $v_0 = 0.2$. The function $\phi(t)$ is defined as follows

$$ \phi(t) = 2\pi \left( (x \cos \theta + y \sin \theta) - (v_x \cos \theta + v_y \sin \theta)t \right). \quad (73) $$

This problem has the exact solution

$$ \rho = 1 + A \sin(\phi(t)), \quad (74) $$

$$ v_x = v_0, \quad (75) $$

$$ v_y = 0, \quad (76) $$

$$ p = 1. \quad (77) $$

The discretization of the computational domain satisfy $N_y = 2N_x$, and periodic boundary conditions are set at all four faces of the domain. In the Table 1 are shown the convergence rates for the RKFV method. The numerical scheme converges to the theoretical order of convergence.

### 4.1.2 One-dimensional Riemann problems

In this section we consider some standard one-dimensional Riemann problems extensively used in relativistic hydrodynamics. The Riemann problem is an initial-value problem for a conservation law defined as

$$ \frac{\partial \boldsymbol{u}}{\partial t} + \frac{\partial f(\boldsymbol{u})}{\partial x} = \mathbf{0}, \quad (78) \\ \boldsymbol{u}(x, 0) = \begin{cases} \boldsymbol{u}_L(x) & \text{for } x < x_m, \\ \boldsymbol{u}_R(x) & \text{for } x > x_m, \end{cases} $$

where $x_m$ is the position of the initial discontinuity. For all examples here discussed, the one-dimensional domain is the closed interval $[0, 1]$. The membrane is localized in the point $x_m = 0.5$. The initial condition for the Riemann problems is characterized for the two constant states written in conservative variables $\boldsymbol{u}_L$ and $\boldsymbol{u}_R$ at left and right sides of the membrane, respectively. After the breakup of
the membrane, different shock patterns are originated. The problems below discussed have been used as benchmark in most of the relativistic codes found in the literature. We refer to Martí & Müller (2003, 2015) for a comprehensive review of relativistic hydrodynamics in special relativity, with a special chapter dedicated to the Riemann problems used here as benchmark, and additionally a detailed literature survey about the numerical methods employed for their solution.

The numerical calculations with high-order schemes are contrasted with the exact solution obtained by using the exact Riemann solver provided by Rezzolla & Zanotti (2001). In the following, we will write all initial states in primitive variables. The left state will be represented by $w_L = (\rho, v_x, v_y, p)_L$ and the right state by $w_R = (\rho, v_x, v_y, p)_R$, where $v_y$ can be interpreted as the tangential velocity in the problem. For all simulations it was necessary to reduce the reconstruction/interpolation order at the vicinity of a strong shock, typically three cells around this point and the adopted reconstruction was the very robust and stable third-order WENO scheme.

**Riemann problem 1 (RP-1).** This test is called also mildly relativistic blast wave. The initial condition is determined by the following left and right states

$$ w_L = (10, 0, 0, 40/3), \\ w_R = (1, 0, 0, 0). \qquad (79) $$

We have employed the ideal gas EOS with adiabatic index $\gamma = 5/3$. The CFL number we have used for the computation is $C_{CFL} = 0.95$. The interval is discretized with 500 computational cells. Transmissive boundary conditions are set at both sides of the interval. The final simulation time is set to $t = 0.4$. In the right state the pressure was set to be $2/3 \times 10^{-6}$ for numerical reasons. This test problem has been considered for many authors in the past (Schneider et al. 1993; Martí & Müller 1996; Donat et al. 1998), but it was in Martí & Müller (1994) that an exact solution was first presented. The waves present in this problem are a left-going transonic rarefaction wave, a contact discontinuity and a right-going shock wave. The fluid behind the shock is moving with a mildly relativistic speed $v = 0.72c$ to the right. In this dense shell behind the shock, the fluid is compressed. Because of that, the fluid is heated to values of the internal energy much larger than the rest-mass energy. This means that the fluid is thermodynamically relativistic, but mildly relativistic dynamically. Plots of the rest-mass density, pressure, and the $x$-component of the velocity are depicted in the Fig. 1 (first row, from top to bottom). Computations with WENO3, WENO5, and WENO7 reconstruction operators are reported. The shock is captured in around 3–4 cells with the WENO3 reconstruction operator.

**Riemann problem 2 (RP-2).** The initial state for the highly relativistic blast wave is given by

$$ w_L = (1, 0, 0, 10^3), \\ w_R = (1, 0, 0, 10^{-2}). \qquad (80) $$

We have employed the ideal gas EOS with adiabatic index $\gamma = 5/3$, and CFL number $C_{CFL} = 0.95$. The interval is discretized with 500 computational cells. Transmissive boundary conditions have been utilized at domain edges. The final simulation time is set to $t = 0.4$. The decay of the initial discontinuity gives rise to a dense intermediate state located between a right-going shock wave and a rarefaction wave propagating to the left. The shock wave and the contact discontinuity are very close to each other. The very thin shell is a very challenging feature for any numerical method solving this problem. Actually, in these simulations, and because of the

smearing at the contact discontinuity, this shell is not well resolved. Plots of the rest-mass density, pressure, and the $x$-component of the velocity are depicted in the Fig. 1 (second row, from top to bottom). Computations with WENO3, WENO5, and WENO7 reconstruction operators are reported.

**Riemann problem 3 (RP-3).** In the transverse blast wave problem, the initial conditions are similar to the RP-2, but with the difference that in this problem a non-zero transverse velocity in the right state is introduced. That is, the initial condition is determined by the following left and right states

$$ w_L = (1, 0, 0, 10^3), \\ w_R = (1, 0, 0.99, 10^{-2}). \qquad (81) $$

We have employed the ideal gas EOS with adiabatic index $\gamma = 5/3$. We have used $C_{CFL} = 0.95$. The interval is discretized with 500 computational cells. The final simulation time is set to $t = 0.4$. Unlike Newtonian hydrodynamics, where the transverse momentum is not coupled with the longitudinal one and it thus is simply advected, the momentum equations in relativistic hydrodynamics are coupled each other through the Lorentz factor. As consequence of that, a non-zero transverse velocity imprints in the solution of the Riemann problem new physical effects that are not present in classical hydrodynamics (Pons et al. 2000; Rezzolla et al. 2003). Plots of the rest-mass density, pressure, and the $x$-component of the velocity are depicted in the Fig. 1 (third row, from top to bottom). The solution was computed by using a finite-volume method with WENO3, WENO5, and WENO7 reconstruction operators. No oscillations are observed in the computed numerical approximations with the WENO3 scheme. A light overshooting is obtained with the WENO5 and WENO7 schemes in the dense intermediate shell.

**Riemann problem 4 (RP-4).** The planar shock reflection consists in an ideal cold fluid colliding a wall, a shock wave propagates backwards, leaving the gas behind at rest. The reflecting wall is located in $x = 1$. At $x = 0$ transmissive boundary conditions are set. The initial state is given by

$$ w = (1, 0.99999, 0, 0.01) \qquad (82) $$

We have employed the ideal gas EOS with adiabatic index $\gamma = 4/3$. We have used $C_{CFL} = 0.95$. The interval is discretized with 500 computational cells. The final simulation time is set to $t = 1.5$. The exact solution of this Riemann problem was first obtained by Blandford & McKee (1976). Plots of the rest-mass density, pressure, and the $x$-component of the velocity computed with the finite volume WENO scheme are depicted in the Fig. 1 (fourth row, from top to bottom). We can see that close to the wall ($x = 1$), the numerical solution of the rest-mass density shows an undershooting, which is due to the wall heating phenomenon (Noh 1987). The shock has been very well captured using around 3–4 cells with the WENO3 reconstruction, but some small oscillations are clearly visible around this discontinuity, although only for the WENO5 and WENO7 schemes.

### 4.1.3 Two-dimensional Riemann problems

In this section we consider two-dimensional Riemann problems, which initial conditions were taken from del Zanna & Bucciantini (2002); Tchekhovskoy et al. (2007), and He & Tang (2012a). The computational domain is the box $[-1, 1] \times [-1, 1]$, divided in four quadrants. The membranes are localized along the lines $x = 0$ and $y = 0$. These Riemann problems simulate practically all essential features found in a two-dimensional flow, for example, shock reflections, vortex, shock interactions, etc. In the two-dimensional
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Solution to four one-dimensional Riemann problems in SRHD. The left and right states, and the description of the flow patterns are given in the text. Plots of the rest-mass density, pressure, and the x-component of the velocity are depicted. The solution was computed with a high order finite-volume WENO scheme, with WENO3, WENO5, and WENO7 reconstruction operators.</figcaption></figure>

Riemann problems, the initial states are defined in the quadrants

Q1, Q2, Q3, and Q4, which are the sets

$$ 
\begin{aligned}
Q1 &:= \{(x, y) \in [-1, 1]^2 | x \ge 0, y \ge 0\}, \\
Q2 &:= \{(x, y) \in [-1, 1]^2 | x < 0, y \ge 0\}, \\
Q3 &:= \{(x, y) \in [-1, 1]^2 | x < 0, y < 0\}, \\
Q4 &:= \{(x, y) \in [-1, 1]^2 | x \ge 0, y < 0\}.
\end{aligned} 
\tag{83} $$

MNRAS 000, 1–24 (2016)
Riemann problem 5 (RP-5). The box [-1, 1] x [-1, 1] is initially filled with a gas with adiabatic index $\gamma = 5/3$. This domain is subdivided into four quadrants, where the states in primitive variables are defined by

$$ 
\begin{aligned}
\mathbf{w}_{Q1} &= (0.035145216124503, 0, 0, 0.162931056509027), \\
\mathbf{w}_{Q2} &= (0.1, 0.7, 0, 1), \\
\mathbf{w}_{Q3} &= (0.5, 0, 0, 1), \\
\mathbf{w}_{Q4} &= (0.1, 0, 0.7, 1).
\end{aligned} 
\quad (84) $$

Transmissive boundary conditions are used in all faces of the computational domain. The final simulation time is $t = 0.8$. Contour plots of the rest-mass density and pressure are depicted in the Fig. 2. After the breakup of the membranes, two contact discontinuities appear on the left and bottom of the domain, and in the first quadrant we can identify two curved front shocks. In the third quadrant it is visible a like-jet structure moving in south-west direction (del Zanna & Bucciantini 2002). The computations are in excellent agreement with those reported in the literature.

Riemann problem 6 (RP-6). The domain is filled with a gas, which initial condition is determined by the following states in primitive variables in the quadrants Q1, Q2, Q3, and Q4

$$ 
\begin{aligned}
\mathbf{w}_{Q1} &= (0.5, 0.5, -0.5, 5), \\
\mathbf{w}_{Q2} &= (1, 0.5, 0.5, 5), \\
\mathbf{w}_{Q3} &= (3, -0.5, 0.5, 5), \\
\mathbf{w}_{Q4} &= (1.5, -0.5, -0.5, 5).
\end{aligned} 
\quad (85) $$

As in the RP-5, the adiabatic index of the gas is $\gamma = 5/3$ and the final simulation time is $t = 0.8$. Contour plots of the rest-mass density and the Lorentz factor are depicted in the Fig. 3. Here we can see the interaction of four vortex sheets, forming a spiral with very low rest-mass density in the centre of the domain (the typical cavitation in gas dynamics). Observe the very well behaviour of the hybrid WENO7/WENO3 scheme in the very low density region, that is the scheme does not produce a negative density or pressure.

Riemann problem 7 (RP-7). In our last two-dimensional Riemann problem, we set the initial condition as the states in primitive variables in the quadrants Q1, Q2, Q3, and Q4 given by

$$ 
\begin{aligned}
\mathbf{w}_{Q1} &= (1, 0, 0, 1), \\
\mathbf{w}_{Q2} &= (0.5771, -0.3529, 0, 0.4), \\
\mathbf{w}_{Q3} &= (1, -0.3529, -0.3529, 1), \\
\mathbf{w}_{Q4} &= (0.5771, 0, -0.3529, 0.4).
\end{aligned} 
\quad (86) $$

The adiabatic index of the gas is $\gamma = 5/3$, and the final simulation time is $t = 0.8$. Contour plots of the rest-mass density and the Mach number are depicted in the Fig. 4. This problem is about the interaction of planar rarefaction waves. The initial discontinuities evolve into rarefaction waves, interacting each other. This interaction generates two symmetric shock waves in the zone where the rarefaction waves have interplayed.

### 4.1.4 Double Mach reflection problem

The double Mach reflection problem was introduced by Woodward & Colella (1984) for Newtonian hydrodynamics. This problem was extended to the equations of relativistic hydrodynamics by Zhang & MacFadyen (2006). The problem consists basically in a shock wave travelling horizontally and collides with an inclined wedge. It is formulated equivalently by sending a shock wave diagonally into a reflecting wall. The computational domain is the box [0, 4] x [0, 1].

We represent the wedge by a reflecting boundary starting at $x = 1/6$ along the x-axis. The strong shock moving initially from the left boundary toward the right makes a 60° angle with the x-axis. Due to the original setup of the problem can not be used in the relativistic case, Zhang & MacFadyen (2006) proposed a new initial state that allows the formation of waves patterns present in the Newtonian version. In fact, they showed that by choosing an ultra-relativistic shock, no Mach reflection is generated. The adiabatic index is set to $\gamma = 7/5$, and the initial state is given by

$$ \mathbf{w}(x, y)|_{t=0} = \begin{cases} \mathbf{w}_L, & \text{for } y > h(x, 0), \\ \mathbf{w}_R, & \text{for } y < h(x, 0), \end{cases} \quad (87) $$

where the left and right states of the shock wave, and the exact position of the shock at the time $t$ are, respectively (see also He & Tang (2012a)),

$$ 
\begin{aligned}
\mathbf{w}_L &= (8.564, 0.4247 \sin 60^\circ, -0.4247 \cos 60^\circ, 0.3808), \\
\mathbf{w}_R &= (1.4, 0, 0, 0.0025),
\end{aligned} 
\quad (88) $$

$$ h(x, t) = \sqrt{3} (x - x_0) - 2v_s t, \quad x_0 = \frac{1}{6}, \quad v_s = 0.4984. $$

Regarding the boundary conditions, we set them as follows: a reflecting wall is placed at $x > x_0$, $y = 0$, as we mentioned before. At the top boundary and depending of the position of the shock wave, we set the boundary conditions to either the post-shock state $w_L$ or the pre-shock state $w_R$. At the right face of the domain the boundary condition is set to the exact pre-shock condition $w_R$. For $x < x_0$, the left and part of the bottom boundaries are set to the exact post-shock state $w_L$. A contour plot of the rest-mass density at simulation time $t = 4.0$ is depicted in the Fig. 5. The computational domain is made of 1600 x 400 cells. The WENO7 reconstruction is used in smooth parts of the flow and the WENO3 in regions with shocks. No carbuncle phenomenon is observed because we employed a Rusanov numerical flux (Woodward & Colella 1984).

### 4.1.5 Cloud-shock interaction

The simulation setup is given after He & Tang (2012a): we employ an EOS with adiabatic index $\gamma = 5/3$. The computational domain is the box [0, 2] x [0, 1] with transmissive boundary conditions at all faces. The centre of the cylindrical cloud is located in the point (1.4, 0.5) and the radius is $r = 0.15$. The initial state for the cloud is $w_c = (3.1538, 0, 0, 0.05)$. The left moving shock wave at time $t = 0$ is located at $x = 1.6$, with left and right states given by

$$ 
\begin{aligned}
\mathbf{w}_L &= (1, 0, 0, 0.05), \\
\mathbf{w}_R &= (1.86522508063118, -0.19678110737829, 0, 0.15).
\end{aligned} 
\quad (89) $$

The computational domain is decomposed into 800 x 400 cells. Contour plots of the rest-mass density and Mach number at simulation time $t = 3.0$ are depicted in the Fig. 6. The combination of WENO3 and WENO7 reconstruction operators makes possible to obtain such profile of the bubble after the interaction with the shock wave. Lower-order schemes require a higher mesh resolution in order to get a similar solution (He & Tang 2012a).

### 4.1.6 Relativistic slab jet propagation

It is widely known that there are three classes of highly collimated and supersonic jets from dense central objects with accretion disks, which depend on the central object, protostars, binary stars, or active galactic nuclei. AGN jets are the largest scale phenomena, and the velocity of the jet beam is highly relativistic, at least close to the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Two-dimensional Riemann problem RP-5. Contour plots of the rest-mass density (left) and pressure (right) at $t = 0.8$. The solution was obtained with a high order finite-volume method with WENO7 reconstruction. The domain $[-1, 1] \times [-1, 1]$ was decomposed into $600 \times 600$ cells.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Two-dimensional Riemann problem RP-6. Contour plots of the rest-mass density (left) and Lorentz factor (right) at $t = 0.8$. The solution was obtained with a high order finite-volume method with WENO7 reconstruction. The domain $[-1, 1] \times [-1, 1]$ was decomposed into $600 \times 600$ cells.</figcaption></figure>

central object (Peterson 1997). The jet, which originates near an accretion disk that surrounds an AGN, can propagate over a long distance, up to a few Mpc, while remaining well collimated. There are two shocks at the end of the jet. One is a bow shock (or a forward shock), which accelerates the ambient gas. The other is a terminal Mach shock (or a reverse shock) at which the beam ends. At the terminal Mach shock, non-thermal particles are accelerated and emit photons through synchrotron radiation and inverse Compton scattering. The gas that crosses the terminal Mach shock into a hot spot is hot and pressurized, and expands laterally, enveloping the beam with the shocked ambient gas, creating a so-called cocoon structure. At the contact discontinuity between the ambient gas and the jet in the cocoon, Kelvin-Helmholtz instabilities develop.

Cygnus A is a suitable object in which to see these features, because it is one of the closest radio galaxies, and the beam propagates

perpendicularly to the line of sight. From observations, its size is about 120 kpc, the beam velocity is ~0.4–1c, and the hot spot’s ram pressure advance speed is 0.03c.

Analytical studies and numerical simulations of the morphology and the dynamics of jets have been performed for the past thirty years. Blandford & Rees (1974) discussed the structure of jets with a theoretical relativistic beam model. The difficulty of numerical relativistic hydrodynamics has delayed the investigation of the relativistic effects on the morphology and the dynamics of jets. Only in the past 20 years stable codes, with or without external magnetic fields, have been developed for the ultra-relativistic regime (Duncan & Hughes 1994; Aloy et al. 1999a).

We solve the SRHD equations (1) with the two-dimensional version of our relativistic hydrodynamic code, based on a high-order RKFV method. We use a two-dimensional Cartesian computational
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Two-dimensional Riemann problem RP-7. Contour plots of the rest-mass density (left) and Mach number (right) at $t = 0.8$. A high-order WENO7 reconstruction was employed for calculating the solution. The domain $[-1, 1] \times [-1, 1]$ is made of $600 \times 600$ cells.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Relativistic double Mach reflection problem. Contour plot of the rest-mass density at $t = 4.0$. A high-order WENO7 reconstruction was used along with a WENO3 operator for the shock capturing. A Rusanov numerical flux was employed for computing the intercell fluxes. The domain was decomposed into $1600 \times 400$ cells. Only the region $[0, 3] \times [0, 1]$ is shown.</figcaption></figure>

region. The grid size is uniform, namely $\Delta x = \Delta y = \text{constant}$. We assume that the ambient gas is homogeneous initially. A relativistic beam flow ($v_b$), is injected at $x = 0$ in the direction of the positive x-axis through a circular nozzle defined by $x^2 + y^2 \le R_b^2$, where $R_b$ is the beam radius, and is in pressure equilibrium with the ambient medium. Outflow boundary conditions are imposed everywhere except at the plane $x = 0$, where injection is assumed through the nozzle and the rest of the plane has a reflecting boundary, that is, at $x = 0$ with $x^2 + y^2 > R_b^2$. The boundary condition at $x = 0$ is crucial for dynamics and the outer shape of the jets. The radius of the injected beam, $R_b$, is used as a scaling unit in this study.

The relativistic jet simulation can be fully specified by only setting the following parameters:

* the beam density, $\rho_b$;

* the beam pressure, $p_b$;

* the flow speed of the beam, $v_b$;

* the ratio of the rest mass densities of the beam fluid and the ambient medium, $\eta$;

* the ratio of the pressures of the beam fluid and the ambient medium, $\kappa$.

For the example considered in this work, the initial parameters are listed in the following. An ideal gas EOS with adiabatic index $\gamma = 5/3$ is used. The ratio of the pressures of the beam fluid and the ambient medium is set to $\kappa = 1$. The density and pressure of the beam are $\rho_b = 1.0$ and $p_b = 1.0 \times 10^{-3}$. The speed of the beam is given by $v_b = 0.999c$. The jet nozzle has a radius $r_b = 0.2$. The computational domain has a size $450r_b \times 150r_b$, and it is discretized with $2400 \times 800$ and $3600 \times 1200$ cells. The simulations were performed with a finite-volume scheme with WENO5 reconstruction operator.

From our simulations, we find that the overall morphology and dynamics of the jets are similar to those discussed in previous works in the literature. The Fig. 7 shows contour plots of the logarithm of the rest-mass density of the jets simulated with the RKFV-WENO5 scheme at two different mesh resolutions. There are some differences in the inner structure of the jets and in the form of the head, mainly because the deceleration phase of the jets, which is caused
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Relativistic cloud-shock interaction. Contour plots of the rest-mass density (top), and Mach number (bottom) at time $t = 3.0$. The computations were done with a finite-volume scheme with WENO7 reconstruction. The computational domain is the box $[0, 2] \times [0, 1]$, discretized by using a mesh of 800 × 400 cells.</figcaption></figure>

by the generation and separation of vortices at the head of the jet, is different for the resolutions used in the simulations. All beams remain collimated from the nozzle, where the beam is injected into the computational region, to the head of the jets. It is also important to note that the beam radius does not increase from the nozzle to the head of the jet. At the end of the beam, a strong ‘terminal Mach shock’ is present. One of the most active points is called a ‘hot spot’, into which shocked beam gas enters through the terminal Mach shock at the head of the jet. The pressure in the hot spot is very high, because of the energy dissipation at the terminal Mach shock, and is matched by ambient gas compressed at the bow shock. It is also observed a back flow that creates a shear flow, and the contact surface becomes unstable because of Kelvin–Helmholtz instabilities. The surface between the back flow and the shocked ambient gas flow also becomes unstable and causes the appearance of vortices in the larger cocoon.

## 4.2 Relativistic magnetohydrodynamics

### 4.2.1 Convergence test

The convergence test is carried out with the computation of the propagation of large amplitude circularly polarized Alfvén waves. This test was first proposed by del Zanna et al. (2007), and it consists of a periodic Alfvén wave which is an exact solution of the SRMHD equations, therefore it is used to check the accuracy of the numerical scheme in smooth flows regimes. As initial condition, the primitive variables are set as follows: $\rho = 1$, $p = 1$, $v_x = 0$, $v_y = -v_A B_y / B_0$, $v_z = -v_A B_z / B_0$, $B_x = 0$, $B_y = \eta B_0 \cos(2\pi x)$, and $B_z = \eta B_0 \sin(2\pi x)$. The two-dimensional computational domain is the box $[0, 1] \times [0, 1]$ (Dumbser et al. 2008), with periodic boundary conditions at all faces of the box. The adiabatic index of the EOS is set to $\gamma = 4/3$. We set $B_0 = 1$, and $\eta = 1$. In del Zanna

<table><caption>**Table 2.** Convergence rates for the SRMHD equations. The propagation of large amplitude circularly polarized Alfvén waves is computed. In these calculations we have employed the Rusanov Riemann solver. Results for $L_2$ norm error of the z-component of magnetic field are given, with reconstruction operators WENO3, WENO5 and WENO7. Simulation time is $t = 2.618033988$.</caption><thead><tr><th>Method</th><th>Cells</th><th><i>L</i><sub>2</sub> error</th><th><i>L</i><sub>2</sub> order</th></tr></thead><tbody><tr><td rowspan="5">WENO3</td><td>20 × 20</td><td>3.291 × 10<sup>-03</sup></td><td rowspan="5">2.93</td></tr><tr><td>40 × 40</td><td>4.303 × 10<sup>-04</sup></td></tr><tr><td>80 × 80</td><td>4.993 × 10<sup>-05</sup></td></tr><tr><td>160 × 160</td><td>5.430 × 10<sup>-06</sup></td></tr><tr><td>320 × 320</td><td>6.129 × 10<sup>-07</sup></td></tr><tr><td rowspan="5">WENO5</td><td>20 × 20</td><td>4.624 × 10<sup>-05</sup></td><td rowspan="5">5.02</td></tr><tr><td>40 × 40</td><td>1.429 × 10<sup>-06</sup></td></tr><tr><td>80 × 80</td><td>4.293 × 10<sup>-08</sup></td></tr><tr><td>160 × 160</td><td>1.228 × 10<sup>-09</sup></td></tr><tr><td>320 × 320</td><td>3.357 × 10<sup>-11</sup></td></tr><tr><td rowspan="5">WENO7</td><td>20 × 20</td><td>4.295 × 10<sup>-06</sup></td><td rowspan="5">7.05</td></tr><tr><td>40 × 40</td><td>3.236 × 10<sup>-08</sup></td></tr><tr><td>80 × 80</td><td>2.359 × 10<sup>-10</sup></td></tr><tr><td>160 × 160</td><td>1.654 × 10<sup>-12</sup></td></tr><tr><td>320 × 320</td><td>1.120 × 10<sup>-14</sup></td></tr></tbody></table>

et al. (2007) is shown that the speed of the Alfvén wave is given by the formula

$$v_A^2 = \frac{B_0^2}{\rho h + B_0^2(1 + \eta^2)} \left[ \frac{1}{2} \left( 1 + \sqrt{1 - \left( \frac{2\eta B_0^2}{\rho h + B_0^2(1 + \eta^2)} \right)^2} \right) \right]^{-1}$$

With the parameters of the initial condition, the speed of the Alfvén wave is $v_A = \sqrt{2/(7 + \sqrt{45})} = 0.38196601125$. The final simulation time is then $t = 1/v_A = 2.618033988$ (one period). In the Table 2 are shown the errors and convergence rates in the $L_2$ norm for the RKFV method. The variable used for this analysis is the z-component of the magnetic field. Observe that the numerical scheme converges to the theoretical order of convergence.

### 4.2.2 One-dimensional Riemann problems

The benchmark for the SRMHD equations involves the solution of one-dimensional Riemann problems. Two fluids with different conditions are separated with a membrane at point $x_m$. After the membrane is removed, the flow is characterized by the appearance of different shock waves and discontinuities. The solution is computed up to certain end time $t$. The considered tests are those discussed by Balsara (2001), and Giacomazzo & Rezzolla (2006). For all these tests, we have used $C_{CFL} = 0.95$, and we set transmissive boundary conditions. The one-dimensional domain is the closed interval $[0, 1]$ was decomposed into 500 computational cells, and the membrane separating the left and right initial states is localized in the point $x_m = 0.5$. The left state will be represented by $w_L = (\rho, v_x, v_y, v_z, p, B_x, B_y, B_z)_L$ and the right state by $w_R = (\rho, v_x, v_y, v_z, p, B_x, B_y, B_z)_R$. Both states are written in terms of the primitive variables. The exact solution is computed using the exact Riemann solver developed by Giacomazzo & Rezzolla (2006).

**Riemann problem 1 (RP-1).** The initial condition is determined by the following left and right states

$$\begin{align}
\boldsymbol{u}_L &= (1, 0, 0, 0, 1, 0.5, 1, 0), \tag{90} \\
\boldsymbol{u}_R &= (0.125, 0, 0, 0, 0.1, 0.5, -1, 0). \nonumber
\end{align}$$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Relativistic astrophysical jet simulation. The logarithm of rest-mass density at simulation time $t = 120$. The RKFV method with WENO5 is used in the simulation. The mesh is made of 1800 × 600 cells (top), and 2400 × 800 (bottom).</figcaption></figure>

This test involves the so-called compound structures, that is, structures that involve a shock and a rarefaction of the same wave family moving together. Their existence was first discussed in Brio & Wu (1988) for the MHD equations. This test involves a left-going slow compound wave. In the analytic solution the slow composed wave is absent. Other structures appearing after the breakup of the membrane are two fast rarefactions waves propagating in opposite directions, two slow shocks also propagating in opposite directions, and a contact discontinuity. For this problem we have assumed $\gamma = 2$. Plots of the rest-mass density, total pressure, and the y-component of the magnetic field are depicted in the Fig. 8 (first row, from top to bottom) at time $t = 0.4$. Computations with WENO3, WENO5, and WENO7 reconstruction operators are rather similar, but with a light undershooting present in the solution computed with the higher-order operators WENO5 and WENO7.

Riemann problem 2 (RP-2). The initial condition is determined by the following left and right states

$$ \mathbf{u}_L = (1, 0, 0, 0, 30, 5, 6, 6), \qquad (91) $$

$$ \mathbf{u}_R = (1, 0, 0, 0, 1, 5, 0.7, 0.7). $$

The solution to this Riemann problem features two left-moving fast and slow rarefactions waves, a contact discontinuity and two right-moving fast and slow shocks. We have employed the ideal gas EOS with adiabatic index $\gamma = 5/3$. Plots of the rest-mass density, total pressure, and the y-component of the magnetic field are depicted in the Fig. 8 (second row, from top to bottom) at time $t = 0.4$. The WENO3, WENO5, and WENO7 reconstruction operators were employed for these computations. No oscillations are present, but the dense intermediate shell is not perfectly resolved.

Riemann problem 3 (RP-3). The initial condition is determined by the following left and right states

$$ \begin{aligned} \mathbf{u}_L &= (1.08, 0.4, 0.3, 0.2, 0.95, 2, 0.3, 0.3), \\ \mathbf{u}_R &= (1, -0.45, -0.2, 0.2, 1, 2, -0.7, 0.5). \end{aligned} \qquad (92) $$

The solution involves a left-going fast shock, a left-going Alfvén discontinuity, a left-going slow rarefaction, a contact discontinuity, a right-going slow shock, a right-going Alfvén discontinuity and a right-going fast shock. We used an ideal gas EOS with adiabatic index $\gamma = 5/3$. Plots of the rest-mass density, total pressure, and y-component of the magnetic field are depicted in the Fig. 8 (third row, from top to bottom) at simulation time $t = 0.55$. Computations were done with different WENO schemes, but the WENO3 was used exclusively in regions with shocks/discontinuities. Observe the light overshoots obtained with the WENO5 and WENO7. Shocks are captured in around 3–4 cells.

Riemann problem 4 (RP-4). The initial condition is determined by the following left and right states

$$ \begin{aligned} \mathbf{u}_L &= (1, 0, 0.3, 0.4, 5, 1, 6, 2), \\ \mathbf{u}_R &= (0.9, 0, 0, 0, 5.3, 1, 5, 2). \end{aligned} \qquad (93) $$

The solution to this Riemann problem features a left-moving fast rarefaction, a left-moving Alfvén discontinuity, a left-moving slow shock, a contact discontinuity, a right-moving slow shock, a right-moving Alfvén discontinuity and a right-moving fast shock. The gas has an adiabatic index $\gamma = 5/3$. The final simulation time is set to $t = 1.5$. Plots of the rest-mass density, total pressure, and the y-component of the magnetic field are depicted in the Fig. 8 (fourth row, from top to bottom). The shock are very well resolved, but some overshooting is observed for the WENO5 and WENO7 schemes,
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Solution to four one-dimensional Riemann problems for the SRMHD equations. The left and right initial states, and the description of the flow patterns are given in the text. From top to bottom, RP-1, RP-2, RP-3, and RP-4. Plots of the rest-mass density, total pressure, and the y-component of the magnetic field are depicted. The solution was computed with a high order finite-volume WENO scheme, with WENO3, WENO5, and WENO7 reconstruction operators. The computational domain is the interval [0, 1], and it was decomposed into 500 cells.</figcaption></figure>
meanwhile for the WENO3 scheme undershooting is present in the region between the dense shell.

### 4.2.3 Magnetic field loop advection

This test has been extensively used for the MHD equations. With this problem one can observe whether the algorithm is able to preserve $\nabla \cdot \mathbf{B} = 0$. In this test a cylindrical current distribution is advected along some direction of the computational domain not aligned with the grid. The magnetic loop remains in magnetostatic balance. For that reason, after some periods, the profile should be preserve its shape. Multidimensional SRMHD codes that do not include algorithms for preserving the constraint $\nabla \cdot \mathbf{B} = 0$, or if they use numerical methods too much diffusive, the magnetic loop will smear over the time (Tóth 2000; Stone et al. 2008; Mignone & Tzeferacos 2010; Beckwith & Stone 2011).

Following the description to this problem given by Beckwith & Stone (2011), the computational domain is the box $[-1.0, +1.0] \times [-0.5, +0.5]$. We set the density to $\rho = 1$, and the pressure $p = 1$ in the whole computational domain. The components of the initial velocity satisfy

$$v_x = v_0 n_x, \quad v_y = v_0 n_y, \quad v_z = 0, \qquad (94)$$

where $v_0$ is the magnitude of the velocity (we use $v_0 = 2/\sqrt{5}$), $n_x$ and $n_y$ are the components of the unit vector in the direction of movement of the loop (we use $n_x = 2/\sqrt{5}$ and $n_y = 1/\sqrt{5}$). The magnetic field is constant everywhere, except for the loop structure of radius $R = 0.3$). For $r \le R$ we have

$$B_x = -B_0 y/r, \quad B_y = +B_0 x/r, \quad B_z = 0, \qquad (95)$$

where $r = \sqrt{x^2 + y^2}$, $B_0$ is the magnitude of the magnetic field. We use a small value of this quantity in order to keep the magnetic pressure smaller than the gas pressure (we use $B_0 = 10^{-3}$). An adiabatic EOS with $\gamma = 5/3$ is considered to close the system. Periodic boundary conditions are set in all four edges of the computational domain. We run the simulation up to time $t = 5.0$. In the Fig. 9 is depicted the magnetic field magnitude at times $t \in \{0.0, 5.0\}$. We discretize the computational domain into $800 \times 400$ cells. The high order finite-volume method with WENO7 reconstruction is used in smooth parts of the flow. A shock capturing method was not necessary for this test. The divergence cleaning of Dedner et al. (2002) was used, with $c_r = 0.18$, and $c_h$ determined by the maximum propagation speed in the system. Observe that the loop profile is only kept when the divergence cleaning is used in order to preserve the solenoidal constraint.

### 4.2.4 Current sheet

The current sheet problem was first discussed by Hawley & Stone (1995) for the MHD equations, and its extension to SRMHD can be found in Beckwith & Stone (2011). This problem comprises a region which is uniformly filled with a gas at rest. The initial configuration for the magnetic field switches signs at the slices $x = +0.25$ and $x = -0.25$. A perturbation to the system with a sinusoidal velocity function in $y$ is added, which generates non-linear, linearly polarized Alfvén waves. These Alfvén waves turn into magnetosonic waves because the magnetic pressure does not remain constant. Magnetic reconnection occurs because of the two current sheets at $x = \pm 0.25$. Moreover, since the parameter $\beta < 1$, the magnetic reconnection drives highly over-pressurised regions, which launch magnetosonic waves transverse to the field, causing

magnetic energy to be transformed into thermal energy (Hawley & Stone 1995). Close to the place where the magnetic reconnection takes place are produced large magnetic field gradients.

The computational domain is the box $[-0.5, +0.5] \times [-0.5, +0.5]$. The density and pressure are set uniform in the whole domain, with $\rho = 1.0$, and $p = 0.5\beta$, where the parameter $\beta$ represents the ratio of gas pressure to magnetic energy density. The velocity has the profile given by $v_y = v_z = 0$, and $v_x = A \sin(2\pi y)$, where the parameter $A$ is used to test the robustness of the algorithm (Hawley & Stone 1995). The components of the magnetic field are given by $B_x = B_z = 0$, and $B_y = 1$ for $|x| > 0.25$ and $B_y = -1$ otherwise. In our simulation we set $A = 0.1$ and $\beta = 0.1$. An adiabatic EOS with $\gamma = 5/3$ is considered. Periodic boundary conditions are used in all faces of the computational domain. The simulation time is $t = 10.0$. The computational domain is decomposed into $600 \times 600$ cells. The WENO7 reconstruction operator is used along with the WENO3 reconstruction as shock capturing strategy. The divergence cleaning of Dedner et al. (2002) was used, with $c_r = 0.18$, and $c_h$ determined by the maximum propagation speed in the system. In the Fig. 10 are depicted the rest-mass density at simulation times $t \in \{5.0, 7.5\}$. For this very complex problem, the hybrid WENO7/WENO3 reconstruction operator yields excellent results, showing the robustness and stability of the considered numerical scheme for problems involving magnetic reconnection.

### 4.2.5 Orszag–Tang vortex

The Orszag–Tang vortex problem was first studied by Orszag & Tang (1979) for the incompressible MHD equations. Many authors have used this problem for the compressible MHD equations in order to know how robust is the used numerical scheme at handling the formation and the interactions of MHD shocks (Zachary et al. 1994; Ryu & Jones 1995; Ryu et al. 1998; Dai & Woodward 1998; Helzel et al. 2011; Jiang & Wu 1999; Tóth 2000; Londrillo & del Zanna 2000). The initial flow profile consists of smooth initial data, and it is obtained by the superposition of a velocity vortex with a magnetic vortex. This initial configuration is highly unstable, which generates a broad range of MHD waves, which interact with each other, making a transition towards turbulence. The extension to SRMHD has been studied by Dumbser et al. (2008), and Beckwith & Stone (2011).

The computational domain is the box $[0, 1] \times [0, 1]$. The density and pressure are set uniform in the whole domain, with $\rho = \gamma^2$ and $p = \gamma$. With this choice of the density and pressure, the sound speed is $c_s = (\gamma - 1)/\gamma$. The velocity has the following profile

$$v_x = -A \sin(2\pi y), \quad v_y = +A \sin(2\pi x), \quad v_z = 0, \qquad (96)$$

where $A = 0.5$. The magnetic field is given by

$$B_x = -\sin(2\pi y), \quad B_y = +\sin(4\pi x), \quad B_z = 0. \qquad (97)$$

An adiabatic EOS with $\gamma = 5/3$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The simulation time is $t = 1.0$, and in the Fig. 11 are depicted the rest-mass density and thermal pressure at time $t = 1.0$. The computational domain is decomposed into $600 \times 600$ cells. A high order finite-volume scheme with WENO7 reconstruction is used in smooth parts of the flow, and the WENO3 in regions with shocks/discontinuities. The divergence cleaning of Dedner et al. (2002) was used, with $c_r = 0.18$. We set the CFL condition to $C_{CFL} = 0.95$. All used schemes (namely, WENO3, WENO5, and WENO7) are very stable and robust when simulating this very challenging problem. The
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Relativistic magnetic field loop advection. Contour plots of the magnetic pressure without divergence cleaning (top) and with divergence cleaning (bottom) at time $t = 0.0$ (left) and at time $t = 5.0$ (right). The solution was calculated with a finite-volume method with WENO7 reconstruction. The computational domain is the rectangle $[-1.0, +1.0] \times [-0.5, +0.5]$, and it was discretized by using a mesh of 800 × 400 cells.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Relativistic current sheet problem. Contour plots of the rest-mass density at times $t = 5.0$ (left) and $t = 7.5$ (right). The calculations were done with a finite-volume scheme with WENO7 reconstruction. The computational domain is the box $[-0.5, +0.5] \times [-0.5, +0.5]$, discretized by using a mesh of 600 × 600 cells.</figcaption></figure>

XTROEM-FV code was able to capture the shock waves and their interactions with the other flow structures emerging in the evolution of this configuration. In the Fig. 12 is plotted the pressure along the slices $y = 0.4277$ and $y = 0.3125$ at time $t = 1.0$ for different mesh resolutions and different WENO reconstruction operators. We stress the absence of spurious oscillations in these slices for all WENO reconstructions.

### 4.2.6 Cylindrical blast wave

The cylindrical blast wave problem comprises a cylindrical region located in the centre of a domain, and filled with a magnetized

over-pressured gas. After the system is allowed to evolve, a strong shock wave moving outwards is formed. Because of the formation of unphysical values in quantities like the rest-mass density, the pressure, and the magnitude of the velocity, a very robust shock capturing has to be used in order to stabilize the simulation. Due to the periodic boundary conditions, the interactions of the shock waves lead to very complex configurations (Zachary et al. 1994; Londrillo & del Zanna 2000; Stone et al. 2008; Mignone et al. 2010).

For this problem, the computational domain is given by the box $[0, 1] \times [0, 1]$. We set the rest-mass density and pressure uniform in the whole domain, with $\rho = 1.0 \times 10^{-4}$ and $p = 5.0 \times 10^{-4}$.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Relativistic Orszag–Tang vortex problem. Contour plots of the rest-mass density (left) and the pressure (right) at time $t = 1.0$. Computations were performed with a finite-volume scheme with WENO7 reconstruction. The computational domain is the box $[0, 1] \times [0, 1]$, discretized by using a mesh of 600 × 600 cells.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Relativistic Orszag–Tang vortex problem. Plots of the pressure along the slices $y = 0.3125$ (top) and $y = 0.4277$ (bottom) at time $t = 1.0$ for different mesh resolutions and different WENO reconstruction operators.</figcaption></figure>

The pressure in the cylindrical region $(x - x_c)^2 + (y - y_c)^2 < R$, with $(x_c, y_c) = (0.5, 0.5)$, and $R = 0.1$, is $p = 1.0$, and the rest-mass density $\rho = 1.0 \times 10^{-2}$. Initially, the velocity is set to zero, that is $v_x = v_y = v_z = 0$. The magnetic field is set to $B_x = B_y = B_z = 0.1$. An adiabatic EOS with $\gamma = 4/3$ is considered. We consider periodic boundary conditions. The simulation time is $t = 0.5$, and in the Fig. 13 are depicted the rest-mass density, the pressure, the Mach number, and Lorentz factor at time $t = 0.3$. The computational domain is made of 600 × 600 cells. A high-order WENO7 reconstruction is used in smooth parts of the flow, and

the WENO3 in regions with shocks. The divergence cleaning of Dedner et al. (2002) was also used. In the Fig. 14 is plotted the pressure along the slices $y = 0.5$ and $x = 0.5$ at time $t = 0.3$ for different mesh resolutions and different WENO reconstruction operators. Observe the lack of spurious oscillations in these slices for all WENO reconstructions.

### 4.2.7 Rotor problem

The rotor problem was first proposed by Balsara & Spicer (1999) for classical MHD. An extension to SRMHD has been considered by del Zanna et al. (2003), and Dumbser et al. (2008). The problem consists of a high-density, rapidly spinning fluid in a low-density fluid. Initially, both fluids are subject to a uniform magnetic field. Torsional Alfvén waves are launched into the fluid at rest because of the rapidly rotating fluid. Then, the rotor decreases its angular momentum.

The computational domain is the square $[-0.5, +0.5] \times [-0.5, +0.5]$. The rest-mass density and the pressure are $\rho = 1.0$ and $p = 1.0$ in the ambient medium. The cylindrical rotor ($0.0 \le (x - x_c)^2 + (y - y_c)^2 \le 0.1$, with $(x_c, y_c) = (0.5, 0.5)$), is filled with a fluid with rest-mass density $\rho = 10.0$. The pressure inside the rotor is the same as in the ambient fluid. The ambient fluid is initially at rest, that is $v_x = v_y = v_z = 0$. The rotor has an angular velocity $\omega$ such that $v = \omega r = 0.995$ at $r = 0.1$. A linear taper is applied to the velocity and rest-mass density field, however only in a very small range $0 \le r \le 1.115$ so that the density and the velocity match those of the ambient fluid at rest at a radius of $R = 1.115$. The magnetic field is given by $B_x = 5.0$, and $B_y = B_z = 0$ in the whole computational domain. An adiabatic EOS with $\gamma = 7/5$ is considered. We apply periodic boundary conditions at the boundaries of the domain. The computational domain is decomposed into 600 × 600 cells. The simulation was carried out with a finite-volume scheme with WENO7 reconstruction along with WENO3 reconstruction operator in regions with shocks/discontinuities. The divergence cleaning of Dedner et al. (2002) was employed, with $c_r = 0.18$, and $c_h$ determined by the maximum propagation speed in the system. In the Fig. 15 are depicted the rest-mass density,
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Relativistic cylindrical blast wave. Contour plots of the rest-mass density (top-left), pressure (top-right), Mach number (bottom-left), and Lorentz factor (bottom-right) at time $t = 0.3$. The solution was obtained with a finite-volume scheme with WENO7 reconstruction. The computational domain is the box $[0, 1] \times [0, 1]$, discretized by using a mesh of 600 × 600 cells.</figcaption></figure>

the gas pressure, the Mach number, and the velocity magnitude at simulation time $t = 0.4$. We observe that the rotor does launch torsional Alfvén waves, as expected. Inside the rotor, the fluid is still in uniform rotation, as it is shown in the Mach number plot.

In the Fig. 16 is plotted the pressure along the slices $y = 0.5$ and $x = 0.5$ at time $t = 0.4$ for different mesh resolutions and different WENO reconstruction operators. Observe the lack of spurious oscillations in these slices for all WENO reconstructions. Observe also that the profile is almost the same for all schemes and resolutions considered in the plot.

### 4.2.8 Cloud-shock interaction

The cloud-shock interaction problem has been used to simulate the disruption of a high density cloud by a strong shock wave. We follow the setup presented by He & Tang (2012b), which is based on the classical MHD version discussed by Dai & Woodward (1994), and Tóth (2000). The computational domain is the box $[0, 2] \times [0, 1]$.

The discontinuity is located at $x = 1.2$ with the left and right states in primitive variables given by

$$ w_L = \begin{pmatrix} 3.86859 \\ 0.68 \\ 0.0 \\ 0.0 \\ 1.251148954517 \\ 0.0 \\ 0.8498108108786 \\ -0.8498108108786 \end{pmatrix}, \quad w_R = \begin{pmatrix} 1.0 \\ 0.0 \\ 0.0 \\ 0.0 \\ 0.05 \\ 0.0 \\ 0.1610642582333 \\ 0.1610642582333 \end{pmatrix}. \quad (98) $$

The cloud is represented by a dense cylinder, and it is located at $x = 1.6$, $y = 0.5$. The radius of the cylinder is $r = 0.15$, and its rest-mass density is $\rho = 30.0$. The cloud is in hydrostatic equilibrium with the ambient gas. An adiabatic EOS with $\gamma = 5/3$ is considered. Transmissive boundary conditions are applied at the domain boundaries. The simulation time is $t = 2.0$, and in the Fig. 17 are depicted the logarithm of the rest-mass density, the gas pressure, and the Lorentz factor computed with a high order finite-volume
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Relativistic cylindrical blast wave. Plots of the pressure along the slices $y = 0.5$ (top) and $x = 0.5$ (bottom) at time $t = 0.3$ for different mesh resolutions and different WENO reconstruction operators.</figcaption></figure>

scheme with WENO7 reconstruction. The physical domain was decomposed into 800 × 400 cells. The divergence cleaning of Dedner et al. (2002) was used. The CFL condition used in this simulation was $C_{CFL} = 0.95$.

### 4.2.9 Kelvin–Helmholtz instability

For this test, we use the configuration proposed by Beckwith & Stone (2011). The Kelvin–Helmholtz instability consists of a perturbation applied to a system with a velocity shear. Here, we run this test problem to demonstrate the algorithm’s ability to evolve a linear perturbation into non-linear magnetohydrodynamic turbulence (Beckwith & Stone 2011; Radice & Rezzolla 2012).

The computational domain is the box $[-0.5, +0.5] \times [-1.0, +1.0]$. The density is set to $\rho = 1.0$ for $|y| \le 0.25$, and $\rho = 2.0$ for $|y| < 0.25$. The pressure is uniform in the whole domain, with $p = 2.5$. The $x$-component of the velocity satisfies $v_x = 0.5$ if $|y| \le 0.25$, and $v_x = -0.5$ if $|y| < 0.25$. The shear velocity is given by

$$v_x(x) = \begin{cases} +v_{\text{shear}} + A_0 \sin(2\pi x), & \text{if } y \ge 0.25; \\ -v_{\text{shear}} - A_0 \sin(2\pi x), & \text{if } y < 0.25; \end{cases} \quad (99)$$

where $v_{\text{shear}} = 0.5$. The instability is seeded by adding a small perturbation in the transverse component of the velocity,

$$v_y(x) = \begin{cases} +A_0 \sin(2\pi x), & \text{if } y \ge 0.25; \\ -A_0 \sin(2\pi x), & \text{if } y < 0.25; \end{cases} \quad (100)$$

where $A_0 = 0.01$ is the perturbation amplitude. The components of the magnetic field are given by $B_x = 0.2$, and $B_y = B_z = 0$. An adiabatic EOS with $\gamma = 4/3$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The final simulation time is $t = 5$, and in the Fig. 18 is depicted the rest-mass density at times $t \in \{3, 4, 5\}$ obtained with the finite-volume

WENO7 scheme. The computation was carried out on a mesh made of 400 × 800 cells.

## 5 Conclusions

In this work we have discussed the extension of the XTROEM-FV code to the SRHD and SRMHD equations. XTROEM-FV is a new code based on very high order finite-volume methods. Its main goal is solving hyperbolic conservation laws found in computational astrophysics, namely, the MHD, SRHD, and SRMHD equations. In order to achieve very high order of accuracy on Cartesian meshes, XTROEM-FV makes use of high-order WENO reconstruction operators. For handling discontinuities and shocks in an efficient and robust manner, XTROEM-FV employs a strategy based on shock detection and further order reduction of the reconstruction operator in those flow regions with discontinuities. Additional building blocks of the XTROEM-FV code for SRHD and SRMHD include the Rusanov and HLL Riemann solvers. These solvers are very robust when a high-order WENO reconstruction is used, allowing to simulate complex problems in relativistic astrophysics. The time discretization is carried out with a fourth-order SSPRK, and for the numerical treatment of the solenoidal constraint XTROEM-FV uses the GLM divergence cleaning. These ingredients are the base of a very efficient and robust numerical framework for solving the MHD, SRHD, and SRMHD equations under very extreme flow conditions.

We have performed a series of numerical simulations to test the ability of XTROEM-FV to tackle complex flow computations in relativistic astrophysics. Among these tests, one- and two-dimensional problems were considered, especially tests to verify the high-order accuracy of the WENO schemes in smooth flows, one-dimensional Riemann problems, and two-dimensional problems with different initial configurations, which after the time evolves, they present several challenging structures even for low-order schemes, like discontinuities, strong rarefactions, shock waves, and the interaction between all of these features. In this work we mainly report results obtained with a WENO7 reconstruction operator, with the only exception of the slab-jet simulation in SRHD, where we have shown calculations with the WENO5 scheme. All computations confirm the capacity of XTROEM-FV to solve complex flow simulations with very high order finite-volume methods.

Currently, in XTROEM-FV we are implementing the equations of resistive relativistic magnetohydrodynamics along with appropriate time discretizations for this kind of equations, like the IMEX schemes. Additional work is being done in the direction of adaptive mesh refinement in order to efficiently simulate complex and highly dynamical flow structures. Applications of XTROEM-FV for accretion flow problems and (relativistic) magnetohydrodynamic turbulence are currently being carried out, emphasizing in the use of high-order numerical schemes.

## Acknowledgements

This research was funded by the Deutscher Akademischer Austauschdienst under grant A0873943. Further support was provided by the Deutsche Forschungsgemeinschaft, grant MU 1319/12-1. We thank Michael Dumbser for the very careful reading of the manuscript and his valuable suggestions that helped to improve the quality of this article. Many thanks to Luciano Rezzolla, Olindo Zanotti, and Bruno Giacomazzo for providing the exact solution of the one-dimensional Riemann problems discussed in this work.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Relativistic rotor problem. Contour plots of the rest-mass density, pressure, Mach number, and the magnitude of the velocity at time $t = 0.4$. These calculations were performed with an RKFV with hybrid WENO7/WENO3 reconstruction. The computational domain is the box $[0, 1] \times [0, 1]$, discretized by using a mesh of 600 × 600 cells.</figcaption></figure>

REFERENCES

Aloy M. A., Pons J. A., Ibáñez J. M., 1999a, *Comput. Phys. Commun.*, 120, 115

Aloy M. A., Ibáñez J. M., Martí J. M., 1999b, *ApJS*, 122, 151

Anderson M., Hirschmann E., Liebling S., Neilsen D., 2006, *Class. Quantum Grav.*, 23, 6503

Anile A., 1989, *Relativistic Fluids and Magnetofluids*. Cambridge University Press, Cambridge

Anile A., Pennisi S., 1987, *Ann. de l'I.H.P. Physique théorique*, 46, 27

Antón L., Zanotti O., Miralles J. A., Martí J. M., Ibáñez J. M., Font J. A., Pons J. A., 2006, *ApJ*, 637, 296

Antón L., Miralles J. A., Martí J. M., Ibáñez J. M., Aloy M. A., Mimica P., 2010, *ApJS*, 188, 1

Balsara D., 2001, *ApJS*, 132, 83

Balsara D., Shu C.-W., 2000, *J. Comput. Phys.*, 160, 405

Balsara D., Spicer D., 1999, *J. Comput. Phys.*, 149, 270

Beckwith K., Stone J., 2011, *ApJS*, 193, 6

Begelman M., Blandford R., Rees M., 1984, *Rev. Mod. Phys.*, 56, 255

Blandford R., McKee C., 1976, *Phys. Fluids*, 19, 1130

Blandford R., Rees M., 1974, *MNRAS*, 169, 395

Blazek J., 2005, *Computational Fluid Dynamics: Principles and Applications*, 2 edn. Elsevier Science, Amsterdam

Brackbill J., Barnes D., 1980, *J. Comput. Phys.*, 35, 426

Bressan A., 2000, *Hyperbolic Systems of Conservation Laws: the One-Dimensional Cauchy Problem*. Oxford University Press, Oxford

Brio M., Wu C.-C., 1988, *J. Comput. Phys.*, 75, 400

Buchmüller P., Helzel C., 2014, *J. Sci. Comput.*, 61, 343

Casper J., Atkins H., 1993, *J. Comput. Phys.*, 106, 62

Choi E., Ryu D., 2005, *New Astron.*, 11, 116

Chorin A., 1967, *J. Comput. Phys.*, 2, 12

Courant R., Friedrichs K., Lewy H., 1928, *Math. Ann.*, 100, 32

Dai W., Woodward P., 1994, *J. Comput. Phys.*, 115, 485

Dai W., Woodward P., 1998, *ApJ*, 494, 317

Dedner A., Kemm F., Kröner D., Munz C.-D., Schnitzer T., Wesenberg M., 2002, *J. Comput. Phys.*, 175, 645

Dolezal A., Wong S., 1995, *J. Comput. Phys.*, 120, 266

Donat R., Font J. A., Ibáñez J. M., Marquina A., 1998, *J. Comput. Phys.*, 146, 58

Dönmez O., 2006, *Appl. Math. Comput.*, 181, 256

Dönmez O., Kayali R., 2006, *Appl. Math. Comput.*, 182, 1296
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Relativistic rotor problem. Plots of the pressure along the slices $y = 0.5$ (top), and $x = 0.5$ (bottom) at time $t = 0.4$ for different mesh resolutions and different WENO reconstruction operators.</figcaption></figure>

Dumbser M., Käser M., 2007, J. Comput. Phys., 221, 693

Dumbser M., Käser M., Titarev V., Toro E., 2007, J. Comput. Phys., 226, 204

Dumbser M., Balsara D., Toro E., Munz C.-D., 2008, J. Comput. Phys., 227, 8209

Dumbser M., Zanotti O., Hidalgo A., Balsara D., 2013, J. Comput. Phys., 248, 257

Duncan C., Hughes P., 1994, ApJ, 436, L119

Evans C., Hawley J., 1988, ApJ, 332, 659

Font J. A., Ibáñez J. M., Marquina A., Martí J. M., 1994, A&A, 282, 304

Giacomazzo B., Rezzolla L., 2006, J. Fluid Mech., 562, 223

Godunov S., 1959, Mat. Sb., 47, 271

Gottlieb S., 2005, J. Sci. Comput., 25, 105

Gottlieb S., Shu C.-W., 1998, Math. Comp., 67, 73

Gottlieb S., Ketcheson D., Shu C.-W., 2009, J. Sci. Comput., 38, 251

Harten A., 1983, J. Comput. Phys., 49, 357

Harten A., Lax P., van Leer B., 1983, SIAM Rev., 25, 35

Harten A., Engquist B., Osher S., Chakravarthy S., 1987, J. Comput. Phys., 71, 231

Hawley J., Stone J., 1995, Comput. Phys. Commun., 89, 127

He P., Tang H., 2012a, Comm. Comput. Phys., 11, 114

He P., Tang H., 2012b, Comput. Fluids, 60, 1

Helzel C., Rossmanith J., Taetz B., 2011, J. Comput. Phys., 230, 3803

Hesthaven J., Warburton T., 2008, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer, Heidelberg

Hundsdorfer W., Ruuth S., Spiteri R., 2003, SIAM J. Numer. Anal., 41, 605

Jameson A., Schmidt W., Turkel E., 1981, in Proceedings of the AIAA 49th Fluid and Plasma Dynamic Conference. No. AIAA 1981-1259. Palo Alto, California

Jiang G.-S., Wu C.-C., 1999, J. Comput. Phys., 150, 561

Keppens R., Meliani Z., van Marle A.-J., Delmont P., Vlasis A., van der Holst B., 2012, J. Comput. Phys., 231, 718

Komissarov S., 1999, MNRAS, 303, 343

Lax P., Wendroff B., 1960, Comm. Pure Appl. Math., 13, 217

LeVeque R., 1992, Numerical Methods for Conservation Laws, 2 edn. Birkhäuser, Bern

LeVeque R., 2002, Finite-Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Relativistic cloud-shock interaction. Contour plots of the logarithm of the rest-mass density (top), pressure (middle), and Lorentz factor (bottom) at time $t = 2$. The computations were done with a finite-volume scheme with WENO7 reconstruction. The computational domain is the box $[0, 2] \times [0, 1]$, discretized by using a mesh of 800 × 400 cells.</figcaption></figure>

Lichnerowicz A., 1967, Relativistic Hydrodynamics and Magnetohydrodynamics: Lectures On the Existence of Solutions. Benjamin, New York

Liu X.-D., Osher S., Chan T., 1994, J. Comput. Phys., 115, 200

Londrillo P., del Zanna L., 2000, ApJ, 530, 508

Martí J. M., Müller E., 1994, J. Fluid Mech., 258, 317

Martí J. M., Müller E., 1996, J. Comput. Phys., 123, 1

Martí J. M., Müller E., 2003, Living Rev. Relat., 6, 1

Martí J. M., Müller E., 2015, Living Rev. Comput. Astrophys., 1, 1

Martí J. M., Font J. A., Ibáñez J. M., Müller E., 1996, in Hardee P., Bridle A., Zensus J., eds, Astronomical Society of the Pacific Conference Series Vol. 100, Energy Transport in Radio Galaxies and Quasars. Boulder, USA

Mignone A., Bodo G., 2005, MNRAS, 364, 126

Mignone A., Bodo G., 2006, MNRAS, 368, 1040

Mignone A., McKinney J., 2007, MNRAS, 378, 1118

Mignone A., Tzeferacos P., 2010, J. Comput. Phys., 229, 2117

Mignone A., Plewa T., Bodo G., 2005, ApJS, 160, 199

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228

Mignone A., Ugliano M., Bodo G., 2009, MNRAS, 393, 1141

Mignone A., Tzeferacos P., Bodo G., 2010, J. Comput. Phys., 229, 5896

Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G.,
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Relativistic Kelvin–Helmholtz instability. From left to right, contour plots of the rest-mass density at simulation times $t \in \{3.0, 4.0, 5.0\}$. The solution was computed with a high order finite-volume WENO7 scheme on a mesh made of 600 $\times$ 600 cells.</figcaption></figure>

2012, ApJS, 198, 7
Munz C.-D., Schneider R., Sonnendrücker E., Voss U., 1999, C. R. Acad. Sci. Math., 328, 431
Noble S., Gammie C., McKinney J., Zanna L. D., 2006, ApJ, 641, 626
Noh W., 1987, J. Comput. Phys., 72, 78
Núñez-de la Rosa J., Munz C.-D., 2016, MNRAS, 455, 3458
Orszag S., Tang C.-M., 1979, J. Fluid Mech., 90, 129
Peterson B., 1997, An Introduction To Active Galactic Nuclei. Cambridge University Press, Cambridge
Pons J. A., Martí J. M., Müller E., 2000, J. Fluid Mech., 422, 125
Powell K., 1994, Technical report, An approximate Riemann solver for magnetohydrodynamics. Institute for Computer Applications in Science and Engineering, NASA Langley Research Center
Radice D., Rezzolla L., 2011, Phys. Rev. D, 84, 024010
Radice D., Rezzolla L., 2012, A&A, 547, A26
Rezzolla L., Zanotti O., 2001, J. Fluid Mech., 449, 395
Rezzolla L., Zanotti O., 2013, Relativistic Hydrodynamics. Oxford University Press, Oxford
Rezzolla L., Zanotti O., Pons J. A., 2003, J. Fluid Mech., 479, 199
Rusanov V., 1961, J. Comput. Math. Phys. (USSR), 1, 267
Ruuth S., Spiteri R., 2002, J. Sci. Comput., 17, 211
Ryu D., Jones T., 1995, ApJ, 442, 228
Ryu D., Miniati F., Jones T., Frank A., 1998, ApJ, 509, 244
Ryu D., Chattopadhyay I., Choi E., 2006, ApJS, 166, 410
Schneider V., Katscher U., Rischke D., Waldhauser B., Maruhn J., Munz C.-D., 1993, J. Comput. Phys., 105, 92
Shu C.-W., 1988, SIAM J. Sci. Comput., 9, 1073
Shu C.-W., 2009, SIAM Rev., 51, 82
Shu C.-W., Osher S., 1988, J. Comput. Phys., 77, 439
Shu C.-W., Osher S., 1989, J. Comput. Phys., 83, 32
Spiteri R., Ruuth S., 2002, SIAM J. Numer. Anal., 40, 469
Spiteri R., Ruuth S., 2003, Math. Comput. Simulat., 62, 125
Stone J., Gardiner T., Teuben P., Hawley J., Simon J., 2008, ApJS, 178, 137
Tchekhovskoy A., McKinney J., Narayan R., 2007, MNRAS, 379, 469
Titarev V., Toro E., 2004, J. Comput. Phys., 201, 238

Titarev V., Toro E., 2005, J. Comput. Phys., 204, 715
Toro E., 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3 edn. Springer, Berlin
Tóth G., 2000, J. Comput. Phys., 161, 605
Wang Z.-J., et al., 2013, Int. J. Numer. Meth. Fluids, 72, 811
Woodward P., Colella P., 1984, J. Comput. Phys., 54, 115
Zachary A., Malagoli A., Colella P., 1994, SIAM J. Sci. Comput., 15, 263
Zhang W., MacFadyen A., 2006, ApJS, 164, 255
del Zanna L., Bucciantini N., 2002, A&A, 390, 1177
del Zanna L., Bucciantini N., Londrillo P., 2003, A&A, 400, 397
del Zanna L., Zanotti O., Bucciantini N., Londrillo P., 2007, A&A, 473, 11
van Leer B., 1979, J. Comput. Phys., 32, 101

This paper has been typeset from a TeX/LATeX file prepared by the author.